Stochastic simulation is an attractive approach for simulating gene-regulatory networks since it readily accounts for noise in transcriptional regulation. Random transitions among the discrete states of operators that control the transcription rates and the finite number of protein and mRNA molecules involved in biochemical reactions are potential noise sources. We consider here a stochastic model for a mutual repressor system: a two gene network in which each protein represses the transcription of the other gene. This type of system has been engineered in E. coli and is referred to as a genetic toggle switch. Our stochastic model of the system predicts the steady-state distribution of protein abundances which may be compared to single cell fluorescence measurements of intercellular variability in protein expression levels. The computational time required to simulate such a system with a stochastic model of even modest complexity is significant.

We present a computer-assisted “equation-free” approach that accelerates computation without requiring an explicit reduction of the underlying stochastic model. Instead the computational speed-up stems from the execution of “intelligently initialized” short bursts of the stochastic simulator with suitable processing of the results. This approach exploits the separation of time scales in the system; synthesis and degradation of new proteins and transcripts typically occur on a slower time scale than processes that change the chemical state of proteins. We extend standard bifurcation analysis (numerical continuation), typically used with systems of ODEs, to our stochastic model to determine the regions of parameter space in which bistability of the genetic toggle switch occurs and compute the mean first passage time between stable steady states using only short-time stochastic simulations. The accuracy of our methods, tested by direct comparison with long-time stochastic simulations, is excellent.

Additionally, we have developed a “variable-free” mode of analysis illustrating that eigenvectors of the weighted graph Laplacian defined on results of stochastic simulation bursts define appropriate “automated" reaction coordinates and may therefore be used to study systems where the correct observables are unknown a priori. We present lifting and restriction procedures for translating between physical system variables and these automated reaction coordinates, enabling all of equation-free analysis previously described to be performed in these new coordinates. This type of equation-free analysis shows promise for the computation of features of the long-time, coarse grained behavior of complex stochastic models of gene regulatory networks, circumventing the need for long Monte Carlo simulations. This latter part of the work is in collaboration with Prof. Coifman and Dr. Nadler at Yale University.