In our algorithm, one method we use follows the partitioning scheme proposed by Cao et al. [13] by identifying the fast and slow processes based on reaction propensities, a_{j}(**x**, t). We identify a fast reaction subset
with m^{f} reactions, and a slow reaction subset
with m^{s} reactions. All species (reactants and products) participating in fast reactions are defined as *fast species*, and the remaining species as *slow species*. The state is given by **X**(t) = (**X**
^{f}(t),**X**
^{s}(t)), where
and
are the states of the fast and slow species, respectively. n^{f} and n^{s} are the number of fast and slow species, respectively, and n = n^{f}+n^{s} is the total number of species. The propensity function of the slow reaction
and of the fast reaction
at state **X**(t) = (**x**
^{
f
},**x**
^{
s
}) can be generally written as
, and
, respectively.

Partitioning is done according to a user-defined threshold of the propensities. The choice of this threshold is influenced by the required accuracy and the computational cost. Using simple model systems [17], it was found that a separation of time scales of at least 2 orders of magnitude is necessary for the MSMC method to be accurate and computationally more efficient; this value is further discussed below using new and more complex examples. This ranked propensity-based partitioning scheme does not always yield a single cut-off point, especially for large networks with multiple gaps in time scales. Furthermore, it does not identify the relaxation time. We propose the use of the computational singular perturbation (CSP) technique to assist with these difficulties.

##### Computational Singular Perturbation (CSP)-assisted stochastic partitioning

The computational singular perturbation (CSP) technique was introduced by Lam, Goussis and co-workers, e.g., [25, 26], as a diagnostic tool for gaining physical insight into the reaction dynamics of large complex networks, and as an accelerated deterministic solver for stiff networks. Our principal interest in the CSP approach is in its diagnostic power, specifically its ability to identify the species present in QE, and the dominating reactions in the fast and slow modes. We refer to our use of CSP techniques in network partitioning as a "CSP-assisted stochastic partitioning".

Consider the deterministic representation of the dynamics of a reaction network

**g(x)** is an n-dimensional vector, such that g

_{i}(t) = dX

_{i}/dt describes the rate of change of species S

_{i}. Differentiating (4), yields

where

is the Jacobian matrix of

**g**. In general, the jacobian matrix,

is not a diagonal matrix. We perform a linear transformation on

**g**,

such that the transformed mode

**f** evolves in a decoupled manner,

The matrix

consists of a set of n linearly independent row vectors,

**b**
_{i}, and is chosen appropriately to decouple the evolution of the transformed modes,

**f**.

is a diagonal matrix with the diagonal elements, |Λ

_{ii}|, arranged in the order of descending magnitudes. Every vector,

**b**
_{
i
}, is associated with a vector,

through the biorthonormality condition, such that,

The matrix
consists of n linearly independent column vectors,
, and is called the basis vectors set. The matrix
, which is the inverse of the basis vector set, is called the dual vector set. Lam defined an ideal basis vector set,
as one that completely decouples the evolution of the transformed modes [25, 26]. The CSP method provides an iterative refinement procedure [25, 27] to obtain the ideal basis vector set from a random trial basis set.

For linear systems, the Jacobian matrix is independent of time, and thus, the right hand eigenvectors of the Jacobian,
, can be conveniently used as the ideal basis vector,
(assuming
is a perfect matrix of linearly independent eigenvectors). In this case,
is the inverse of the eigenvector matrix of
. For nonlinear systems, the eigenvectors of
are time-dependent, and do not constitute an ideal basis vector set. Lam [25] proved that using the eigenvectors as an ideal basis set for nonlinear systems provides a leading order accuracy. Since we employ the CSP technique as a partitioning guide, leading order accuracy suffices as our examples indicate.

The elements Λ_{ii} represent the time scales of the system, with a possible gap in the magnitude of Λ_{ii} reflecting a time-scale separation. Assuming that such a gap can be identified, i.e., the system is stiff, we get a set of n^{f} fast modes, **f**
^{
f
}, and n^{s} slow modes, **f**
^{
s
}. If the separation is large enough, i.e.,
, then the fast modes relax rapidly, and the reduced system consisting of only the slow modes can be evolved over the low-dimensional manifold using larger time increments. While such a transformation in deterministic systems is helpful in removing stiffness from the system and reducing the system dimension, here we are more interested in extracting the physical meaning of the modes, and identifying reactions and/or species that contribute to the fast modes.

Lam [

25] defined the mode participation index (PI) of reaction j in mode i as

where
is given by the dot product **b**
_{i} ⊙ **υ**
_{j}. Δy* and Δt* are user-defined values of acceptable error in temporal solution provided by the CSP method. To simplify our analysis, we neglect the second term, O(|b_{i}·Δy*/Δt*|), that sets the threshold for magnitude of relaxed f. The participation index,
, for the fast modes (i = 1, 2, ..., n^{f}) gives the fractional contribution of reaction j to the fast mode, i.

Herein, the CSP method is used at any stochastic state to understand the instantaneous time scales and appropriately partition the network at that state. This constitutes our second partitioning method. The evolution of the reaction network is still carried out in a stochastic manner, using the HyMSMC method. The Jacobian at a stochastic state can be approximated using the corresponding deterministic model or its variant where the mean-field propensity functions are used instead of the deterministic reaction rates. We propose the latter approach. In addition, we recommend using an analytical Jacobian, if possible, due to its higher accuracy and lower computational cost. Aside from partitioning, the CSP method provides an estimate of the relaxation time of the system that could be used for converging the fast network. Since there are multiple eigenvalues, and thus multiple time scales, in our implementation we choose the smallest eigenvalue (in magnitude) from the n^{f} fast modes, to relax the fastest dynamics. We discuss this concept further below and in the results section.

Obviously, there is a computational overhead associated with the evaluation of the Jacobian and the eigenvalue analysis. This overhead depends on whether the system is linear (eigenvalue analysis is done only once) or not, the system size and the frequency of carrying out the CSP analysis. Ideally, network partitioning, either CSP-assisted or ranked propensity-based, should be performed prior to every call to the microscopic solver. Practically however, when the trajectory evolves fairly smoothly and gradually, a few calls of the CSP routine during a simulation are sufficient. For the networks considered in this work, the overhead of CSP is negligible in comparison to the cost of a stochastic simulation (see results section). Below, we demonstrate the potential of the CSP technique using several examples and discuss the CPU of the CSP analysis for the most complicated example.