A Dominated Coupling From The Past algorithm for the stochastic simulation of networks of biochemical reactions

<p>Abstract</p> <p>Background</p> <p>In recent years, stochastic descriptions of biochemical reactions based on the Master Equation (ME) have become widespread. These are especially relevant for models involving gene regulation. Gillespie’s Stochastic Simulation Algorithm (SSA) is the most widely used method for the numerical evaluation of these models. The SSA produces exact samples from the distribution of the ME for finite times. However, if the stationary distribution is of interest, the SSA provides no information about convergence or how long the algorithm needs to be run to sample from the stationary distribution with given accuracy. </p> <p>Results</p> <p>We present a proof and numerical characterization of a Perfect Sampling algorithm for the ME of networks of biochemical reactions prevalent in gene regulation and enzymatic catalysis. Our algorithm combines the SSA with Dominated Coupling From The Past (DCFTP) techniques to provide guaranteed sampling from the stationary distribution. The resulting DCFTP-SSA is applicable to networks of reactions with uni-molecular stoichiometries and sub-linear, (anti-) monotone propensity functions. We showcase its applicability studying steady-state properties of stochastic regulatory networks of relevance in synthetic and systems biology.</p> <p>Conclusion</p> <p>The DCFTP-SSA provides an extension to Gillespie’s SSA with guaranteed sampling from the stationary solution of the ME for a broad class of stochastic biochemical networks.</p

Stochastic analysis of nonlinear dynamics and feedback control for gene regulatory networks with applications to synthetic biology

The focus of the thesis is the investigation of the generalized repressilator model (repressing genes ordered in a ring structure). Using nonlinear bifurcation analysis stable and quasi-stable periodic orbits in this genetic network are characterized and a design for a switchable and controllable genetic oscillator is proposed. The oscillator operates around a quasi-stable periodic orbit using the classical engineering idea of read-out based control. Previous genetic oscillators have been designed around stable periodic orbits, however we explore the possibility of quasi-stable periodic orbit expecting better controllability. The ring topology of the generalized repressilator model has spatio-temporal symmetries that can be understood as propagating perturbations in discrete lattices. Network topology is a universal cross-discipline transferable concept and based on it analytical conditions for the emergence of stable and quasi-stable periodic orbits are derived. Also the length and distribution of quasi-stable oscillations are obtained. The findings suggest that long-lived transient dynamics due to feedback loops can dominate gene network dynamics. Taking the stochastic nature of gene expression into account a master equation for the generalized repressilator is derived. The stochasticity is shown to influence the onset of bifurcations and quality of oscillations. Internal noise is shown to have an overall stabilizing effect on the oscillating transients emerging from the quasi-stable periodic orbits. The insights from the read-out based control scheme for the genetic oscillator lead us to the idea to implement an algorithmic controller, which would direct any genetic circuit to a desired state. The algorithm operates model-free, i.e. in principle it is applicable to any genetic network and the input information is a data matrix of measured time series from the network dynamics. The application areas for readout-based control in genetic networks range from classical tissue engineering to stem cells specification, whenever a quantitatively and temporarily targeted intervention is required

A Dominated Coupling From The Past algorithm for the stochastic simulation of networks of biochemical reactions-2

He DCFTP-SSA (+), as in Fig. 2 (); the SSA with = 1000(○), as in Fig. 2 (); the SSA started from the two modes (*), as in Fig. 2 (); the SSA started from uniform initial conditions (∇), as in Fig. 2 (); and the SSA uniformly sampled from a long run (□), as in Fig. 2 (). For each scheme, we produced = 100, 316, 1000, 3162 and 10000 samples to show how the error improves as the number of samples increases. The DCFTP-SSA converges to the stationary distribution at the expected rate, whereas the approximate estimates obtained using the SSA level off in a similar manner as in Fig. 1. () The distribution of coalescence times for the DCFTP-SSA for this network is bimodal with a very long tail for the second mode, indicating the likelihood of long coalescence times. The data presented corresponds to 6000 runs.<p><b>Copyright information:</b></p><p>Taken from "A Dominated Coupling From The Past algorithm for the stochastic simulation of networks of biochemical reactions"</p><p>http://www.biomedcentral.com/1752-0509/2/42</p><p>BMC Systems Biology 2008;2():42-42.</p><p>Published online 8 May 2008</p><p>PMCID:PMC2529164.</p><p></p

A Dominated Coupling From The Past algorithm for the stochastic simulation of networks of biochemical reactions-3

. The escape time from the mode located away from the origin is 2 × 10. () The mean first passage time from the origin to the other mode is 3 × 10.<p><b>Copyright information:</b></p><p>Taken from "A Dominated Coupling From The Past algorithm for the stochastic simulation of networks of biochemical reactions"</p><p>http://www.biomedcentral.com/1752-0509/2/42</p><p>BMC Systems Biology 2008;2():42-42.</p><p>Published online 8 May 2008</p><p>PMCID:PMC2529164.</p><p></p

A Dominated Coupling From The Past algorithm for the stochastic simulation of networks of biochemical reactions-5

5 sampled with the DCFTP-SSA (+) and with the standard SSA with stopping times = 2(○), 4(□), 6(◇). For this simple ME, the limiting value of the Euclidean error of the finite-time SSA is , where = 1 - exp(-) and () is the modified Bessel function of the first kind []. This means that SSA simulations that are run for a time will converge to a systematic sampling error, indicated by the dotted lines. This source of error is eliminated when using the DCFTP-SSA, which shows no flooring for and the expected scaling with the number of Monte Carlo samples []. The guarantees provided by the DCFTP-SSA come at a modest computational cost, which is comparable to that of long SSA runs. () The distribution of coalescence times for the DCFTP-SSA is relatively symmetric and concentrated around the mean with a rapid decay for long times. The data presented corresponds to 6000 runs. This distribution reflects the benign structure of the unimodal stationary distribution of this particular ME, which makes long coalescence times unlikely.<p><b>Copyright information:</b></p><p>Taken from "A Dominated Coupling From The Past algorithm for the stochastic simulation of networks of biochemical reactions"</p><p>http://www.biomedcentral.com/1752-0509/2/42</p><p>BMC Systems Biology 2008;2():42-42.</p><p>Published online 8 May 2008</p><p>PMCID:PMC2529164.</p><p></p

A Dominated Coupling From The Past algorithm for the stochastic simulation of networks of biochemical reactions-0

5 sampled with the DCFTP-SSA (+) and with the standard SSA with stopping times = 2(○), 4(□), 6(◇). For this simple ME, the limiting value of the Euclidean error of the finite-time SSA is , where = 1 - exp(-) and () is the modified Bessel function of the first kind []. This means that SSA simulations that are run for a time will converge to a systematic sampling error, indicated by the dotted lines. This source of error is eliminated when using the DCFTP-SSA, which shows no flooring for and the expected scaling with the number of Monte Carlo samples []. The guarantees provided by the DCFTP-SSA come at a modest computational cost, which is comparable to that of long SSA runs. () The distribution of coalescence times for the DCFTP-SSA is relatively symmetric and concentrated around the mean with a rapid decay for long times. The data presented corresponds to 6000 runs. This distribution reflects the benign structure of the unimodal stationary distribution of this particular ME, which makes long coalescence times unlikely.<p><b>Copyright information:</b></p><p>Taken from "A Dominated Coupling From The Past algorithm for the stochastic simulation of networks of biochemical reactions"</p><p>http://www.biomedcentral.com/1752-0509/2/42</p><p>BMC Systems Biology 2008;2():42-42.</p><p>Published online 8 May 2008</p><p>PMCID:PMC2529164.</p><p></p