BACKGROUND
at a molecular level, the functioning of cellular processes is unavoidably stochastic.
recent advances in real-time single cell imaging, micro-fluidic manipulation and synthetic biology have shown that gene expression and protein abundance in single cells are dynamic random variables submitted to significant variability  <cit> .
markov processes approaches to gene networks dynamics, originating from the pioneering ideas of delbrück  <cit> , capture diverse features of the experimentally observed variability, such as transcription bursts  <cit> , various types of steady-state distributions of rna and protein numbers  <cit> , noise amplification or reduction  <cit> .
in the case of molecular switches, random transitions between two or several metastable states with different transcriptional activities, limit the possibility of the corresponding genetic circuits to store cellular memory and generate variability  <cit> .
even the simplest stochastic model, such as the production module of a protein  <cit>  contains tens of variables and biochemical reactions, and at least the same number of parameters.
direct simulation of such networks by stochastic simulation algorithm   <cit>  is extremely time consuming.
network reverse engineering, model parameter identification, and design principles studies suffer from the same drawbacks when the full models are attacked with traditional tools.
various approximation methods were proposed, based on time discretization such as binomial, poisson or gaussian time leaping  <cit> , or based on fast/slow partitions  <cit> .
the time complexity of exact simulation algorithms scales with the number of jumps  per unit time.
currently available hybrid algorithms treat fast reactions as continuous variables, which significantly reduces the number of jumps  <cit> .
although computationally efficient, the use of slow reactions as discrete variables, and of rapid reactions as continuous variables, does not always provide a convenient description of the hybrid stochastic process.
indeed, the discrete/continuous structure of the state space of a hybrid molecular system is better expressed in terms of species components.
activity of some promoters and production of the corresponding proteins occurs in bursts  <cit> .
in a bursting event, a steep increase of the number of transcripts is followed by relatively smoother exponential degradation.
in this example, the natural candidates for continuous variables are the gene transcripts and products concentrations, while the states of the dna promoters are the discrete variables.
this picture also emphasizes the possibility for noise transfer from discrete, microscopic components, to continuous, mesoscopic, and closer to the phenotype, variables.
in such processes, both microscopic and mesoscopic variables are noisy.
at the mesoscopic scale, the stochastic fluctuations result from the intermittency of the dynamics and do not have a simple dependence on the numbers of molecules of the continuous components.
counter-intuitively, the noise amplitude could increase with the number of transcripts or proteins in a burst  <cit> .
stochastic intermittence or bursting occurs not only in transcription, but also in many other biochemical processes, such as action potential firing  <cit> , blips in calcium signaling  <cit> , possibly happloinsufficiency  <cit> , etc.
thus, a large class of molecular stochastic models can be approximated by stochastic hybrid processes in which some state components are discrete, while others are continuous.
processes of this kind have been intensively studied in relation to applications in technology   <cit> .
several classes of stochastic hybrid processes were studied, among which the most important are the switching diffusion processes and the piecewise deterministic processes.
for switching diffusions, the continuous state component has continuous trajectories governed by stochastic differential equations, punctuated by jumps controlled by the discrete component.
for piecewise deterministic processes, the continuous component follows deterministic ordinary differential equations between two jumps.
the jumps can be changes of the parameters in the  differential equations , or instantaneous changes of the continuous variables , or both.
piecewise deterministic processes have been recently proposed as models for various biological systems  <cit> .
diffusion approximations  of markov processes justify approximate versions of the gillespie algorithm, such as the τ-leap method  <cit> .
a general approach allowing to obtain and classify various hybrid approximations that can be obtained from molecular markov jump processes is still needed and represents a first result of this paper.
hybrid stochastic processes can be obtained by applying the law of large numbers  to the continuous components.
this will be performed here by using a partial kramers-moyal expansion of the chemical master equation.
in many cases, this procedure provides only a moderate reduction of the number of stochastic jumps to be simulated.
indeed, remaining discrete components can jump rapidly, in which case the simulation by a hybrid process performs ineffectively.
the second result of this paper allows us to cope with frequent jumps of the discrete components by using averaging.
averaging techniques are widely used in physics and chemistry to simplify models by eliminating fast microscopic variables  <cit> .
in our case, the microscopic variables are the discrete state components  and the resulting approximations are averaged switched diffusion or averaged  deterministic processes.
standard averaging techniques for discrete markov chains  <cit>  use state aggregation.
however, state aggregation algorithms are effective only for markov chains with a finite, small number of states.
a more effective averaging method, that we present here, is cycle averaging in chemical reaction networks.
our method exploits the formal analogy between the quasi-stationarity assumption for fast deterministic linear cycles  <cit>  and the stochastic averaging assumption for the same type of submodels.
stochastic quasi-steady-state approximations presented in  <cit>  combine singular perturbations for the master equations and van kampen's Ω-expansion  <cit>   and are very close to our approach.
however, we provide here a much more general setting and completely characterize the hybrid limits that result from this setting.
the stochastic dynamics of our simplified hybrid models provide good approximations of the un-reduced chemical master equation.
we show rigorously elsewhere  <cit>  that the simplified hybrid processes are weak approximations of the un-reduced markov process.
this means that all the statistical properties of the trajectories are approximated in distribution.
in particular, steady state distributions of molecular species concentrations should be accurately reproduced.
also, distributions of intermittence times like the times between bursts in the production of mrna or proteins, or the commutation time for stochastic switches are expected to be similar for the reduced model and for the initial model.
the structure of this paper is the following.
in the methods section we obtain hybrid simplifications for pure jump markov processes and also propose a new averaging algorithm for these models.
in the "results and discussion" section we present several examples of hybrid models obtained from stochastic chemical kinetics models, including a medium scale protein production model for which we demonstrate averaging.
methods
introductory notes on stochastic hybrid systems
stochastic hybrid processes are couples of the type , x) where θ is a discrete component , x is a vector with piecewise continuous components, and t is time.
the discrete component θ can be considered as a controlled markov chain, whose transition matrix may depend on the continuous variable x. the discrete state space Θ  can be finite or infinite.
for instance, supposing that the discrete variables are r "rare" molecular species, whose numbers never go over a small value n, then Θ ⊂ { <dig> ..., n}r has at most r states.
different values of θ may also correspond to various states of a single molecule.
for proteins with multiple phosphorylation sites the number of states is 2r where r is the number of sites.
in the case of competitive transcriptional regulation, when n transcription factors are competing on r binding sites, the number of states is at most r.
there are several possible ways to specify the transitions of the discrete variables.
generally, these can be given by the stochastic matrix λi, j  where λi, j dt + o is the probability of a transition from the state i to the state j, between t and t + dt.
however, this representation is not handy when the matrices have very large dimension.
in such cases, the usual molecular markov jump process description is more appropriate.
then, possible transitions are specified by the stoichiometry vectors γi.
the intensity  function λi  represents the probability per unit time that a transition from θ to θ + γi takes place between t and t + dt  <cit> .
it follows that the probability for a transition  between t and t + dt is λ dt + o, where λ = ∑i λi is the total intensity.
the inverse of the total intensity represents the mean waiting time between two successive transitions.
the transition θ → θ + λi is chosen with probability λi/λ.
jumps can also occur in the continuous variables, in which case the size of the jump can be continuously distributed.
there are several possibilities for the evolution of the continuous variables leading to various classes of hybrid processes.
piecewise deterministic processes
piecewise deterministic processes were first formalized by  <cit>  and found many applications in various areas ranging from industry to mathematical finance and biology.
the state of a piecewise deterministic process  is a couple ζ =  where θ ∈ Θ denotes the discrete variable and x ∈ ℝn denotes the continuous variable.
a pdp is given by three local characteristics namely:
• a vector field χθ  ,
• a jump intensity λ and
• a transition measure qt such that qt is the probability distribution of the post-jump positions  ∈ Θ × ℝn.
related to this is the jump distribution qj giving the probability distribution of the jumps: qj = qt .
such a process involves deterministic motion punctuated by jumps.
the lengths of the time intervals between successive jumps are random variables called waiting times τi.
on the time interval  the discrete variable remains constant θ = θi, while the continuous variable evolves according to the differential equation:
let ϕ) be the solution of the differential equation .
we suppose that this solution exists and is unique for all t ≥ ti.
the waiting time is a random variable whose distribution reads  <cit> :
the reader can easily obtain eq.. the probability not to jump in the time interval  is  <dig> - λds + o. thus p = p + o. by taking the logarithm we get log p - log p = -λds + o. summing over ds leads to eq.. the more familiar expression used in the gillespie algorithm p = exp corresponds to the particular case when λ does not depend on the continuous variable.
after a deterministic evolution on the -th interval, the discrete variable and the continuous variable can change according to:
where ,  are sampled from the jump distribution qj.
the discrete variables are parameters of the differential equations describing the dynamics of the continuous variables.
notice that if  =  <dig> the only effect of a jump on the continuous variables is a change of regime .
the trajectories of the continuous variables are globally continuous.
however, the velocities of the continuous variables are discontinuous at jumps of the discrete variables.
we call this possibility "switching".
if  ≠  <dig>  then continuous variables can have instantaneous jumps and their trajectories are only piecewise continuous.
we call this possibility "breakage".
it is possible to have both "switching" and "breakage".
a pdp can be also defined by its generator  <cit> , acting on functions f defined on the phase space:
the adjoint of the generator is used to obtain the fokker-planck equation  <cit>  describing the time evolution of the probability distribution p of the process:
the structure of the generic pdp corresponds to the following simulation algorithm:
set the initial state condition x = x <dig>  θ = θ <dig>
generate a random variable u uniformly distributed in  <cit> ,
integrate the system of differential equations
between ti and ti + τi with the stopping condition f  = u,
generate a second uniform variable v and use it to randomly choose the jumps  ,
change the system state ) into  + ),
reiterate the system from 2) with the new state until a time tmax previously defined is reached.
hybrid diffusions
a rather general class of hybrid diffusions has been proposed in  <cit> .
however, in that setting, jumps of continuous variables are commanded by the hitting of some predefined phase space sets, which is not our case.
we propose a simpler and different setting, which is natural for biochemical systems.
in our setting, a hybrid diffusion is defined by:
• a vector field χ ,
• a diffusion matrix σ,
• a jump intensity and a transition measure defined like for pdps.
between two consecutive jumps at ti and ti + τi, the continuous variables x satisfy the itô stochastic differential equations  <cit> :
where wj  are independent one-dimensional wiener processes.
the waiting time τ is determined like for pdps.
one should integrate the system:
with the initial condition f =  <dig>  x = xi and the stopping condition f = u, where u is a random variable, independent from , uniformly distributed on  <cit> .
the generator for such a process is:
where "."
stands for the scalar product and ":" stands for the double contracted tensor product, i.e.
.
notice that hybrid diffusions contain pdps as the particular case when diffusion is zero .
hybrid simplification of markov pure jump processes
markovian stochastic chemical systems
most of the existing numerical methods for stochastic chemical systems are based on the representation of the chemical reaction system as a markov pure jump process , that corresponds to gillespie's  <cit>  stochastic simulation algorithm .
the state of a system of n chemical species a <dig> ..., an is a vector x whose components are the numbers of molecules from each species.
the state space of the process is e = ℤn.
from now on, we shall use upper case letters x for numbers of molecules, and lower case letters x = x/ for concentrations, where  is the reaction volume, typically the volume of the cell's compartment.
for each reaction,
we have the jump γi = βi - αi ∈ ℤn, i =  <dig> ..., nr .
we consider equally the reverse reaction, corresponding to the jump -γi.
the intensities and the transition measures for the markov jump process are:
vi, v-i are the reaction rates and probabilities  satisfy .
here δx is the dirac measure centered on x, and x' is the post-jump position.
several choices are possible for the reaction rates  <cit> .
the most popular is the mass action law, when rates are polynomial functions of concentrations:
as for any homogeneous markov process, we can characterize this process by its generator  <cit> , whose action on test functions f defined on the state space e is:
the adjoint of the generator defines the chemical master equation which describes the time evolution of the probability distribution p in phase space.
notice that hybrid diffusions or pdps with no continuous components are markov pure jump processes.
partial kramers-moyal expansion
diffusion approximations  for stochastic chemical kinetics was proposed by  <cit>  and was rigorously justified by  <cit> .
this approximation is valid when the numbers of molecules of all species are much larger than one.
it concerns for instance the thermodynamic limit  → ∞ where  is the volume of a well stirred reactor.
then, the diffusion approximation applies to intermediate scales  in phase space.
formally, this approximation can be obtained by the kramers-moyal expansion  <cit> , which is the second order taylor expansion of the master equation with respect to jumps divided by the volume; in the first order one gets the deterministic dynamics, and in the second order the langevin dynamics.
in  <cit> , we have used partial kramers-moyal expansions to obtain hybrid approximations of mpjp.
the possibilities of this method are examined here in more generality.
first, let us partition the species of the system into two sets x = , where xd represents the species in small numbers and xc the species in large numbers.
note that this partition is different from reaction based partitions used in other works  <cit> .
all these works have as objective the reduction of simulation time.
with respect to this objective, many schemes that avoid individual simulation of rapid reactions by regrouping them and by treating the corresponding channels as continuous variables in "fluid" limits, are more or less equivalent.
we fix ourself also another objective, which is to find a natural representation of the simplified processes.
a species based partition is more suitable for this end.
let m be the set of reactions in which every reversible reaction contributes with two reactions, one for each direction.
the species partition determines a partition of m into three sets .
the reaction partition distinguishes between the reactions whose reactants and products are in xc (), reactions whose reactants and products are in xd () and reactions with at least a reactant or a product in xd and at least a reactant or product in xc ().
we should emphasize that although  reactions are rapid, some other rapid reactions can be contained in  or .
the second step towards the hybrid processes is to rescale the continuous variables by dividing them with the volume and obtaining xc = xc/.
the xd variables will not be rescaled and will be considered discrete.
their evolution is given by discrete transitions.
the third step is to consider a second order taylor expansion of the master equation, for the continuous variables xc and with respect to the small parameters γi/.
this expansion can be performed equivalently on the generator  or on the master equation .
for simplicity, we use the generator.
the test functions depend now on two variables  and are considered to be twice differentiable with respect to the continuous variables xc.
let us consider that rates of reactions in  are proportional to the volume , , i ∈ .
the second order taylor expansion of  with respect to xc reads:
where ⊗ is for tensor product kl = γk γl).
the approximated generator corresponds to hybrid diffusions with drift function  and diffusion matrix .
as it can be easily seen, the drift and diffusion coefficients for the variables xc do not depend on the discrete variables xd.
thus, switching is not present in this approximation and discrete variables can be just forgotten if not observed.
furthermore, jumps  vanish in the limit  → ∞, implying that continuous variables have no instantaneous jumps .
however, switching and breakage is present in many examples in molecular biology.
how is this possible?
in order to cope with the possibility of switching and breakage we need to consider that some reactions from  exhibit particular properties.
we can consider two types of special reactions:
<dig>  reactions in  that are frequent enough to change the continuous variables and induce switching,
<dig>  reactions in  with a large stoichiometric vector, that induce breakage.
these reactions will be referred to as "super-reactions".
the first type of super reactions are denoted by the set  and the second type of super-reactions are denoted .
as a consequence, we distinguish two types of hybrid approximations following which one of the two "super-reactions" are present in the mechanism.
a. hybrid approximations with switching
this first case of approximate process is induced by "super-reactions"  with stoichiometric vectors of order one, but with rates that scale with the volume: , i ∈ .
reactions from  induce rapid transitions of the continuous variables.
the set  is considered to be empty.
for this approximation, even if the rates of the "super-reactions" depend on the discrete variables, we suppose that the "super-reactions" do not change the discrete variables.
in other words
if this condition fails, another type of approximation arises: some discrete variables change rapidly.
the resulting process will be an averaged hybrid process.
the evolution of the discrete variables is given by slow jumps .
the component xc has rapid jumps but their amplitude is relatively small.
thus, between two consecutive transitions of the discrete variables, the continuous variables can be approximated by a smooth deterministic trajectory.
the continuous variables have continuous trajectory, however, discontinuous velocities.
the first order of the kramers-moyal expansion corresponds to the piecewise deterministic approximation defined by the flow function
the jump intensity
and the jump probability
where ,  are the projections of γi on continuous and on discrete components, respectively.
the relation  clearly shows that the jump intensity depends both on continuous and on discrete variables.
nevertheless, only the discrete variables are changed by jumps.
the second order of the kramers-moyal expansion corresponds to the hybrid diffusion approximation whereas the diffusion matrix reads:
b. hybrid approximations with breakage
this type of approximation is induced by the "super-reactions" with a large stoichiometric vector .
we suppose also, that the set  is empty.
in other words, to have this approximation, we suppose that at least one reaction i in the set  has a large stoichiometric vector j/ is comparable with |xj| for i ∈  and for some species j ∈ c. a reaction  can change both continuous and discrete components.
the main difference consists in the fact that a reaction  significantly change, in one step, the concentration of the xc component.
the reactions of the type  do not contribute to the deterministic flow.
they only appear as jumps.
as before, we write down the main characteristics of the approximated pdp process , namely the flow function
the jump intensity
the jump probabilities for the discrete variables
and the jump probability for the continuous variables
if the super-reactions i ∈  have  ≠  <dig>  jumps occur simultaneously in the continuous and in the discrete variables.
the deterministic flow is defined by the flow function χ ) which does not depend on the discrete variable xd.
in this case there is no switching.
in the second order of the kramers-moyal expansion, we include phase space diffusion, which contains contributions only from reactions of :
if both sets ,  are not empty, then both types of discontinuities can appear in the trajectory of xc:
• switching: discontinuity of velocity , with no discontinuity of trajectory,
• breakage: discontinuity of trajectory .
the first discontinuities are induced by reactions from , while the latter are induced by reactions from .
c. hybrid approximations with singular switching
although mathematically distinct, breakage and switching could be physically indistinguishable in certain cases.
indeed, the jump of the continuous variable can result from the repeated application of one or several very fast reactions from .
we can obtain processes with breakage as singular limits of processes with switching when very short lasting steep variations alternate with long lasting slower variations of the continuous variables.
let us consider that there is a subset of extremely rapid super-reactions  such that:
where  <dig> <ϵ <<  <dig> is a small positive parameter, and where  are discrete species, substrates of the reaction i ∈ .
like for reactions , we consider that  =  <dig>  for all i ∈ .
for each reaction i ∈  we consider that there are two subsets of reactions in :  contains reactions that produce the species  and  contains reactions that consume the species .
we define .
we consider that all reactions in  are rapid:
first, we apply the kramers-moyal expansion for large  and obtain a switched pdp.
the flow function of this pdp reads:
where
the jump intensity of the pdp reads
we show in the additional file  <dig> that, in the limit ϵ →  <dig>  the switched pdp converges to a pdp having the following generator:
where  is a probability density satisfying  and Φ is the solution of the differential equation .
we recognize above the generator of a pdp with breakage.
contrary to the previous case of breakage resulting from super-reactions , the breakage size is now continuously distributed.
the switch transitions  disappear in the singular limit .
practical criteria for identifying small parameters and super-reactions leading to piecewise deterministic approximations
the law of large numbers is applicable in the limit  → ∞.
certainly, in cell biology, the idea of infinite volumes should be considered with care.
for this reason we will replace this condition by a set of easy to check criteria concerning relative orders of parameters of the models.
these criteria will concern piecewise deterministic approximations.
criteria for hybrid diffusion approximation, involving central limit theorem, will not be discussed in this paper.
our criteria are relative to two reference quantities.
the first reference is a large number n representing a lower bound for the numbers of molecules of continuous species.
the second reference is a time τ, which is a lower bound both for the characteristic times of the deterministic dynamics of the continuous variables τc and for the average time between two successive jumps of the discrete variables τd, namely τ = min.
τ can be related to kinetic parameters by methods exposed in  <cit> .
the applicability of the law of large numbers to continuous variables implies two conditions.
the first condition is that jump sizes should be small with respect to the numbers of molecules, i.e.
<<n for all i ∈ .
the second condition is that the number of jumps of the continuous variables on the timescale τ should be big, i.e.
vi τ >>  <dig> for all i ∈ .
the two conditions are fulfilled simultaneously in the limit  → ∞ because n, vi scale with  and , τ do not depend on .
a similar condition must be satisfied by super-reactions of the type 1: vi >> τ- <dig> also for all i ∈ .
the super-reactions of the type  <dig> do not satisfy the small jump condition, i.e.
must be comparable to n  for all i ∈ .
in order to have singular switching we need a new small parameter ϵ <<  <dig> and the following set of conditions:
a) there are super-reactions of the type  <dig>  these are just very quick super-reactions of the type  <dig>
b) super-reactions of the type  <dig> act only on a laps of time τϵ = ϵτ i.e.
vi = - <dig> >> τ- <dig> for i ∈ .
reactions in  that inactivate super-reactions of type  <dig> could thus be as frequent as reactions in  and in  .
c) super-reactions of the type  <dig> are quick enough to produce breakages comparable to n during the time τϵ, i.e.
vi >> - <dig> for i ∈ .
the conditions of the type vi >> τ- <dig> can be simplified if the reactions are first order with respect to discrete species.
in this case vi = ki xd with xd ≈  <dig>  therefore the condition is ki >> τ- <dig>  all these criteria are summarized in the table  <dig>
see text for definitions of various reaction sets.
τc, n are typical timescale and number of molecules of the continuous components, respectively.
τ = min, where τd is the mean time between two successive jumps of the discrete components.
vi are the average numbers of reactions per unit time , k are first order kinetic constants and γi are the stoichiometric vectors.
averaging for stochastic chemical kinetics
the performance of the hybrid algorithm can be very bad when the discrete mechanism contains rapid cycles which effectuate many reactions on the deterministic timescale τc .
indeed, in this situation the deterministic solver is artificially sampling the interval between two discrete cycle reactions.
first, this leads to unreasonable increase of the simulation time.
second, the condition on the number of jumps of continuous variables is not satisfied and the hybrid approximation is not accurate.
in this case, the hybrid scheme performs worse than ssa.
it is therefore important to eliminate rapid cycling from the system before implementing numerical schemes for the hybrid approximations.
this can be done by averaging.
averaging principles are widely used for deterministic  and stochastic  systems.
the classical averaging idea is to identify fast ergodic  variables and to average the dynamics of the slow variables with respect to the quasi-stationary distribution of the fast variables.
an important difficulty is to identify the right slow and fast variables.
for perturbed hamiltonian dynamical systems the action-angle variables provide the natural framework for averaging  <cit> .
in the case of chemical kinetics, there are two kinds of slow variables that should be taken into account for averaging.
the first kind are just slow components .
the second kind are linear combinations of components that are conserved by the fast cycling dynamics.
these new slow variables, that play a role similar to the action variables in hamiltonian systems, provide new "aggregated" or "lumped" species in the averaged system.
notice that aggregation of states has been used for averaging markov chains  <cit> .
the aggregation of species that we propose here adapts this type of argument to the case of chemical kinetics.
variable aggregation for simple deterministic reactive models has been used in relation to applications in ecology .
in  <cit>  we proposed a general solution for simplification of reaction networks, that works for any linear mechanism with well separated constants.
in this algorithm, species aggregation is systematically applied to hierarchies of rapid cycles.
we show here that the same algorithm can be effectively used also in the stochastic case.
averaging principles for reaction mechanisms
averaging can be applied to multi-scale reaction mechanisms.
this procedure leads to averaged hybrid models.
in order to apply averaging, one should first identify a sub-model of discrete components, satisfying the following conditions:
c1) the dynamics of the sub-model is much faster than the dynamics of continuous and of other discrete, slower modes.
c2) the dynamics of the sub-model is ergodic: starting with any state, the system can reach any other state in finite time.
the general procedure to obtain averaged hybrid simplifications is described in the additional file  <dig>  some cases are particularly interesting.
the first case corresponds to fast discrete cycles producing continuous species.
in this case rapid super-reactions  change both discrete and continuous components  is not valid).
the discrete components that are affected by such reactions are fast discrete variables.
in the resulting hybrid model, both the continuous dynamics and the slow discrete reaction rates should be averaged with respect to the fast discrete variables.
the second case corresponds to fast, purely discrete cycles.
in this case, some of the cycles of the sub-mechanism  are at least as fast or faster than reactions .
the other reactions of  are much slower.
in the resulting model, rates of the slow reactions of  should be averaged with respect to the fast discrete variables.
in both cases, one needs the steady probability distribution for the fast discrete sub-model.
all the slow processes will be averaged with respect to this distribution.
the calculation of the steady probability distribution can be easily done only when the sub-model is linear.
considering linear sub-models has also another advantage.
in this case, averaging and aggregation lead to a coarse-grained reaction mechanisms.
for these reasons, we have developed an algorithm for linear sub-models.
of course, gene network models are generally nonlinear.
however, this does not mean that all the parts of such models behave nonlinearly.
many sub-models, in particular monomolecular cycles, can be simplified by averaging methods that we designed for linear networks.
cycle averaging for linear sub-models
linear sub-models
the are two types of linear reaction networks: monomolecular networks and first order networks.
the structure of monomolecular reaction networks can be completely defined by a simple digraph, in which vertices correspond to chemical species ai, edges correspond to reactions ai → aj with rate constants kji >  <dig>  for each vertex, ai, a positive real variable ci  is defined.
the deterministic kinetic equation is
first order reaction networks include monomolecular networks as a particular case, and are characterized by a single substrate and by reaction rates that are proportional to the concentration of the substrate.
first order reaction networks can contain reactions that are not monomolecular, such as a → a + b, or a → b + c. we shall restrict ourselves to pseudo-conservative first order reactions, ie reactions that do not change the total number of molecules in a given submechanism .
with such constraints, the total number of molecules in the sub-mechanism is conserved and the kinetic equations are the same as .
degradation reactions can be studied in this framework by considering a special component , that collects degraded molecules.
further release of the constraints is possible.
for instance, the system can be opened by allowing constant  production terms in eq.. these terms will change the steady state, but will not influence the relaxation times of the system.
the linear sub-mechanisms can be considered as part of a nonlinear network, given fixed  values of external inputs .
for instance, even in systems of binary reactions, one can define pseudo-monomolecular reactions when one of the substrates of the binary reaction is not changing .
this condition can be fulfilled if the substrate is in excess, for instance.
the stochastic dynamics of a unique molecule in such linear reaction network is given by the probability p that the molecule is in aj at the time t. we can easily show that the master equation for p is the same as the deterministic kinetic equation .
considering only one molecule does not restrict generality because when several molecules are present in a linear network, these behave independently.
thus, the simplification algorithm proposed for deterministic networks  <cit>  can be also applied to stochastic networks  <cit> .
the algorithm is based on a set of operations transforming the reaction graph into an acyclic digraph without branching .
for such type of graphs the eigenvectors and eigenvalues of the kinetic matrix can be straightforwardly calculated, which completely solves the problem of deterministic dynamics.
we could follow precisely the same approach to simplify and then solve the stochastic dynamics.
however, we argue here that applying only a few steps of the algorithm is enough for effectively reducing computational time of the ssa method.
cycle averaging algorithm
let us recall that the reduction method presented in  <cit>  uses two types of operations acting on the reaction graph.
i. construction of an auxiliary network .
for each node ai of a linear sub-mechanism, let us define κi as the maximal kinetic constant for reactions ai → aj: ki = maxj{kji}.
for the corresponding j we use the notation ϕ: ϕ = arg maxj{kji}.
an auxiliary reaction network  is the set of reactions ai → aϕ with kinetic constants κi.
in such a network there is no branching: if several reactions consume the same component ai, only the quickest one is kept and all the other discarded.
ii.
glueing cycles .
rapid cycles are replaced by a single node.
constants of reactions leaving these cycles are renormalized according to an averaging principle .
in order to present the simplification algorithm let us use two simple examples.
first, let us consider a chain of molecular reactions a <dig> → a <dig> → ... am.
the reaction rate constant for ai → ai+ <dig> is ki.
all rate constants are considered well separated, i.e.
either ki <<kj or ki >> kj for any i ≠ j.
the smallest rate constant in the chain is called limiting, and denoted by klim.
if 1/klim <<τc , then all molecules a <dig> are transformed into molecules am on a timescale much quicker than the deterministic time τc.
we can thus ignore the chain reactions and consider that the entire mass of the chain is practically always in am.
this is equivalent to considering the chain at quasi-stationarity because the steady state probability distribution of a chain is a dirac delta measure localized at the end of the chain.
however, if we do not simplify chains, then simulating them by gillespie's ssa will not be computationally expensive because the mass of the chain is transferred to the end of the chain amin a number of steps that is relatively small .
as a second example, let us consider the cycle c be the following sequence of mono-molecular reactions a <dig> → a <dig> → ... am → a <dig>  let all rate constants be well separated and the smallest one be klim like before.
we add to the cycle one branching reaction; this transforms aj a component of the cycle into b a component exterior to the cycle.
we consider the following distinct situations:
the branching reaction is aj → b of rate constant k and k >> kj,
the branching reaction is aj → b and k <<kj,
the branching reaction is aj → aj + b, or
the branching reaction is aj → aj+ <dig> + b of rate constant kj.
in the situation  the exit reaction is faster and dominates the cycling reaction aj → aj+ <dig>  according to the rule for auxiliary networks in this case  the cycle can be opened  and the resulting multiscale dynamics is the one of a chain; we recover the previous example and in this case we can safely decide to do nothing.
in the situation  the exit reaction is much slower than the cycling reaction.
in this case the molecules inside the cycle have rapid transformations and the mass distribution inside the cycle can be considered to reach quasi-stationary distribution.
we call this cycle "unbroken".
as discussed in  <cit> , the relaxation time of a cycle with separated constants is the inverse of the second slowest rate constant k >> k = klim.
to understand this, one should consider the two possible paths to equilibrate a cycle, one passing by the slowest step and the quicker one passing by the second slowest step: the quicker short-cuts the first one.
thus, a cycle can be considered quasi-stationary if k >> 1/τc.
a non-averaged fast cycle is computationally expensive in ssa, if a molecule can perform a huge number of steps along the cycle on the timescale τc.
the corresponding condition involves the quasi-stationary flux  and reads k = klim >> 1/τc.
from a quasi-stationary cycle, the mass is lost stochastically, but slowly by the branching reaction.
the intensity of the loss process can be calculated by replacing xj by its average with respect to the quasi-stationary distribution of the cycle.
the average of xj is  = n klim/kj, where n is the total mass inside the cycle .
we obtain the average intensity  = n kklim/kj.
in the situations  or  the average intensities of the branching reactions are  = n kklim/kj and  = n klim, respectively.
all these operations can be presented as a:
simplification algorithm
input:
a first order reaction mechanism g with separated kinetic constants.
output:
a simplified first order reaction mechanism.
while there are fast unbroken cycles
for each cycle in g not containing reactions of the type  having a sufficiently intense flux  do
1: "glue" the cycle into a single node c having the total mass n;
2: replace the exit reaction of the type  aj → b of rate constant k by a reaction c → b of effective constant k' = kklim/kj;
3: replace the reaction of the type  aj → aj + b or rate constant k by a reaction c → c + b of effective constant k' = kklim/kj;
4: replace the reaction of the type  aj → aj+ <dig> + b of rate constant kj by a reaction c → c + b of effective constant k' = klim;
imposing the "glueing" condition on the cycling flux and not on the cycle relaxation time allows to iterate the cycle averaging algorithm until no more fast unbroken cycles remain.
this would not be possible when adopting a relaxation time criterion to perform averaging.
indeed, the relaxation time of a cycle is given by the second slowest constant.
this leaves place for a counter-intuitive possibility: one can have slow cycles embedded into rapid cycles.
for instance, a fast cycle whose nodes are "glued" cycles can have a slower node as the beginning of its limiting step.
indeed, the internal constants of this node are necessarily more rapid than the limiting step of the large cycle, but can be slower than the second constant which gives the relaxation time of the large cycle.
the slow scales are lost by glueing.
in order to recover the full multi-scale dynamics, a restoration operation has been considered at the end of the algorithm presented in  <cit> .
in this operation, all "glued" cycles are restored without their limiting steps.
thus, slower cycle sub-dynamics can be recovered.
restoration is not needed for glued cycles satisfying the flux condition klim >> , because these are automatically faster than the timescale τc.
our averaging method works for mechanisms with well separated constants.
generalization for partially separated constants are possible.
for instance, one can consider cycles containing reactions with equivalent constants, but such that for each node of the cycle, the constant of the branching reaction is much smaller than the constant of the cycling reaction.
noisy networks and design rules for hierarchies of cycles
as we have shown in the previous sections, intrinsic noise in biochemical networks is generated at a "microscopic level" by the discrete variables and can be observed at the "mesoscopic level" of the continuous variables either as switching, or as breaking events.
thus, when there are no discrete variables , there is no intrinsic noise.
also, if the switching events are much faster than the deterministic time scale, averaging principles apply and noise is not transmitted to the continuous variables: the deterministic approximation is again recovered.
the only way to transmit noise by switching to the mesoscopic level is by intermittency and this needs particular combinations of slow reactions that change the values of the discrete species and frequent reactions that change the continuous species.
intrinsic noise transmission from micro- to meso-scopic level results from certain  combinations of low numbers, fast and rapid timescales.
some general design rules relating topology to dynamic properties relative to noise can be derived from our approach.
by hierarchies of cycles we understand reaction mechanisms containing cycles connected by branching reactions forming higher level cycles.
the nodes of higher level cycles are "glued" cycles from the lower levels.
hierarchies of cycles in discrete variables have interesting properties with respect to noise production.
generally, in cycle hierarchies, effective constants of branching reactions are at least as slow or slower than the limiting step  of the node  from which they start.
coupling cycles into hierarchies is a way to produce slower and slower reactions from initially rapid reactions and generate thus intermittency.
as a possible design rule, we could state: exit reactions of the type  or  ) generate intermittency when the exit node is not the beginning of the limiting step in some unbroken fast cycle.
indeed, unless kj is the limiting step in the cycle, one has klim/kj <<  <dig>  then, the average intensity of the exit reaction of the type  or  is weak and could represent a source of intermittency in the system.
this situation should be avoided for less noise in the system, or created when noise is wanted.
an example of how this rule applies to the biochemistry of bacteria will be given below.
a systematic investigation of the possibilities of this class of design rules will be presented elsewhere in relation to synthetic biology.
RESULTS
methodology to obtain hybrid simplifications
we demonstrate how hybrid simplifications can be obtained from markov pure jump models for gene networks.
the simplified models preserve the main stochastic properties of the initial complex models and can replace these models in simulations.
the simplification procedure consists of four successive steps:
i identification of discrete and continuous variables, partition of the reactions.
ii cycle averaging.
iii identification of super-reactions.
iv construction of the hybrid simplifications.
in order to justify the utility of our approach we show by examples that all types of hybrid processes that we have discussed are represented in gene networks models.
to introduce the examples we employ the following notation.
reactions are numbered by integers.
if the i-th reaction is reversible, then  is the i-th reaction in the forward direction and  the i-th reaction in the reverse direction.
irreversible reactions are denoted simply ri.
the rate constants units are s- <dig> for monomolecular and s-1- <dig> for binary reactions.
in order to obtain pseudo-monomolecular rate constants from binary reaction constants with slowly varying substrates x we have used the following formula:
where na is the avogadro number and  is the cell volume.
for a bacterium, na  ≈ 1- <dig>
similarly, the reaction rates are calculated as
we discuss two simple hybrid models, one with switching the second one with breakage, then two more complex models that need averaging.
all our simulations were performed using matlab  <dig> in a windows xp <dig> environment with a dual core intel  <dig>   <dig>  ghz processor.
hybrid model with switching: cook's model
the simplest model with switching has been introduced by cook  <cit>  as a model for haploinsufficiency phenomena.
this model can be described by the following system of biochemical reactions:
in this model, g, g* represent inactive, respectively active states of the gene.
in the active state, the gene produces some protein p. let g <dig> be the gene copy number, which is a conserved quantity of the dynamics g + g* = g <dig>  the haploinsufficiency regime corresponds to a small value of g <dig>  for the simplicity of the argument we consider g <dig> =  <dig>
let x =  be the state vector.
we notice that, the partitions of the species is the following: the species in large number xc = {p}, and the species in small number xd = {g, g*}.
this partition of the species defines a partition of reactions.
with the above notations we get:
the deterministic timescale is τ = - <dig>  the stoichiometric vectors have all lengths of order one, much smaller than n .
for the parameters values used in  <cit>  we have k2/k <dig> >>  <dig> and r <dig> is a super-reaction of the type .
the reaction r <dig> satisfies vi = n k <dig> >> k <dig> = τ- <dig>  which means that the law of large numbers can be applied to the continuous variable.
in the first order of the kramers-moyal development we obtain a pdp approximation with switching.
in the following, we will denote by x =  the continuous variable and by θ = g* the discrete variable.
this model is a piecewise deterministic process with state space  ∈ { <dig>  1} × ℝ.
the flow function χ is given by
and the jump intensity λ is defined by
the resulting pdp is the same as on-off systems studied in operational research  <cit> .
the model has been proposed as abstract model for stochastic protein production in several other places  <cit> .
with the pdp description of the system and using the hybrid algorithm introduced in the section methods, we draw a possible time evolution of this system .
the same initial condition was used when the gillespie algorithm was implemented for the model described above.
the simulation time to generate a trajectory using a pdp approximation, was  <dig>  seconds, while with the gillespie algorithm the simulation time was  <dig>  seconds in the average.
the trajectories look qualitatively the same: production intervals are followed by degradation intervals, with no discontinuity in the continuous variable.
in order to quantitatively test the accuracy of the pdp algorithm, we have computed the stationary distribution of the number of proteins using a piecewise deterministic simulation and compared it to the estimated distribution from gillespie trajectories.
the numerical methods to estimate stationary distributions are described in the additional file  <dig>  the theoretical curve for the pdp model is a beta distribution.
the variable x = p/ follows the beta distribution , b is the beta function  <cit> .
the result of various comparisons is represented in figure 1b.
cook's model operates in the region k <dig> >> km <dig>  which corresponds to a broken cycle in the discrete variables.
if the opposite inequality is valid k <dig> <<km <dig>  then cycle averaging is needed.
the cycle  should be glued to a point  of mass  = g <dig> =  <dig> and the branching reaction r <dig> becomes , where .
the resulting model is a birth and death model with effective birth rate  and death rate constant k <dig>  provided that k2/k <dig> >>  <dig> , the model can be approximated by a completely deterministic process with flow given by the averaged flow function χ = -k <dig> x + .
hybrid models with breakage: neuroscience and bacterium operator sites
neural systems exhibit stochastic behavior.
stein  <cit>  proposed a simple markovian model for the evolution of a neuron membrane potential.
in this model, discontinuous random jumps of the potential are followed by exponential decay.
we do no discuss here how to obtain this hybrid model from a microscopic pure jump model.
stein's model is a phenomenological representation of very complex electric and biochemical processes.
we demonstrate its properties in terms of trajectories and stationary distribution.
the subthreshold behavior of the membrane potential of a neuron cell is described in this model by a hybrid markov process of continuous variable v with jumps of constant intensity λ and constant amplitude in the continuous variable.
although in the original stein model there are jumps of positive and negative sign, corresponding to excitatory and inhibitory synaptic activations, here we consider only positive sign jumps.
if t is the moment of jump, then v - v = a, where a >  <dig> is the amplitude.
between two jumps v decays according to  where α is a constant.
the generator of such process is:
the steady probability distribution p for such a process satisfies the following delay differential equation:
and p =  <dig>  for all x <  <dig>
analytical solutions are not known for this delay differential equation.
in the additional file  <dig> we describe the numerical scheme to calculate the steady probability distribution.
the comparison between the simulated and calculated distribution is shown in figure 2c-d. in figure 2a-b we show trajectories of the system.
stein's model could also apply  to gene networks.
generalizations of this model, allowing for continuous distributions of the jumps, could be used as an archetype of intermittent activity of a promoter.
hybrid models with breakage can result as singular limits of hybrid models with switching as discussed in the methods section hybrid approximations with singular switching.
the typical case is a bacterium operon.
this model can be found in many places in literature.
a detailed molecular example  <cit>  will be studied as our last example .
a simpler version of this model is proposed by  <cit> .
it consists of the following reactions:
the first reaction is zero order, all the other reactions are first order.
the parameters satisfy k <dig> = ak <dig>  k <dig> = bk <dig>  k <dig> = ϵk <dig>  we consider that b >>  <dig>  ϵ <<  <dig>
from the aspect of the trajectories , the authors of  <cit>  hypothetize that a pdp approximation with breakage is the natural simplification and solve the corresponding stationary hybrid fokker-planck equation for the protein component.
we show here that the hybrid fokker-planck equation given in  <cit>  can be found as an application of our approach.
first, we notice that a pdp limit is applicable to the markov jump process.
the species partition is xd = mrna, xc = protein.
the mrna component follows a birth and death process with birth intensity k <dig> and death intensity  mrna.
using the master equation for birth and death processes  <cit>  and considering that k1/ = aϵ/ <<  <dig>  it follows that the probability to have zero or one molecule of mrna is close to one .
this justifies that mrna is a discrete species.
the partition of the reactions is  = {r <dig>  r3},  = {r2},  = {r4}.
the timescale of the continuous variables is τ = - <dig>  v <dig> = k <dig> x >> τ, where x is the number of proteins, provided that x >>  <dig>  this condition, allowing application of the pdp approximation, is satisfied because b >>  <dig> .
let us show that we have singular switching, equivalent to breakage.
the reaction r <dig> is a super-reaction of the type  because v <dig> = k <dig> = b- <dig> the discrete variable mrna can be considered to remain zero almost all the time, except during the negligible duration of the bursts when mrna =  <dig>  the reaction  corresponds to the transition  <dig> →  <dig>  while the reaction  corresponds to the transition  <dig> →  <dig> of the discrete component mrna.
the intensity of the reaction r <dig> satisfies v <dig> = - <dig>  thus, the mean duration of the burst is ϵτ.
according to the section methods, we are in the case of a pdp with singular switching.
the evolution of the unique continuous variable x  is well approximated by the following hybrid generator ):
which shows that b is the mean number of proteins produced in a burst .
we can notice a continuous distribution of the breakage size , situation different from the stein model.
the fokker-planck equation can be solved in this case.
the steady distribution of the continuous variable x is the gamma distribution  <cit> .
first complex example: λ-phage toggle switch
in this section we revisit a classical example of toggle molecular switch.
the model of cro-repressor  switch in λ-phage, a temperate bacteriophage of escherichia coli, was investigated by many authors  <cit> .
the life cycle of phages has two alternative pathways, lysogenic when the phage duplicates synchronously with the host genome and lytic when the phage produces large amounts of its own mrna.
the switch between these two pathways is controlled by the level of ci protein.
the lysogenic pathway corresponds to large levels of ci, while the lytic pathway corresponds to low levels of ci.
a ci dimer may bind to any of the three operators: or <dig>  or <dig>  and or <dig>  in this order.
by cooperativity, or <dig> and or <dig> are almost simultaneously occupied by ci.
the third site or <dig> will be occupied only when the ci concentration is high enough.
a simpler model, in which or <dig> is absent, may be considered  <cit> .
let ci, ci <dig> and d denote the repressor, repressor dimer, and dna promoter site.
the dynamics of the operator sites is described by the following reaction system:
where the dci <dig> and  complexes denote binding to the or <dig> or or <dig> sites, respectively binding to both sites.
the other reactions concern production, degradation of molecules ci, ci2:
where n is the number of proteins per mrna transcript.
the state vector is x =  where xd = {d, dci <dig>  , dci2ci2}, xc = {ci, ci2}.
note that the choice of discrete variables is dictated here, like in our first example, by the conservation law d + dci <dig> +  + dci2ci <dig> = d <dig> that restricts the numbers of d, dci <dig>  , dci2ci <dig> to small values.
the partition of the species induces the following partition of the reactions:
the interesting property of the λ-phage model is its bistability.
the naive calculation to find steady states uses quasi-stationarity of the deterministic dynamics.
then, we can use the equilibrium equations for reactions ri, i ∈  <cit> , the steady state equation for ci and the conservation law to compute the concentration x of ci at steady-state.
x =  <dig> is one of the steady state, namely the lytic attractor.
it also corresponds to an absorbing state of the markov process .
the lytic state can be rendered non-absorbing by including, like in  <cit> , a small basal production rate of ci protein in the model.
however, this is not important for the illustration of our approach.
if k2/km <dig> = k3/km <dig> it follows that steady states  must satisfy:
where .
let us notice that failure of naive quasi-stationarity calculations was demonstrated for nonlinear subsystems  <cit> .
however, these calculations can be justified here by linear cycle averaging .
bistability occurs when the quartic equation  has real roots, which is the case when σ is small and α is big.
more precisely, the model is bistable if .
like in  <cit>  we have used σ =  <dig>  in order to ensure bistability we have chosen α =  <dig>
the sub-mechanism rdc contains rapid unbroken cycles that need to be averaged before any further approximation of the process.
these rapid cycles  correspond to rapid binding-unbinding of the dimer ci <dig> on the dna.
three steps lead from the unreduced model  <dig> to the averaged model  <dig> :
<dig>  the cycle dci <dig>  dci2ci <dig> is unbroken.
it is glued to the node  whose total mass is equal to the mass of dci <dig> and dci2ci <dig>
<dig>  the limiting step of the cycle is klim = km <dig> <<k <dig> ci <dig>
<dig>  the branching reaction dci <dig> → n x + dci <dig> is replaced by  → n x +  of effective constant .
the reaction dci <dig> → d + ci <dig> is replaced by  → d + ci <dig> with the reaction constant .
after averaging, two more cycles remain in the resulting model  <dig>  however, the rates of the remaining reactions d + dci <dig> → dci <dig> and d + dci <dig> →  are equivalent, which does not allow further application of our algorithm.
furthermore, the slow reactions  → d + ci <dig>  and  → d + ci <dig> produce intermittence and should by no means be averaged.
the next approximation is a first order partial kramers-moyal expansion.
the reactions in rc and the super-reaction r <dig> ∈  contribute to the deterministic dynamics of the continous species defined by the following differential equations:
where x, y are the concentrations of ci and ci <dig>  respectivelly.
the remaining reactions induce the jump mechanism.
the time evolution towards the lysogenic attractor is represented in figure 3b for the un-reduced model  <dig> and in figure 3c for the pdp model  <dig> .
steady probability distribution close to this attractor is represented in figure 3d for all models in this study.
we can notice the intermittent behavior of the ci, ci <dig> components that is well captured in the switching pdp approximations .
second complex example: stochastic bursting of a repressed bacterium operon
under strong repression, protein production from a bacterium operator site undergoes stochastic bursting.
a stochastic model for protein production in prokaryotes has been introduced in  <cit> .
the behavior of the operator site under the regulation of a repressor molecule that prevents protein production has been considered in  <cit> .
the corresponding model is represented in the table  <dig> and in figure  <dig>
in this model, the bacterium is considered to be in exponential growth phase, increasing size and dividing normally.
cell growth is simulated by a linear increase of the volume in time.
during replication, the nuclear material doubles .
at fission, the nuclear material is halved and all other components are divided among daughter cells according to a binomial distribution.
i. species and reaction partition
some of the components  are confined to small numbers by the conservation law d + d.r + d.rn ap + tr rn ap = const.
two more components rbs and ribrbs have small lifetimes ≈ 2s and can not accumulate to significant levels.
thus, the discrete variables are:
the remaining components are in large numbers and form the continuous variables:
we notice that rib, rnap and r are in relatively large numbers and practically constant, which justifies a first order reaction approximation for the mechanism rd.
the sets of discrete species d, d.r, d.rn ap and rbs, ribrbs form rapid unbroken cycles and can be averaged.
ii.
cycle averaging
the cycle averaging procedure can be applied three times:
<dig>  the cycle d, d.r is unbroken.
it is glued to the node d* whose total mass is equal to the mass of d and d.r.
<dig>  the limiting step of the cycle is klim = km <dig> <<k <dig>
<dig>  the branching reaction d → d.rn ap is replaced by d* → d.rn ap of effective constant .
<dig>  the cycle d*, d.rn ap is unbroken.
it is glued to the node d** whose total mass is equal to the mass of d and d.r and d.rn ap.
<dig>  we have  <<km <dig> hence the limiting step of the cycle is .
<dig>  the branching reaction d.rn ap → trrn ap is replaced by d** → trrn ap of effective constant .
<dig>  the cycle rbs, rib.rbs is unbroken.
it is glued to the node rbs* whose total mass is the one of rbs and of rib.rbs.
<dig>  the limiting step is km <dig> <<k <dig> rib.
<dig>  the branching reaction rib.rbs → elrib + rbs is replaced by the reaction rbs* → elrib + rbs* of effective constant k7' = k <dig>
<dig>  the branching reaction rbs → ∅ is replaced by the reaction rbs* → ∅ of effective constant .
notice that a loss of accuracy should be expected from the application of the third averaging step.
the separation of the branching and cycling reaction rate constants is not that good.
indeed, k7/km <dig> ≈  <dig>  while in theory we need k7/km <dig> <<  <dig>
iii.
identification of super-reactions
notice that after cycle averaging, a low intensity reaction d** → trn ap results, producing intermittency .
the reduced mechanism is represented in figure 4b.
the discrete/continuous partition of the species in the reduced mechanism is inherited from the initial model.
in the reduced mechanism, the reaction rbs* → rbs* + elrib is very quick, even quicker than the continuous reactions elrib → prot, prot → foldedprot, foldedprot → ∅, therefore it is a super-reaction of the type .
iv.
hybrid approximation
first order kramers-moyal expansion applied to the averaged markov jump process leads to a pdp with switching.
the continuous variables obey to the following differential equations:
the remaining three discrete components  form a markov jump process.
inside the reaction mechanism there is a rapid chain leading from trrn ap to rbs* and to elrib production which is fed by a extremely slow reaction producing trrn ap.
thus rbs* presents unfrequent bursts of activity leading to rapid production of elrib.
the increasing part of the burst does not last long, because the discrete component rbs* rapidly switches back to zero by the reaction rbs* → ∅.
thus, in the continuous variables, after a steep  increase of elrib one observes a slower decrease.
we have the case of a pdp with switching, where the steep increase event could be assimilated to a breakage .
however, the timescales of the decreasing and increasing parts of the peak are not well separated.
thus, the breakage approximation could be used only to obtain qualitative results.
to summarize, we have considered five models: the unreduced markov jump model, two averaged reduced jump markov models  and two pdp models with switching obtained from each of the previous two averaged jump markov models.
the three jump markov models are simulated using the gillespie algorithm, and the two hybrid models are simulated with the pdp algorithm.
of course, the process is markovian only between two cell cycle events that are under external command.
considering the external command, the process is semi-markovian  <cit> .
we could restore the markovian framework, by considering a cell cycle variable that has periodic autonomous dynamics and triggers the transitions between the cell cycle stages.
in order to compare the performance of the models  we have represented the total jump intensities for the first three models  as functions of time on a trajectory.
the model that demands the least computer time is the one with the smallest jump intensity.
in figure 5b, we notice a decrease of several orders of magnitudes of the total intensity from the exact ssa model to the models obtained after the second and third averaging steps.
we can not use the same method to estimate the computation time of the pdp algorithm.
indeed, although it is true generally that the computation time increases with the number of discrete transitions, the number of operations to compute the deterministic parts depends on the deterministic solver.
in order to reduce the number of operations per time point we have favored one-step schemes .
we also favored implicit stiff solvers that can function with large time steps, reducing computation time.
all the calculations were performed by using ode23s solver of matlab .
a comparison of the times to generate a  <dig> cycles long trajectories with the different methods is given in table  <dig>  this table is rather illustrative of the advantages of various approximations.
averaging allows a tremendous reduction of the execution time at least  <dig> times for weak repression and  <dig> times for strong repression .
the kramers-moyal expansion leading to pdps produce an extra 2-fold decrease of the execution time, but this is true only for weak repression and after averaging of all rapid cycles.
the pdp model m <dig> with incomplete averaging of fast cycles has a very bad performance, that can be even worse than the non-averaged ssa.
this can be explained by the strong stiff character of the deterministic dynamics with steep elrib peaks.
the time, expressed in  correspond to generating a  <dig> division cycles long trajectory  for models m1-m <dig>  the values represent log-averages for  <dig> launches.
model m <dig> is obtained from m <dig> and model m <dig> is obtained from m <dig>  via the kramers-moyal expansion.
to test the accuracies of various approximations we have computed the steady probability distribution of the protein and of the folded protein using trajectories generated by the five models.
the distribution obtained by ssa for the un-reduced model is used as the "exact" reference.
the observed errors are the consequences of less good separation between systems constants.
for instance, in figure  <dig>  k <dig> =  <dig>  and km <dig> =  <dig> , while in theory we need k <dig> <<km <dig>  however, the approximate models are qualitatively correct.
all models render correctly the bursting behavior of the system.
another advantage of a simplified model is the reduced number of parameters.
the full ssa model has  <dig> parameters.
after the first two averaging steps only  <dig> parameters remain, and after the three averaging steps only  <dig>  the pdps have the same number of parameters as the corresponding averaged markov jump processes, namely  <dig> and  <dig> parameters.
furthermore, the parameters of the simplified models are monomials of parameters of the unreduced model.
this can be used for sensitivity analysis.
after identification of the monomials that are critical for a given property the backtracked in uence of the initial parameters is given by the power of the corresponding parameter in the critical monomial.
the details of the method can be found in  <cit> .
CONCLUSIONS
we have presented, in a unified framework, various hybrid simplifications of stochastic biochemical models.
these simplifications are based on partial kramers-moyal expansions and on averaging.
the use of simplified models in stochastic studies of cellular processes has several advantages.
the first advantage of simplified models is the reduction of computational time.
we have shown that cycle averaging leads to the most drastic reduction of the computation time.
this averaging algorithm, that can be used independently of the kramers-moyal expansion, represents a general preconditioning method for both stochastic and deterministic simulations.
in stochastic studies, it reduces the number of discrete events to be simulated.
in deterministic studies, it produces less stiff systems.
the preconditioned models can be the starting point for other reduction methods.
mathematically, our simplified models are weak approximations of the fully detailed jump markov processes.
this means that all statistical properties of the trajectories of the full model including steady distributions, transition times, etc.
should be rendered without significant loss of accuracy by the simplified models.
of course, after the reduction procedure, some variables and reactions may disappear and some resulting parameters are synthetic.
because the correspondence with the original variables and parameters is known, it is always possible to go back to the details of the full model.
in particular, the identification of the critical parameters of the reduced model allows to backtrack the critical parameters of the full model.
the kramers-moyal expansion, leading to hybrid simplifications, produced only a moderate decrease of the computation time.
this limitation was partly due to our choice of low to medium size models with rather small numbers of molecules and with simple dynamics of the continuous variables.
more obvious computational advantage can be obtained for more complex models.
particularly, this could be the case for models with oscillating dynamics of the continuous variables, such as molecular clocks.
the second advantage of simplified models lies in the understanding that these bring with respect to the stochastic properties of the system.
averaging and hybrid simplifications unravel the origin of noise in multiscale biochemical systems.
noise is generated at the microscopic level of discrete variables and transferred to the mesoscopic level of the the continuous variables.
in this transfer, the topology of the reaction network plays a role, but also other properties are important such as the hierarchy of timescales of the system.
our simplification algorithm explains how and when hierarchies of cycles lead to intermittence and noise in gene networks.
in many cases, important properties such as the stationary distribution, noise amplitude, waiting times between successive bursts can be analytically calculated for the hybrid simplifications.
we have discussed several types of hybrid approximations: piecewise deterministic and hybrid diffusions with switching, with breakage, and singular switching that can be assimilated to breakage.
taken together, these results offer a rather complete panorama of various intermittence phenomena in biochemical systems.
authors' contributions
or proposed the methodology to obtain hybrid simplifications.
ac and or chose the examples and implemented the algorithms.
ac, ad and or contributed to the mathematical aspects of the methods part.
all authors drafted, read and approved the final manuscript.
supplementary material
additional file 1
additional material.
this file contains mathematical proofs of some results concerning averaging and singular switching.
it also contains a description of the methods to estimate steady probability density function for pdps.
click here for file
acknowledgements
ac acknowledges support from the region britanny council .
we thank alexander n. gorban for valuable discussions, the cycle averaging algorithm is inspired by previous joint works with him  <cit> .

 at a molecular level, the functioning of cellular processes is unavoidably stochastic. recent advances in real-time single cell imaging, micro-fluidic manipulation and synthetic biology have shown that gene expression and protein abundance in single cells are dynamic random variables submitted to significant variability. markov processes approaches to gene networks dynamics, originating from the pioneering ideas of delbrück, capture diverse features of the experimentally observed variability, such as transcription bursts, various types of steady-state distributions of rna and protein numbers, noise amplification or reduction. in the case of molecular switches, random transitions between two or several metastable states with different transcriptional activities, limit the possibility of the corresponding genetic circuits to store cellular memory and generate variability. even the simplest stochastic model, such as the production module of a protein contains tens of variables and biochemical reactions, and at least the same number of parameters. direct simulation of such networks by stochastic simulation algorithm is extremely time consuming. network reverse engineering, model parameter identification, and design principles studies suffer from the same drawbacks when the full models are attacked with traditional tools. various approximation methods were proposed, based on time discretization such as binomial, poisson or gaussian time leaping, or based on fast/slow partitions. the time complexity of exact simulation algorithms scales with the number of jumps per unit time. currently available hybrid algorithms treat fast reactions as continuous variables, which significantly reduces the number of jumps. although computationally efficient, the use of slow reactions as discrete variables, and of rapid reactions as continuous variables, does not always provide a convenient description of the hybrid stochastic process. indeed, the discrete/continuous structure of the state space of a hybrid molecular system is better expressed in terms of species components. activity of some promoters and production of the corresponding proteins occurs in bursts. in a bursting event, a steep increase of the number of transcripts is followed by relatively smoother exponential degradation. in this example, the natural candidates for continuous variables are the gene transcripts and products concentrations, while the states of the dna promoters are the discrete variables. this picture also emphasizes the possibility for noise transfer from discrete, microscopic components, to continuous, mesoscopic, and closer to the phenotype, variables. in such processes, both microscopic and mesoscopic variables are noisy. at the mesoscopic scale, the stochastic fluctuations result from the intermittency of the dynamics and do not have a simple dependence on the numbers of molecules of the continuous components. counter-intuitively, the noise amplitude could increase with the number of transcripts or proteins in a burst. stochastic intermittence or bursting occurs not only in transcription, but also in many other biochemical processes, such as action potential firing, blips in calcium signaling, possibly happloinsufficiency, etc. thus, a large class of molecular stochastic models can be approximated by stochastic hybrid processes in which some state components are discrete, while others are continuous. processes of this kind have been intensively studied in relation to applications in technology. several classes of stochastic hybrid processes were studied, among which the most important are the switching diffusion processes and the piecewise deterministic processes. for switching diffusions, the continuous state component has continuous trajectories governed by stochastic differential equations, punctuated by jumps controlled by the discrete component. for piecewise deterministic processes, the continuous component follows deterministic ordinary differential equations between two jumps. the jumps can be changes of the parameters in the differential equations, or instantaneous changes of the continuous variables, or both. piecewise deterministic processes have been recently proposed as models for various biological systems. diffusion approximations of markov processes justify approximate versions of the gillespie algorithm, such as the τ-leap method. a general approach allowing to obtain and classify various hybrid approximations that can be obtained from molecular markov jump processes is still needed and represents a first result of this paper. hybrid stochastic processes can be obtained by applying the law of large numbers to the continuous components. this will be performed here by using a partial kramers-moyal expansion of the chemical master equation. in many cases, this procedure provides only a moderate reduction of the number of stochastic jumps to be simulated. indeed, remaining discrete components can jump rapidly, in which case the simulation by a hybrid process performs ineffectively. the second result of this paper allows us to cope with frequent jumps of the discrete components by using averaging. averaging techniques are widely used in physics and chemistry to simplify models by eliminating fast microscopic variables. in our case, the microscopic variables are the discrete state components and the resulting approximations are averaged switched diffusion or averaged deterministic processes. standard averaging techniques for discrete markov chains use state aggregation. however, state aggregation algorithms are effective only for markov chains with a finite, small number of states. a more effective averaging method, that we present here, is cycle averaging in chemical reaction networks. our method exploits the formal analogy between the quasi-stationarity assumption for fast deterministic linear cycles and the stochastic averaging assumption for the same type of submodels. stochastic quasi-steady-state approximations presented in combine singular perturbations for the master equations and van kampen's Ω-expansion and are very close to our approach. however, we provide here a much more general setting and completely characterize the hybrid limits that result from this setting. the stochastic dynamics of our simplified hybrid models provide good approximations of the un-reduced chemical master equation. we show rigorously elsewhere that the simplified hybrid processes are weak approximations of the un-reduced markov process. this means that all the statistical properties of the trajectories are approximated in distribution. in particular, steady state distributions of molecular species concentrations should be accurately reproduced. also, distributions of intermittence times like the times between bursts in the production of mrna or proteins, or the commutation time for stochastic switches are expected to be similar for the reduced model and for the initial model. the structure of this paper is the following. in the methods section we obtain hybrid simplifications for pure jump markov processes and also propose a new averaging algorithm for these models. in the "results and discussion" section we present several examples of hybrid models obtained from stochastic chemical kinetics models, including a medium scale protein production model for which we demonstrate averaging. methods introductory notes on stochastic hybrid systems stochastic hybrid processes are couples of the type, x) where θ is a discrete component, x is a vector with piecewise continuous components, and t is time. the discrete component θ can be considered as a controlled markov chain, whose transition matrix may depend on the continuous variable x. the discrete state space Θ can be finite or infinite. for instance, supposing that the discrete variables are r "rare" molecular species, whose numbers never go over a small value n, then Θ ⊂ {..., n}r has at most r states. different values of θ may also correspond to various states of a single molecule. for proteins with multiple phosphorylation sites the number of states is 2r where r is the number of sites. in the case of competitive transcriptional regulation, when n transcription factors are competing on r binding sites, the number of states is at most r. there are several possible ways to specify the transitions of the discrete variables. generally, these can be given by the stochastic matrix λi, j where λi, j dt + o is the probability of a transition from the state i to the state j, between t and t + dt. however, this representation is not handy when the matrices have very large dimension. in such cases, the usual molecular markov jump process description is more appropriate. then, possible transitions are specified by the stoichiometry vectors γi. the intensity function λi represents the probability per unit time that a transition from θ to θ + γi takes place between t and t + dt. it follows that the probability for a transition between t and t + dt is λ dt + o, where λ = ∑i λi is the total intensity. the inverse of the total intensity represents the mean waiting time between two successive transitions. the transition θ → θ + λi is chosen with probability λi/λ. jumps can also occur in the continuous variables, in which case the size of the jump can be continuously distributed. there are several possibilities for the evolution of the continuous variables leading to various classes of hybrid processes. piecewise deterministic processes piecewise deterministic processes were first formalized by and found many applications in various areas ranging from industry to mathematical finance and biology. the state of a piecewise deterministic process is a couple ζ = where θ ∈ Θ denotes the discrete variable and x ∈ ℝn denotes the continuous variable. a pdp is given by three local characteristics namely: • a vector field χθ, • a jump intensity λ and • a transition measure qt such that qt is the probability distribution of the post-jump positions ∈ Θ × ℝn. related to this is the jump distribution qj giving the probability distribution of the jumps: qj = qt. such a process involves deterministic motion punctuated by jumps. the lengths of the time intervals between successive jumps are random variables called waiting times τi. on the time interval the discrete variable remains constant θ = θi, while the continuous variable evolves according to the differential equation: let ϕ) be the solution of the differential equation. we suppose that this solution exists and is unique for all t ≥ ti. the waiting time is a random variable whose distribution reads: the reader can easily obtain eq.. the probability not to jump in the time interval is - λds + o. thus p = p + o. by taking the logarithm we get log p - log p = -λds + o. summing over ds leads to eq.. the more familiar expression used in the gillespie algorithm p = exp corresponds to the particular case when λ does not depend on the continuous variable. after a deterministic evolution on the -th interval, the discrete variable and the continuous variable can change according to: where, are sampled from the jump distribution qj. the discrete variables are parameters of the differential equations describing the dynamics of the continuous variables. notice that if = the only effect of a jump on the continuous variables is a change of regime. the trajectories of the continuous variables are globally continuous. however, the velocities of the continuous variables are discontinuous at jumps of the discrete variables. we call this possibility "switching". if ≠ then continuous variables can have instantaneous jumps and their trajectories are only piecewise continuous. we call this possibility "breakage". it is possible to have both "switching" and "breakage". a pdp can be also defined by its generator, acting on functions f defined on the phase space: the adjoint of the generator is used to obtain the fokker-planck equation describing the time evolution of the probability distribution p of the process: the structure of the generic pdp corresponds to the following simulation algorithm: set the initial state condition x = x θ = θ generate a random variable u uniformly distributed in, integrate the system of differential equations between ti and ti + τi with the stopping condition f = u, generate a second uniform variable v and use it to randomly choose the jumps, change the system state ) into + ), reiterate the system from 2) with the new state until a time tmax previously defined is reached. hybrid diffusions a rather general class of hybrid diffusions has been proposed in. however, in that setting, jumps of continuous variables are commanded by the hitting of some predefined phase space sets, which is not our case. we propose a simpler and different setting, which is natural for biochemical systems. in our setting, a hybrid diffusion is defined by: • a vector field χ, • a diffusion matrix σ, • a jump intensity and a transition measure defined like for pdps. between two consecutive jumps at ti and ti + τi, the continuous variables x satisfy the itô stochastic differential equations: where wj are independent one-dimensional wiener processes. the waiting time τ is determined like for pdps. one should integrate the system: with the initial condition f = x = xi and the stopping condition f = u, where u is a random variable, independent from, uniformly distributed on. the generator for such a process is: where "." stands for the scalar product and ":" stands for the double contracted tensor product, i.e.. notice that hybrid diffusions contain pdps as the particular case when diffusion is zero. hybrid simplification of markov pure jump processes markovian stochastic chemical systems most of the existing numerical methods for stochastic chemical systems are based on the representation of the chemical reaction system as a markov pure jump process, that corresponds to gillespie's stochastic simulation algorithm. the state of a system of n chemical species a..., an is a vector x whose components are the numbers of molecules from each species. the state space of the process is e = ℤn. from now on, we shall use upper case letters x for numbers of molecules, and lower case letters x = x/ for concentrations, where is the reaction volume, typically the volume of the cell's compartment. for each reaction, we have the jump γi = βi - αi ∈ ℤn, i =..., nr. we consider equally the reverse reaction, corresponding to the jump -γi. the intensities and the transition measures for the markov jump process are: vi, v-i are the reaction rates and probabilities satisfy. here δx is the dirac measure centered on x, and x' is the post-jump position. several choices are possible for the reaction rates. the most popular is the mass action law, when rates are polynomial functions of concentrations: as for any homogeneous markov process, we can characterize this process by its generator, whose action on test functions f defined on the state space e is: the adjoint of the generator defines the chemical master equation which describes the time evolution of the probability distribution p in phase space. notice that hybrid diffusions or pdps with no continuous components are markov pure jump processes. partial kramers-moyal expansion diffusion approximations for stochastic chemical kinetics was proposed by and was rigorously justified by. this approximation is valid when the numbers of molecules of all species are much larger than one. it concerns for instance the thermodynamic limit → ∞ where is the volume of a well stirred reactor. then, the diffusion approximation applies to intermediate scales in phase space. formally, this approximation can be obtained by the kramers-moyal expansion, which is the second order taylor expansion of the master equation with respect to jumps divided by the volume; in the first order one gets the deterministic dynamics, and in the second order the langevin dynamics. in, we have used partial kramers-moyal expansions to obtain hybrid approximations of mpjp. the possibilities of this method are examined here in more generality. first, let us partition the species of the system into two sets x =, where xd represents the species in small numbers and xc the species in large numbers. note that this partition is different from reaction based partitions used in other works. all these works have as objective the reduction of simulation time. with respect to this objective, many schemes that avoid individual simulation of rapid reactions by regrouping them and by treating the corresponding channels as continuous variables in "fluid" limits, are more or less equivalent. we fix ourself also another objective, which is to find a natural representation of the simplified processes. a species based partition is more suitable for this end. let m be the set of reactions in which every reversible reaction contributes with two reactions, one for each direction. the species partition determines a partition of m into three sets. the reaction partition distinguishes between the reactions whose reactants and products are in xc (), reactions whose reactants and products are in xd () and reactions with at least a reactant or a product in xd and at least a reactant or product in xc (). we should emphasize that although reactions are rapid, some other rapid reactions can be contained in or. the second step towards the hybrid processes is to rescale the continuous variables by dividing them with the volume and obtaining xc = xc/. the xd variables will not be rescaled and will be considered discrete. their evolution is given by discrete transitions. the third step is to consider a second order taylor expansion of the master equation, for the continuous variables xc and with respect to the small parameters γi/. this expansion can be performed equivalently on the generator or on the master equation. for simplicity, we use the generator. the test functions depend now on two variables and are considered to be twice differentiable with respect to the continuous variables xc. let us consider that rates of reactions in are proportional to the volume,, i ∈. the second order taylor expansion of with respect to xc reads: where ⊗ is for tensor product kl = γk γl). the approximated generator corresponds to hybrid diffusions with drift function and diffusion matrix. as it can be easily seen, the drift and diffusion coefficients for the variables xc do not depend on the discrete variables xd. thus, switching is not present in this approximation and discrete variables can be just forgotten if not observed. furthermore, jumps vanish in the limit → ∞, implying that continuous variables have no instantaneous jumps. however, switching and breakage is present in many examples in molecular biology. how is this possible? in order to cope with the possibility of switching and breakage we need to consider that some reactions from exhibit particular properties. we can consider two types of special reactions: reactions in that are frequent enough to change the continuous variables and induce switching, reactions in with a large stoichiometric vector, that induce breakage. these reactions will be referred to as "super-reactions". the first type of super reactions are denoted by the set and the second type of super-reactions are denoted. as a consequence, we distinguish two types of hybrid approximations following which one of the two "super-reactions" are present in the mechanism. a. hybrid approximations with switching this first case of approximate process is induced by "super-reactions" with stoichiometric vectors of order one, but with rates that scale with the volume:, i ∈. reactions from induce rapid transitions of the continuous variables. the set is considered to be empty. for this approximation, even if the rates of the "super-reactions" depend on the discrete variables, we suppose that the "super-reactions" do not change the discrete variables. in other words if this condition fails, another type of approximation arises: some discrete variables change rapidly. the resulting process will be an averaged hybrid process. the evolution of the discrete variables is given by slow jumps. the component xc has rapid jumps but their amplitude is relatively small. thus, between two consecutive transitions of the discrete variables, the continuous variables can be approximated by a smooth deterministic trajectory. the continuous variables have continuous trajectory, however, discontinuous velocities. the first order of the kramers-moyal expansion corresponds to the piecewise deterministic approximation defined by the flow function the jump intensity and the jump probability where, are the projections of γi on continuous and on discrete components, respectively. the relation clearly shows that the jump intensity depends both on continuous and on discrete variables. nevertheless, only the discrete variables are changed by jumps. the second order of the kramers-moyal expansion corresponds to the hybrid diffusion approximation whereas the diffusion matrix reads: b. hybrid approximations with breakage this type of approximation is induced by the "super-reactions" with a large stoichiometric vector. we suppose also, that the set is empty. in other words, to have this approximation, we suppose that at least one reaction i in the set has a large stoichiometric vector j/ is comparable with |xj| for i ∈ and for some species j ∈ c. a reaction can change both continuous and discrete components. the main difference consists in the fact that a reaction significantly change, in one step, the concentration of the xc component. the reactions of the type do not contribute to the deterministic flow. they only appear as jumps. as before, we write down the main characteristics of the approximated pdp process, namely the flow function the jump intensity the jump probabilities for the discrete variables and the jump probability for the continuous variables if the super-reactions i ∈ have ≠ jumps occur simultaneously in the continuous and in the discrete variables. the deterministic flow is defined by the flow function χ ) which does not depend on the discrete variable xd. in this case there is no switching. in the second order of the kramers-moyal expansion, we include phase space diffusion, which contains contributions only from reactions of: if both sets, are not empty, then both types of discontinuities can appear in the trajectory of xc: • switching: discontinuity of velocity, with no discontinuity of trajectory, • breakage: discontinuity of trajectory. the first discontinuities are induced by reactions from, while the latter are induced by reactions from. c. hybrid approximations with singular switching although mathematically distinct, breakage and switching could be physically indistinguishable in certain cases. indeed, the jump of the continuous variable can result from the repeated application of one or several very fast reactions from. we can obtain processes with breakage as singular limits of processes with switching when very short lasting steep variations alternate with long lasting slower variations of the continuous variables. let us consider that there is a subset of extremely rapid super-reactions such that: where <ϵ << is a small positive parameter, and where are discrete species, substrates of the reaction i ∈. like for reactions, we consider that = for all i ∈. for each reaction i ∈ we consider that there are two subsets of reactions in: contains reactions that produce the species and contains reactions that consume the species. we define. we consider that all reactions in are rapid: first, we apply the kramers-moyal expansion for large and obtain a switched pdp. the flow function of this pdp reads: where the jump intensity of the pdp reads we show in the additional file that, in the limit ϵ → the switched pdp converges to a pdp having the following generator: where is a probability density satisfying and Φ is the solution of the differential equation. we recognize above the generator of a pdp with breakage. contrary to the previous case of breakage resulting from super-reactions, the breakage size is now continuously distributed. the switch transitions disappear in the singular limit. practical criteria for identifying small parameters and super-reactions leading to piecewise deterministic approximations the law of large numbers is applicable in the limit → ∞. certainly, in cell biology, the idea of infinite volumes should be considered with care. for this reason we will replace this condition by a set of easy to check criteria concerning relative orders of parameters of the models. these criteria will concern piecewise deterministic approximations. criteria for hybrid diffusion approximation, involving central limit theorem, will not be discussed in this paper. our criteria are relative to two reference quantities. the first reference is a large number n representing a lower bound for the numbers of molecules of continuous species. the second reference is a time τ, which is a lower bound both for the characteristic times of the deterministic dynamics of the continuous variables τc and for the average time between two successive jumps of the discrete variables τd, namely τ = min. τ can be related to kinetic parameters by methods exposed in. the applicability of the law of large numbers to continuous variables implies two conditions. the first condition is that jump sizes should be small with respect to the numbers of molecules, i.e. <<n for all i ∈. the second condition is that the number of jumps of the continuous variables on the timescale τ should be big, i.e. vi τ >> for all i ∈. the two conditions are fulfilled simultaneously in the limit → ∞ because n, vi scale with and, τ do not depend on. a similar condition must be satisfied by super-reactions of the type 1: vi >> τ- also for all i ∈. the super-reactions of the type do not satisfy the small jump condition, i.e. must be comparable to n for all i ∈. in order to have singular switching we need a new small parameter ϵ << and the following set of conditions: a) there are super-reactions of the type these are just very quick super-reactions of the type b) super-reactions of the type act only on a laps of time τϵ = ϵτ i.e. vi = - >> τ- for i ∈. reactions in that inactivate super-reactions of type could thus be as frequent as reactions in and in. c) super-reactions of the type are quick enough to produce breakages comparable to n during the time τϵ, i.e. vi >> - for i ∈. the conditions of the type vi >> τ- can be simplified if the reactions are first order with respect to discrete species. in this case vi = ki xd with xd ≈ therefore the condition is ki >> τ- all these criteria are summarized in the table see text for definitions of various reaction sets. τc, n are typical timescale and number of molecules of the continuous components, respectively. τ = min, where τd is the mean time between two successive jumps of the discrete components. vi are the average numbers of reactions per unit time, k are first order kinetic constants and γi are the stoichiometric vectors. averaging for stochastic chemical kinetics the performance of the hybrid algorithm can be very bad when the discrete mechanism contains rapid cycles which effectuate many reactions on the deterministic timescale τc. indeed, in this situation the deterministic solver is artificially sampling the interval between two discrete cycle reactions. first, this leads to unreasonable increase of the simulation time. second, the condition on the number of jumps of continuous variables is not satisfied and the hybrid approximation is not accurate. in this case, the hybrid scheme performs worse than ssa. it is therefore important to eliminate rapid cycling from the system before implementing numerical schemes for the hybrid approximations. this can be done by averaging. averaging principles are widely used for deterministic and stochastic systems. the classical averaging idea is to identify fast ergodic variables and to average the dynamics of the slow variables with respect to the quasi-stationary distribution of the fast variables. an important difficulty is to identify the right slow and fast variables. for perturbed hamiltonian dynamical systems the action-angle variables provide the natural framework for averaging. in the case of chemical kinetics, there are two kinds of slow variables that should be taken into account for averaging. the first kind are just slow components. the second kind are linear combinations of components that are conserved by the fast cycling dynamics. these new slow variables, that play a role similar to the action variables in hamiltonian systems, provide new "aggregated" or "lumped" species in the averaged system. notice that aggregation of states has been used for averaging markov chains. the aggregation of species that we propose here adapts this type of argument to the case of chemical kinetics. variable aggregation for simple deterministic reactive models has been used in relation to applications in ecology. in we proposed a general solution for simplification of reaction networks, that works for any linear mechanism with well separated constants. in this algorithm, species aggregation is systematically applied to hierarchies of rapid cycles. we show here that the same algorithm can be effectively used also in the stochastic case. averaging principles for reaction mechanisms averaging can be applied to multi-scale reaction mechanisms. this procedure leads to averaged hybrid models. in order to apply averaging, one should first identify a sub-model of discrete components, satisfying the following conditions: c1) the dynamics of the sub-model is much faster than the dynamics of continuous and of other discrete, slower modes. c2) the dynamics of the sub-model is ergodic: starting with any state, the system can reach any other state in finite time. the general procedure to obtain averaged hybrid simplifications is described in the additional file some cases are particularly interesting. the first case corresponds to fast discrete cycles producing continuous species. in this case rapid super-reactions change both discrete and continuous components is not valid). the discrete components that are affected by such reactions are fast discrete variables. in the resulting hybrid model, both the continuous dynamics and the slow discrete reaction rates should be averaged with respect to the fast discrete variables. the second case corresponds to fast, purely discrete cycles. in this case, some of the cycles of the sub-mechanism are at least as fast or faster than reactions. the other reactions of are much slower. in the resulting model, rates of the slow reactions of should be averaged with respect to the fast discrete variables. in both cases, one needs the steady probability distribution for the fast discrete sub-model. all the slow processes will be averaged with respect to this distribution. the calculation of the steady probability distribution can be easily done only when the sub-model is linear. considering linear sub-models has also another advantage. in this case, averaging and aggregation lead to a coarse-grained reaction mechanisms. for these reasons, we have developed an algorithm for linear sub-models. of course, gene network models are generally nonlinear. however, this does not mean that all the parts of such models behave nonlinearly. many sub-models, in particular monomolecular cycles, can be simplified by averaging methods that we designed for linear networks. cycle averaging for linear sub-models linear sub-models the are two types of linear reaction networks: monomolecular networks and first order networks. the structure of monomolecular reaction networks can be completely defined by a simple digraph, in which vertices correspond to chemical species ai, edges correspond to reactions ai → aj with rate constants kji > for each vertex, ai, a positive real variable ci is defined. the deterministic kinetic equation is first order reaction networks include monomolecular networks as a particular case, and are characterized by a single substrate and by reaction rates that are proportional to the concentration of the substrate. first order reaction networks can contain reactions that are not monomolecular, such as a → a + b, or a → b + c. we shall restrict ourselves to pseudo-conservative first order reactions, ie reactions that do not change the total number of molecules in a given submechanism. with such constraints, the total number of molecules in the sub-mechanism is conserved and the kinetic equations are the same as. degradation reactions can be studied in this framework by considering a special component, that collects degraded molecules. further release of the constraints is possible. for instance, the system can be opened by allowing constant production terms in eq.. these terms will change the steady state, but will not influence the relaxation times of the system. the linear sub-mechanisms can be considered as part of a nonlinear network, given fixed values of external inputs. for instance, even in systems of binary reactions, one can define pseudo-monomolecular reactions when one of the substrates of the binary reaction is not changing. this condition can be fulfilled if the substrate is in excess, for instance. the stochastic dynamics of a unique molecule in such linear reaction network is given by the probability p that the molecule is in aj at the time t. we can easily show that the master equation for p is the same as the deterministic kinetic equation. considering only one molecule does not restrict generality because when several molecules are present in a linear network, these behave independently. thus, the simplification algorithm proposed for deterministic networks can be also applied to stochastic networks. the algorithm is based on a set of operations transforming the reaction graph into an acyclic digraph without branching. for such type of graphs the eigenvectors and eigenvalues of the kinetic matrix can be straightforwardly calculated, which completely solves the problem of deterministic dynamics. we could follow precisely the same approach to simplify and then solve the stochastic dynamics. however, we argue here that applying only a few steps of the algorithm is enough for effectively reducing computational time of the ssa method. cycle averaging algorithm let us recall that the reduction method presented in uses two types of operations acting on the reaction graph. i. construction of an auxiliary network. for each node ai of a linear sub-mechanism, let us define κi as the maximal kinetic constant for reactions ai → aj: ki = maxj{kji}. for the corresponding j we use the notation ϕ: ϕ = arg maxj{kji}. an auxiliary reaction network is the set of reactions ai → aϕ with kinetic constants κi. in such a network there is no branching: if several reactions consume the same component ai, only the quickest one is kept and all the other discarded. ii. glueing cycles. rapid cycles are replaced by a single node. constants of reactions leaving these cycles are renormalized according to an averaging principle. in order to present the simplification algorithm let us use two simple examples. first, let us consider a chain of molecular reactions a → a →... am. the reaction rate constant for ai → ai+ is ki. all rate constants are considered well separated, i.e. either ki <<kj or ki >> kj for any i ≠ j. the smallest rate constant in the chain is called limiting, and denoted by klim. if 1/klim <<τc, then all molecules a are transformed into molecules am on a timescale much quicker than the deterministic time τc. we can thus ignore the chain reactions and consider that the entire mass of the chain is practically always in am. this is equivalent to considering the chain at quasi-stationarity because the steady state probability distribution of a chain is a dirac delta measure localized at the end of the chain. however, if we do not simplify chains, then simulating them by gillespie's ssa will not be computationally expensive because the mass of the chain is transferred to the end of the chain amin a number of steps that is relatively small. as a second example, let us consider the cycle c be the following sequence of mono-molecular reactions a → a →... am → a let all rate constants be well separated and the smallest one be klim like before. we add to the cycle one branching reaction; this transforms aj a component of the cycle into b a component exterior to the cycle. we consider the following distinct situations: the branching reaction is aj → b of rate constant k and k >> kj, the branching reaction is aj → b and k <<kj, the branching reaction is aj → aj + b, or the branching reaction is aj → aj+ + b of rate constant kj. in the situation the exit reaction is faster and dominates the cycling reaction aj → aj+ according to the rule for auxiliary networks in this case the cycle can be opened and the resulting multiscale dynamics is the one of a chain; we recover the previous example and in this case we can safely decide to do nothing. in the situation the exit reaction is much slower than the cycling reaction. in this case the molecules inside the cycle have rapid transformations and the mass distribution inside the cycle can be considered to reach quasi-stationary distribution. we call this cycle "unbroken". as discussed in, the relaxation time of a cycle with separated constants is the inverse of the second slowest rate constant k >> k = klim. to understand this, one should consider the two possible paths to equilibrate a cycle, one passing by the slowest step and the quicker one passing by the second slowest step: the quicker short-cuts the first one. thus, a cycle can be considered quasi-stationary if k >> 1/τc. a non-averaged fast cycle is computationally expensive in ssa, if a molecule can perform a huge number of steps along the cycle on the timescale τc. the corresponding condition involves the quasi-stationary flux and reads k = klim >> 1/τc. from a quasi-stationary cycle, the mass is lost stochastically, but slowly by the branching reaction. the intensity of the loss process can be calculated by replacing xj by its average with respect to the quasi-stationary distribution of the cycle. the average of xj is = n klim/kj, where n is the total mass inside the cycle. we obtain the average intensity = n kklim/kj. in the situations or the average intensities of the branching reactions are = n kklim/kj and = n klim, respectively. all these operations can be presented as a: simplification algorithm input: a first order reaction mechanism g with separated kinetic constants. output: a simplified first order reaction mechanism. while there are fast unbroken cycles for each cycle in g not containing reactions of the type having a sufficiently intense flux do 1: "glue" the cycle into a single node c having the total mass n; 2: replace the exit reaction of the type aj → b of rate constant k by a reaction c → b of effective constant k' = kklim/kj; 3: replace the reaction of the type aj → aj + b or rate constant k by a reaction c → c + b of effective constant k' = kklim/kj; 4: replace the reaction of the type aj → aj+ + b of rate constant kj by a reaction c → c + b of effective constant k' = klim; imposing the "glueing" condition on the cycling flux and not on the cycle relaxation time allows to iterate the cycle averaging algorithm until no more fast unbroken cycles remain. this would not be possible when adopting a relaxation time criterion to perform averaging. indeed, the relaxation time of a cycle is given by the second slowest constant. this leaves place for a counter-intuitive possibility: one can have slow cycles embedded into rapid cycles. for instance, a fast cycle whose nodes are "glued" cycles can have a slower node as the beginning of its limiting step. indeed, the internal constants of this node are necessarily more rapid than the limiting step of the large cycle, but can be slower than the second constant which gives the relaxation time of the large cycle. the slow scales are lost by glueing. in order to recover the full multi-scale dynamics, a restoration operation has been considered at the end of the algorithm presented in. in this operation, all "glued" cycles are restored without their limiting steps. thus, slower cycle sub-dynamics can be recovered. restoration is not needed for glued cycles satisfying the flux condition klim >>, because these are automatically faster than the timescale τc. our averaging method works for mechanisms with well separated constants. generalization for partially separated constants are possible. for instance, one can consider cycles containing reactions with equivalent constants, but such that for each node of the cycle, the constant of the branching reaction is much smaller than the constant of the cycling reaction. noisy networks and design rules for hierarchies of cycles as we have shown in the previous sections, intrinsic noise in biochemical networks is generated at a "microscopic level" by the discrete variables and can be observed at the "mesoscopic level" of the continuous variables either as switching, or as breaking events. thus, when there are no discrete variables, there is no intrinsic noise. also, if the switching events are much faster than the deterministic time scale, averaging principles apply and noise is not transmitted to the continuous variables: the deterministic approximation is again recovered. the only way to transmit noise by switching to the mesoscopic level is by intermittency and this needs particular combinations of slow reactions that change the values of the discrete species and frequent reactions that change the continuous species. intrinsic noise transmission from micro- to meso-scopic level results from certain combinations of low numbers, fast and rapid timescales. some general design rules relating topology to dynamic properties relative to noise can be derived from our approach. by hierarchies of cycles we understand reaction mechanisms containing cycles connected by branching reactions forming higher level cycles. the nodes of higher level cycles are "glued" cycles from the lower levels. hierarchies of cycles in discrete variables have interesting properties with respect to noise production. generally, in cycle hierarchies, effective constants of branching reactions are at least as slow or slower than the limiting step of the node from which they start. coupling cycles into hierarchies is a way to produce slower and slower reactions from initially rapid reactions and generate thus intermittency. as a possible design rule, we could state: exit reactions of the type or ) generate intermittency when the exit node is not the beginning of the limiting step in some unbroken fast cycle. indeed, unless kj is the limiting step in the cycle, one has klim/kj << then, the average intensity of the exit reaction of the type or is weak and could represent a source of intermittency in the system. this situation should be avoided for less noise in the system, or created when noise is wanted. an example of how this rule applies to the biochemistry of bacteria will be given below. a systematic investigation of the possibilities of this class of design rules will be presented elsewhere in relation to synthetic biology.  methodology to obtain hybrid simplifications we demonstrate how hybrid simplifications can be obtained from markov pure jump models for gene networks. the simplified models preserve the main stochastic properties of the initial complex models and can replace these models in simulations. the simplification procedure consists of four successive steps: i identification of discrete and continuous variables, partition of the reactions. ii cycle averaging. iii identification of super-reactions. iv construction of the hybrid simplifications. in order to justify the utility of our approach we show by examples that all types of hybrid processes that we have discussed are represented in gene networks models. to introduce the examples we employ the following notation. reactions are numbered by integers. if the i-th reaction is reversible, then is the i-th reaction in the forward direction and the i-th reaction in the reverse direction. irreversible reactions are denoted simply ri. the rate constants units are s- for monomolecular and s-1- for binary reactions. in order to obtain pseudo-monomolecular rate constants from binary reaction constants with slowly varying substrates x we have used the following formula: where na is the avogadro number and is the cell volume. for a bacterium, na ≈ 1- similarly, the reaction rates are calculated as we discuss two simple hybrid models, one with switching the second one with breakage, then two more complex models that need averaging. all our simulations were performed using matlab in a windows xp environment with a dual core intel ghz processor. hybrid model with switching: cook's model the simplest model with switching has been introduced by cook as a model for haploinsufficiency phenomena. this model can be described by the following system of biochemical reactions: in this model, g, g* represent inactive, respectively active states of the gene. in the active state, the gene produces some protein p. let g be the gene copy number, which is a conserved quantity of the dynamics g + g* = g the haploinsufficiency regime corresponds to a small value of g for the simplicity of the argument we consider g = let x = be the state vector. we notice that, the partitions of the species is the following: the species in large number xc = {p}, and the species in small number xd = {g, g*}. this partition of the species defines a partition of reactions. with the above notations we get: the deterministic timescale is τ = - the stoichiometric vectors have all lengths of order one, much smaller than n. for the parameters values used in we have k2/k >> and r is a super-reaction of the type. the reaction r satisfies vi = n k >> k = τ- which means that the law of large numbers can be applied to the continuous variable. in the first order of the kramers-moyal development we obtain a pdp approximation with switching. in the following, we will denote by x = the continuous variable and by θ = g* the discrete variable. this model is a piecewise deterministic process with state space ∈ { 1} × ℝ. the flow function χ is given by and the jump intensity λ is defined by the resulting pdp is the same as on-off systems studied in operational research. the model has been proposed as abstract model for stochastic protein production in several other places. with the pdp description of the system and using the hybrid algorithm introduced in the section methods, we draw a possible time evolution of this system. the same initial condition was used when the gillespie algorithm was implemented for the model described above. the simulation time to generate a trajectory using a pdp approximation, was seconds, while with the gillespie algorithm the simulation time was seconds in the average. the trajectories look qualitatively the same: production intervals are followed by degradation intervals, with no discontinuity in the continuous variable. in order to quantitatively test the accuracy of the pdp algorithm, we have computed the stationary distribution of the number of proteins using a piecewise deterministic simulation and compared it to the estimated distribution from gillespie trajectories. the numerical methods to estimate stationary distributions are described in the additional file the theoretical curve for the pdp model is a beta distribution. the variable x = p/ follows the beta distribution, b is the beta function. the result of various comparisons is represented in figure 1b. cook's model operates in the region k >> km which corresponds to a broken cycle in the discrete variables. if the opposite inequality is valid k <<km then cycle averaging is needed. the cycle should be glued to a point of mass = g = and the branching reaction r becomes, where. the resulting model is a birth and death model with effective birth rate and death rate constant k provided that k2/k >>, the model can be approximated by a completely deterministic process with flow given by the averaged flow function χ = -k x +. hybrid models with breakage: neuroscience and bacterium operator sites neural systems exhibit stochastic behavior. stein proposed a simple markovian model for the evolution of a neuron membrane potential. in this model, discontinuous random jumps of the potential are followed by exponential decay. we do no discuss here how to obtain this hybrid model from a microscopic pure jump model. stein's model is a phenomenological representation of very complex electric and biochemical processes. we demonstrate its properties in terms of trajectories and stationary distribution. the subthreshold behavior of the membrane potential of a neuron cell is described in this model by a hybrid markov process of continuous variable v with jumps of constant intensity λ and constant amplitude in the continuous variable. although in the original stein model there are jumps of positive and negative sign, corresponding to excitatory and inhibitory synaptic activations, here we consider only positive sign jumps. if t is the moment of jump, then v - v = a, where a > is the amplitude. between two jumps v decays according to where α is a constant. the generator of such process is: the steady probability distribution p for such a process satisfies the following delay differential equation: and p = for all x < analytical solutions are not known for this delay differential equation. in the additional file we describe the numerical scheme to calculate the steady probability distribution. the comparison between the simulated and calculated distribution is shown in figure 2c-d. in figure 2a-b we show trajectories of the system. stein's model could also apply to gene networks. generalizations of this model, allowing for continuous distributions of the jumps, could be used as an archetype of intermittent activity of a promoter. hybrid models with breakage can result as singular limits of hybrid models with switching as discussed in the methods section hybrid approximations with singular switching. the typical case is a bacterium operon. this model can be found in many places in literature. a detailed molecular example will be studied as our last example. a simpler version of this model is proposed by. it consists of the following reactions: the first reaction is zero order, all the other reactions are first order. the parameters satisfy k = ak k = bk k = ϵk we consider that b >> ϵ << from the aspect of the trajectories, the authors of hypothetize that a pdp approximation with breakage is the natural simplification and solve the corresponding stationary hybrid fokker-planck equation for the protein component. we show here that the hybrid fokker-planck equation given in can be found as an application of our approach. first, we notice that a pdp limit is applicable to the markov jump process. the species partition is xd = mrna, xc = protein. the mrna component follows a birth and death process with birth intensity k and death intensity mrna. using the master equation for birth and death processes and considering that k1/ = aϵ/ << it follows that the probability to have zero or one molecule of mrna is close to one. this justifies that mrna is a discrete species. the partition of the reactions is = {r r3}, = {r2}, = {r4}. the timescale of the continuous variables is τ = - v = k x >> τ, where x is the number of proteins, provided that x >> this condition, allowing application of the pdp approximation, is satisfied because b >>. let us show that we have singular switching, equivalent to breakage. the reaction r is a super-reaction of the type because v = k = b- the discrete variable mrna can be considered to remain zero almost all the time, except during the negligible duration of the bursts when mrna = the reaction corresponds to the transition → while the reaction corresponds to the transition → of the discrete component mrna. the intensity of the reaction r satisfies v = - thus, the mean duration of the burst is ϵτ. according to the section methods, we are in the case of a pdp with singular switching. the evolution of the unique continuous variable x is well approximated by the following hybrid generator ): which shows that b is the mean number of proteins produced in a burst. we can notice a continuous distribution of the breakage size, situation different from the stein model. the fokker-planck equation can be solved in this case. the steady distribution of the continuous variable x is the gamma distribution. first complex example: λ-phage toggle switch in this section we revisit a classical example of toggle molecular switch. the model of cro-repressor switch in λ-phage, a temperate bacteriophage of escherichia coli, was investigated by many authors. the life cycle of phages has two alternative pathways, lysogenic when the phage duplicates synchronously with the host genome and lytic when the phage produces large amounts of its own mrna. the switch between these two pathways is controlled by the level of ci protein. the lysogenic pathway corresponds to large levels of ci, while the lytic pathway corresponds to low levels of ci. a ci dimer may bind to any of the three operators: or or and or in this order. by cooperativity, or and or are almost simultaneously occupied by ci. the third site or will be occupied only when the ci concentration is high enough. a simpler model, in which or is absent, may be considered. let ci, ci and d denote the repressor, repressor dimer, and dna promoter site. the dynamics of the operator sites is described by the following reaction system: where the dci and complexes denote binding to the or or or sites, respectively binding to both sites. the other reactions concern production, degradation of molecules ci, ci2: where n is the number of proteins per mrna transcript. the state vector is x = where xd = {d, dci, dci2ci2}, xc = {ci, ci2}. note that the choice of discrete variables is dictated here, like in our first example, by the conservation law d + dci + + dci2ci = d that restricts the numbers of d, dci, dci2ci to small values. the partition of the species induces the following partition of the reactions: the interesting property of the λ-phage model is its bistability. the naive calculation to find steady states uses quasi-stationarity of the deterministic dynamics. then, we can use the equilibrium equations for reactions ri, i ∈, the steady state equation for ci and the conservation law to compute the concentration x of ci at steady-state. x = is one of the steady state, namely the lytic attractor. it also corresponds to an absorbing state of the markov process. the lytic state can be rendered non-absorbing by including, like in, a small basal production rate of ci protein in the model. however, this is not important for the illustration of our approach. if k2/km = k3/km it follows that steady states must satisfy: where. let us notice that failure of naive quasi-stationarity calculations was demonstrated for nonlinear subsystems. however, these calculations can be justified here by linear cycle averaging. bistability occurs when the quartic equation has real roots, which is the case when σ is small and α is big. more precisely, the model is bistable if. like in we have used σ = in order to ensure bistability we have chosen α = the sub-mechanism rdc contains rapid unbroken cycles that need to be averaged before any further approximation of the process. these rapid cycles correspond to rapid binding-unbinding of the dimer ci on the dna. three steps lead from the unreduced model to the averaged model: the cycle dci dci2ci is unbroken. it is glued to the node whose total mass is equal to the mass of dci and dci2ci the limiting step of the cycle is klim = km <<k ci the branching reaction dci → n x + dci is replaced by → n x + of effective constant. the reaction dci → d + ci is replaced by → d + ci with the reaction constant. after averaging, two more cycles remain in the resulting model however, the rates of the remaining reactions d + dci → dci and d + dci → are equivalent, which does not allow further application of our algorithm. furthermore, the slow reactions → d + ci and → d + ci produce intermittence and should by no means be averaged. the next approximation is a first order partial kramers-moyal expansion. the reactions in rc and the super-reaction r ∈ contribute to the deterministic dynamics of the continous species defined by the following differential equations: where x, y are the concentrations of ci and ci respectivelly. the remaining reactions induce the jump mechanism. the time evolution towards the lysogenic attractor is represented in figure 3b for the un-reduced model and in figure 3c for the pdp model. steady probability distribution close to this attractor is represented in figure 3d for all models in this study. we can notice the intermittent behavior of the ci, ci components that is well captured in the switching pdp approximations. second complex example: stochastic bursting of a repressed bacterium operon under strong repression, protein production from a bacterium operator site undergoes stochastic bursting. a stochastic model for protein production in prokaryotes has been introduced in. the behavior of the operator site under the regulation of a repressor molecule that prevents protein production has been considered in. the corresponding model is represented in the table and in figure in this model, the bacterium is considered to be in exponential growth phase, increasing size and dividing normally. cell growth is simulated by a linear increase of the volume in time. during replication, the nuclear material doubles. at fission, the nuclear material is halved and all other components are divided among daughter cells according to a binomial distribution. i. species and reaction partition some of the components are confined to small numbers by the conservation law d + d.r + d.rn ap + tr rn ap = const. two more components rbs and ribrbs have small lifetimes ≈ 2s and can not accumulate to significant levels. thus, the discrete variables are: the remaining components are in large numbers and form the continuous variables: we notice that rib, rnap and r are in relatively large numbers and practically constant, which justifies a first order reaction approximation for the mechanism rd. the sets of discrete species d, d.r, d.rn ap and rbs, ribrbs form rapid unbroken cycles and can be averaged. ii. cycle averaging the cycle averaging procedure can be applied three times: the cycle d, d.r is unbroken. it is glued to the node d* whose total mass is equal to the mass of d and d.r. the limiting step of the cycle is klim = km <<k the branching reaction d → d.rn ap is replaced by d* → d.rn ap of effective constant. the cycle d*, d.rn ap is unbroken. it is glued to the node d** whose total mass is equal to the mass of d and d.r and d.rn ap. we have <<km hence the limiting step of the cycle is. the branching reaction d.rn ap → trrn ap is replaced by d** → trrn ap of effective constant. the cycle rbs, rib.rbs is unbroken. it is glued to the node rbs* whose total mass is the one of rbs and of rib.rbs. the limiting step is km <<k rib. the branching reaction rib.rbs → elrib + rbs is replaced by the reaction rbs* → elrib + rbs* of effective constant k7' = k the branching reaction rbs → ∅ is replaced by the reaction rbs* → ∅ of effective constant. notice that a loss of accuracy should be expected from the application of the third averaging step. the separation of the branching and cycling reaction rate constants is not that good. indeed, k7/km ≈ while in theory we need k7/km << iii. identification of super-reactions notice that after cycle averaging, a low intensity reaction d** → trn ap results, producing intermittency. the reduced mechanism is represented in figure 4b. the discrete/continuous partition of the species in the reduced mechanism is inherited from the initial model. in the reduced mechanism, the reaction rbs* → rbs* + elrib is very quick, even quicker than the continuous reactions elrib → prot, prot → foldedprot, foldedprot → ∅, therefore it is a super-reaction of the type. iv. hybrid approximation first order kramers-moyal expansion applied to the averaged markov jump process leads to a pdp with switching. the continuous variables obey to the following differential equations: the remaining three discrete components form a markov jump process. inside the reaction mechanism there is a rapid chain leading from trrn ap to rbs* and to elrib production which is fed by a extremely slow reaction producing trrn ap. thus rbs* presents unfrequent bursts of activity leading to rapid production of elrib. the increasing part of the burst does not last long, because the discrete component rbs* rapidly switches back to zero by the reaction rbs* → ∅. thus, in the continuous variables, after a steep increase of elrib one observes a slower decrease. we have the case of a pdp with switching, where the steep increase event could be assimilated to a breakage. however, the timescales of the decreasing and increasing parts of the peak are not well separated. thus, the breakage approximation could be used only to obtain qualitative results. to summarize, we have considered five models: the unreduced markov jump model, two averaged reduced jump markov models and two pdp models with switching obtained from each of the previous two averaged jump markov models. the three jump markov models are simulated using the gillespie algorithm, and the two hybrid models are simulated with the pdp algorithm. of course, the process is markovian only between two cell cycle events that are under external command. considering the external command, the process is semi-markovian. we could restore the markovian framework, by considering a cell cycle variable that has periodic autonomous dynamics and triggers the transitions between the cell cycle stages. in order to compare the performance of the models we have represented the total jump intensities for the first three models as functions of time on a trajectory. the model that demands the least computer time is the one with the smallest jump intensity. in figure 5b, we notice a decrease of several orders of magnitudes of the total intensity from the exact ssa model to the models obtained after the second and third averaging steps. we can not use the same method to estimate the computation time of the pdp algorithm. indeed, although it is true generally that the computation time increases with the number of discrete transitions, the number of operations to compute the deterministic parts depends on the deterministic solver. in order to reduce the number of operations per time point we have favored one-step schemes. we also favored implicit stiff solvers that can function with large time steps, reducing computation time. all the calculations were performed by using ode23s solver of matlab. a comparison of the times to generate a cycles long trajectories with the different methods is given in table this table is rather illustrative of the advantages of various approximations. averaging allows a tremendous reduction of the execution time at least times for weak repression and times for strong repression. the kramers-moyal expansion leading to pdps produce an extra 2-fold decrease of the execution time, but this is true only for weak repression and after averaging of all rapid cycles. the pdp model m with incomplete averaging of fast cycles has a very bad performance, that can be even worse than the non-averaged ssa. this can be explained by the strong stiff character of the deterministic dynamics with steep elrib peaks. the time, expressed in correspond to generating a division cycles long trajectory for models m1-m the values represent log-averages for launches. model m is obtained from m and model m is obtained from m via the kramers-moyal expansion. to test the accuracies of various approximations we have computed the steady probability distribution of the protein and of the folded protein using trajectories generated by the five models. the distribution obtained by ssa for the un-reduced model is used as the "exact" reference. the observed errors are the consequences of less good separation between systems constants. for instance, in figure k = and km =, while in theory we need k <<km however, the approximate models are qualitatively correct. all models render correctly the bursting behavior of the system. another advantage of a simplified model is the reduced number of parameters. the full ssa model has parameters. after the first two averaging steps only parameters remain, and after the three averaging steps only the pdps have the same number of parameters as the corresponding averaged markov jump processes, namely and parameters. furthermore, the parameters of the simplified models are monomials of parameters of the unreduced model. this can be used for sensitivity analysis. after identification of the monomials that are critical for a given property the backtracked in uence of the initial parameters is given by the power of the corresponding parameter in the critical monomial. the details of the method can be found in.  we have presented, in a unified framework, various hybrid simplifications of stochastic biochemical models. these simplifications are based on partial kramers-moyal expansions and on averaging. the use of simplified models in stochastic studies of cellular processes has several advantages. the first advantage of simplified models is the reduction of computational time. we have shown that cycle averaging leads to the most drastic reduction of the computation time. this averaging algorithm, that can be used independently of the kramers-moyal expansion, represents a general preconditioning method for both stochastic and deterministic simulations. in stochastic studies, it reduces the number of discrete events to be simulated. in deterministic studies, it produces less stiff systems. the preconditioned models can be the starting point for other reduction methods. mathematically, our simplified models are weak approximations of the fully detailed jump markov processes. this means that all statistical properties of the trajectories of the full model including steady distributions, transition times, etc. should be rendered without significant loss of accuracy by the simplified models. of course, after the reduction procedure, some variables and reactions may disappear and some resulting parameters are synthetic. because the correspondence with the original variables and parameters is known, it is always possible to go back to the details of the full model. in particular, the identification of the critical parameters of the reduced model allows to backtrack the critical parameters of the full model. the kramers-moyal expansion, leading to hybrid simplifications, produced only a moderate decrease of the computation time. this limitation was partly due to our choice of low to medium size models with rather small numbers of molecules and with simple dynamics of the continuous variables. more obvious computational advantage can be obtained for more complex models. particularly, this could be the case for models with oscillating dynamics of the continuous variables, such as molecular clocks. the second advantage of simplified models lies in the understanding that these bring with respect to the stochastic properties of the system. averaging and hybrid simplifications unravel the origin of noise in multiscale biochemical systems. noise is generated at the microscopic level of discrete variables and transferred to the mesoscopic level of the the continuous variables. in this transfer, the topology of the reaction network plays a role, but also other properties are important such as the hierarchy of timescales of the system. our simplification algorithm explains how and when hierarchies of cycles lead to intermittence and noise in gene networks. in many cases, important properties such as the stationary distribution, noise amplitude, waiting times between successive bursts can be analytically calculated for the hybrid simplifications. we have discussed several types of hybrid approximations: piecewise deterministic and hybrid diffusions with switching, with breakage, and singular switching that can be assimilated to breakage. taken together, these results offer a rather complete panorama of various intermittence phenomena in biochemical systems. authors' contributions or proposed the methodology to obtain hybrid simplifications. ac and or chose the examples and implemented the algorithms. ac, ad and or contributed to the mathematical aspects of the methods part. all authors drafted, read and approved the final manuscript. supplementary material additional file 1 additional material. this file contains mathematical proofs of some results concerning averaging and singular switching. it also contains a description of the methods to estimate steady probability density function for pdps. click here for file acknowledgements ac acknowledges support from the region britanny council. we thank alexander n. gorban for valuable discussions, the cycle averaging algorithm is inspired by previous joint works with him.
sets of discrete and continuous variables are treated with different methods and are coupled together in a physically justified approach. stochastic simulation of gene networks by markov processes has important applications in molecular biology. we propose a unified framework for hybrid simplifications of markov models of multiscale stochastic gene networks dynamics. the simplified models reproduce with good accuracy the stochastic properties of the gene networks, including waiting times in intermittence phenomena, fluctuation amplitudes and stationary distributions. by averaging and variable aggregation we drastically reduce simulation time and eliminate non-critical reactions. also, answering important questions about the relation between network topology and intrinsic noise generation and propagation should be based on general mathematical results. these general results are difficult to obtain for exact models. approximate schemes reduce the computational time by reducing the number of simulated discrete events. the complexity of exact simulation algorithms scales with the number of discrete jumps to be performed. we discuss several possible hybrid simplifications, and provide algorithms to obtain them from pure jump processes. in hybrid simplifications, some components are discrete and evolve by jumps, while other components are continuous. hybrid simplifications are obtained by partial kramers-moyal expansion which is equivalent to the application of the central limit theorem to a sub-model. hybrid simplifications can be used for onion-like approaches to multi-scale biochemical systems, in which various descriptions are used at various scales.
 stochastic simulation of gene networks by markov processes has important applications in molecular biology. the complexity of exact simulation algorithms scales with the number of discrete jumps to be performed. approximate schemes reduce the computational time by reducing the number of simulated discrete events. also, answering important questions about the relation between network topology and intrinsic noise generation and propagation should be based on general mathematical results. these general results are difficult to obtain for exact models.  we propose a unified framework for hybrid simplifications of markov models of multiscale stochastic gene networks dynamics. we discuss several possible hybrid simplifications, and provide algorithms to obtain them from pure jump processes. in hybrid simplifications, some components are discrete and evolve by jumps, while other components are continuous. hybrid simplifications are obtained by partial kramers-moyal expansion which is equivalent to the application of the central limit theorem to a sub-model. by averaging and variable aggregation we drastically reduce simulation time and eliminate non-critical reactions. hybrid and averaged simplifications can be used for more effective simulation algorithms and for obtaining general design principles relating noise to topology and time scales. the simplified models reproduce with good accuracy the stochastic properties of the gene networks, including waiting times in intermittence phenomena, fluctuation amplitudes and stationary distributions. the methods are illustrated on several gene network examples.  hybrid simplifications can be used for onion-like approaches to multi-scale biochemical systems, in which various descriptions are used at various scales. sets of discrete and continuous variables are treated with different methods and are coupled together in a physically justified approach. 