# A hybrid multiscale Monte Carlo algorithm (HyMSMC) to cope with disparity in time scales and species populations in intracellular networks

- Asawari Samant
^{1}, - Babatunde A Ogunnaike
^{1}and - Dionisios G Vlachos
^{1}Email author

**8**:175

**DOI: **10.1186/1471-2105-8-175

© Samant et al; licensee BioMed Central Ltd. 2007

**Received: **13 February 2007

**Accepted: **24 May 2007

**Published: **24 May 2007

## Abstract

### Background

The fundamental role that intrinsic stochasticity plays in cellular functions has been shown via numerous computational and experimental studies. In the face of such evidence, it is important that intracellular networks are simulated with stochastic algorithms that can capture molecular fluctuations. However, separation of time scales and disparity in species population, two common features of intracellular networks, make stochastic simulation of such networks computationally prohibitive. While recent work has addressed each of these challenges separately, a generic algorithm that can *simultaneously* tackle disparity in time scales *and* population scales in stochastic systems is currently lacking. In this paper, we propose the hybrid, multiscale Monte Carlo (HyMSMC) method that fills in this void.

### Results

The proposed HyMSMC method blends stochastic singular perturbation concepts, to deal with potential stiffness, with a hybrid of exact and coarse-grained stochastic algorithms, to cope with separation in population sizes. In addition, we introduce the computational singular perturbation (CSP) method as a means of systematically partitioning fast and slow networks and computing relaxation times for convergence. We also propose a new criteria of convergence of fast networks to stochastic low-dimensional manifolds, which further accelerates the algorithm.

### Conclusion

We use several prototype and biological examples, including a gene expression model displaying bistability, to demonstrate the efficiency, accuracy and applicability of the HyMSMC method. Bistable models serve as stringent tests for the success of multiscale MC methods and illustrate limitations of some literature methods.

## Background

Stochastic behavior, resulting from the small population of species in an intracellular medium, can stimulate interesting molecular phenomena, such as noise-induced bifurcations, phenotypic heterogeneity, etc. Through such phenomena, intracellular noise can have an impact on the physiology of the cell [1, 2]. It is therefore important that *in silico* analysis of intracellular networks, aimed at relating cell function to molecular events, is capable of capturing fluctuations. Exact stochastic algorithms, such as the stochastic simulation algorithm (SSA) [3, 4] and its variants [5], can capture the molecular noise by accounting for the random nature of biochemical processes.

The SSA [3, 4] provides the exact solution but can be computationally intensive, especially when large populations and/or reaction networks are involved. τ-leap methods [6–9], which accelerate the simulation by firing multiple reactions in a coarse time step, provide an excellent alternative when large populations are involved. However, the original Poisson-based τ-leap method [6] and its modifications, the binomial τ-leap [7, 8] and the implicit τ-leap [10] methods, do not perform too well at low populations. Therefore, an intracellular network comprised of mixed population scales is not amenable to a τ-leap treatment. A modified τ-leap algorithm [11] was recently developed to avoid negative populations that can occur in the Poisson τ-leap method. Implicitly, this modification results in a hybrid, stochastic algorithm that handles mixed population levels by seamlessly switching between the SSA for low population species and the Poisson τ-leap method for large population species [11]. Similarly, the partitioned leaping method [12] combines the exact next reaction method [5] with the Poisson τ-leap [6] method to enable stochastic simulations over a disparate range of populations.

Another aspect of biological networks plaguing stochastic simulation is the presence of processes with disparate reaction rates. Fast reactions, which are sampled more frequently by the SSA, reduce the size of the time step. As a result, it is difficult to simulate macroscopic real-time. τ-leap methods are also inefficient when a large separation of time scales is encountered. Most multiscale, stochastic algorithms [13–21] accelerate simulation of stiff networks by reducing the time spent in simulating the fast network. One way to achieve this is to use an approximate, accelerated algorithm, such as the Langevin method [22], to speed up the simulation of fast reactions [15, 19]. This approach tacitly assumes that the faster kinetics arise purely from the large populations, and is not applicable to cases involving fast reactions with small populations. An alternative multiscale approach [13, 16–18, 20] based on the slow-scale SSA (SS-SSA) [13, 14] concept, uses information from the quasi-equilibrium (QE) description of the fast network to evolve the slow network (see Appendix A). In this multiscale approach, individual networks are modeled with the SSA, and thus, large populations in either the fast or the slow network cannot be handled effectively.

*a priori*assumptions about the scales and collapses to the exact SSA and/or non-multiscale stochastic solver, when a scale separation is absent. Additional new elements in our work include: (1) the introduction of the

*computational singular perturbation (CSP)*framework, as an auxiliary tool, to aid network partitioning and determine relaxation times of the fast network for its convergence and (2) of a

*new statistical relaxation criterion*that eliminates the need for the complete description of the QE probability distribution function (PDF), and (3) the applicability of the HyMSMC method to networks displaying bistability in the fast or the slow networks. To the best of our knowledge, this is the first successful application of a multiscale, stochastic algorithm to a bistable network, reported in the literature.

The paper is organized as follows. We begin with a brief introduction of the SSA notation, and then provide a detailed explanation of the CSP-assisted partitioning technique and the new statistical relaxation criterion. Next, we introduce the new hybrid multiscale Monte Carlo (HyMSMC) algorithm. Finally, we demonstrate the efficiency, accuracy and generality of the new algorithm with a few prototype and real biological examples.

## Results and Discussion

### Stochastic simulation algorithm (SSA) Notation

The notation we follow is identical to that of Gillespie [3, 4]. Consider a well-mixed, isothermal system of n species S≡ {S_{1}, S_{2}, ..., S_{n}} reacting through m reactions R≡ {R_{1}, R_{2}, ..., R_{m}}. Let the state of the system at any time t be denoted by a n-dimensional vector, **X**(t) = (X_{1}(t), X_{2}(t), ..., X_{n}(t)), where X_{i}(t) is the number of molecules of species S_{i} at time t. Let the n-dimensional vector **υ**_{
j
}correspond to the stoichiometry vector of reaction **R**_{j}, such that υ_{ij} is the stoichiometric coefficient of species S_{i} in reaction R_{j}. Given that the system is in a state **X**(t) = **x** at time t, we define a propensity function, a_{j}(**x**, t) such that a_{j}(**x**,t)dt gives the transition probability of the j^{th} reaction, Rj, occurring in an infinitesimally small time interval (t, t+dt). This function is typically dependent on the state of the species and the reaction conditions of the network, such as temperature and pressure.

### Stochastic Quasi-Equilibrium, Network Partitioning, Relaxation, and Evolution

#### a. Quasi-equilibrium (QE) in stochastic systems

where $\widehat{P}$(**x**^{f}, t'|**x**, t) is the probability of observing ${\widehat{X}}^{f}$(t') = **x**^{f}, given that the system is in a state **x** at time t. ${\widehat{X}}^{f}$(t) is a stochastic process that is identical to **X**^{f}(t), but describes the evolution of the fast species purely via the fast network.

_{∞}(

**x**

^{ f }|

**x**), at t'→∞ (asymptotic limit); in practice, this happens over a relaxation time, ${\tau}_{Rel}^{f}$, of the fast network. Accuracy of the QE approximation requires that the relaxation time, ${\tau}_{Rel}^{f}$, be much smaller than the smallest time scale, $\underset{j\in {R}^{s}}{\mathrm{min}}({\tau}_{j}^{s})$, of the slow network

Obtaining this probabilistic description and then using it in the evolution of the slow network forms the core of multiscale, stochastic MC methods [13, 16, 17, 20] (see also Appendix A). The solver used to compute it is called microscopic solver. A difficulty is that the relaxation time, ${\tau}_{Rel}^{f}$, is unknown in complex systems. We have indirectly addressed this issue in Refs. [17, 24] and revisit it below.

#### b. Partitioning of a stiff stochastic system

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 ${R}^{f}\equiv \left\{{R}_{1}^{f},{R}_{2}^{f},\mathrm{...},{R}_{{m}^{f}}^{f}\right\}$ with m^{f} reactions, and a slow reaction subset ${R}^{s}\equiv \left\{{R}_{1}^{s},{R}_{2}^{s},\mathrm{...},{R}_{{m}^{s}}^{s}\right\}$ 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 ${X}^{f}(t)=\left({X}_{1}^{f},{X}_{2}^{f},\mathrm{..},{X}_{{n}^{f}}^{f}\right)$ and ${X}^{s}(t)=\left({X}_{1}^{s},{X}_{2}^{s},\mathrm{..},{X}_{{n}^{s}}^{s}\right)$ 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 ${R}_{j}^{s}$ and of the fast reaction ${R}_{j}^{f}$ at state **X**(t) = (**x**^{
f
},**x**^{
s
}) can be generally written as ${a}_{j}^{s}\left(x;t\right)={a}_{j}^{s}\left({x}^{f},{x}^{s};t\right)$, and ${a}_{j}^{f}\left(x;t\right)={a}_{j}^{f}\left({x}^{f};t\right)$, 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".

**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

**g**. In general, the jacobian matrix, $\underset{\xaf}{\underset{\xaf}{J}}$ is not a diagonal matrix. We perform a linear transformation on

**g**,

**f**evolves in a decoupled manner,

**b**

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

**f**. $\underset{\xaf}{\underset{\xaf}{\Lambda}}$ 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, ${b}_{i}^{*}$ through the biorthonormality condition, such that,

The matrix $\underset{\xaf}{\underset{\xaf}{{b}^{*}}}$ consists of n linearly independent column vectors, ${b}_{i}^{*}$, and is called the basis vectors set. The matrix $\underset{\xaf}{\underset{\xaf}{b}}$, which is the inverse of the basis vector set, is called the dual vector set. Lam defined an ideal basis vector set, $\underset{\xaf}{\underset{\xaf}{{b}^{*}}}$ 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, $\underset{\xaf}{\underset{\xaf}{J}}$, can be conveniently used as the ideal basis vector, $\underset{\xaf}{\underset{\xaf}{{b}^{*}}}$(assuming $\underset{\xaf}{\underset{\xaf}{J}}$ is a perfect matrix of linearly independent eigenvectors). In this case, $\underset{\xaf}{\underset{\xaf}{b}}$ is the inverse of the eigenvector matrix of $\underset{\xaf}{\underset{\xaf}{J}}$. For nonlinear systems, the eigenvectors of $\underset{\xaf}{\underset{\xaf}{J}}$ 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., $\mathrm{min}\left(\left|{\Lambda}_{ii}^{f}\right|\right)\gg \mathrm{max}\left(\left|{\Lambda}_{ii}^{s}\right|\right)$, 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.

where ${B}_{j}^{i}$ 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, ${P}_{j}^{i}$, 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.

#### c. Evolution of the slow network

*slow-scale approximation*[13] evolves the slow network using the slow-scale propensities as the effective transition probabilities of the slow reactions. The slow-scale propensity of reaction ${R}_{j}^{s}$ is the expectation of the propensity function, ${a}_{j}^{s}$(

**x**

^{ f },

**x**

^{ s }), over the temporal QE PDF given by

where the summation is performed over all the equilibrium states **x**^{
f'
}visited by the fast species [13]. Slow reactions are picked and the elapsed time is computed the same way as in the SSA.

After a slow reaction has been picked, the populations of the species need to be updated. This is not as straightforward, and while critical, has been largely overlooked. A rigorous way of picking slow reactions and subsequently updating the slow network is based on a joint PDF [17] (see Appendix A). In this work, we use the last state present after the relaxation criterion (see section d) has been met. In doing this, we minimize the computational effort and also ensure that the fast state belongs to stationary PDF P_{∞}(**x**^{
f
}| **x**). In some cases, the selected slow reaction cannot be executed from some statesof the fast species, e.g., a reactant species with zero population, even though these states might be visited frequently [17]. Care is taken that the chosen slow reaction can be executed from that state by checking the feasibility of firing the chosen slow reaction from the selected state. Failure to satisfy this condition calls for additional simulation of the relaxed fast network until the chosen slow reaction can occur from the current state of fast species or selection from a database of states from P_{∞}(**x**^{
f
}| **x**). This eliminates the problem of having negative populations.

#### d. Relaxation of the Fast Network

The next question is how does one gauge the relaxation of the fast network and decide the time, τ_{Eqm}, required to obtain an accurate approximation of the stationary PDF, P_{→}(**x**^{
f
}| **x**), or the slow-scale propensities, $\overline{{a}_{j}^{s}}$**(x)** ? E et al. [20] decide τ_{Eqm} via an implicit relation to the deviation of the numerical estimate of the slow-scale propensity from its true value. Salis and Kaznessis [16] relax the fast network for a user-defined number of MC events, $MC{E}_{Eqm}^{f}$. A predetermined choice of τ_{Eqm} or $MC{E}_{Eqm}^{f}$ overlooks the likelihood of the relaxation time varying with time, and may make the multiscale simulation inaccurate and potentially inefficient. Thus, we need a more generic relaxation criterion that decides the relaxation time (or number of MC events) of the fast network on-the-fly. This issue is addressed next.

_{∞}(x

^{ f }|

**x**). In this work, we test convergence of the slow-scale propensities rather than of the entire PDF. We propose the 2 sample t-test to check for convergence of the slow-scale propensities. The 2 sample t-test is a commonly used statistical tool to verify with a certain degree of confidence, if two independently drawn samples of a random variable come from the same probability distribution. In the current context, the state of the fast species (and also the propensity, ${a}_{j}^{s}$(

**x**

^{ f },

**x**

^{ s }), of the slow reaction) is the random variable, such that a sample of a size

*p*corresponds to

*p*events of the fast network, with the state at every fast event corresponding to one data point in the sample. The t-statistic for a 2 sample t-test, assuming unequal sample sizes, N

_{1}and N

_{2}, and unequal variances, is given by

Here ${\overline{{a}_{j}^{s}}|}_{1}$ and ${{\sigma}_{j}^{2}|}_{1}$ are, respectively, the mean and variance of the propensity ${a}_{j}^{s}$ evaluated from N_{1} = (w - 1)·MCE_{win} MC events in the first (w-1) simulation windows. ${\overline{{a}_{j}^{s}}|}_{2}$ and ${{\sigma}_{j}^{2}|}_{2}$ are the mean and the variance of the propensity ${a}_{j}^{s}$ from N_{2} = MCE_{Win} events in the w^{th} (last) window. We define a window as a short stochastic simulation of the fast network, consisting of MCE_{win} fast MC events. The slow-scale propensities are converged to a stationary value when

-t_{α/2,υ} <t_{
j
}, <t_{α/2,υ,} for all j∈R^{
s
},

_{α/2,υ}is the value of the t-statistic for υ degrees of freedom and a significance level of α/2. We use a significance level, α = 0.05, in all the simulations shown in the results section. υ is evaluated as

To simplify the implementation of the above test and to avoid the recurring evaluation of the t_{α/2,υ}, we assume that υ > 40, and hence fix the critical value at t_{α/2,υ} ≅ 1.96. Since t_{α/2,υ} increases with decreasing υ, using t_{α/2,υ} ≅ 1.96 in cases where υ < 40 imposes a more stringent criterion on the relaxation, and hence can be used safely.

The accuracy of the multiscale method depends on the convergence of the slow-scale propensities. For instance, if the QE PDF is bimodal, then it is crucial that both the branches be sufficiently sampled in evaluating the slow-scale propensities. Also, it is possible that by pure chance, the MCE_{win} fast events in the first relaxation window do not alter the populations of any fast species participating in the slow network. So the t-test might falsely indicate convergence of slow-scale propensities as a result of inadequate sampling. To avoid such a situation, we suggest that the window length accounts for the size of, and the separation of scales within, the fast network, to ensure that all the fast reactions are adequately sampled.

We propose that the relaxation time, computed via CSP, should be used in conjunction with the relaxation time estimated via the statistical one to ensure sufficient sampling of QE state space. In general, we find good agreement between the statistical and CSP relaxation methods as detailed below in several numerical examples. In general, one should relax the system for the longest of the two times estimated by the CSP and the statistical criterion. Convergence of bistable systems is interesting (and more complex) as discussed in the results section.

### The Hybrid Multiscale Monte Carlo (HyMSMC) Method

The new relaxation criterion discussed above enhances the efficiency of the original MSMC method. There are two other areas that can be exploited to further accelerate the method: one is to execute multiple firings at each microscopic time step and the other is to reduce the number of evaluations of the QE PDF. One could reuse the same QE PDF for a few consecutive coarse (macroscopic) time steps, assuming that the occurrence of the slow events does not significantly alter this PDF. This is reminiscent of the τ-leap (TL) methods [6–10], where the leap condition allows one to perform multiple firings of the reactions in a pre-selected leap time. The inherent assumption of the TL methods is that the propensities remain unchanged in the leap time, and hence one can approximate the number of firings in the reaction channels using a Poisson [6, 10] or a binomial random variable [7, 8]. While maintaining the same QE PDF for a few coarse (macroscopic) time steps can potentially result in substantial computational savings, we need a criterion to indicate how long one could maintain the same PDF without sacrificing accuracy. We address this question by using the framework of the hybrid SSA-TL solver [11] to systematically coarse grain the macroscopic solver in time. Since leaping the slow reactions significantly moves the fast network away from its QE, we also propose the hybrid SSA-TL method [11] as the microsolver to speedup the relaxation of the fast network at every coarse time step.

The τ-leap solver [6–8], which temporally coarse-grains the exact SSA, works well only in the limit of large population [7, 8, 11]. At low population of the reactant species, there is a risk of observing negative populations due to the unbounded nature of the Poisson random variable. The binomial τ-leap methods [7, 8] eliminate this problem, but are not as accurate and efficient at low populations. Recently, the hybrid SSA-TL method was introduced by Cao et al. [11], primarily, to reduce the possibility of negative population of the Poisson τ-leap solver. A similar hybrid solver, the partitioned leaping method [12], which combines the next reaction method [5] and the Poisson τ-leap method [6], was introduced by Harris and Clancy. We incorporate the hybrid SSA-TL solver [11] into the MSMC method as the macrosolver and the microsolver. This practically eliminates the likelihood of negative populations while using the explicit Poisson TL method. We call this method the hybrid multiscale Monte Carlo (HyMSMC) method.

The overall concept of HyMSMC is depicted in Figure 1. Integration of the hybrid solver with the multiscale framework begins with classification of the fast (and the slow) reactions into SSA reactions and TL reactions subset, ${R}_{SSA}^{f}$ and ${R}_{TL}^{f}$ and (${R}_{SSA}^{s}$ and ${R}_{TL}^{s}$), respectively. All reactions whose reactant(s) population is less than some critical population, X_{crit}, are defined as SSA reactions [11]. We have successfully used X_{crit} = 10 in all simulations to obtain accurate results. In principle, the cost and accuracy of the hybrid solver increases with X_{crit}, with the method reducing to SSA for X_{crit} → ∞, and to a Poisson TL for X_{crit} = 0. Generating Poisson random numbers is more expensive than generating uniform random numbers. As a result, there is a X_{crit} value below which the hybrid solver is more expensive than the SSA. This computational overhead of the hybrid solver should be considered in deciding the optimum X_{crit} that maximizes the efficiency of the solution.

#### a. Hybrid Solver as the Microscopic Solver

^{f}is chosen as

_{1}and ξ

_{2}are uniform random numbers, sampled from the unit interval (0,1]. The number of firings, ${k}_{j}^{f}$ of a TL reaction ${R}_{j}^{f}$ ∈ ${R}_{TL}^{f}$is sampled from a Poisson distribution with mean ${a}_{j}^{f}{\tau}^{f}$

If ${\tau}_{TL}^{f}<{\tau}_{SSA}^{f}$ then only the TL reactions are fired according to Eq. (18).

#### b. Hybrid Solver as the Macroscopic Solver

*τ*

^{s}

If ${\tau}_{SSA}^{s}>{\tau}_{TL}^{s}$, only the TL reactions are executed according to Eq. (23).

### Algorithm Implementation

- 1.
Initialize the system state X(t = 0) at time t = 0 and the simulation parameters, such as the critical population X

_{crit}, the number of fast MC events, MCE_{win}, for each relaxation window, the coarse graining factor, r, in the leap criterion, the partitioning criterion, and the significance level a of the statistical test used for the relaxation. - 2.
At the state

**X**(t) = (**x**^{ f },**x**^{ s }), evaluate the propensities of all the reactions in the

network. Partition the reaction network into a fast and a slow network, and the species into fast and slow species.

#### a. Microscopic Solver – Fast Network

- 3.
Perform MCE

_{win}MC events of the fast network using the hybrid SSA-TL solver. In each event, - a.
Identify the SSA and the TL subsets in the fast network, ${R}_{SSA}^{f}$ and ${R}_{TL}^{f}$, respectively, based on the populations of the participating reactants;

- b.
Evaluate the time increment τ

^{f}using (14)-(16); - c.
Execute ${k}_{j}^{f}$ reactions, ${R}_{j}^{f}$ ∈ R

^{ f }using (17) and (18); - d.
Update the state of the fast species, ${x}_{i}^{f}\leftarrow {x}_{i}^{f}+{\displaystyle \sum _{j\in {R}^{f}}{k}_{j}^{f}{\nu}_{ij}^{f}}fori\in {S}^{f}$

- 4.
Evaluate the slow-scale propensities $\overline{{a}_{j}^{s}}$(see Appendix, Eq. (4)).

- 5.
If the convergence criterion (12) is satisfied and the time is longer than the relaxation time predicted by CSP, go to 6. Else, go to 3 and perform another MCE

_{win}fast events using the hybrid solver.

#### b. Macroscopic Solver – Slow Network

- 6.
Evolve the slow network using the hybrid SSA-TL solver:

- a.
Identify the TL and the SSA reactions subsets ${R}_{TL}^{s}$ and ${R}_{TL}^{s}$ at the state (

**x**^{ f },**x**^{ s }); - b.
- c.
- d.
If any reaction ${R}_{j}^{s}$ ∈ R

^{ s }cannot be fired ${k}_{j}^{s}$ times, execute steps 3a-3d until state**X**^{ f }=**x**^{ f }allows ${k}_{j}^{s}$ executions of all reactions, ${R}_{j}^{s}$ ∈ R^{ s }or pick a state from a database that can be fired. Else go to 6e; - e.
Update the state of all species, ${x}_{i}\leftarrow {x}_{i}+{\displaystyle \sum _{i=1}^{n}{k}_{j}^{s}{\nu}_{ij}^{s}}fori\in S$.

- f.
Update the time t = t + τ

^{s}. - 7.
If the desired time is reached, stop. Else, go to step 2.

Step 6d is performed to avoid negative populations that may arise from using the stationary PDF to select the representative state of the fast network (see discussion in sub-section c above and Appendix A). Next we demonstrate the strength of the HyMSMC algorithm with the help of simple prototype examples. Real biochemical networks are also considered to demonstrate the generality of the approach to complex networks.

### Algorithm Testing

In each of the following examples, we demonstrate the efficiency and accuracy of the proposed HyMSMC algorithm. Simulations were performed on 2.40 GHz Intel^{®} Xeon(TM) processors. The relaxation of the fast network is assessed with a relaxation tolerance, t_{α/2,υ} = 1.96 (α = 0.05,υ > 40) and using windows of 25 MC events. The hybrid solver uses a critical population of X_{crit} = 10 and a leap condition tolerance of r = 0.05. In all the examples shown below, the method is accurate with this parameter set. Each numerical example uniquely serves to highlight other novelties of the new HyMSMC scheme.

#### a. Coupled Isomerization Reactions

The above network is one of the simplest networks conducive to QE treatment. The simplicity of this network allows us to clearly demonstrate the key features of the HyMSMC algorithm. It also clarifies the applicability and adaptability of the algorithm to different levels of time-scale and population-scale separation. The rate constants are chosen such that the first reaction pair is much faster than the second reaction pair, and thus it constitutes the fast network. A four orders of magnitude time scale separation exists between the partitioned networks. Based on our definition, species A and B are the fast species and C is the slow species.

_{A}(t = 0) = 1500, X

_{B}(t = 0) = X

_{c}(t = 0) = 0, we generate stochastic trajectories of network (Al) using the exact SSA, the MSMC method, and the HyMSMC method. As can be seen in Figure 3 and Figure 4, the multiscale methods correctly predict the temporal evolution of the network. The HyMSMC trajectories in Figure 3c and Figure 4c have been generated from states recorded after the occurrence of every slow event, i.e., at every macroscopic step, compared to every millionth MC events in the SSA (Figure 3a and Figure 4a) trajectories. Approximately, 1000 data points have been used to generate all the trajectories seen in Figure 3 and Figure 4. SSA takes around 8 minutes for a single stochastic run from t = 0 to t = 10000 s. The MSMC run is completed in half a minute, and the HyMSMC run, in a fraction of a second. The speedup over the exact SSA is 16 and 320 using the MSMC technique and the HyMSMC technique, respectively. Using number of MC events as a metric, the speedup of the MSMC and HyMSMC methods over the SSA is ~20 and 3600, respectively.

_{1}/(k

_{1}+ k

_{-1}) and N = (X

_{A}+ X

_{B}). The slow-scale propensity of the linear reaction, B → C, is analytically given by k2Np (see dashed line in Figure 6).

^{6}), and in each case we ran the simulation until a certain time and monitored the CPU time and the number of MC events required to complete the run. The parameters in 5 cases are presented in Table 1. The time scale separation is maintained at around 10

^{4}in all cases.

Rate constants and initial populations in system (Al) used to maintain the time scale separation at ~O(10^{4}) as the population X_{0} increases.

Initial Population X | Fast Network k | Slow Network k |
---|---|---|

5 | 10 | 1.0 |

50 | 10 | 10 |

500 | 10 | 10 |

5000 | 10 | 10 |

50000 | 1.0 | 10 |

In contrast to the MSMC method, the speedup achieved with the HyMSMC method improves with increasing population. At low populations, the HyMSMC speedup approaches that of the MSMC technique, because the hybrid solver reduces to the SSA at these population levels that are below the critical population limit. At the other extreme, as the population increases, the number of slow events required to reach a certain time approaches 1 (see Figure 7). As a result, the CPU time and the number of MC events required to reach this end time reduce to the requirements for one slow event. Thus, in the limit of high population, the speedup of the HyMSMC over the exact SSA tends to plateau off. By increasing the end time used for the speedup measurement, say from 10^{4} s to 10^{5} s, the leveling of the HyMSMC speedup curve occurs at a higher population (see Figure 7a). In theory, the speedup of the HyMSMC method has a power law dependence on the population scale.

_{1}/k

_{2}), given that the populations of all species are in the same range (~500). Figure 8 shows that the acceleration achieved with both multiscale schemes has a power law dependence on stiffness. The speedup of the HyMSMC method lies above that of the MSMC method, highlighting the additional speedup achieved through temporal coarse- graining. As the network stiffness reduces to below 10

^{3}, the speedup of the MSMC drops below unity, implying that the MSMC method is more expensive than the SSA, and hence should not be used. This is caused by the cost of the MSMC method associated with relaxing the fast network. However, even with this small time-scale separation in the system, the HyMSMC still gives a speedup of around 30 over the SSA due to temporal coarse-graining.

Rate constants versus network stiffness in system (A1).

Network Stiffness ~O(k | Slow Network k |
---|---|

10 | 1.0 |

10 | 10 |

10 | 10 |

10 | 10 |

10 | 10 |

**x**, the Jacobian is time-invariant, and thus its eigenvalues and eigenvectors do not evolve with time. The eigenvalues of the Jacobian matrix are Λ

_{11}≈ -200, Λ

_{22}≈ -0.015, and Λ

_{33}= 0. The relaxation time of the fast network approximated from the largest in magnitude eigenvalue, (~|5/Λ

_{11}|~2.5 × 10

^{-2}), is consistent with the time estimated from the statistical convergence criterion over the course of a simulation (~ 2.5 × 10

^{-2}). The corresponding eigenvectors are

_{B}builds up, the contribution of reaction 2 to the fast mode increases, and eventually reactions 1 and 2 contribute almost equally to the fast mode. Thus, the partitioning of the networks with the CSP analysis agrees with that based on the propensities.

#### b. A System with Bistability in the Fast Network

_{1}and X

_{2}are in excess and act as buffering species. Consequently, the population of species X

_{1}and X

_{2}can be combined with the rate constants k

_{1}and k

_{2}, to give lumped rate constants ${{k}^{\prime}}_{1}$ and ${{k}^{\prime}}_{2}$, respectively. The network can be written as

The reaction rates in (B3) are based on the discrete nature of the reacting species. Solving (B3) for X gives ${X}_{Eqm}^{1}$ = 5, ${X}_{Eqm}^{2}$ = 20 and ${X}_{Eqm}^{3}$ = 50 as the steady state solutions of (B2) with ${X}_{Eqm}^{1}$ and ${X}_{Eqm}^{3}$ being stable and ${X}_{Eqm}^{2}$ being unstable. The equilibrium solution of the above steady state rate equation, using a generic algebraic solver, such as the Newton's method, depends on the initial guess.

To assess the accuracy of the HyMSMC technique, we compared the PDFs at 1000 s from 10000 trajectories, with those of the SSA. As was mentioned earlier, by selecting the last state at convergence as the representative state at the macroscopic time step, we are implicitly selecting the fast state from the frequency PDF (see Appendix A) present at that macroscopic time. For an accurate histogram of the fast species at any time, their states need to be sampled from the temporal PDF (see Appendix A) at any state of the slow network. In most cases, sampling from the frequency PDF does not cause any noticeable errors in the histogram of the fast species. However, these sampling-induced discrepancies surface in circumstances involving low populations and/or multimodality. This network exemplifies such a situation. Even though the statistics of the slow species Y is correctly captured (open triangles in Figure 11b), there is a distinct mismatch with the SSA results for the fast species (see filled triangles in Figure 11a). It is possible to obtain the correct probabilistic information of the fast network at any time, provided the slow network has been evolved correctly. To demonstrate this, we performed additional sampling of the fast network to obtain the temporal PDF. We then compute the average of the temporal PDFs at 1000 s using 10000 trajectories to account for the joint probability of observing a slow state, **x**^{s}, and the relaxed PDF, PDF (**x**^{s}), of the fast species at that slow state. Using this simple procedure, the PDFs of both the species match with those obtained using the SSA (see Figure 11 a, open triangles). This corroborates the fact that QE information of the fast network is embedded in the evolution of the slow species, and can be extracted by sufficient time-average sampling when needed. Usually, we do not have an *a priori* knowledge of the differences between the temporal and the frequency PDFs. Therefore, we propose that only the temporal PDF should be used.

_{Eqm}, of the fast network are the same at every slow step. Assuming that the same initial guess is given to the solver at every macroscopic time step, the bistability of the fast species will not be identified by the SS-SSA. This error creeps into the slow species evolution too. This can be seen in Figure 11b (filled symbols), where using X

_{Eqm}= 5 or X

_{Eqm}= 50 to evaluate the slow-scale propensities in the SS-SSA results in a wrong PDF of the slow species. This example emphasizes the need for accurate transfer of information (stochastic closure) between scales, and the failure of the mean-field approximation of the slow-scale propensities for bistable systems.

_{22}| is ~10

^{-4}and |Λ

_{11}| fluctuates between negative and positive values (this is characteristic of sampling states from the unstable attractor as the system transitions from one branch to the other) but is of the order of 10

^{2}-10

^{3}. The temporal profile of the participation indices in the fast mode is presented in Figure 12. The combined participation of reactions 5 and 6 is negligible at all times. The individual contribution of reactions 1 to 4 fluctuates a lot, but their combined participation to the fast mode equals ~1, see inset in Figure 12, implying that the fast mode is comprised of the first four reactions. Thus, the partitioning of the two methods (propensity and CSP based) is consistent.

Assessment of convergence in this bistable system is more complicated. The relaxation time estimated via CSP provides the time needed to sample a branch. For bistable systems, aside from two relaxation times, there is also a time scale related to the (residence) time the system stays in a branch. In order to adequately sample both branches, one needs to simulate the system for times longer than all three time scales. The mean residence time of the fast network was estimated to be ~0.01 s on the higher branch (X ~50) and ~0.04 s on the lower branch (X ~5). The relaxation time of the fast network based on magnitude of the eigenvalues lies in the range O(10^{-2})-O(10^{-3}). Thus, in this case, due to the comparable magnitude of residence times and relaxation times, the relaxation time based on the largest in magnitude eigenvalue gives a reasonable estimate of the sampling time of the fast network. In our simulation, we simulated the fast network for an additional 0.1 s, which was decided as a multiple of the sum of the residence times on each of the bifurcation branches.

#### c. Gene Expression Model: Bistability in the Slow Network

Having established the benefits of the algorithms with simple networks, we focus on real biological networks. The positive feedback-regulated gene expression model proposed in the work of Kepler and Elston [1] mathematically proved that fluctuations between operator states are a potential source of stochasticity in the level of translated protein, even at high mean populations. The gene expression model used by Kepler and Elston is an excellent biological example of "multiscales plaguing stochastic simulation". Unlike the previous example where the bistability was present in the fast network, in this gene expression model, the bistability is rooted in the slow network, and propagates to the fast network.

_{samp}= 0.001 s) the fast network at 10 s to collect the temporal PDF, followed by the ensemble average eliminates this discrepancy (see inset in Figure 15b).

_{0}= 0 and θ

_{0}= 1) obtained from a single equilibrium run of the HyMSMC method and the SSA. Figure 16 shows that the HyMSMC method correctly captures the PDF of life time of both states.

_{11}| ~10

^{6}) at all times, implying the presence of one fast mode. The microsolver relaxes the fast network in ~10

^{-4}time units based on the statistical criterion, which is greater than the eigenvalue-based relaxation time of 5 × l0

^{-6}, implying that our statistical relaxation criterion is more conservative and dominates convergence in this example. The contribution of the reactions to the fast mode is shown in Figure 17. The fast mode almost completely consists of reactions 1 and 2 (see inset in Figure 17). However, the contribution of reaction 2 drops to zero when the D = 0 state is encountered. In such a situation, the CSP-identified fast network will consist of only reaction 1. Even the rank based partitioning faces similar problems in such a situation. One way to tackle the problem is to base the partitioning of the network on average propensity values evaluated from a few (depending on network size) MC events. These averaged propensities could effectively be used to perform the CSP analysis.

#### d. Heat Shock Response (HSR) Model

Reaction network of the heat shock response model.

RxnNo. | Reactions | Rate const | Reaction Type |
---|---|---|---|

1 | Protein → UnfProt | 0.2 | Unfolding |

2 | UnfProt + DnaJ → UnfProt:DnaJ | 0.0108 | Association |

3 | UnfProt:DnaJ → Protein + DnaJ | 0.2 | Dissociation/Folding |

4 | DNA(σ32) → mRNA(σ32) | 1.40E-03 | Transcription |

5 | mRNA(σ32) → σ32 | 0.07 | Translation |

6 | mRNA(σ32) → degrades | 1.40E-06 | Degradation |

7 | σ32 → σ32 | 0.7 | Conformational |

8 | σ32 | 1.3E-01 | Conformational |

9 | DNA(DnaJ) + σ32 | 4.88E-03 | Transcription combined with Translation |

10 | DNA(FtsH) + σ32 | 4.88E-03 | |

11 | DNA(GroEL) + σ32 | 6.29E-03 | |

12 | DnaJ → degrades | 6.40E-10 | Degradation |

13 | FtsH → degrades | 7.40E-11 | |

14 | GroEL → degrades | 1.80E-08 | |

15 | σ32 + DnaJ → σ32:DnaJ | 3.62E-04 | Association |

16 | σ32:DnaJ → σ32 + DnaJ | 4.40E-04 | Dissociation |

17 | FtsH + σ32:DnaJ → DnaJ + FtsH | 1.42E-06 |

^{32}and DnaJ were chosen from the fast and slow network, respectively, to demonstrate the accuracy of the algorithm at all scales. Using the HyMSMC method requires 6 s (7 × l0

^{5}MC events) to reach t = 300, compared to 1200 s (9 × 10

^{8}MC events) required by SSA, thus accelerating the simulation by a factor of 200 (1000 based on MC events). As a test of accuracy, the normalized histograms of states at t = 1 s from 1000 trajectories were generated. The results are in good agreement, as shown in Figure 20. The algorithm accurately captures the noise of species σ

^{32}despite the low population involved (see Figure 20b).

^{-3}) is in good agreement with the one based on the statistical convergence criterion (order of ~1.2 × l0

^{-3}). The temporal profile of the participation indices, shown in Figure 21, identifies the first 3 reactions in the network as the fast reactions that contribute the most to the fastest mode. The combined participation index of reactions 4 to 17 is negligible compared to the participation index of reactions 1 to 3. Thus, reactions 4–17 are the slow reactions. Figure 21 also shows that for the time considered, the network partitioning is maintained. Thus, even for a moderately large network, the CSP method offers a systematic method to partition the network. For the HSR model, the CSP analysis increased the CPU cost of the HyMSMC method by less than 4%. Thus, for the largest network simulated in this paper, the CPU overhead of the CSP-assisted partitioning method is negligible.

## Conclusion

In this paper, we have presented a hybrid, multiscale Monte Carlo (HyMSMC) algorithm that can *simultaneously* handle disparity of timescales and species population in well-mixed stochastic networks. The HyMSMC method uses a hybrid SSA-i-leap method [28] as the microscopic and the macroscopic solvers, and employs stochastic singular perturbation concepts when time scale separation exists. As a result, it is able to systematically exploit large populations in the fast and/or slow networks to accelerate stochastic simulations without losing accuracy. Importantly, the absence of *a priori* assumptions about the scales in the network makes the HyMSMC method applicable over a wide range of the population and time-scale separation. Consequently, the HyMSMC method can switch between the exact SSA [3, 4], the τ-leap [6], or the MSMC [17] method in a dynamic manner, depending on the instantaneous timescales and populations in the network.

The superiority of the HyMSMC method over the original MSMC method [17] and the SSA [3, 4] was demonstrated using various examples, including two networks displaying bistability. The core of most multiscale algorithms is the closures between scales, which strongly influence accuracy. Using the Schlögl model and a simple gene expression model, we showed that the closures used in the HyMSMC method accurately capture the bimodality in the fast and the slow networks.

A second contribution of this paper is the stochastic partitioning and convergence that have not extensively been discussed in previous work. We introduced the computational singular perturbation (CSP) technique [25, 26] as a diagnostic tool to systematically identify the fast and the slow networks. In all examples, CSP was able to correctly partition the network and provide estimates of relaxation time for converging the fast network. This latter use of eigenvalues in conjunction with statistical tests is strongly recommended in order to eliminate numerical artifacts of the latter. The computational overhead associated with the evaluation of the Jacobian matrix and the eigenvalue analysis has been small for the examples studied here.

## Appendix A: Evaluation of PDFs and Slow Scale Propensities

### Calculation of PDFs

*(microscopic solver)*until the stochastic low-dimensional manifold has been reached, then compute the QE PDF of the fast network, and subsequently use it to evolve the slow network. In our original MSMC method [17], this was accomplished using the SSA but this can also be done with the τ-leap method. We compute the time-averaged PDF of the fast states

where Δt^{f} is the total life time of the state **x**^{
f
}in an equilibrium simulation of the fast network for time ${\tau}_{Eqm}={\displaystyle \sum _{{x}^{{f}^{\prime}}}\Delta {t}^{{f}^{\prime}}}$, given that the slow system is at state **x**^{
s
}. We obtain Δt^{f} by monitoring the life-time of states accessed once at equilibrium.

Here, MCE^{f} is the number of times a state **x**^{
f
}is observed in an equilibrium simulation of MCE_{Eqm} MC events of the fast network. We refer to the PDFs based on the life-time of states (Eq. (1)) and frequency of states (Eq. (2)) as the temporal and frequency PDF, respectively. In most cases, the frequency and temporal PDFs are practically the same. However, differences arise when the QE PDFs are multimodal and/or when low populations are encountered, as was shown in Ref. [17]. Numerical examples of this paper further indicate that the frequency PDF method can result in erroneous results in such cases, so we strongly recommend the use of the temporal PDF.

### Calculation of Slow-Scale Propensities

The propensity function of the slow reaction ${R}_{j}^{s}$, given by ${a}_{j}^{s}$(**x**^{
f
},**x**^{
s
}), is a function of the states of the slow and the fast species. Thus, the distribution of **X**^{
f
}results in a distribution of ${a}_{j}^{s}$(**x**^{
f
},**x**^{
s
}) for all j ∈ ${R}_{j}^{s}$. The transition probabilities used to evolve the slow network should be projected onto the QE PDF of the fast network. In Ref. [17], we underscored that the most rigorous way of evolving the slow processes should be based on a joint PDF that accounts for the life time of various fast states along with the rate at which these states are being changed by slow processes.

Evaluating the first moment of the propensity, ${a}_{j}^{s}$(**x**^{
f
},**x**^{
s
}), over the temporal QE PDF, P_{∞}(**x**^{
f
}| **x**), is a simpler way to assimilate the stochastic QE description into the dynamics of the slow network [13]. The SS-SSA [13] uses an analytical expression of QE-PDF to evaluate Eq. (10), which restricts the applicability of the algorithm to simple networks. Alternatively, it uses a deterministic, steady-state solution of the fast network to evaluate the slow-scale propensities. As correctly pointed out by Cao et al. [13], this mean-field approach neglects the correlation in the stochastic QE and is accurate at large populations only. In a numerical example, we show that the mean-field approach also fails when the fast network has multiple steady states. Other stochastic multiscale methods [16, 17, 20], which employ the slow-scale approximation, obtain the slow-scale propensities via a numerical approximation, and are, therefore, more general than the SS-SSA, at the expense, of course, of increased computational cost.

where $MC{E}_{Eqm}^{f}$ is the number of MC events executed by the *microscopic solver* to simulate time (after relaxation has been achieved) ${\tau}_{Eqm}={\displaystyle \sum _{p=1}^{MC{E}_{Eqm}^{f}}{\tau}_{p}^{f}}$ of the quasi- equilibrated fast network. ${{\tau}^{f}|}_{p}$ is the time increment in the p^{th} MC event, and ${{a}_{j}^{s}|}_{p}$ is the propensity of the j^{th} slow reaction before the occurrence of the p^{th} MC event.

## Declarations

### Acknowledgements

The research was partially supported by the Department of Energy (DE-FG02-04ER25626 and DE-FG02-05ER25702). However, any opinions, findings, conclusions, or recommendations expressed herein are those of the authors and do not necessarily reflect the views of the DOE.

## Authors’ Affiliations

## References

- Kepler TB, Elston TC: Stochasticity in transcriptional regulation: Origins, consequences and mathematical representations. Biophys J. 2001, 81: 3116-3136.PubMed CentralView ArticlePubMed
- Elowitz MB, Levine AJ, Siggia ED, Swain PS: Stochastic gene expression in a single cell. Science. 2002, 297 (5584): 1183-1186. 10.1126/science.1070919.View ArticlePubMed
- Gillespie DT: A general method for numerically simulating the stochastic evolution of coupled chemical reactions. Journal of Computational Physics. 1976, 22: 403-434. 10.1016/0021-9991(76)90041-3.View Article
- Gillespie DT: Exact stochastic simulation of coupled chemical reactions. Journal of Physical Chemistry. 1977, 81: 2340-2361. 10.1021/j100540a008.View Article
- Gibson MA, Bruck J: Efficient exact stochastic simulation of chemical systems with many species and many channels. Journal of Physical Chemistry A. 2000, 104: 1876-1889. 10.1021/jp993732q.View Article
- Gillespie DT: Approximate accelerated stochastic simulation of chemically reacting systems. Journal of Chemical Physics. 2001, 115 (4): 1716-1733. 10.1063/1.1378322.View Article
- Chatterjee A, Vlachos DG, Katsoulakis MA: Binomial distribution based τ-leap accelerated stochastic simulation. Journal of Chemical Physics. 2005, 122 (2): 024112-10.1063/1.1833357.View ArticlePubMed
- Tian T, Burrage K: Binomial leap methods for simulating stochastic chemical kinetics. Journal of Chemical Physics. 2004, 121 (21): 10356-10364. 10.1063/1.1810475.View ArticlePubMed
- Auger A, Chatelain P, Koumoutsakos P: R-leaping: Accelerating the stochastic simulation algorithm by reaction leaps. Journal of Chemical Physics. 2006, 125 (8):
- Rathinam M, Petzold LR, Cao Y, Gillespie DT: Stiffness in stochastically reacting systems: The implicit τ-leaping method. Journal of Chemical Physics. 2003, 119 (24): 12784-12794. 10.1063/1.1627296.View Article
- Cao Y, Gillespie DT, Petzold LR: Avoiding negative populations in explicit Poisson tau-leaping. Journal of Chemical Physics. 2005, 123 (5): 054104-10.1063/1.1992473.View ArticlePubMed
- Harris LA, Clancy P: A "partitioned leaping" approach for the multiscale modeling of chemical reaction dynamics. Journal of Chemical Physics. 2006, 125: 144107-10.1063/1.2354085.View ArticlePubMed
- Cao Y, Gillespie DT, Petzold LR: The slow-scale stochastic simulation algorithm. Journal of Chemical Physics. 2005, 122 (1): 14116-10.1063/1.1824902.View ArticlePubMed
- Cao T, Gillespie DT, Petzold LR: Multiscale stochastic simulation algorithm with stochastic partial equilibrium assumption for chemically reacting system. Journal of Computational Physics. 2005, 206: 395-411. 10.1016/j.jcp.2004.12.014.View Article
- Sails H, Kaznessis Y: Accurate hybrid stochastic simulation of a system of coupled chemical or biochemical reactions. Journal of Chemical Physics. 2005, 122 (5): 54103-10.1063/1.1835951.View Article
- Sails H, Kaznessis Y: An equation-free probabilistic steady-state approximation: Dynamics application to the stochastic simulation of biochemical reaction networks. Journal of Chemical Physics. 2005, 123: 214106-10.1063/1.2131050.View Article
- Samant A, Vlachos DG: Overcoming stiffness from stochastic simulation stemming from partial equilibrium: a multiscale Monte Carlo algorithm. Journal of Chemical Physics. 2005, 123 (14): 144114-10.1063/1.2046628.View ArticlePubMed
- Rao CV, Arkin AP: Stochastic chemical kinetics and the quasi-steady-state assumption: Application to the Gillespie algorithm. Journal of Chemical Physics. 2003, 118 (11): 4999-5010. 10.1063/1.1545446.View Article
- Haseltine EL, Rawlings JB: Approximate simulation of coupled fast and slow reactions for stochastic chemical kinetics. Journal of Chemical Physics. 2002, 117 (15): 6959-6969. 10.1063/1.1505860.View Article
- E W, Liu D, Vanden Eijnden E: Nested stochastic simulation algorithm for chemical kinetic systems with disparate rates. Journal of Chemical Physics. 2005, 123: 194107-10.1063/1.2109987.View ArticlePubMed
- Goutsias J: Quasi-equilibrium approximation of fast reaction kinetics in stochastic biochemical systems. Journal of Chemical Physics. 2005, 122 (18): 184102-10.1063/1.1889434.View ArticlePubMed
- Gillespie DT: The chemical Langevin equation. Journal of Chemical Physics. 2000, 113 (1): 297-306. 10.1063/1.481811.View Article
- Vanden-Eijnden E: Numerical techniques for multi-scale dynamical systems with stochastic effects. Comm Math Sci. 2003, 1 (2): 385-391.View Article
- Chatterjee A, Vlachos DG: Multiscale spatial Monte Carlo simulations: Multigriding, computational singular perturbation, and hierarchical stochastic closures. Journal of Chemical Physics. 2006, 124 (6): 0641101-06411016. 10.1063/1.2166380.View Article
- Lam SH: Using CSP to understand complex chemical kinetics. Combustion Science and Technology. 1993, 89 (5–6): 375-404.View Article
- Lam SH: The CSP method of simplifying kinetics. International Journal of Chemical Kinetics. 1994, 26: 461-486. 10.1002/kin.550260408.View Article
- Valorani M, Goussis DA, Creta F, Najm HN: Higher order corrections in the approximation of low-dimensional manifolds and the construction of simplified problems with the CSP method. Journal of Computational Physics. 2005, 209: 754-786. 10.1016/j.jcp.2005.03.033.View Article
- Cao Y, Gillespie DT, Petzold LR: Efficient step size selection for the τ-leaping simulation method. Journal of Chemical Physics. 2006, 124: 044109-10.1063/1.2159468.View ArticlePubMed
- Takahashi K, Kaizu K, Bin H, Tomita M: A multi-algorithm, multi-timescale method for cell simulation. Bioinformatics. 2004, 20 (4): 538-546. 10.1093/bioinformatics/btg442.View ArticlePubMed

## Copyright

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.