Hybrid stochastic simplifications for multiscale gene networks

Background 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. Results 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 [1-3] 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. Conclusion Hybrid simplifications can be used for onion-like (multi-layered) 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.

We get the following system for all 0 ≤ i ≤ n-1. We let X i = g(ip) and get the following linear system for all i ≤ n − 1 (where j is chosen to be the positive integer closest to The truncated trajectory of a pure jump process is a series of times (jump instants) t0 = T M IN < t 2 < . . . < t N = T M AX and a corresponding series of states x 0 , . . . , x N . Considering that the process is constant between t i and t i+1 , the steady probability density function can be estimated as a histogram of the values x i with weights t i+1 − t i .
For piece-wise deterministic processes, the assumption that the process is constant between two successive jump instants is no longer satisfied. In this case, the interval between two jump instants is sampled by the deterministic solver. Considering that the time step δt = t i+1 − t i is small enough, a linear approximation is suitable for the histogram, corresponding to the trapezium rule: x i should be replaced by (x i + x i+1 )/2. Sometimes, we choose to resample the deterministic pieces of the trajectories, using uniform discretisation and spline interpolation.

Justification of averaging results
In order to obtain averaged hybrid approximations from the master equation we generally need two parameters. A large parameter V (representing volume) is used to rescale large molecule numbers as in Kramers-Moyal or Ω expansions. A small parameter (representing ratio of fast to slow timescales) is used to separate fast dynamics within cycles and slow dynamics between cycles. In the most general case, Kramers-Moyal and averaging can be performed in one step starting with the master equation and developing with respect to 1/V and . This is a singular perturbation problem because the dynamics at > 0, 1/V > 0 is not an uniform approximation of the dynamics at = 0, 1/V = 0. Arbitrarily large differences between the two dynamics occur if we wait long enough. For instance, weakly coupled cycles have a different long time dynamics and reach a different steady state than totally uncoupled cycles.
We can consider the following cases: Case 1 : Fast discrete cycles producing continuous species. There are fast super-reactions of the type 1 that act on discrete species γ D i = 0 for some i ∈ S 1 . In this case there is only one small parameter η = 1/V .
Case 2 : Fast discrete cycles. In this case there are some fast transitions between discrete species. The small parameter is .
Case 3 : A combination of the above two cases. There are two small parameters, η and .
We discuss only the first two cases.
Case 1 The conditions defining the Case 1 mean that some rapid cycles change in the same time discrete and continuous variables. The discrete species change within a finite set of values, and remain in small numbers. The continuous species can be produced continuously, in large numbers. This case is similar to the QSSA SPA-Ω expansion in [1]. Our discrete species are the quasi-stationary species in [1]. The difference is that here we consider rapid cycling (implying production and consumption) of discrete species, while in [1] only rapid consumption is considered (meaning that γ D i < 0 for i ∈ S 1 . To simplify, we consider that all reactions in R DC are of this type (R DC = S 1 and γ D i = 0 for all reactions in R DC ).
Let x c = ηX c be the concentrations of the continuous species. Let X 2 D be the discrete species that are affected by super-reactions of the type 1 and X 1 D the remaining discrete species.
We start from the following master equation: We develop p as a power series in η: This series will provide the large time solution of the master equation. Using a Taylor expansion of (1) we obtain at various orders: This means that at X 1 D , x c , t fixed, p 0 is the steady-state distribution of the process defined by the fast reactions R DC . Now, we consider that this process is ergodic, meaning that starting from any state we can reach any other state in finite time. The reaction mechanism R DC should be made of one or several interconnected cycles. This condition ensures the uniqueness of the steady state distribution and we can write that: where ρ(X 2 D ) is the unique steady state distribution of the fast process and ψ(X 1 D , x c , t) is the time-dependent distribution of the remaining variables.
At order η 0 where By using (2) and by summing (3) with respect to all possible values of X 2 D we obtain: Eq.(4) shows that the zeroth order approximation is an averaged PDP.
Case 2 We do not need to distinguish here between discrete and continuous species. The reader could imagine the situation when there are only discrete species. We consider that there are fast reactions R 2 and slow reactions R 1 . We consider that V i = −1 v i for i ∈ R 2 . Like for Case 2 we consider two types of species. X 2 D are the species that are affected by fast reactions and X 1 D the remaining species.
This case can be treated exactly in the same way as the preceding case. We start from the following master equation: With the same assumption on the ergodicity of the process defined by fast reactions, we have In the zero-th approximation, ρ(X 2 D ) is the steady state distribution of the fast mechanism and ψ(X 1 D , t) represents the dynamics of an averaged mechanism: When the mechanism R 2 is made of a single cycle A 2 1 → A 2 2 → ...A 2 m → A 2 1 , the steady state distribution can be easily calculated. This distribution is parametrised by the value of the total mass of the cycle

Justification of the singular switching results
Consider the following PDP with singular switching. The flow function is: and the jump intensity is We shall obtain the limit process by two methods. The first method employs a scaling argument, and the second method employs the expansion of the hybrid Fokker-Planck equation.
The state of the system is (x, X D , X D ) where x are continuous species, X D are discrete species, and X D are discrete species substrates of of rapid reactions of the type S 3 that are responsible for the breakage. To simplify the calculations we consider that variables X D take only two values X D ∈ {0, 1}. If X D = 1 species x are rapidly produced by reactions S 3 , and if X D = 0 production stops. Once in the state X D = 1, the system rapidly switches back to X D = 0 by fast reactions R − D .
to the state (X D + γ D i , X D = 1, x), where it stays only for a very short time. In this state, the continuous variables are submitted to the fast dynamics: during the short random time τ that satisfies then jumps to the state X D + γ D i + γ D j , X D = 0 by a reaction j ∈ R − D .
Let Φ(s; x, X D ) be the solution of the equation Passing to short timescales s = t/ , we can easily show that the variation of the continuous variable starting from x and during the time τ is: where the random times satisfies Keeping the random variable ∆x as jump (breakage) of the continuous variable and noticing from (7) that τ → 0 (almost surely), we obtain the following generator of the limit hybrid process: is a probability density satisfying ∞ 0 ρ ij (s)ds = 1. From the generator we can obtain the hybrid Fokker-Planck equation (using the adjoint of the generator).
In the adjoint, we should replace trajectories starting from x to trajectories arriving in x, thus changing s into −s and renormalizing the density to take into account the phase space volume transformation. Thus, the hybrid Fokker-Planck equation reads: Fokker-Planck equation argument To simplify formulae, we consider with no loss of generality that The hybrid Fokker-Planck equation for the singular switching process reads: where Consider that p(X D , X D , x, t) = p 0 (X D , X D , x, t) + p 1 (X D , X D , x, t) + . . .
Then we Taylor expand the Fokker-Planck equations (10) and we obtain: At order −1 p 0 (X D , 1, x, t) = 0 From the first of the Eqs.11 we obtain The linear differential equation (12) has the solution: By replacing (13) in the second of the Eqs.11 we obtain again Eq.(9).