Depending on the frequency of reactions, the dynamics of biomolecular networks are accurately represented by discrete and stochastic models, the stochastic simulation algorithm (SSA) and the chemical Langevin equation (CLE), respectively. At higher concentrations, the CLE is approximated by deterministic ordinary differential equations. In such networks, large number of species are typically involved in reactions characterized by a wide range of time scales. In this context, reduced order models capturing the essential dynamical features of the overall dynamics are necessary to design a self-consistent and efficient simulation algorithm integrating the different descriptions of a system of reactions.

In this work, we focus on the CLE in the presence of reactions whose propensities span different orders of magnitude, yet remain in the regime where continuous stochastic descriptions are valid. The objective is to develop a systematic framework to derive a non-stiff description of the dominant slow dynamics of this multi-scale CLE. A first step consists of seeking coordinate changes that lead to an explicit identification of fast and slow variables. To this end, we present sufficient and necessary conditions to derive a diffeomorphism to decouple fast and slow variables. This decoupling step allows us to extend the method of adiabatic elimination (see Gardiner, Handbook of stochastic methods) to the systems under consideration. In this approach, the fast variables are assumed to relax to a pseudo-stationary density given by assuming that the slow variables are constant. After determination of this pseudo-stationary density, projection formalism is applied to obtain an approximate Fokker-Planck equation for the probability density of the slow variables. An approximate solution of the original Fokker-Planck equation is then computed by multiplying the fast variable pseudo-stationary density to the slow variable density.

As an illustration, this analysis framework is applied to a representative biochemical network modeled by a CLE system exhibiting time-scale multiplicity due to the presence of fast and slow reactions. An approximate distribution is derived using the proposed method. Through simulations, we show how well this approximate distribution compares to the distribution obtained by solving the original Fokker-Planck equation.