Exact derivation and practical application of a hybrid stochastic simulation algorithm for large gene regulatory networks
EExact derivation and practical application of a hybrid stochasticsimulation algorithm for large gene regulatory networks
Jaroslav [email protected]
Abstract
We present a highly efficient and accurate hybrid stochastic simulation algorithm (HSSA) for thepurpose of simulating a subset of biochemical reactions of large gene regulatory networks (GRN).The algorithm relies on the separability of a GRN into two groups of reactions, A and B, suchthat the reactions in A can be simulated via a stochastic simulation algorithm (SSA), while thosein group B can yield to a deterministic description via ordinary differential equations. First, wederive exact expressions needed to sample the next reaction time and reaction type, and then givetwo examples of how a GRN can be partitioned. Although the methods presented here can beapplied to a variety of different stochastic systems within GRN, we focus on simulating mRNAsin particular. To demonstrate the accuracy and efficiency of this algorithm, we apply it to athree-gene oscillator, first in one cell, and then in an array of cells (up to 64 cells) interacting viamolecular diffusion, and compare its performance to the Gillespie algorithm (GA). Depending onthe particular numerical values of the system parameters, and the partitioning itself, we show thatour algorithm is between 11 and 445 times faster than the GA. a r X i v : . [ q - b i o . M N ] S e p NTRODUCTION
The business of modeling chemical reaction networks generally falls into two camps: 1)chemical master equation (CME); and 2) stochastic simulation algorithms (SSA). Histori-cally, research in the former has centered on solving the CME using analytic methods, bothexact and approximate [1–6], or by developing efficient numerical techniques [7–15]. Thelatter, popularly known and the Gillespie algorithm (GA) [16], has emerged as an alternativeto the CME and is generally more preferred for modeling systems with many species and/orreaction channels. Although the GA is simple in construction and guarantees a solution, ittends to be computationally expensive. Hence, research on faster versions of the GA hasbeen ongoing [17–21]. Although the choice between the CME and an SSA depends on sev-eral factors, such as the complexity of the reaction network, efficiency of each approach etc.,the general rule is this: use the CME if it can be solved either analytically or numerically;otherwise employ an SSA. In gene regulatory networks (GRN), which will be the focus ofthis paper, most of the systems of interest require the use of an SSA.In recent years, a new camp has emerged in which the CME and the GA are combinedinto a hybrid [22–30]. The purpose of hybrid models is to avoid stochastically simulatingevery reaction and instead to use ordinary differential equations (ODE), e. g. the CME,to describe a subset of reactions, which can lead to a significant increase in computationalspeed. If the ODE for such a subset can be solved efficiently (preferably analytically),then the numerical cost would be incurred only by the remaining reactions. The first hybridalgorithm was proposed by Haseltine and Rawlings [22], and later improved upon by Burrage et. al. [23] and Salis and Kaznessis [24]. The principle behind both works is the segregationof reactions into fast and slow reactions, wherein the fast reactions are simulated by theLengevin equation, while the slow reactions are handled by an SSA. The applicability of thisalgorithm, of course, is limited to systems that contain both fast and slow reactions. Jahnkeand Altntan [25] developed a hybrid method that splits the system into two (or more) subsetsof reactions in such a way that the CME for one subset can be solver analytically, while thereactions in the other subset are simulated in a τ -leaping fashion (see reference [18]). Lateron, other hybrid models, ones that could be implemented for any set of reactions, fast orotherwise, were developed [26–30]. In these later models, the CME is used to describe asubset of reactions, while a modified SSA is applied to the remaining set. The main obstacle2o this approach is to work out the probability distributions for the next reaction time andthe reaction index of the latter subset, as they both depend on the state history of thespecies described by the CME. One way to solve this issue is to treat the first subset asdeterministic, thereby collapsing the system’s many possible histories into one, as was donein [30]. However, this approach restricts the ways a system can be partitioned into thosethat lend themselves to such approximations. Duso and Zechner [29] developed a hybridmodel that coupled statistical moments of the chemical species in the first subset to theprobabilities of the reaction time and reaction type of the latter subset.In this paper, we build on our earlier work [26, 27] in which we derived the probabilityfor the next reaction time and the probabilities for the reaction types for the letter reactionsubset. However, the method was demonstrated only on systems without any reactionsin which two molecules combine into another molecular species, e.g. dimerization. In thepresent work, we demonstrate how any reaction network can be partitioned into a subsetof reactions describable by a CME and a subset that is simulated with an SSA by derivingexact: probability distribution for the next reaction time; probability for the reaction type;and the conditional joint probability for all the species in the first subset, given the samplednext reaction time and reaction type. Although the formalism of our approach is exact andworks for any reaction network, the application of it varies in difficulty depending on themanner in which a network is partitioned. Therefore, to demonstrate our method in practice,we show two examples of how a system of three interacting genes can be partitioned in orderto render all calculations tractable. The dynamics of the system we chose is oscillatory, whichis suitable for testing the effects of systematic errors, as these types of error would manifestthrough phase shifts. Our results show excellent accuracy for all parameter sets and anincrease in efficiency by factors ranging from 11 to 445. MASTER EQUATION AND THE GILLESPIE ALGORITHM
Given a set of molecular copy numbers x = { x , ..., x V } , and a set of reactions V (cid:88) i =1 α iµ x i a µ ( x , t ) −−−−→ V (cid:88) i =1 β iµ x i µ = 1 , , ..., J (1)where the integers α ik and β ik are the stoichiometric coefficients and a µ ( x , t ) are the reactionpropensities, the joint probability distribution, P ( x , t ), for x is a solution of the chemical3aster equation (CME) dP ( x , t ) dt = J (cid:88) µ =1 a µ ( x − f µ , t ) P ( x − f µ , t ) − P ( x , t ) J (cid:88) µ =1 a µ ( x , t ) . (2)The state-change vector f µ specifies the change in x due to the µ th reaction: x → x + f µ .Except for a very specific class of systems, Eq. (2) cannot be solved exactly. Attempts tosolve it approximately, using either analytic or numerical methods, are usually thwarted bythe curse of dimensionality , which tends to cast itself on many-variable systems.An alternative approach is to simulate the reaction events in Eq. (1) by sampling the timebetween reactions, τ , and the reaction index µ = 1 , ..., J from a joint probability distribution P ( µ ; τ ). The time is then advanced by τ and the system variables are updated based onwhich reaction took place. This procedure is called the Gillespie algorithm (GA) and it canbe derived as follows:The probability that no reaction occurs within an infinitesimal time interval dt is givenby P ( dt ) = (1 − R ( x , dt ) , (3)where R ( x , t ) = (cid:88) µ a µ ( x , t ) . (4)The probability that no reaction occurs within a finite time t = dtN is simply P ( t ) = N (cid:89) n =1 (1 − R ( x , t n ) dt ) = exp (cid:20) − (cid:90) t dt (cid:48) R ( x , t (cid:48) ) dt (cid:21) , (5)where t n = ndt . The probability that no reaction occurs up to t and that reaction µ occursbetween t and t + dt is given by P ( µ ; t ) = P ( t ) a µ ( x , t ) dt. (6)Since t can be sampled independently of µ , via Eq. (5), µ must be sampled from a conditionalprobability P ( µ | τ ) that µ occurs provided t has been observed to be some time τ . Accordingto the Bayes relation, P ( µ | τ ) can be expressed as P ( µ | τ ) = P ( µ ; τ ) (cid:80) ν P ( ν ; τ ) = a µ ( x , τ ) R ( x , τ ) . (7)4or systems in which the reaction propensities are time-independent, the two probabilitiesreduce to the well known expressions: P ( τ ) = e − R ( x ) τ probability to observe t = τ (8) P ( µ | τ ) = a µ ( x ) R ( x ) probability to observe reaction µ , given τ (9)With these simple relations, the GA can be implemented as follows:1. At t = 0 choose an initial state x and compute the propensities a µ ( x ).2. Select two random numbers ξ and ξ .3. Compute τ by solving P ( τ ) = ξ .4. Find the smallest integer µ that satisfies µ (cid:88) ν =1 P ( ν | τ ) > ξ . (10)5. Update the variables by letting x → x + f µ and set t to t + τ .6. Return to step 1.A solution to Eq. (2) can be obtained by running the GA enough times to generate astatistically significant ensemble. While this process guarantees a solution, it can be verytime consuming. A HYBRID APPROACH
Let us consider a reaction network that can be partitioned into two sets of reactions, A = { a ( x ) , ..., a κ ( x ) } and B = { a κ +1 ( x ) , ..., a J ( x ) } , such that a subset of variables y = { y = x , ..., y M = x M } is affected only by reactions in A , while the remaining set z = { z = x M +1 , ..., z V − M = x V } can be affected by reactions in both A and B . We define g µ asthe state-change vector that affects y only, and h µ as the state-change vector that affectsexclusively z .Let us now imagine that we know the exact stochastic evolution of the variable set z during some time t ; in other words, we know the path of each variable in z in the t − z plane.We could then ask: what are the probabilities that 1) no reaction in A occurs for during t ;5nd 2) reaction µ occurs between t and t + dt ? The former can be constructed the same wayas in Eq. (5): P ( t | Z ) = N (cid:89) n =1 [1 − R A ( z ( t n ) , y ) dt ] = exp (cid:20) − (cid:90) t R A ( z ( t (cid:48) ) , y ) dt (cid:48) (cid:21) , (11)where R A ( z ( t ) , y ) = κ (cid:88) ν =1 a ν ( z ( t ) , y ) (12)and Z is the set of paths taken by the variables in z , i.e. Z = { z ( t ) , ..., z ( t N ) (cid:124) (cid:123)(cid:122) (cid:125) ; z ( t ) , ..., z ( t N ) (cid:124) (cid:123)(cid:122) (cid:125) ; ... ; z V − M ( t ) , ..., x V − M ( t N ) } path of z path of z ... (13)Since we do not actually know the paths, i. e. the values of all the entries in Z , we mustmultiply Eq. (11) by the probability for Z , P ( Z ), and sum over all the variables z i ( t n ),except for the endpoints of each path, z (0) = z and z ( t ) = z t , which we keep fixed: Q ( z t , t | z ) = (cid:88) Z (cid:54) = z , z t P ( Z ) P ( τ | Z ) . (14)If the initial set z is known only with some probability P in ( z ), we must multiply Eq. (14)by P in ( z ) and sum over z : Q ( z t , t ) = (cid:88) z P in ( z ) Q ( z t , t | z ) . (15)We will refer to Eq. (15) as the Q -distribution; it is the joined probability that no reactionoccurs until t and that the variables z will have the values z t at t . Finally, the probabilitythat no reaction occurs until t is given by Q ( t ) = (cid:88) z t Q ( z t , t ) . (16)Next, we would like to compute the conditional probability, Q ( µ | t ), that reaction µ occursbetween t and t + dt , provided no reaction occured until t . This can be obtained from thejoined probability, Q ( z t , µ, t ), that no reaction occurs until t , that z = z t at t , and thatreaction µ occurs between t and t + dt , which is simply this: Q ( z t , µ, t ) = Q ( z t , t ) a µ ( z t , y ) dτ. (17)6o, the conditional probability that reaction µ occurs, provided no reaction occured up to t , reads Q ( µ | t ) = (cid:80) z t Q ( z t , µ, t ) (cid:80) z t ,µ Q ( z t , µ, t ) . (18)To update the initial probability, P in ( z ), after t and µ are sampled, we must compute theconditional probability to observe z t , provided no reaction occured until t and reaction µ occured at t , which is given by Q ( z t | , µ, t ) = Q ( z t , µ, t ) (cid:80) z t Q ( z t , µ, t ) = Q ( z t , t ) a µ ( z t ) (cid:80) z t Q ( z t , t ) a µ ( z t ) , (19)and set z t = z . Note that this “update” only considers how P in ( z ) changes given newinformation, namely that reaction µ has occurred at t = τ . However, we must also takeinto account how the reaction itself affects P in ( z ). For example, if a transcription factor(TF) bonded to a promoter, there would be one less TF available to bind to promoters.The probability of there being n number of available TFs would be shifted to the leftalong the n -axis. This shift is simple when the distribution is far enough from zero, i. e.when P ( n = 0) (cid:28) P ( n ) → P ( n + 1). For a general case, one must shift P ( n ) to theleft by one, but add P ( n = 1) to P ( n = 0): P ( n = 0) → P ( n = 0) + P ( n = 1) and P ( n = i ) → P ( n = i + 1) for i >
0. On the other hand, when a TF dissociates form apromoter, thereby adding a TF to the system, the shift is always P ( n ) → P ( n − S µ , whose action will update appropriately P in ( z ) asa function of µ .The last piece of the puzzle is to figure out how to compute Eq. (15). Fortunately, thishas already been done in a previous work of ours [26], in which we showed that Eq. (15) isthe solution of this equation: dQ ( z , t ) dt = J (cid:88) ν = κ +1 a ν ( z − h ν , t ) Q ( z − h ν , t ) − Q ( z , t ) J (cid:88) ν = κ +1 a ν ( z , t ) − R A ( z , y ) Q ( z , t ) (20)with the initial conditions Q ( z ,
0) = P in ( z ). Note that if we set R A to zero, Eq. (20) wouldreduce to the CME for B . Hence, we will refer to Eq. (20) as modified chemical masterequation (MCME).We now have all the ingredients to implement the hybrid stochastic simulation algorithm(HSSA). The steps are as follows:1. At t = 0 choose an initial state y and initial probability P in ( z ).7. Select two random numbers ξ and ξ .3. Solve Eq. (20) for Q ( z , t ), and compute Q ( t ), Q ( ν | t ) and Q ( z t | , µ, τ ) from Eqs. (16),(18) and (19), respectively.4. Compute τ by solving Q ( τ ) = ξ .5. Find the smallest integer µ that satisfies µ (cid:88) ν =1 Q ( ν | τ ) > ξ . (21)6. Update the variables y by letting y → y + g µ and set t to t + τ .7. Update the initial probability: P in ( z ) = ˆ S µ Q ( z | , µ, τ ).8. Return to step 1.The above steps give a recipe on how to implement an HSSA. However, step 3 relies onone important assumption, namely, that Eq. (20) is tractable. This is where the mannerin which a reaction network is partitioned becomes important for practical reasons. If thenumber and nature of reactions in B are such that Eq. (20) can be solved either analyticallyor numerically, then the HSSA can be implemented in practice. Furthermore, since thepurpose of deriving the HSSA is to outperform conventional SSAs, such as the GA, step3 must not only be feasible, but needs to be done efficiently. These considerations willessentially dictate how a reaction networks is partitioned. In the next section, we give twoexamples of how a reaction system can be partitioned. The last section will be reserved fordiscussing the performances of said examples. PARTITIONING OF THE SYSTEM: PRACTICAL APPLICATIONS
To make our examples concrete, we will focus exclusively on gene regulatory networks(GRN). 8he system we want to partition, shown in Fig. 1, is made up of the following variables: w kj promoter of gene k in state j range: 0 or 1 m k mRNA copy number of gene k range: 0 to ∞ n k protein copy number of gene k range: 0 to ∞ s k dimer copy number composed of 2 n k range: 0 to ∞ ; (22)and its evolution is driven by these reactions:1 . w kj + s q −→ w kj (cid:48) a jj (cid:48) k w kj s q . w kj (cid:48) −→ w kj + s q b jj (cid:48) k w kj (cid:48) . ∅ −→ m k r k w k . m k −→ ∅ d mk m k . m k −→ m k + n k K mk m k . n k −→ ∅ d nk n k . n k −→ s k c + k n k ( n k − / . s k −→ n k c − k s k . s k −→ ∅ d sk s k . (23)Reaction 1 takes a promoter of gene k from the state labeled by j to the state labeled by j (cid:48) ; while reaction 2 is the reverse of reaction 1 (in general, the state label j (cid:48) in a jj (cid:48) k and b jj (cid:48) k depends on the index q ). Reactions 3 and 4 synthesize and degrade m k , respectively.Reactions 5 and 6 synthesize and degrade n k , respectively. Reactions 7 and 8 bind two copiesof n k to form one copy of s k and dissociate s k into two copies of n k , respectively. Finally,reaction 9 degrades s k . The list of reactions (23) is very standard in systems biology, but isnot exhaustive. We could have added the formation of (homo/hetero) trimers, quadrimersand higher oligamers, post-transcription and post-translation processes, transport across thenuclear membrane and others. For simplicity however, we will stick to the reactions in (23).In the following examples, we will focus on obtaining the stochastic behavior of thevariable set y only; in the process, the information about the set z will be relegated to itsaverage: (cid:104) z (cid:105) . 9 IG. 1. A system of three interacting genes. The lower index of all parameters refers to the genein question. The upper indices of a ijk correspond to the transitions between promoter states (leftto right); same notation applies to b ijk . The parameters r k and K mk give the transcription rate andtranslation rate, respectively, while in that same order, c + k and c − k represent the rate of homodimerformation and dissociation. Finally, d mk , d nk and d sk are the degradation rates of mRNA, proteinand homodimer, respectively. xample 1 Let us partition the system as follows:group A group B . ∅ −→ m k . w k + s k −→ w k . m k −→ ∅ . w k + s k −→ w k . w k −→ w k + s k . w k −→ w k + s k . m k −→ m k + n k . n k −→ ∅ . n k −→ s k . s k −→ n k . s k −→ ∅ , (24)The propensities for A , a µ ( w , w , w ), are given by a = r δ w a = r δ w a = r δ w a = d m m a = d m m a = d m m . (25)Before we deal with the modified CME (MCME) for B , let us first write down the CME: d P dt = W ( s ) P + (cid:88) k K k m k (cid:2) P ( n k − − P (cid:3) + d nk (cid:2) ( n k + 1) P ( n k + 1) − n k P (cid:3) + (cid:88) k c + k (cid:2) ( n k + 2)( n k + 1) P ( n k + 2 , s k − − n k ( n k + 1) P (cid:3) + (cid:88) k c − k (cid:2) ( s k + 1) P ( n k − , s k + 1) − s k P (cid:3) + (cid:88) k d sk (cid:2) ( s k + 1) P ( s k + 1) − s k P (cid:3) , (26)11here P is a vector whose dimension is equal to the number of all unique combinations ofpromoter states ( w , w , w ) and the elements of the matrix W ( s ) give the propensities forthe transitions between the combinations. We have employed a short hand notation in which P is short for P ( n , n , n , s , s , s , t ), P ( n + 2 , s −
1) is short for P ( n + 2 , n , n , s − , s , s , t ), etc.Note that reactions 1 to 4 in B can be approximately decoupled from reactions 5 to 9 onaccount of the fact that reactions 1 and 2 change the copy number of s k cyclically and atmost by two: they can happen one after the other before the only available reaction channelinvolving the promoter are reactions 3 and 4, which return one or both copies of s k to thesystem. Hence, if s k (cid:29) k = 1 , ,
3, the CME for reactions 5 to 9 reads: dPdt = (cid:88) k K k m k (cid:2) P ( n k − − P (cid:3) + d nk (cid:2) ( n k + 1) P ( n k + 1) − n k P (cid:3) + (cid:88) k c + k (cid:2) ( n k + 2)( n k + 1) P ( n k + 2 , s k − − n k ( n k + 1) P (cid:3) + (cid:88) k c − k (cid:2) ( s k + 1) P ( n k − , s k + 1) − s k P (cid:3) + (cid:88) k d sk (cid:2) ( s k + 1) P ( s k + 1) − s k P (cid:3) . (27)Although this equation cannot be solved exactly, it can be shown that the variances of n k and s k are very close to Poisson [31]. Thus, provided the averages (cid:104) n k (cid:105) and (cid:104) s k (cid:105) are muchgreater than 1, we can make the approximation (cid:88) n , s f ( n , s ) P ( n , s ) = f ( (cid:104) n (cid:105) , (cid:104) s (cid:105) ) , (28)for any smooth function f . Thus, if we write the joint probability P ( l, n , s ), where the index l labels the elements of P , as P ( l, n , s ) = P ( l | n , s ) P ( n , s ) , (29)then, it follows from relation (28) that (cid:88) n , s (cid:88) q W lq ( s ) P ( q, n , s ) = (cid:88) n , s (cid:34)(cid:88) q W lq ( s ) P ( q | n , s ) (cid:35) P ( n , s )= (cid:88) q W lq ( (cid:104) s (cid:105) ) P ( q |(cid:104) n (cid:105) , (cid:104) s (cid:105) ) = ( ¯WP ) l , (30)12here ¯W = W ( (cid:104) s (cid:105) ). Hence, summing both sides of Eq. (26) over n and s , we obtain d P dt = ¯WP . (31)The average (cid:104) s k (cid:105) ( (cid:104) n k (cid:105) ) can be obtained by multiplying Eq. (27) by s k ( n k ), summing over n and s , and remembering the relation (28): ddt (cid:104) s k (cid:105) = c + k (cid:104) n k (cid:105) ( (cid:104) n k (cid:105) − − ( c − k + d sk ) (cid:104) s k (cid:105) ddt (cid:104) n k (cid:105) = K k m k + 2 c − k (cid:104) s k (cid:105) − c + k (cid:104) n k (cid:105) ( (cid:104) n k (cid:105) − − d nk (cid:104) n k (cid:105) . (32)Since we are treating n and s as deterministic, the CME for the entire system B , is Eq.(31). For convenience, let us write it component by component: ddt P ( w , w , w ) = (cid:88) w (cid:48) ,w (cid:48) ,w (cid:48) ¯ W ( w , w , w ; w (cid:48) , w (cid:48) , w (cid:48) ) P ( w (cid:48) , w (cid:48) , w (cid:48) ) . (33)Since we know that a promoter state of one gene does not influence a promoter state ofanother gene, P ( w , w , w ) must be a product, P ( w ) P ( w ) P ( w ), such that P k ( w k )satisfies d P k dt = M k P k , (34)where M k = − a k (cid:104) s k (cid:105) b k a k (cid:104) s k (cid:105) − a k (cid:104) s k (cid:105) − b k b k a k (cid:104) s k (cid:105) − b k . Here, a k = a k = a k and b k = b k = b k . Taking a derivative of P ( w ) P ( w ) P ( w ), weobtain ddt P ( w , w , w ) = (cid:88) w (cid:48) M ( w , w (cid:48) ) P ( w (cid:48) , w , w ) + (cid:88) w (cid:48) M ( w , w (cid:48) ) P ( w , w (cid:48) , w )+ (cid:88) w (cid:48) M ( w , w (cid:48) ) P ( w , w , w (cid:48) ) . (35)Following Eq. (20), the MCME acquires the form ddt Q ( w , w , w ) = (cid:88) w (cid:48) M ( w , w (cid:48) ) Q ( w (cid:48) , w , w ) + (cid:88) w (cid:48) M ( w , w (cid:48) ) Q ( w , w (cid:48) , w )+ (cid:88) w (cid:48) M ( w , w (cid:48) ) Q ( w , w , w (cid:48) ) − (cid:88) i ( r i δ w i , + d mi m i ) Q ( w , w , w ) . (36)13e can simplify this equation by defining a new function ˜ Q ( w , w , w ), such that Q ( w , w , w ) = exp (cid:34) − (cid:88) i d mi m i t (cid:35) ˜ Q ( w , w , w ) . (37)Hence, ddt ˜ Q ( w , w , w ) = (cid:88) w (cid:48) M ( w , w (cid:48) ) ˜ Q ( w (cid:48) , w , w ) + (cid:88) w (cid:48) M ( w , w (cid:48) ) ˜ Q ( w , w (cid:48) , w )+ (cid:88) w (cid:48) M ( w , w (cid:48) ) ˜ Q ( w , w , w (cid:48) ) − (cid:88) i r i δ w i , ˜ Q ( w , w , w ) . (38)Since the initial conditions demand that Q ( w , w , w , t = 0) = ˜ Q ( w , w , w , t = 0) = P in, ( w , P in, ( w , P in, ( w , Q ( w , w , w ) can also be written as a product ˜ Q ( w ) ˜ Q ( w ) ˜ Q ( w ),which, when inserted into Eq. (38) yields˜ Q ( w ) ˜ Q ( w ) (cid:34) ddt ˜ Q ( w ) − (cid:88) w (cid:48) M w ,w (cid:48) ˜ Q ( w (cid:48) ) + r δ w , ˜ Q ( w ) (cid:35) +˜ Q ( w ) ˜ Q ( w ) (cid:34) ddt ˜ Q ( w ) − (cid:88) w (cid:48) M w ,w (cid:48) ˜ Q ( w (cid:48) ) + r δ w , ˜ Q ( w ) (cid:35) +˜ Q ( w ) ˜ Q ( w ) (cid:34) ddt ˜ Q ( w ) − (cid:88) w (cid:48) M w ,w (cid:48) ˜ Q ( w (cid:48) ) + r δ w , ˜ Q ( w ) (cid:35) = 0 . (39)Thus, our Q -distribution is given by Q ( w , w , w , t ) = exp (cid:34) − (cid:88) i d mi m i t (cid:35) ˜ Q ( w , t ) ˜ Q ( w , t ) ˜ Q ( w , t ) , (40)with ˜ Q k ( w k ) satisfying ddt ˜Q k = [ M k − R k ] ˜Q k , (41)where R k = r k . We are now ready to compute Q ( t ), Q ( ν | t ) and Q ( z t | , µ, τ ) (step 3 of the HSSA). Thefirst one is simply Q ( t ) = exp (cid:34) − (cid:88) i d mi m i t (cid:35) ˜ Q ( t ) ˜ Q ( t ) ˜ Q ( t ) , (42)14here ˜ Q k ( t ) = (cid:88) w k ˜ Q k ( w k , t ); (43)the second one reads Q ( µ | τ ) = (cid:88) w w w ˜ Q ( w , τ ) ˜ Q ( w , τ ) ˜ Q ( w , τ ) a µ ( w , w , w ) / Q ( τ ) , (44)where Q ( τ ) = (cid:88) µ (cid:88) w w w ˜ Q ( w , τ ) ˜ Q ( w , τ ) a µ ( w , w , w ); (45)while the third one is given by Q ( w , w , w | µ, τ ) = ˜ Q ( w , τ ) ˜ Q ( w , τ ) a µ ( w , w , w ) / Q ( τ ) . (46)Inserting the propensities listed in Eq. (25) into (44), we obtain Q (1 | τ ) = r ˜ Q (1 , τ ) ˜ Q ( τ ) ˜ Q ( τ ) Q ( τ ) Q (2 | τ ) = r ˜ Q ( τ ) ˜ Q (1 , τ ) ˜ Q ( τ ) Q ( τ ) Q (3 | τ ) = r ˜ Q ( τ ) ˜ Q ( τ ) ˜ Q (1 , τ ) Q ( τ ) Q (4 | τ ) = d m m ˜ Q ( τ ) ˜ Q ( τ ) ˜ Q ( τ ) Q ( τ ) Q (5 | τ ) = d m m ˜ Q ( τ ) ˜ Q ( τ ) ˜ Q ( τ ) Q ( τ ) Q (6 | τ ) = d m m ˜ Q ( τ ) ˜ Q ( τ ) ˜ Q ( τ ) Q ( τ ) (47)Similarly, if we insert the explicit reaction propensities into (46), we obtain the updates forthe initial probability distributions, P in ,k ( w k ,
0) = Q k ( w k , τ ): µ = 1 : Q ( w | , τ ) = δ w , Q ( w | , τ ) = ˜ Q ( w , τ )˜ Q ( τ ) , Q ( w | , τ ) = ˜ Q ( w , τ )˜ Q ( τ ) ,µ = 2 : Q ( w | , τ ) = δ w , Q ( w | , τ ) = ˜ Q ( w , τ )˜ Q ( τ ) , Q ( w | , τ ) = ˜ Q ( w , τ )˜ Q ( τ ) ,µ = 3 : Q ( w | , τ ) = δ w , Q ( w | , τ ) = ˜ Q ( w , τ )˜ Q ( τ ) , Q ( w | , τ ) = ˜ Q ( w , τ )˜ Q ( τ ) ,µ = 4 , , Q ( w | µ, τ ) = ˜ Q ( w , τ )˜ Q ( τ ) , Q ( w | µ, τ ) = ˜ Q ( w , τ )˜ Q ( τ ) , Q ( w | µ, τ ) = ˜ Q ( w , τ )˜ Q ( τ ) (48)15ince none of the reactions in A change either of the copy numbers n or s , the shiftingoperator ˆ S µ is just unity.In principle, we could now run the HSSA as prescribed in the previous section. However,in order to make it as efficient as possible, we must find a way to deal with Eqs. (32) and(41). Although in general Eqs. (32) cannot be solved analytically, we can exploit the factthat τ will always be small compared to some characteristic time during which (cid:104) n k (cid:105) and (cid:104) s k (cid:105) will change significantly. For such τ we can write (cid:104) n k (cid:105) and (cid:104) s k (cid:105) as polynomials on τ : (cid:104) n k (cid:105) = (cid:88) q =0 u kq τ q (cid:104) s k (cid:105) = (cid:88) q =0 v kq τ q . (49)Plugging these into Eqs. (32), we obtain algebraic equations for u kq and v kq , which, up toorder τ , read v k = c + k u k ( u k − − ( c − k + d sk ) v k u k = K k m k + 2 c − k v k − c + k u k ( u k − − d nk u k v k = c + k u k u k − u k ) − ( c − k + d sk ) v k u k = 2 c − k v k − c + k (2 u k u k − u k ) − d nk u k . (50)The smallness of τ can also be useful for simplifying Eq. (41). If we write the matrix M k as M k = M (0) k + ∆ k ( t ), where M (0) k = − a k u k b k a k u k − a k u k − b k b k a k u k − b k , and ∆ k ( t ) = (cid:88) q =1 a k u kq t q − − , Eq. (41) becomes ddt ˜Q k = [( M (0) k − R k ) + ∆ k ( t )] ˜Q k . (51)Since ∆ k ( t ) is a correction to the matrix M (0) k − R k , we can write ˜Q k = ˜Q (0) k + ˜Q (1) k + ... ,where the superscript indicates the order on ∆ k ( t ), i. e. ˜Q ( n ) k ∼ O ( ∆ k ( t ) n ). Plugging ˜Q k ∆ k ( t ), we obtainthis series of equations: ddt ˜Q (0) k = ( M (0) k − R k ) ˜Q (0) k ddt ˜Q ( n ) k = ( M (0) k − R k ) ˜Q ( n ) k + ∆ k ( t ) ˜Q ( n − k ∀ n > . (52)These equations can be solved analytically for an arbitrary n in terms of the solution to ˜Q (0) k , which is ˜Q (0) k ( t ) = (cid:88) p e E kp t (cid:0) U k B p U − k (cid:1) P in ,k , (53)where E kp and U k are the p th eigenvalue and the column eigenvectors of M (0) k − R k , respec-tivelly, and B p = δ p δ p
00 0 δ p . For the purpose of testing, we will consider Eq. (53) to be our Q -distribution. Plugging Eq.(53) into Eq. (42), we obtain Q ( t ) = exp (cid:34) − (cid:88) i d mi t (cid:35) (cid:16) tr ˜Q (0)1 ( t ) (cid:17) (cid:16) tr ˜Q (0)2 ( t ) (cid:17) (cid:16) tr ˜Q (0)3 ( t ) (cid:17) . (54) τ can now be computed by choosing a random real number ξ = (0 ,
1] and solving Q ( τ ) = ξ for τ . Once we have a value of τ , expressions in (47) and (48) can be computed from thisrelation: ˜ Q k ( w k , τ ) = ( e w k ) T ˜Q (0) k ( τ ) , (55)where e i = δ i δ i δ i . In order to speed up the numerical search for τ we first computed the average τ , τ Av = (cid:90) ∞ dtQ ( t ) , (56)and then evaluated Q ( τ ) − ξ for τ = τ Av n/
10 for n = 0 , , , ... and stopped when Q ( τ Av n/ − ξ became negative. Finally, we passed a straight line y = mt + b through thepoints Q ( τ Av ( n − /
10) and Q ( τ Av n/ τ = ( ξ − b ) /m . The accuracy and speedgain relative to the GA are discussed in the results section.17 xample 2 In this example we choose a less ambitious partitioning of the system:group A group B ∅ −→ m k m k −→ m k + n k m k −→ ∅ n k −→ ∅ w k + s k −→ w k n k −→ s k w k + s k −→ w k s k −→ n k w k −→ w k + s k s k −→ ∅ w k −→ w k + s k . (57)The propensities for A are the same as in the previous example, but in addition we have a = a w s a = a w s a = b w a = b w a = a w s a = a w s a = b w a = b w a = a w s a = a w s a = b w a = b w . (58)The CME for B is the same as Eq. (27), while the MCME reads dQdt = H ( Q ) − (cid:34) (cid:88) i =1 ( r i δ w i , + d mi m i ) + (cid:88) i =1 3 (cid:88) j =1 ( a ji s i − + b ji ) w ji (cid:35) Q, (59)18here, for simplicity of notation, we made the definition H ( Q ) = (cid:88) k K k m k (cid:2) Q ( n k − − Q (cid:3) + d nk (cid:2) ( n k + 1) Q ( n k + 1) − n k Q (cid:3) + (cid:88) k c + k (cid:2) ( n k + 2)( n k + 1) Q ( n k + 2 , s k − − n k ( n k + 1) Q (cid:3) + (cid:88) k c − k (cid:2) ( s k + 1) Q ( n k − , s k + 1) − s k Q (cid:3) + (cid:88) k d sk (cid:2) ( n k + 1) Q ( n k + 1) − n k Q (cid:3) (60)Defining a new variable ˜ Q , such that Q ( n , s , t ) = ˜ Q ( n , s , t ) × exp (cid:40) − (cid:34) (cid:88) i =1 ( r i δ w i , + d mi m i ) + (cid:88) j =1 3 (cid:88) j =1 b ji w ji ) (cid:35) t (cid:41) × exp (cid:40) − (cid:88) j =1 3 (cid:88) j =1 (cid:90) t dt (cid:48) a ji (cid:104) s i − (cid:105) ( t (cid:48) ) (cid:41) , (61)where (cid:104) s i (cid:105) ( t ) = (cid:88) n , s s i P ( n , s , t ) , (62)and inserting (61) into Eq. (60), we obtain d ˜ Qdt = H ( ˜ Q ) − ∆( t ) ˜ Q, (63)where ∆( t ) = (cid:88) i =1 3 (cid:88) j =1 a ji ( s i − − (cid:104) s i − (cid:105) ) . (64)As before, we will treat ∆( t ) as small and write ˜ Q = ˜ Q (0) + ˜ Q (1) + ... where ˜ Q ( n ) ∼ O (∆( t ) n ),and insert it into Eq. (63) to obtain d ˜ Q (0) dt = H ( ˜ Q (0) ) d ˜ Q ( n ) dt = H ( ˜ Q ( n ) ) − ∆( t ) ˜ Q ( n − ∀ n > . (65)The first equation is identical to the CME in (27), and since the initial conditions mustsatisfy the relation ˜ Q (0) ( n , s ,
0) = P in ( n , s , Q (0) is the probability P ( n , s , t ). If we sumthe second equation for n = 1 over n and s , we obtain d ˜ Q (1) dt = 0 , (66)19hich implies that ˜ Q (1) ( n , s , t ) = 0. Hence, we can write Q ( t ) = (cid:88) n , s Q ( n , s , t )= (cid:104) Q (2) ( t ) + ... (cid:105) × exp (cid:40) − (cid:34) (cid:88) i =1 ( r i δ w i , + d mi m i ) + (cid:88) j =1 3 (cid:88) j =1 b ji w ji ) (cid:35) t (cid:41) × exp (cid:40) − (cid:88) j =1 3 (cid:88) j =1 (cid:90) t dt (cid:48) a ji (cid:104) s i − (cid:105) ( t (cid:48) ) (cid:41) . (67)For the purpose of testing, we set Q ( n> ( t ) = 0. As in the previous example, we expand theaverages (cid:104) n k (cid:105) and (cid:104) s k (cid:105) in powers of τ and keep only the lowest terms, which yields Q ( τ ) = e − Rτ , (68)where R = (cid:88) i =1 ( r i δ w i , + d mi m i ) + (cid:88) j =1 3 (cid:88) j =1 ( a ji v (0) k + b ji w ji ) . (69)In this case, computing τ from Q ( τ ) − ξ = 0 is particularly simple: τ = − R ln ξ . (70)To compute the probability for reaction µ to occur once we sample τ , we simply follow theprescription in Eq. (18): Q ( µ | τ ) = (cid:80) n , s ˜ Q (0) ( n , s , τ ) a µ ( n , s ) (cid:80) µ (cid:80) n , s ˜ Q (0) ( n , s , τ ) a µ ( n , s ) = a µ ( ¯n ( τ ) , ¯s ( τ )) (cid:80) µ a µ ( ¯n ( τ ) , ¯s ( τ )) , (71)where ¯n ( τ ) = (cid:104) n (cid:105) = (cid:80) q u ( q ) τ q and ¯s ( τ ) = (cid:104) s (cid:105) = (cid:80) q v ( q ) τ q . To update the probabilities, wefollow Eq. (19): Q ( n , s | µ, τ ) = ˜ Q (0) ( n , s , τ ) a µ ( n , s ) (cid:80) n , s ˜ Q (0) ( n , s , τ ) a µ ( n , s ) = ˜ Q (0) ( n , s , τ ) a µ ( n , s ) a µ ( (cid:104) n (cid:105) , (cid:104) s (cid:105) ) . (72)By virtue of the relation (28), the updated averages (cid:104) n (cid:105) and (cid:104) s (cid:105) become n = ˆ S µ (cid:80) n , s ˜ Q (0) ( n , s , τ ) n a µ ( n , s ) a µ ( (cid:104) n (cid:105) , (cid:104) s (cid:105) ) = ˆ S µ ns = ˆ S µ (cid:80) n , s ˜ Q (0) ( n , s , τ ) s a µ ( n , s ) a µ ( (cid:104) n (cid:105) , (cid:104) s (cid:105) ) = ˆ S µ s . (73)20e can approximate the effect of ˆ S µ by letting ¯n (0) ¯s (0) → ¯n ( τ ) ¯s ( τ ) + h µ , for those ¯ n k and ¯ s k that update to a value greater than or equal to zero; those that updateto a negative value are left unchanged. RESULTS
FIG. 2. Computing times for HSSAex1 and HSSAex2 relative to the computing time of GA(normalized to one), for different control parameters f and h . The computing times were averagedover 100 realizations running from 0 to 33 hours. IG. 3. Average and variance of m as a function of time and different control parameters f and h for the GA (blue), HSSAex1 (green) and HSSAex2 (red). The ensemble size was 500. IG. 4. A) Speed gain relative to the GA for HSSAex2 as a function of the control parameters f m and f λ and number of cells. B) Absolute computing time as a function of the control parameters f m and f λ and number of cells. C) Average and variance of m belonging to the first cell as afunction of time f m = f λ = 1. The ensemble size was 500. For A) and B), the simulations wereperformed up to 166 hours. We have simulated the system described in Fig. 1 with the GA, and the two variationsof the HSSA detailed in examples 1 and 2, which we label as HSSAex1 and HSSAex2. Theparameters were chosen so as to make the system oscillatory:23 , = a , = a , = a , = h × − min − a , = a , = h × × − min − b , = b , = b , = b , = b , = b , = h − min − r = r = r = 10 × f − min − K m = K m = K m = f min − d m = d m = d m = 5 × − min − d n = d n = d n = 1 × − min − d s = d s = d s = 5 × − min − c +1 = c +2 = c +3 = 5 × − min − c − = c − = c − = 5 × − min − . (74)The parameters h and f allow us to change the system without loosing oscillations. Fig. 2shows the efficiency of HSSAex1 and HSSAex2 relative to the GA, which has been normalizedto one for convenience. The green column represents the computing time of HSSAex1 relativeto the computing time of the GA; the red column gives the relative computing time ofHSSAex2. The HSSAex2 is clearly more efficient compared to HSSAex2 in all but the lastcase ( f = h = 10), but even in this case they are very close.Regarding accuracy, Fig. 3 shows the average and variance of the mRNA copy number ofthe first gene for the GA (blue), HSSAex1 (green) and HSSAex2 (red) for the 6 cases shownin Fig. 2. Although the efficiency between the HSSAex1 and HSSAex2 may differ greatly,their accuracy is equally good.Given the impressive speed gain of HSSAex2, we decided to test it on a 2-dimensionalarray of coupled identical cells. The coupling was introduced via a simple diffusion processof the homodimer of gene 3: s j → s j (cid:48) , where j and j (cid:48) label two adjacent cells. The diffusioncoefficient was chosen to be λγ λ , where λ = 10 − min − and γ λ is a dimensionless controlparameter. Another control parameter f m was introduced: K m = K m = K m = f f m . Theearlier parameters f and h were set to 1. Figs. 4 A and B show the speed gain of HSSAex2relative to the GA and the absolute computing time of HSSAex2, respectively. Fig. 4 Cshows the average and variance of the mRNA copy number of the first gene in the first cellfor the GA (blue) and HSSAex2 (red). 24 ISCUSSION AND CONCLUSION
We have presented exact derivation and practical applications of a hybrid stochastic simu-lation algorithm (HSSA) that can be, depending on system parameters, orders of magnitudefaster than the Gillespie algorithm (GA), and highly accurate. The principal behind theHSSA is the partitioning of a system of reactions into two groups, A and B ; the reactions in A are simulated using a Gillespie-type algorithm, while the reactions in B are described bythe chemical master equation (CME). We have derived exact formulas and equations whichallow, in principal, any reaction network to be partitioned in an arbitrary way, up to thecondition that there exists a subset of variables that is affected only by reactions in A . Forbiological systems such as gene regulatory networks (GRN), this condition is nearly alwayssatisfied. One way to violate it would be to have a fully connected reaction network in whichevery species of molecule interacts with every other species of molecules – which is rare atbest.Although the prescribed steps of the HSSA are straight forward, carrying them out mayrange from tractable, to difficult to impossible, depending on the specific partitioning of thesystem. To demonstrate how a reaction network can be partitioned, and how the HSSAmay be implemented in practice, we chose a GRN comprised of three interacting genesand gave two detailed examples of how this particular system lends itself to partitioning.In the first example, the reactions in A were merely the transcription and degradationof the mRNA of all three genes (6 reactions in total); while in the second example, wealso included the reactions that change promoter states of all three genes (18 reactions intotal). In carrying out the steps of the HSSA in both examples, we made the assumptionthat the reactions consisting of translation, forward and backward homodimerization, anddegradation of the monomers and homodimes for all three genes, lead to fluctuations thatare close to Poisson. This allowed us to set these fluctuations to zero, thereby reducing thecomplexity present in the steps of the HSSA. Consequently, the information about the copynumbers for those species that were described via the CME was reduced to their averages. Toobtain information about the fluctuations, one can write down equations for the statisticalmoments for these species which would lead only to a marginal loss of efficiency. This willbe demonstrated in a future work.The comparison in speed and accuracy presented in the “Results” section establishes the25SSA as an extremely useful tool for studying stochasticity in GRN, especially as imple-mented in example 2. Depending on the parameters, we found that the HSSAex2 was atleast 11 times faster than the GA (Fig. 2, h = 10, f = 1) and at most 96 times faster (Fig.2, h = 1, f = 10) for simulations of a single cell. For an array of cells, the HSSAex2 was upto 445 times faster (Fig. 4 A, f m = 5, γ λ = 10, in the case of 4 cells). The HSSAex1 didconsiderably worse compared to example 2. The best case scenario in single-cell simulationswas a speed gain factor of 14 (Fig. 2, h = 10, f = 10). Given the superior performance ofthe implementation in example 2, we did not simulate a multi-cell array using HSSAex1.It may seem counterintuitive that HSSAex1 be slower than HSSAex2, given that theformer has fewer reactions to simulate compared to the latter. However, when we considerthe set of tasks that HSSAex1 has to perform in B , it becomes understandable. In partic-ular, computing eigenvalues and eigenvectors in Eq. (53) and searching for the solution to Q ( t ) = ξ are computationally more expensive than evaluating the relatively simple alge-braic expressions in HSSAex2. However, as we saw for the last case in Fig. 2, the HSSAex1performed slightly better than HSSAex2 due to the faster promoter dynamics engendered byan increase in h . This emphasizes the dependence of the network topology and parameterson its partitioning.Out of many types of oligomers, the system we chose to work with contained only ho-modimers. Addition of heterodimers and higher ologomers is trivial in both HSSAex1 andHSSAex2: one only needs to modify Eqs. (32) to include the formation, dissociation anddegradation of these species. Although these reactions would couple the hitherto separateequations for the proteins and homodimers belonging to different genes, the polynomialexpansion on τ would still be applicable. For this type of system, the HSSA would likelyperform even better relative to the GA, given that the number of additional reactions sucha coupling would introduce is proportional to the square of the number of genes (for dimersand a fully connected protein interaction network). This would lead to a significant loss ofefficiency for the GA, but only a moderate one for the HSSA.In conclusion, we have derived exact expressions needed to partition an arbitrary reactionnetwork for the purpose of implementing an HSSA. We showed on a three-gene network howthe system may be partitioned. The two ways of partitioning lead to similar accuracybut significant overall difference in efficiency. The largest speed gain, compared to the GA,reached a factor of 445 for an array of 4 identical cells. Given that the reactions of our system26f choice are ubiquitous in systems biology, we believe that the methodologies advanced inthis paper will not only serve as preferred tools for discovering stochastic properties oflarge GRN, but will also open doors to further research in the technical aspects of systempartitioning. [1] Jahnke, T Wilhelm Huisinga W, (2007) Solving the chemical master equation for monomolec-ular reaction systems analytically J. Math. Biol. 54, 126[2] Pendar H, Platini T, Kulkarni RV, (2013) Exact protein distributions for stochastic modelsof gene expression using partitioning of Poisson processes Phys. Rev. E, 87, 042720[3] Shahrezaei V, Swain PS, (2008) Analytical distributions for stochastic gene expression PNAS,105 (45) 17256-17261[4] Walczak AM, Mugler A, Wiggins CH, (2012) Analytic methods for modeling stochastic regu-latory networks. Methods Mol Biol. 880, 273322[5] Bokes P, King JR, Wood ATA, Loose M, (2012) Exact and approximate distributions ofprotein and mRNA levels in the low-copy regime of gene expression J. Math. Biol. 64, 5,829854[6] Popovi´c N, Marr C, Swain PS (2016) A geometric analysis of fast-slow models for stochasticgene expression J. Math. Biol. 72, 12, 87122[7] Mugler A, Walczak AM, Wiggins CH, (2011) Spectral solutions to stochastic models of geneexpression with bursts and regulation. Phys. Rev. E, 80[8] Wolf V, Goel R, Mateescu M, Henzinger TA, (2010) Solving the chemical master equationusing sliding win- dows. BMC Systems Biology, 4, 42[9] Albert J, Rooman M, (2016) Probability distributions for multimeric systems. J. Math. Biol.72, 157169[10] Albert J, (2020) Dimensionality reduction via path integration for computing mRNA distri-butions arXiv:2006.08192[11] Albert J, (2019) Path integral approach to generating functions for multistep post-transcription and post-translation processes and arbitrary initial conditions Authors J. Math.Biol. 79(6-7): 2211-2236[12] Bokes P, King JR, Wood ATA, Loose M, (2012) Multiscale stochastic modelling of gene xpression J. Math. Biol. 65, 3, 493520[13] Veerman F, Marr C, Popovi´c N (2018) Time-dependent propagators for stochastic models ofgene expression: an analytical method J. Math. Biol. 77, 2, 261312[14] Munsky B, Khammash M, (2006) The finite state projection algorithm for the solution of thechemical master equation J. Chem. Phys. 124, 044104[15] Gupta A, Mikelson J, Khammash M, (2017) A finite state projection algorithm for the sta-tionary solution of the chemical master equation J. Chem. Phys. 147(15)[16] Gillespie DT, (1977) Exact Stochastic Simulation of Coupled Chemical Reactions. J. Phys.Chem. 81(25), 2340-2361[17] Gibson MA, Bruck J, (2000) Efficient Exact Stochastic Simulation of Chemical Systems withMany Species and Many Channels. J. Phys. Chem. 104(9), 18761889[18] Gillespie DT, (2001) Approximate accelerated stochastic simulation of chemically reactingsystems. J. Chem. Phys. 115(4), 1716[19] Cao Y, Li H, Petzold L, (2004) Efficient formulation of the stochastic simulation algorithmfor chemically reacting systems. J. Chem. Phys. 121, 4059[20] Cao Y, Gillespie DT, Petzold LR, (2005) Avoiding negative populations in explicit Poissontau-leaping. J. Chem. Phys. 123(5), 054104[21] Cao Y, Gillespie DT, Petzold LR, (2005) Efficient step size selection for the tau-leaping sim-ulation method. J. Chem. Phys. 124(4), 044109[22] Haseltine EL, Rawlings JB, (2002) Approximate simulation of coupled fast and slow reactionsfor stochastic chemical kinetics J. Chem. Phys. 117, 6959[23] Burrage K, Tian T, Burrage P, (2004) A multi-scaled approach for simulating chemical reactionsystems. Progress in Biophysics & Molecular Biology, 85, 217-234[24] Howard Salis H, Kaznessis Y, (2005) Accurate hybrid stochastic simulation of a system ofcoupled chemical or biochemical reactions J. Chem. Phys. 122, 054103[25] Jahnke T, Altntan D, (2010) Efficient simulation of discrete stochastic reaction systems witha splitting method. BIT Num Math 50(4), 797-822[26] Albert J, (2016) A hybrid of the chemical master equation and the Gillespie algorithm forefficient stochastic simulations of sub-networks. PloS one 11 (3), e0149909[27] Albert J, (2016) Stochastic simulation of reaction subnetworks: Exploiting synergy betweenthe chemical master equation and the Gillespie algorithm AIP Conference Proceedings 17901), 150026[28] Zechner C, Koeppl H, (2014) Uncoupled analysis of stochastic reaction networks in fluctuatingenvironments Plos Comp Biol, doi:10.1371/journal.pcbi.1003942.[29] Duso L, Zechner C, (2018) Selected-node stochastic simulation algorithm J. Chem. Phys, 148,164108[30] Kurasov P, L¨uck A, Mugnolo D, Wolf V, (2018) Stochastic Hybrid Models of Gene RegulatoryNetworks Mathematical Biosciences, 305, 170-177[31] Pucci F, Rooman M, (2018) Deciphering noise amplification and reduction in open chemicalreaction networks J. R. Soc. Interface 15: 20180805 https://doi.org/10.1098/rsif.2018.0805[32] Work in progress (paper will follow in near future)