Constrained Approximation of Effective Generators for Multiscale Stochastic Reaction Networks and Application to Conditioned Path Sampling

Efficient analysis and simulation of multiscale stochastic systems of chemical kinetics is an ongoing area for research, and is the source of many theoretical and computational challenges. In this paper, we present a significant improvement to the constrained approach, which is a method for computing effective dynamics of slowly changing quantities in these systems, but which does not rely on the quasi-steady-state assumption (QSSA). The QSSA can cause errors in the estimation of effective dynamics for systems where the difference in timescales between the"fast"and"slow"variables is not so pronounced. This new application of the constrained approach allows us to compute the effective generator of the slow variables, without the need for expensive stochastic simulations. This is achieved by finding the null space of the generator of the constrained system. For complex systems where this is not possible, or where the constrained subsystem is itself multiscale, the constrained approach can then be applied iteratively. This results in breaking the problem down into finding the solutions to many small eigenvalue problems, which can be efficiently solved using standard methods. Since this methodology does not rely on the quasi steady-state assumption, the effective dynamics that are approximated are highly accurate, and in the case of systems with only monomolecular reactions, are exact. We will demonstrate this with some numerics, and also use the effective generators to sample paths of the slow variables which are conditioned on their endpoints, a task which would be computationally intractable for the generator of the full system.


Introduction
[1] Define a dominating process to have transition rates given by the matrix M = 1 ρ G + I.
[2] This process has uniformly distributed reaction events on the time interval [t 0 , t 1 ]. The number r of such events is given by (1).
[3] Once r =r has been sampled, the type of each event must be decided, by sampling from the distribution (2), starting with the first event. An event which corresponds to rate m i,i indicates that no reaction event has occurred at this event.
[4] Once all event types have been sampled, we have formed a sample from the conditioned path space. We will briefly review the method presented in [16] for the exact sampling of conditioned paths in stochastic chemical networks. Suppose that we have a Markov jump process, possibly constructed from such a network, with a generator G. The generator of such a process is the operator G such that the master equation of the system can be expressed as where p is the (often infinite dimensional) vector of probabilities of being in 92 a particular state in the system. We wish to sample a path, conditioned on 93 X(t 0 ) = x 0 and X(t 1 ) = x 1 . Such a path can be found by creating a domi-94 nating process (i.e. a process whose rate is greater than the fastest rate of any 95 transitions of the original system) with a uniform rate. 96 We define the rate to be greater than the fastest rate of the process with generator G, so that ρ > max Then we define the transition operator of the dominant process by: We can then derive the number of reaction events N U of the dominating process  Table 2. Here the α i (X(t)) denote the propensity of the reaction R i when the 142 system is in state X(t) = [X 1 (t), X 2 (t), . . .], where ∆tα i (X(t) is the probability 143 that this reaction will fire in the infinitesimally small time interval (t, t + ∆t) 144 with 1 ∆t > 0. The stoichiometric vectors ν i denote the change in the state 145 vector X(t) due to reaction R i firing. 146 In order to describe the constrained approach, we first introduce some defi-147 nitions that will be helpful. tem. This is equivalent to applying the constrained stoichiometric projector P 164 to each of the stoichiometric vectors in the system. This may leave some reac- 165 tions with a null stoichiometric vector, and so these reactions can be removed 166 from the system. This projection can lead to aphysical systems where one or 167 more variables may become negative; in these cases we set the propensities of 168 the offending reactions at states where a move to a negative rate is possible, to 169 zero.

170
Let us illustrate this using an example which we shall be using later in the 171 paper. 172 −→ X 1 .
[1] Define a basis of the state space in terms of slow and fast variables.
[2] Initialise the value of the state, X(t 0 ) = x.
[3] Calculate propensity functions at the current state α i (X(t)).
[4] Sample the waiting time to the next reaction in the system [5] Choose one j ∈ {1, . . . , M }, with probability α j /α 0 , and perform reaction R j , with stoichiometry which has been projected using the constrained stoichiometric projector: [6] Repeat from step [3]. in this case, we choose F = X 2 . Following the projection of the stoichiometric 182 vectors using the constrained projector, the constrained system can be written 183 in the following way: analysed.

211
Suppose that we have a constrained system with N f fast variables, F 1 , F 2 , . . . , F N f . The generator for the constrained system with S = s is given by G F (s). Since the system is ergodic, there is a one-dimensional null space for this generator. This can be found by using standard methods for identifying eigenvectors, by searching for the eigenvector corresponding to the eigenvalue equal to zero. Krylov subspace methods allow us to find these eigenvectors with very few iterations. Suppose we have found such a vector v = [v 1 , v 2 , . . .], such that Then our approximation to the invariant distribution of this system is given by the discrete probability distribution represented by the vector Our aim is now to use this distribution to find the effective propensities of the 212 slow reactions of the original system.

213
Suppose that we have N s slow reactions in the original system. Each has an 214 associated propensity function α 1 (S, F ), α 2 (S, F ), . . . , α Ns (S, F ). We now sim-215 ply want to find the expectation of each of these propensity functions with 216 respect to the probability distribution p(s): [1] For each value of the slow variable S = s ∈ Ω, compute the generator G s of the constrained subsystem.
[4] Construct an effective generator G of the slow processes of the system using these effective propensities.  this error can be significant.

232
The CMA does not rely on the QSSA, and is able to take into account 233 the effect that the slow reactions have on the invariant distribution of the fast 234 variables, conditioned on a value of the slow variables. In a true fast-slow 235 system, this will yield the same results as the QSSA, but for most systems of 236 interest, the constrained approach will have a significant increase in accuracy.

237
If we follow the approach outlined in Table 3, we don't even need to conduct 238 any stochastic simulations to approximate the effective dynamics.

239
The assumptions for the CMA are weaker than the QSSA, namely that variables. This means that we have made the assumption that the current value 243 of the fast variables has no effect on the transition rates of the slow variables 244 once a slow reaction has occurred. This is subtly weaker than the QSSA, and  There will be many systems for which the constrained subsystem is itself a 252 highly complex and multiscale system. In this event, it will not be feasible to find 253 the null space of a sensibly truncated generator for the constrained subsystem.

254
Therefore, we need to consider how we might go about approximating this.

255
Fortunately, we already have the tools to do this, since we can iteratively apply

264
In this section we will present some analytical and numerical results produced We will consider this system in the following parameter setting: Here V denotes the volume of the well-mixed thermally-equilibrated reactor. The QSSA tells us that the fast subsystem (made up of reactions R 3 and R 4 ) 282 reaches probabilistic equilibrium on a timescale which is negligible in comparison 283 with the timescale on which the slow reactions are occurring. Therefore we may 284 treat this subsystem in isolation with fixed S: This is a very simple autocatalytic reaction system, for which a great deal of analytical results are available. For instance, we can compute the invariant distribution for this system [22], which gives us that X 2 is a binomial random variable Therefore, we can compute the conditional expectation E(X 2 |S) = k3S k3+k4 in this 286 fast subsystem, and use this to approximate the effective rate of reaction R 2 .

287
Therefore, the effective slow system is given by the reactions: where We can compute the invariant distribution for this effective system [22], which 290 in this instance is a Poisson distribution: We can quantify the error we have made in using the quasi-steady state as-292 sumption by, for example, comparing this distribution with the true invariant 293 distribution. Once again, using the results of [22], we can compute the invariant 294 distribution of the system (5), which is a multivariate Poisson distribution: , andλ 2 = k1V k2 . Trivially one can compute the marginal 296 distribution on the slow variable S: Therefore S is also a Poisson variable with intensity λ =λ 1 +λ 2 = k1V (k2+k3+k4) k2k3 , 298 which differs from the intensity approximated invariant density (10) by k1V k3 .
Note that k 3 is one of the fast rates, and k 1 V is one of the slow rates, and therefore as the difference in timescales of the fast and slow reactions increases, 301 this error decreases to zero, so that the QSSA gives us an asymptotically exact 302 approximation of the slow dynamics. For systems with a finite timescale gap, 303 the QSSA approximation will incur error over and above the error incurred in 304 any approximation of the marginalised slow process by a Markov process. For comparison, let us compute approximations of the effective slow rates 307 by using the CMA. The CMA for this system tells us that we need to analyse 308 the constrained system (6).

309
The constrained system in this example only contains monomolecular reactions, and as such can be analysed using the results of [22]. The invariant distribution for this system is a binomial, such that Using this, we can compute the effective propensity of reaction R 2 , giving us the effective ratek 2 = k2k3 k2+k3+k4 . The invariant distribution of (9) with this effective rate fork 2 is once again a Poisson distribution with intensity which is identical to the intensity of the true distribution on the slow vari-310 ables. In other words, for this example, the CMA produces an approximation 311 of the effective dynamics of the slow variables for this system, whose invari- approach gives us a system whose invariant distribution is exactly equal to the 317 marginal distribution on the slow variables for the full system. Another example 318 is presented in Section 5.3, for which the constrained system is itself too large to 319 easily compute expectations directly through its generator, and requires another 320 iteration of the CMA to be applied.

321
For this example, we did not even need to compute the invariant distri-322 butions of the constrained systems numerically. In Section 5.2, we will come 323 across a system for which it is necessary to numerically compute the invariant 324 distribution of the constrained system.  The approaches described in Section 1.1 hit problems when the system for 335 which we are trying to generate a conditioned path is multiscale. In a multiscale 336 system, the rate ρ of the dominating process will be very large, and as such 337 the number of reaction events will be large, even if the path we are trying to 338 sample is short. Therefore M r is likely to be a full matrix, and the number of 339 calculations of (2) will be large. Moreover, the size of M is also likely to be Suppose that we truncate the domain for this system to This truncated system has 201 2 = 40401 different states, and therefore the gen-354 erator G ∈ R 40401×40401 . Although this matrix is sparse, the matrix exponential 355 required in (1)  This means far fewer calculations of (2). What is more, as we saw in Section 366 5.1.2, for some systems the CMA exactly computes the effective dynamics of 367 the slow variables, with no errors.

368
The system (5), in order to highlight more effectively the differences between 369 the CMA and a QSSA-based approach, is in a parameter range (8), for which 370 the difference in time scales between the "fast" and "slow" variables is relatively 371 small, and of course for systems with larger timescale difference, the difference 372 in ρ between the full and effective generators would be far larger. = 44. Therefore, we are attempting to sample paths which start 390 and finish at the the mean of the invariant distribution, which in itself is not a 391 particularly interesting thing to do, but it will allow us to highlight again the 392 advantages of using the CMA over QSSA-based approaches.

393
Since the system is open, we are required to truncate the domain in order to be able to store and manipulate the effective generator. We truncate the domain to Ω = {[X 1 , X 2 ]|S = X 1 + X 2 ≤ 400}. Therefore we aim to sample paths   Let us consider the following chemical system, which in certain parameter 441 regimes exhibits bistable behaviour.
R 5 , R 6 : In particular, we consider parameter regimes where the occurrence of reactions R 5 and R 6 are on a relatively faster timescale than the other reactions. The following is just such a parameter regime: k 4 = 92.8, k 5 V = 10, k 6 = 500, k 7 = 6.
As before, V denotes the volume of the well-mixed thermally-equilibrated reac- is attributed to the fast reactions, R 5 and R 5 , given by: .
(13) This proportion, which is a measure of the gap in timescales between the "fast" reactions R 5 and R 6 , and the rest of the reactions, varies across the domain. We can approximate the expected proportion of propensity attributed to the fast reactions: where π Ω is the approximate invariant density of the full generator on the trun- R 5 , R 6 : Lines denoted by C i in this and what follows denotes a constraint. It is impor-474 tant to keep a track of these constraints, since each one reduces the dimension 475 of the state space by one.

476
For a fixed value of S = X 1 + 2X 2 ∈ {0, 1, . . . , S max }, we wish to find the 477 generator for the process X 2 (or equivalently X 1 = S − 2X 2 ) within this fast 478 subsystem. The generator can be found by considering the master equation for 479 each state X 2 = i: where p i (t) is the probability of X 2 (t) = i. Putting this set of differential equations into vector form gives us: where G is the generator of the fast subsystem (14). Note that since we are restricted to states such that X 1 + X 2 = S for some value of S, there are only the system for the state S as in (7), and in turn be entered into the (truncated) 488 effective generator for the slow variable. −→ X 1 + X 1 . This leads to 509 the following constrained system: 510 C 1 : Note that since reactions R 6 and R 7 have the same stoichiometry, this system 511 can be simplified by removing R 7 and adding its rate to R 6 : 512 C 1 :  being the altered rate for reaction R 6 . Following this methodology, an effective 515 generator G can be computed. have been chosen to demonstrate how far apart these two approximations can 532 be, but since the CMA has no additional costs associated with it, the advantages 533 of this approach are significant. The relative l 2 errors of these two approaches, 534 when compared with the approximate density π Ω , are given in Table 4, along 535 with the position and heights of the two local maxima in the densities.    In this section, we will illustrate how the nested approach outlined in Section 569 4 can be applied. We will consider an example for which we know the invariant 570 distribution of the slow variables. This gives us a way of quantifying any errors 571 that we incur by applying the nested CMA and QSSA approaches.
We will consider this system in the following parameter regime: As before, V denotes the volume of the well-mixed thermally-equilibrated reac- are occurring. This is demonstrated in Figure 6, where there is a clear gap in 576 the frequency of reactions R 1 and R 2 (the slowest), R 5 and R 6 (fast reactions) 577 and R 3 and R 4 (fastest reactions).

578
This system was chosen as we are able to, using the results in [22], find the exact invariant distribution of the slow variable S 1 = X 1 + X 2 + X 3 . In this instance, it is a Poisson distribution with mean parameter λ = k 1 V k 2 γκ (γk 2 + 3γκ + 2k 2 κ) = 64.2.
Total number of reactions Note that the quantity S 1 = X 1 + X 2 + X 3 is a conserved quantity with 585 respect to these reactions, and as such is the slow variable in this system. This 586 is in itself also a system with more than one timescale, and as such, we may 587 want to iterate again and apply a second QSSA assumption, based on the fact 588 that reactions R 3 and R 4 are fast reactions in comparison with reactions R 5 589 and R 6 . This leads to a second fast subsystem: 590 C 1 : Note that the quantity S 2 = X 1 +X 2 = S 1 −X 3 is a conserved quantity with 591 respect to these reactions, and as such is the slow variable in this system. At this 592 point in [11], the authors simulate the system using the Gillespie SSA. We could 593 adopt the approach that we used in Section 5.2, in which we find the invariant 594 distribution of the system by constructing its generator and then finding the 595 normalised eigenvector corresponding to the null space of that operator. This 596 would not be expensive since there are only S 2 different states. However, as 597 in Section 5.1, as this system only contains monomolecular reactions, we can 598 exactly find its invariant distribution. In this case, X 1 and X 2 follow a binomial 599 distribution with meanX 1 ,X 2 = S2 2 . This can then be used to compute the 600 effective rate of reaction R 5 in the first subsystem (18), α 5 (X 1 , X 2 ) ≈ γX 2 = 601 γ 2 S. This fast subsystem is then reduced to the following: 602 C 1 : Note that we have completely eliminated the fast variables X 1 and X 2 , and 603 instead consider the slower variable S 2 = X 1 + X 2 , with effective rate for R 5 604 given by the analysis above. This system is exactly solvable, and its invariant 605 distribution is a gamma distribution with means given byX 3 = S1 3 andS 2 = 606 2S1 3 , found by computing the steady states of the mean field ODEs [22]. This 607 in turn can be used to compute the effective rate of reaction R 2 in the full 608 system, where we now lose all of the fast variables X 1 , X 2 , X 3 and instead wish 609 to understand the dynamics of the slow variable S 1 = X 1 + X 2 + X 3 , which is 610 only altered by reactions R 1 and R 2 . This system is given by the following: Here the effective rate for R 2 has been found by using the approximation of the 612 effective rate α 2 (S 1 ) = k 2X3 = k2 3 S. We will now go through the same procedure, but this time using the con-615 strained subsystems instead of the fast subsystems as used in the previous sec-616 tion. There are 3 choices for the fast reactions, each involving two out of X 1 , 617 X 2 and X 3 . Since X 1 is the product of a zeroth order reaction, it is preferable 618 not to include this as one of the fast variables, and so we pick F 1 = [X 2 , X 3 ]. variables: Note that R 1 is removed, since it does not change the fast variables. R 2 is the 622 only other reaction which has changes to its stoichiometry due the constrained 623 stoichiometric projector. We have reduced the dimension of the system (due 624 to the constraint X 1 + X 2 + X 3 = σ for some σ ∈ N), but we are still left 625 with a multiscale system, which in theory could be computationally intractable 626 for us to find the invariant distribution for, through funding the null space of 627 its generator. Therefore, we can apply another iteration of the CMA to this 628 constrained system.

629
Reactions R 3 and R 4 are the fastest reactions in the system, and therefore 630 we pick our next slow variable that we wish to constrain to be S 2 = X 1 + X 2 , 631 which is invariant with respect to these reactions. Due to the previous constraint 632 S 1 = X 1 + X 2 + X 3 , we are only required to define one fast variable for this 633 system. Both choices F 2 = X 1 , X 2 , are essentially equivalent, and so we pick 634 F 2 = X 1 . These choices lead us to the following second constrained system: 635 C 1 : Here ν i denotes the stoichiometric vector associated with reaction R i , i.e. the 636 vector which is added to the state X(t) if reaction R i fires at time t. Notice 637 that we now have two separate constraints, and as such reactions R 5 and R 6 638 now have zero stoichiometric vectors. Moreover, these constraints lead us to 639 a somewhat unphysical reaction for R 2 . The reactant for this reaction is X 3 , 640 but only X 2 and X 1 are affected by this altered reaction. In system (19) when 641 reaction R 2 fires, we lose one X 3 , and gain X 1 . Therefore, both constraints 642 within (20) have been violated. In order to reset these constraints, without 643 changing the fast variable F = X 3 , we arrive at the stoichiometry presented This is a closed system, with a very limited number of different states. Therefore, it is computationally cheap to construct its generator, and to find that 648 generator's null space. Our aim with this system, is to find the invariant distri-649 bution of the fast variable given particular values for the constraints C 1 and C 2 .

650
This distribution will then allow us to compute the expectation of the reaction 651 R 4 within the constrained system (6), which is the only reaction which is depen-652 dent on the results of the second constrained system (since X 3 = S 1 −S 2 ). Once 653 the invariant distribution has been found, this can be used to find the effective 654 propensity of reaction R 5 given values of S 1 = X 1 + X 2 + X 3 and S 2 = X 1 + X 2 .

655
In turn, the constrained system (19)    to the invariant measure which was accurate up to 12 digits. 746 We showed how these effective generators can be used in the sampling of 747 paths conditioned on their endpoints. Such an approach could be employed as 748 a method to sample missing data within a Gibb's sampler when attempting to 749 find the structure of a network that was observed [16]. This approach could also 750 be used simply to simulate trajectories of the slow variables, in the same vein as

751
[6] or [11]. In this instance, it would only be necessary to compute the column of 752 the effective generator corresponding to the current value of the slow variables.

753
The constrained approach consistently significantly outperforms approxima-754 tions computed using the more standard QSSA-based approach, and at negligi-755 ble additional cost. Furthermore, in the limit of large separation of timescales, 756 the constrained approach asymptotically approaches the QSSA approximation.

757
The computational savings that we make in using the CMA depends on the 758 application with which we wish to use the effective generators. Similarly, if we 759 wish to approximate the invariant distribution of the slow variables, then the 760 CMA will always be less costly than exhaustive stochastic simulation. This is because we are able to directly compute the invariant distribution, whereas in 762 the simulation setting, to obtain the same statistics we would be required to 763 compute a very long simulations.

764
If, on the other hand, we simply wish to use the CMA to compute a tra-765 jectory of the slow variables, then the savings will vary, based on the size of 766 the chosen domain, and the relative differences in propensity of the fast and 767 slow reactions in the relevant regions. If our aim is only to produce one rel-768 atively short trajectory, then it is possible that stochastic simulation will be 769 more efficient than using the CMA. However this is such a trivial task, that any