Optimal explicit stabilized postprocessed $\tau$-leap method for the simulation of chemical kinetics

The simulation of chemical kinetics involving multiple scales constitutes a modeling challenge (from ordinary differential equations to Markov chain) and a computational challenge (multiple scales, large dynamical systems, time step restrictions). In this paper we propose a new discrete stochastic simulation algorithm: the postprocessed second kind stabilized orthogonal $\tau$-leap Runge-Kutta method (PSK-$\tau$-ROCK). In the context of chemical kinetics this method can be seen as a stabilization of Gillespie's explicit $\tau$-leap combined with a postprocessor. The stabilized procedure allows to simulate problems with multiple scales (stiff), while the postprocessing procedure allows to approximate the invariant measure (e.g. mean and variance) of ergodic stochastic dynamical systems. We prove stability and accuracy of the PSK-$\tau$-ROCK. Numerical experiments illustrate the high reliability and efficiency of the scheme when compared to other $\tau$-leap methods.


Introduction
The modeling of kinetic chemical processes involves multiple chemical species with different population's size and reacting time-scales.A typical ordinary differential equation (ODE) model for the simulation of such problems are the reaction rate equations (RRE), but this model is accurate only at the thermodynamic limit (i.e. when the populations size and system volume tend to infinity, but the concentrations remain constant).In contrast, for systems with small populations, as living cells, discrete and stochastic modeling is required to capture the correct kinetics.Assuming proper mixing and thermal equilibrium a discrete dynamical system in the form of a Markov process as well as its accompanying master equation, the chemical master equation (CME), can be derived for the evolution of the probability density function of a chemical system [21,35].The stochastic simulation algorithm (SSA) [19,20] gives an exact method to compute samples from the distribution of the CME.However, while very easy to implement, the SSA can become overwhelmingly slow due to the presence of multiple scales in the system (stiffness) and presence of large and small populations size, leading to reactions that fire extremely often.
By fixing a step size (or leap time) τ encompassing several reactions, Gillespie proposes the τleap method [22].This approximate procedure lumps together the reactions that would occur in a time lapse τ and fires them simultaneously.If the number of reactions fired in a time τ is large then the τ -leap scheme approximates the Euler-Maruyama method for the chemical Langevin equation (CLE) and in the thermodynamic limit this latter scheme approaches the explicit Euler method for the RRE [22,29].
In this paper we focus on τ -leap methods.As the reactions fire at disparate time-scales and the systems are typically stiff, the standard explicit τ -leap method [22] faces stability issues [16].Furthermore, even when the stability conditions are met, amplification properties due to explicitness of the scheme prevent to capture the correct statistics of the process.Implicit τ -leap schemes [14,16,25,39] in contrast usually do not have stability constraints.However, for ergodic dynamical systems, these schemes generally fail in capturing the exact statistics of the system.Hence, unless the fast processes are resolved, both explicit and implicit τ -leap methods fail to correctly integrate ergodic dynamical systems (see [33] for a similar discussion in the context of stochastic differential equations (SDEs) driven by diffusion processes).We mention further the trapezoidal τ -leap method, that is accurate in sampling the invariant measure for linear equations, but might fail for nonlinear problems [33].Very recently, in [40], a split step scheme generalizing the θ-method is introduced.The method is accurate in sampling the invariant measure of the process, however at each step it requires the solution of two nonlinear systems and an optimization problem for the scheme's coefficients.
Apart from the aforementioned implicit or explicit Runge-Kutta-like methods, several hybrid schemes making use of different models and levels of granularity exist in the literature.Such multirate (or multiscale) methods exploit the multiscale nature of chemical reaction systems, which often consist of multiple reactions firing at disparate time-scales.Roughly speaking, most of such schemes divide the reactions (or species) into fast and slow ones.Then, the fast dynamics are resolved by making use of a quasiequilibrium assumption and the slow terms are integrated employing larger step sizes -see [13,18,26,27,28,30,37].In this paper we do not assume that the system is clearly separable into fast and slow dynamics, therefore multirate methods are not discussed in what follows.
We now briefly describe explicit stabilized methods, that are the basis building blocks of our new scheme.In the ODE context, explicit stabilized methods are a compromise between explicit and implicit schemes.No linear algebra solutions are needed as for implicit methods, while quadratic growth (with the number of stages) of their stability domains allows for much better stability properties than classical explicit methods.Well known explicit stabilized methods are the Runge-Kutta-Chebyshev (RKC) methods [41,45,47], the DUMKA methods [31,32,36] and the Runge-Kutta orthogonal Chebyshev (ROCK) [1,7] methods (note that the first-order RKC and ROCK scheme coincide).More recently, the first-order RKC (or ROCK) scheme has been extended to SDEs, yielding the S-ROCK family [3,4,6,10] and higher order extensions in [8].For mean-square stable problems the S-ROCK scheme introduced in [6] represents an important improvement over the Euler-Maruyama method thanks to its improved stability properties (it does however not preserve the optimal stability domain of the first-order RKC method).However, for non mean-square stable problems and for problems with too large variance the efficiency of the S-ROCK scheme deteriorates.Starting from the same S-ROCK scheme as in [6], the authors in [5] derive the τ -ROCK method for equations driven by discrete noise; however, this method inherits the same issues as the S-ROCK method of [6].The SK-ROCK scheme [2] is an improvement over the previous S-ROCK method for SDEs, this scheme has an optimal stability domain's size growing quadratically with the number of stages and is second-order accurate in sampling the invariant measure of a class of ergodic SDEs.
The main contribution of this paper is the design and analysis of a new τ -leap method for stiff chemical systems.Inspired by the SK-ROCK scheme we propose here the PSK-τ -ROCK leap method for chemical kinetics.This method has several desirable properties: • it is fully explicit avoiding any linear algebra computations and is as easy to implement as the explicit τ -leap method; • it has an extended and optimal stability domain growing quadratically with the number of function evaluations, avoiding any step size restriction as the standard explicit τ -leap method; • thanks to a postprocessing technique adapted from [2,48] it shows remarkable properties in sampling correct statistics of non mean-square stable chemical systems, even when fast reactions are not resolved.
We analyze the accuracy and stability properties of the scheme and its long-time dynamics for ergodic dynamical systems.The efficiency and accuracy of the new scheme are illustrated through a sequence of numerical experiments, where we also compare the method against other τ -leap schemes as the implicit τ -leap and the trapezoidal τ -leap method [16,39].The rest of the paper is organized as follows.In Section 2 we give an introduction to the SSA and τ -leap method, in Section 3 we introduce the PSK-τ -ROCK method and provide a detailed pseudocode.The accuracy and stability analysis of the scheme is given in Section 4, while in Section 5 we provide the numerical examples.Conclusions are found in Section 6.

The SSA and the τ -leap method
In this section we briefly recall the modeling of a well stirred chemical reaction system at thermal equilibrium and introduce the SSA and τ -leap scheme.
A model for chemical reaction systems.Consider a chemical system composed by N species (of molecules) S 1 , . . ., S N which interact in M reactions, denoted R 1 , . . ., R M .We are interested in the number of molecules of each specie in an instant of time t.We denote by X(t) = (X 1 (t), . . ., X N (t)) ⊤ the state vector, where X j (t) ∈ N is the number of molecules of specie S j at time t.Each reaction R j is characterized by the propensity function a j (x) and the state-change vector ν j .Given a state x ∈ N N and an infinitesimal time dt, the quantity a j (x) dt is the probability that reaction R j fires within dt units of time.The state-change vector ν j describes the change in state x when reaction R j is fired, i.e. reaction R j has the effect of changing the state vector from x to x + ν j .We will denote by a(x) = (a 1 (x), . . ., a M (x)) ⊤ the vector of propensity functions and by ν = (ν 1 , . . ., ν M ) the stoichiometric matrix.
Example 2.1.We provide here an illustrative example of the above description.To do so, we consider the famous Michaelis-Menten system describing the mechanism of enzymatic catalysis.The model consists in four species: a substrate S 1 , an enzyme S 2 , a complex enzyme-substrate S 3 and the product S 4 .The three reactions may be written as (2.1) The state vector is X(t) = (X 1 (t), X 2 (t), X 3 (t), X 4 (t)) ⊤ and represents the number of molecules of each specie S 1 , S 2 , S 3 , S 4 at time t.If the first reaction fires, the value of X 3 (t) is increased by one molecule and X 1 (t), X 2 (t) are decreased by one molecule each, hence the state vector is updated as X(t) + ν 1 , where ν 1 = (−1, −1, 1, 0) ⊤ .In the same manner we define ν 2 = (1, 1, −1, 0) ⊤ and ν 3 = (0, 1, −1, 1) ⊤ .The propensity function a 1 (x) is the probability that the first reaction takes place within one unit of time and is given by a 1 (x) = c 1 x 1 x 2 , where the product x 1 x 2 is the number of possible distinct combinations of S 1 , S 2 molecules and c 1 is the probability that given two reactants S 1 , S 2 the reaction actually fires.Similarly, a 2 (x) = c 2 x 3 and a 3 (x) = c 3 x 3 .
There are also other types of reaction which are not listed in (2.1), for instance 2S n cj → S m and 3S n cj → S m , whose propensity functions are a j (x) = c j x n (x n − 1)/2! and a j (x) = c j x n (x n − 1)(x n − 2)/3!, respectively, and their structure follows from a combinatoric argument.
Given the vector of propensity functions a(x) and the stoichiometric matrix ν the system evolves following two simple rules [19].i) Given a state vector X(t) at time t, in an infinitesimal time dt the reactions R j are independent and the probability that reaction R j fires is given by a j (X(t)) dt.
ii) If R j fires the system is updated as X(t + dt) = X(t) + ν j .
The Stochastic Simulation Algorithm.From i), Gillespie [19] derived a probability density function from which we can sample a random pair (τ, j), where τ is the waiting time until the next reaction and j is the index of the next reaction.This is the core of the stochastic simulation algorithm (SSA), given by: 1) Sample the waiting time τ from an exponentially distributed random variable with rate a 0 (X(t)) = M j=1 a j (X(t)).2) Sample the next reaction j from an M point random variable, where index j has probability a j (X(t))/a 0 (X(t)).
3) Update the state vector as X(t + τ ) = X(t) + ν j and the time as t ← t + τ .
4) Return to 1), unless a stopping criteria is satisfied.
The most important property of the SSA is that it is exact in sampling the statistics of the system.However, if there is at least one reaction with high probability of firing, then a 0 (X(t)) will be large and the waiting time τ will likely be very short.Therefore, the SSA will use an excessively large number of time steps and become practically unreasonably expensive.
The τ -leap method.The τ -leap method [22] speeds the simulation by fixing a step size τ and firing all the reactions that occur within time τ simultaneously.This leaping strategy leads to a good approximation of the SSA if the so-called leap condition is satisfied: the propensity functions a j (x) must not change appreciably in the time interval [t, t + τ ].
First, suppose that in the time interval [t, t + τ ] the propensity functions a j (x) are constant and thus the reaction events are independent.Under this assumption, the number of times that reaction R j fires in the time interval [t, t + τ ] is described by a Poisson random variable with rate a j (x)τ , that we denote as P j (a j (x)τ ).Hence, under the leap condition, the τ -leap scheme is a good approximation to the SSA, where X n is an approximation of X(t n ) with t n = nτ .We note that in order to satisfy the leap condition the reactants population cannot be too small, otherwise a few reactions change considerably the number of reactants and thus the propensity functions change substantially as well.
Since the mean (and variance) of P j (a j (x)τ ) is a j (x)τ , it is useful to decompose the right-hand side of (2.2) in a drift term and a zero-mean noise term: where We note that (2.3) is very similar to the Euler-Maruyama scheme for SDEs, where the diffusion is replaced by the zero-mean discrete noise Q(x, τ ).
For stiff chemical systems, the approximation (2.3) can face severe step size τ restrictions to be stable [16].Using implicit time-stepping can cure stability issues at the expense of solving nonlinear problems.But implicit methods might fail to capture the correct statistics of a chemical system, in form of mean and variance, due to damping introduced by implicitness.

The PSK-τ -ROCK method
In this section we introduce the PSK-τ -ROCK scheme.This explicit stabilized τ -leap method is composed of: i) a time-marching scheme (denoted SK-τ -ROCK) for the computation of approximate solutions X n ; ii) a postprocessing procedure (denoted P) used to improve the accuracy of X n whenever needed, usually only at the very last time step.
In Section 3.1 we define the time-marching scheme SK-τ -ROCK while in Section 3.2 we motivate and introduce the postprocessing procedure P. The combination of the SK-τ -ROCK time-marching scheme with the postprocessor P yields the PSK-τ -ROCK scheme.In Section 3.3 we provide a detailed pseudocode for the PSK-τ -ROCK scheme and discuss some implementation details.
Considering a test problem, we will show in Section 4 that the PSK-τ -ROCK scheme has an optimal stability domain growing quadratically with the number of stages and thanks to the postprocessing procedure accurate sampling of the process' statistics is achieved.The properties shown on the test problem in Section 4 are verified numerically in Section 5 on more involved problems.

The SK-τ -ROCK step
Let τ be the step size, ε ≥ 0 be the damping parameter and β = 2 − 4ε/3; typically ε = 0.05.We denote by ρ the spectral radius of the Jacobian of f evaluated in X n , with f as in (2.4), and let the number of stages s ∈ N satisfy τ ρ ≤ βs 2 .The SK-τ -ROCK step, of size τ , is given by where f, Q are given in (2.4).The coefficients µ j , ν j , κ j , for j = 1, . . ., s, are as follows.We let where T s (x) is the Chebyshev polynomial of the first kind of degree s, defined by Finally, we define µ 1 = ω 1 /ω 0 , ν 1 = sω 1 /(2ω 0 )1 , κ 1 = sω 1 /ω 0 and, for j = 2, . . ., s, In (3.1),only one evaluation of the drift term f is required for accuracy, while the additional s− 1 evaluations are used to increase stability.Indeed, as we will see in Section 4, the SK-τ -ROCK step involves the first and second kind shifted Chebyshev polynomials, that are instrumental to obtain optimal stability domains.The parameter ε in (3.2) is called damping parameter.For ε = 0 the stability domain of the method (3.1) will have a finite number of points z i along the negative real axis for which the absolute value of the stability function is exactly one.This is avoided setting ε > 0. Also, introduction of damping is essential to study the ergodic properties of the numerical scheme (see Section 4 for details).
The main difference with respect to the previous τ -ROCK scheme [5] is that here the noise term is put at the beginning of the iteration and therefore it is stabilized by the drift.In the reversed τ -ROCK scheme, also introduced in [5], the noise is as well put at the beginning of the iteration but with different parameters ν 1 = 1 and κ 1 = 0, yielding in an overly damped noise.

The postprocessing procedure
In chemical reactions, one is often interested in the stationary state of a given system.Hence, an algorithm capable of capturing the invariant measure of the system is of considerable interest.Therefore, we propose here a postprocessing procedure for the SK-τ -ROCK time-marching scheme introduced in Section 3.1, which allows to considerably improve its accuracy when applied to non mean-square stable problems.We stress that the postprocessor is applied only when higher accuracy is required and it is not needed to advance the solution in time.However, before introducing the postprocessing procedure for chemical kinetics we briefly motivate it recalling the postprocessors' theory for linear SDEs.
Postprocessors for linear SDEs.Postprocessors are since long employed to increase the accuracy of numerical solutions to ODEs [11].However, a postprocessors framework for ergodic SDEs has been only recently proposed in [48].We recall here the ideas developed in [48] but restricted to the very particular case of the Ornstein-Uhlenbeck process.Consider the SDE where Applying a Runge-Kutta method to (3.4) yields X n+1 = A(z)X n + B(z) √ τ σξ n , where z = τ λ and Therefore, the numerical method has order r 1 for the invariant measure (i.e.| lim n→∞ Var( ) as z → 0. However, higher order is easily achieved applying a postprocessing procedure.Indeed, applying the postprocessor and therefore higher order The postprocessing procedure.Based on the ideas developed in [48] for the Ornstein-Uhlenbeck process, we define here the postprocessor for the SK-τ -ROCK time-marching scheme (3.1).Up to our knowledge, in the literature of chemical kinetics no such postprocessors have been used.In order to obtain higher accuracy for the invariant measure of the system at a certain time t n = nτ , the postprocessor with is employed.We stress that the PSK-τ -ROCK scheme does not need to compute (3.5) at each time step but only whenever higher accuracy is required.Due to the damping properties of the SK-τ -ROCK steps (see Theorem 4.3 below), the variance of the numerical solution X n in (3.1) is smaller than the exact variance.Adding the random variable αQ(X n , τ ) in (3.5) allows to increase the variance of the numerical solution, yielding in a better approximation.

The PSK-τ -ROCK method: the algorithm
The PSK-τ -ROCK method, thus, advances the solution in time using the SK-τ -ROCK scheme (3.1) to (3.3) and applies the cheap postprocessing step (3.5) and (3.6), whenever higher accuracy for the invariant measure is needed.We summarize in this section such method by providing a detailed pseudocode in Algorithm 1.
The input parameters of Algorithm 1 are the initial value X 0 , the step size τ , the end time T and the drift and compensated Poisson noise terms f (x) and Q(x, τ ), respectively, which are defined in (2.4).The output is the postprocessed numerical solution X N , which is an approximation to the exact solution X(T ), with T = N τ .The procedure for computing the method's coefficients at Line 6 of Algorithm 1 is given in Function Coefficients(s, ε) below.We conclude the section with a few comments on Algorithm 1.
• Approximation of the spectral radius at Line 4 is very cheap if performed with nonlinear power methods [34,46].In our experience those methods usually converge with at most two function evaluations, see Tables 2, 3 and 5.It is good practice to store the eigenvector associated to the largest eigenvalue and use it as starting guess for the next call to the nonlinear power method.
• We emphasize that Algorithm 1 has low memory requirement as it needs three stage vectors K −1,0,1 only, disregarding the size of s.Moreover, Line 10 has zero cost if performed by simply swapping memory addresses.
• It is common to replace s by s + 1 after Line 5, this enlarges the stability domain and ensures stability of the method even if the spectral radius ρ increases within one time step.
• The call to Coefficients(s, ε) at Line 6 is needed only if the number of stages s changes from one time step to the next.This does not happen too frequently.
• Algorithm 1 and Function Coefficients(s, ε) can be merged.Indeed the computation of coefficients µ j , ν j , κ j for j = 2, . . ., s can be done inside the for loop beginning at Line 9 of Algorithm 1, avoiding the execution of a for loop exclusively for the coefficients' definition.However, this does not improve significantly the performance unless the number of stages changes frequently.
• Finally, Line 1 in Function Coefficients(s, ε) has negligible cost if the values of ω 1 are precomputed and stored in table.

Accuracy and stability analysis
We analyze here the long-time accuracy and stability of the PSK-τ -ROCK scheme on a model problem: the reversible isomerization reaction which was first introduced in [16] and plays the role of test equation.This is a reversible system and therefore the number of molecules is constant, i.e.X 1 (t) + X 2 (t) = X T .As a consequence specie S 2 Function Coefficients(s, ε) Input : s and ε Output: Method's coefficients µ j , ν j , κ j for j = 1, . . ., s can be neglected, we consider only specie S 1 and denote X(t) = X 1 (t).The system is described by the two reactions Note that, for this particular model, λ = −(c 1 + c 2 ) is the only eigenvalue of the Jacobian of f , with f as in (2.4).Hence, the spectral radius of the Jacobian of f is ρ = |λ|.In what follows, λ represents the stiffness of the system.Problem (4.1) has a stationary state X ∞ with a binomial distribution B(n, p), where n = X T and p = c 2 /|λ| [23].Hence, Note that if c 2 = 0, i.e. (4.1) is not reversible, then E(X ∞ ) = Var(X ∞ ) = 0 and the problem is considered to be mean-square stable.
Definition 4.1.A numerical method applied to (4.1) is said to have absolutely stable mean and variance if, and only if, E(X n ) and Var(X n ) remain bounded as n → ∞.Moreover, if c 2 = 0 and lim n→∞ E(X n ) = 0 and lim n→∞ Var(X n ) = 0 then the method is said to be mean-square stable.
Preliminary results.In the foregoing analysis we will need the polynomials where U n (x) is the Chebyshev polynomial of the second kind of degree n defined recursively by and ω 0 , ω 1 are defined in (3.2); therefore, A s , B s depend implicitly on ε.In the next lemma we collect known results about the stability polynomial A s (z) [45].In particular, observe as the stability domain of A s (z) grows quadratically with s.
Proof.From (2.4) and (4.2), we deduce that for the test equation (4.1) it holds where Replacing f as in (4.6) into (3.1)yields and thus, considering the change of variables where z = τ λ.Using ν j + κ j = 1 for j = 2, . . ., s we also obtain It is shown recursively [2,47] that For the variance we use the law of total variation Var(X) = E(Var(X|Y )) + Var(E(X|Y )), where X, Y are two random variables, and we obtain    see for instance [9,42,43,48].Therefore, the postprocessing strategy for the SK-ROCK scheme cannot be extended straightaway to SDEs driven by Poisson noise, as those considered here.We therefore analyze in this section the PSK-τ -ROCK scheme.A full analysis for the linear problem (4.1) is provided, while we verify numerically in Section 5 that the postprocessing techniques proposed in Section 3.2 also successfully apply to nonlinear problems.
The next theorem provides an expression for the mean and variance of the postprocessed step (3.5), for a general parameter α. ) and c s (z) as in (4.5).
We discuss here our choice for α in (3.6), which aims at providing an amplifying factor cs (z) ≤ 1 as close to 1 as possible.From (4.11) we note that it is possible to obtain cs (z) = 1 for all z ∈ (−ℓ ε s , 0] only if c s (z) is an affine function of the form d s (z) = 1 + 2α2 z.From Figure 2(a) we deduce that this is untrue, unless ε = 0.However, even though c s (z) oscillates, we observe in Figure 2(a) that it often approaches the affine function d s (z) passing through the end points (0, c s (0)) and (−ℓ ε s , c s (−ℓ ε s )), i.e.
where we used c s (0) = 1, c s (−ℓ ε s ) = 0 2 and ℓ ε s = 2ω 0 /ω 1 .Therefore, we choose α so that The resulting amplifying factor cs (z) (4.11) for the postprocessed scheme PSK-τ -ROCK is displayed in Figure 3(a) for different values of ε and s = 5.In Figure 3(b) we display cs (z) as a function of s for fixed z = −200, −20, −2 and ε = 0.05.We see that the variance is not amplified with this definition of cs (z).Furthermore its damping remains bounded.This is in sharp contrast with the τ -ROCK or Rev-τ -ROCK methods (or some standard explicit or implicit τ -leap methods), where the variance is either amplified or strongly damped, respectively.We conclude this section showing that for ε = 0 if (4.10) holds then the PSK-τ -ROCK scheme captures the variance exactly, that is cs (z) = 1 for all z, as we see in Figure 3(a).We recall, however, that ε = 0 may lead to instabilities.Proof.For ε = 0 it holds ω 0 = 1, ω 1 = 1/s 2 and α = 1/(2s).Moreover, The identity

Numerical experiments
In this section we consider different numerical experiments in order to verify the stability and accuracy properties of the PSK-τ -ROCK scheme of Section 3 and also asses its efficiency compared to other τ -leap methods for stiff problems.To do so, we first evaluate its accuracy on a nonstiff bistable problem, where we investigate the effect of the postprocessing step (3.5) on nonlinear problems.Then, we consider a mean-square stable problem containing fast variables.Next, we tackle a nonlinear problem where the fast variables have very large variance and thus the equation is non mean-square stable.Finally we consider an application to a more involved problem, namely a genetic feedback loop.In Sections 5.2 to 5.4 we compare the efficiency of the new scheme against classical state of the art τ -leap methods for chemical kinetics.
Numerical methods and implementation details.Let us provide some general implementation details concerning the next experiments.In what follows, we often compare the PSK-τ -ROCK method (3.1) and (3.5) (also Algorithm 1 in Section 3.3) against other τ -leap schemes: the SK-τ -ROCK scheme defined by (3.1) but without postprocessing (3.5), the explicit stabilized τ -ROCK and Rev-τ -ROCK method [5], the implicit τ -leap (Imp-τ -leap) method and its postprocessed version (PImp-τ -leap) [39], the trapezoidal τ -leap (Trap-τ -leap) method [16] and the more recent split-step implicit τ -leap (SSI-τ -leap) method [25].All the aforementioned methods are implemented in C++ using the Eigen library [24] for the linear algebra routines.The recently introduced split step method [40] is not considered in the following experiments; nevertheless, notice that its cost is roughly twice the cost of the Imp-τ -leap method.
For the PSK-τ -ROCK and SK-τ -ROCK method we always use ε = 0.05 as damping parameter, i.e. the standard choice for explicit stabilized methods, and α as in (3.6).The number of stages s is chosen before each step according to the condition τ ρ ≤ βs 2 with β = 2 − 4ε/3, where the spectral radius ρ of the Jacobian of f (see (2.4)) is approximated employing a nonlinear power method [34,46].If relevant, we report the number of iterations of this nonlinear power method (PM).Here the condition τ ρ ≤ βs 2 guarantees stability and as ε is small βs 2 is a good approximation of the true stability domains' size ℓ ε s = 2ω 0 /ω 1 .For the τ -ROCK and Rev-τ -ROCK method the damping parameter ε = ε(s) depends on s and might be large, therefore βs 2 does not approximate the stability domain's size accurately.Hence, for these methods, the number of stages s is the smallest integer verifying τ ρ ≤ ℓ ε s , where ℓ ε s = 2ω 0 /ω 1 is the exact size of the stability domain of a scheme with s stages and damping ε.The damping parameter is chosen according to the strategy described in [6] for the τ -ROCK method and according to [5, eq. 4.32] for the Rev-τ -ROCK method.
Finally, for the postprocessing step of the PImp-τ -leap scheme we use ten steps of size δτ = 0.2/ρ of the Imp-τ -leap method (note that the relaxation time of the fast variables is proportional to 1/ρ).
We conclude this paragraph commenting on negative populations.As the Poisson random variables are unbounded, it is a common issue of τ -leap methods that negative populations arise when too many reactions fire during one step, compared to the available number of reactants.In the literature, numerous strategies have been developed in order to guarantee positive populations [12,17,38,44,49]; however, developing or adapting such a strategy for the SK-τ -ROCK methods is not the focus of this paper.Therefore, in the following numerical experiments we employ the usual trick of considering the absolute value of the components of x whenever Q(x, τ ) is evaluated.

Accuracy experiment on the nonstiff Schlögl model
In this first numerical experiment we want to investigate the effect of the postprocessing procedure (3.5) on the accuracy of the solution.To do so, we compare the PSK-τ -ROCK method with the non postprocessed SK-τ -ROCK method on the nonstiff nonlinear Schlögl model: where B 1 and B 2 are buffered species whose populations are kept constant at N 1 = 10 5 and N 2 = 2 • 10 5 molecules, respectively.The state vector X(t) represents the number of molecules of S, we set X(0) = 250 and the final time T = 50.The state-change vectors are ν 1 = ν 3 = 1, ν 2 = ν 4 = −1, the propensity functions are and we set c 1 = 3 • 10 −7 , c 2 = 10 −4 , c 3 = 10 −3 and c 4 = 3.5.For this choice of parameters the solution X has a bistable distribution.
In this experiment we fix τ = 0.5 and as the problem is nonstiff the PSK-τ -ROCK and SK-τ -ROCK methods are stable already for s = 1; however, note that even though s = 1 the PSK-τ -ROCK and SK-τ -ROCK methods do not correspond to the standard τ -leap scheme (2.2).Using 10 6 samples, we estimate the probability density function (pdf) of X(T ) approximated by the PSK-τ -ROCK or the SK-τ -ROCK method with s = 1.We display the results in Figure 4(a) and compare them against a reference pdf computed with the SSA.We observe that the PSK-τ -ROCK method matches very well the reference solution, while the SK-τ -ROCK scheme tends to cluster its solutions too close to the two stable points; both results are in line with the results of Section 4.  Now we want to investigate the accuracy of the schemes with respect to the number of stages s.
To do so, we define the density distance area (DDA) [15] where p(x) is the probability density function of X(T ) computed by the SSA and p(x) the one computed by the PSK-τ -ROCK or SK-τ -ROCK method.We display in Figure 4(b) the DDA of the two schemes for different stages s and fixed τ = 0.5, where the pdf p(x) is estimated using 10 6 samples.We note that the SK-τ -ROCK method becomes more accurate as s increases.In contrast, the PSK-τ -ROCK method is accurate already for low values of s.

Efficiency experiment on the Michaelis-Menten system
Here we consider the Michaelis-Menten system already described in Example 2.1, where we set X(0) = (3000, 120, 0, 0) ⊤ and the reaction rate constants as c 1 = 1.66 • 10 −3 , c 2 = 10 −4 and c 3 = 10 3 .With this set of coefficients, the variables X 1 , X 4 are slow and X 2 , X 3 are fast; however, this is a mean-square stable problem and therefore all variances tend to zero.A typical solution of the Michaelis-Menten system is displayed in Figure 5.For this model, the quantities of interest are the slow variables X 1 , X 4 ; therefore the byproducts X 2 , X 3 are neglected in what follows.
In this experiment, we want to compare the accuracy and efficiency of the PSK-τ -ROCK method against the other τ -leap methods listed at the beginning of Section 5. To do so, we fix τ = 0.25, integrate the equations with the different τ -leap methods and compare the expectations and standard deviations, computed over 10 6 samples, against a reference solution computed with the SSA.Also, we verify accuracy at the transient time T = 5 and at the equilibrium time T = 50.The results are reported in Table 1.We observe that all the schemes approximate relatively well the slow variables X 1 , X 4 , as the step size τ = 0.25 is small enough to resolve them and the fast variables have too small variance to perturb the accuracy of slow dynamics.Comparing the PSK-τ -ROCK and SK-τ -ROCK schemes, we observe that for this mean-square stable problem the postprocessing procedure (3.5) has no effect on the solution's accuracy.X 1 (5) X 4 (5) X Table 2 displays the computational times of the different methods, together with the average number of stages s and the average damping parameter ε.Moreover, for the explicit stabilized schemes we also display the average number of iterations, per time step, needed to approximate the spectral radius ρ with the nonlinear power method (PM) and for the implicit methods we display the average number of iterations needed by the Newton (N) method.The PSK-τ -ROCK and SK-τ -ROCK methods are the most efficient schemes.

Stiff nonlinear reversible reaction
In this experiment we consider the stiff nonlinear reversible reaction with c 1 = 50, c 2 = 10 3 and X(0) = (400, 3990) ⊤ .In this setting, (5.1) is at equilibrium and we can illustrate the accuracy of the schemes in capturing the invariant measure of the system, moreover the variances are large and thus the problem is not mean-square stable.As the quantity X C = X 1 + 2X 2 is constant over time we can eliminate X 2 from the system and let X = X 1 .The propensity functions and state-change vectors are We fix T = 0.2, τ = 0.01 and integrate the system with the τ -leap methods listed at the beginning of the section; the results are reported in Table 3, where we use 10 6 samples to approximate the statistics.For the τ -ROCK scheme we could not use the standard strategy to define the number of stages, as for this nonlinear problem with large variance we found that it does not guarantee stability.Searching for a set of parameters ε, s providing at least 1% of stable paths we found ε = 3500, s = 800, which is an unusable number of stages due to round-off errors [47].The same holds for Rev-τ -ROCK.We therefore did not include these methods in our numerical comparisons.All the other schemes were successful in 100% of the Monte-Carlo iterations.We observe in Table 3 that the SSI-τ -leap scheme is completely off.The other methods approximate well the expectation, while only PSK-τ -ROCK and PImp-τ -leap provide good approximations to the standard deviation, with PSK-τ -ROCK being the most accurate at smaller cost.For this very small problem, the Impτ -leap method is slightly faster than the PSK-τ -ROCK scheme; in contrast, it is significantly less precise.In [16] it is shown that the trapezoidal rule captures the exact invariant measure of ( Table 3. Nonlinear reversible reaction.For different methods, we report the mean and standard deviation of X, the per step average of number of stages, the average number of nonlinear power method (PM) or Newton (N) iterations and the total computational time.
as (5.1).The same phenomenon is observed in [33] for the trapezoidal rule for stochastic differential equations driven by diffusion processes.

Numerical experiment on a genetic positive feedback loop
We consider here a stiff biological system modeling a genetic positive feedback loop.The system describes the production of a protein, which auto regulates its production rate by binding to its gene promoter [38,49].This system is described by the set of reactions We consider the same methods as in the previous experiments and apply them to (5.2) using a fixed step size τ = 0.05 and computing 10 5 samples.We start by observing that the system contains the nonlinear reversible reaction (5.1), which induces significant fluctuations in X 1 , X 2 .The same reaction has been considered in Section 5.3, where we found that at equilibrium the variance of X 1 , X 2 is relatively large and the τ -ROCK and Rev-τ -ROCK methods need an excessively large number of stages and damping parameter in order to be stable.For (5.2) the same considerations hold and indeed we found that simulation of (5.2) with the τ -ROCK and Rev-τ -ROCK schemes is too often unstable; therefore, we will discard these methods from the rest of the experiment.For this example, the Trap-τ -leap and the SSI-τ -leap schemes had severe convergence issues in the Newton method, therefore their results are not reported neither.
In Table 4(a) we note that all methods are accurate in sampling the means.In Table 4(b) we observe that the PSK-τ -ROCK is the most precise in estimating the variance of all the species.The PImp-τ -leap scheme seems to be almost as precise as PSK-τ -ROCK, however we see in Table 5 that it is slower.Indeed, in Table 5 we see that PSK-τ -ROCK is not only the most accurate scheme but also the fastest one, together with SK-τ -ROCK.The speed-up is significantly larger than for the implicit methods.Table 5. Genetic positive feedback loop.For different methods, we report the per step average of the number of stages, the average number of nonlinear power method (PM) or Newton (N) iterations, the computational time and the speed-up compared to the SSA.

Conclusion
Based on stabilized methods, second kind Chebyshev polynomials and a postprocessing procedure, we have proposed a new explicit τ -leap method for discrete ergodic stochastic systems with multiple scales (Algorithm 1).Robustness (accuracy of the scheme) and extended stability domains growing quadratically with the stage number have been established (Theorem 4.3).Accurate approximation of the invariant measure of ergodic systems has been shown (Theorem 4.5) thanks to the cheap postprocessing procedure.Compared to other existing methods, the PSK-τ -ROCK method is shown to be: • faster, • more accurate, • easier to implement.
Numerical experiments confirmed the theoretical stability and accuracy results illustrating the efficiency of the PSK-τ -ROCK scheme when compared to other τ -leap schemes for stiff and ergodic discrete stochastic systems: in all considered cases the PSK-τ -ROCK method was the fastest and most accurate scheme.
Displaying the damping parameter ε, with respect to s, of the SK-τ -ROCK and τ -ROCK methods.Displaying the stability domain size ℓ ε s , with respect to s, of the SK-τ -ROCK and τ -ROCK methods.

Figure 1 .
Figure 1.Comparison between the damping parameter ε and the stability domain size ℓ ε s of the SK-τ -ROCK and the τ -ROCK method.For SK-τ -ROCK ε = 0.05, for τ -ROCK ε depends on s.

Figure 3 .
Figure 3. Illustration of cs(z) for different values of ε, s and z.
-ROCK, s = 1 SK-τ -ROCK, s = 1 SSA (a) Comparing the reference pdf computed by the SSA against the pdf computed by the two methods, with s = 1.Errors committed on the pdf by the two methods, as functions of the stages s.

Figure 4 .
Figure 4. Schlögl reaction.Approximation of the probability density function (pdf ) of X(T ) computed by the PSK-τ -ROCK and SK-τ -ROCK methods, with τ = 0.5 and different stages s.

Table 1 .
Michaelis-Menten.Empirical means and standard deviations of X 1 , X 4 at times T = 5 and T = 50.

Table 2 .
Michaelis-Menten.Per step average of number of stages, of damping parameter, of nonlinear power method (PM) or Newton (N) iterations and total computational times taken by the different methods.

Table 4 .
Genetic positive feedback loop.Empirical means and standard deviations.