Biochemical Network Stochastic Simulator (BioNetS): software for stochastic modeling of biochemical networks
 David Adalsteinsson^{1}Email author,
 David McMillen^{2} and
 Timothy C Elston^{1}
DOI: 10.1186/14712105524
© Adalsteinsson et al 2004
Received: 07 November 2003
Accepted: 08 March 2004
Published: 08 March 2004
Abstract
Background
Intrinsic fluctuations due to the stochastic nature of biochemical reactions can have large effects on the response of biochemical networks. This is particularly true for pathways that involve transcriptional regulation, where generally there are two copies of each gene and the number of messenger RNA (mRNA) molecules can be small. Therefore, there is a need for computational tools for developing and investigating stochastic models of biochemical networks.
Results
We have developed the software package Biochemical Network Stochastic Simulator (BioNetS) for efficientlyand accurately simulating stochastic models of biochemical networks. BioNetS has a graphical user interface that allows models to be entered in a straightforward manner, and allows the user to specify the type of random variable (discrete or continuous) for each chemical species in the network. The discrete variables are simulated using an efficient implementation of the Gillespie algorithm. For the continuous random variables, BioNetS constructs and numerically solvesthe appropriate chemical Langevin equations. The software package has been developed to scale efficiently with network size, thereby allowing large systems to be studied. BioNetS runs as a BioSpice agent and can be downloaded from http://www.biospice.org. BioNetS also can be run as a stand alone package. All the required files are accessible from http://x.amath.unc.edu/BioNetS.
Conclusions
We have developed BioNetS to be a reliable tool for studying the stochastic dynamics of large biochemical networks. Important features of BioNetS are its ability to handle hybrid models that consist of both continuous and discrete random variables and its ability to model cell growth and division. We have verified the accuracy and efficiency of the numerical methods by considering several test systems.
Background
Mathematical modeling of complex biological networks has a lengthy history [1–5]. In the past, the standard approach for modeling these systems has been to derive ordinary differential equations (ODEs) based on the law of mass action for the concentrations of the biochemical species involved in the network [6–16]. Experimental studies [17–19] have demonstrated, however, that stochastic effects can be significant in cellular reactions, particularly in the case of transcriptional regulation, where generally there are two copies of each gene and the number of messenger RNA (mRNA) molecules can be small. A number of recent experimental and modeling studies have addressed the role of fluctuations in gene expression [20–31]. Many modeling studies have employed the wellestablished Gillespie Monte Carlo algorithm [32] or one of its more recent variants [33, 34]. These algorithms offer an exact solution to the stochastic evolution of chemical systems, but they are computationally very expensive. A much more efficient approach is to approximate the species as continuous variables and formulate the problem in terms of stochastic differential equations (SDEs), often referred to as chemical Langevin equations [24, 28, 35]. This approximation works remarkably well for many cases, even when the number of particles involved is as small as ten, and the resulting simulations can run orders of magnitude more quickly than the discrete Monte Carlo approach. In other cases, when some or all of the particle numbers are very small, the system may need to be modeled using the discrete approach, or a hybrid method in which some species are treated discretely while others are evolved using the continuum approximation. With the increasing interest in formulating accurate models of large biochemical networks, there is a need for reliable software packages that correctly incorporate stochastic effects, yet are fast enough to simulate large interconnected sets of reacting species (as found, for example, in signaling cascades or genetic regulatory networks). We have developed the BIOchemical NETwork Stochastic Simulator, "BioNetS," to meet this need. BioNetS is capable of performing full discrete simulations using an efficient implementation of the Gillespie algorithm. It is also able to set up and solve the chemical Langevin equations, which are a good approximation to the discrete dynamics in the limit of large abundances. Finally, BioNetS can handle hybrid models in which chemical species that are present in low abundances are treated discretely, whereas those present at high abundances are handled continuously. Thus, the user can pick the simulation method that is best suited to their needs. All aspects of the software are highly optimized for efficiency.
The remainder of this manuscript is arranged in the following way. In the Implementation section, the mathematical background for the Gillespie method, chemical Langevin equations and hybrid models is presented, along with a discussion of the numerical algorithms used in BioNetS. Under Results and Discussion we provide a brief introduction to BioNetS along with several examples. The examples serve two purposes: 1) to illustrate how to use the software and 2) to verify its efficiency and accuracy. More complete documentation can be found at http://x.amath.unc.edu/BioNetS, and in the documentation included with the package.
Implementation
We first develop the mathematical methodology on which BioNetS is built. Readers interested in using BioNetS without going into its underlying structure can proceed directly to the Results and discussion section.
Discrete reactions and the gillespie algorithm
BioNetS makes use of elementary reactions (zeroth, first and second order). The following examples illustrates each type of reaction:
In the above reactions, the calligraphic letters denote a single molecule of a chemical species. The number of molecules of a particular species in the system at time t is denoted with uppercase letters (e.g., A(t), B(t), A_B(t), and V(t)). All the rate constants, γ, δ, and k_{1}k_{6}, have units of per time. Eq. 1 represents a process in which a molecule A is produced when the reaction proceeds in the forward direction and is degraded in the reverse direction. In the forward direction the reaction is zeroth order and proceeds with an average rate of γ. In the backward direction, the reaction is first order, and the average rate of degradation is δA(t). The forward reaction in Eq. 2 represents a process in which chemical species A is converted to species B. In this case A and B might represent two different conformations of the same molecule. In Eq. 2 both the forward and backward reactions are first order because the reaction rates are proportional to the respective concentrations. The forward reaction given in Eq. 3 is a second order reaction in which an A molecule and a B molecule come together to form the complex A_B. The average rate for the reaction is k_{1}A(t)B(t). The backward reaction is a first order reaction in which A_B dissociates at an average rate of k_{2}A_B(t). In Eq. 4 the forward reaction produces a molecule V. The difference between this reaction and the forward reaction in Eq. 1 is that the average rate is k_{3}V(t). This leads to exponential growth of V(t). This reaction is particularly useful if V(t) is interpreted as the cell volume. In the backward reaction, two V molecules come together and degrade one of the V molecules. The average rate for this reaction is k_{4}V(t)(V(t)  1). The V(t)  1 term arises because two of V(t) molecules must be chosen to react. This type of term also arises in reactions that produce homodimers. This reaction eventually stops the exponential growth of V. The net effect of these two reactions is to produce logistic growth. The total average reaction rate for the set of reactions given in Eqs. 1–4 is
where F_{ i }and B_{ i }are the average forward and backward rates, respectively, for the i th reaction.
For the rest of this section, we assume that the volume of the cell is not changing and only consider Eqs. 1–3. In the Examples we consider a case in which the volume is changing. If A(t), B(t) and A_B(t) are present in large numbers, then the law of mass action can be applied to derive equations that govern the concentrations [A]= A(t)/V, [B]= B(t)/V and [AB]= A_B(t)/V, where V is the cell volume. These equations are
The primed rate constants indicate that they have been appropriately scaled by the volume (i.e, k'_{3}= k_{3}V and γ' = γ/V), and, therefore, have units of either per time per concentration or concentration per time. Note that to convert to units of molar, we also have to appropriately scale the rate constants by Avagadro's number. Eqs. 6–8 represent a macroscopic description of the process, because they ignore fluctuations in the concentration that arise from the stochastic nature of chemical reactions.
In general, A(t), B(t) and A_B(t) are random variables that take on any nonnegative integer value. The Gillespie algorithm [32] can be used to generate sample paths of the process. This algorithm assumes that the random time ΔT_{ i }, between the i th and i + 1 reaction, is exponentially distributed. For the simple example given by Eqs. 1–3, the mean waiting time between reactions, which characterizes the exponential distribution, is μ_{ΔTi}= γ + δ A(t_{ i }) + k_{1}A(t_{ i }) + k_{2}B(t_{ i }) + k_{3}A(t_{ i })B(t_{ i }) + k_{4}A_B(t_{ i }), where t_{ i }is the time at which the ith reaction occurred. Therefore, t_{ i }^{+1} = t_{ i }+ ΔT_{ i }. Once the time at which the next reaction occurs has been determined, the following probabilities are used to determine which reaction occurred:
Once the reaction has been determined, the chemical species are updated accordingly. As discussed in the Numerical methods section, BioNetS uses an efficient implementation of the Gillespie algorithm [33].
Another description of discrete stochastic processes is achieved through use of the master equation that governs how the probabilities of the various random variables in the process evolve in time. Let p_{a, b,a_b}(t) = Pr [A(t) = a, B(t) = b, A_B(t) = a_b], then P_{ a,b,a_b }(t) satisfies the master equation
The master equation is the starting point for deriving various approximate schemes for describing the system [28]. In the next section, we discuss an approximate scheme that is valid in the limit of large, but finite molecule numbers. The simplest approximation scheme is achieved by considering the first moments of the process. We will use over bars to denote averaging. For example, . Eq. 15 can be used to derive equations that govern the time evolution of all the first moments. Because of the second order reaction in Eq. 3, the equations for the means are coupled to the second moments. In fact, the n th moment equations contain terms that involve the n+ l moments. Thus, there is no closure to the system. The simplest closure scheme is to assume that all moments factorize (e.g., ). This represents the macroscopic limit in which fluctuations are ignored. In this limit, we recover Eqs. 6–8 from the master equation.
The diffusion limit and the chemical langevin equations
The general form of the master equation for a system consisting of N chemical species and M reactions is
where n is a Ndimensional vector of species numbers, F_{ i }and B_{ i }are the backward and forward rates for the i th reaction, and the vectors δ_{ i }contain the stoichiometric constants for the ith reaction. For the simple model given by Eqs. 1–3, N = 3, M = 3, and p_{ n }(t) = Pr[A(t) = n_{1}, B(t) = n_{2}, and A_B(t) = n_{3}]. The forward and backward rates are F_{1} = γ, B_{1} = δn_{1}, F_{2} = k_{1}n_{1}, B_{2} = k_{2}n_{2}, F_{3} = k_{3}n_{1}n_{2}, and B_{3} = k_{4}n_{3}. The δ_{ i }vectors are the rows of the stoichiometric matrix
The (i,j) element in the above matrix represents the change in the j th chemical species when the i th reaction proceeds in the forward direction.
If the molecule numbers are large as compared to 1, then the master equation Eq. 16 can be approximated by the continuous process [28, 35]
where
This result can be derived in several ways. One method is to note that Eq. 15 represents a second order finite differencing of Eq. 18, with a grid size of 1. Another method is to make use of the shift operator
where f(n) is an arbitrary smooth function and for our purposes k is an integer. If the shift operator is used in Eq. 15, the diffusion limit is achieved when the Taylor series expansion given in Eq. 21 is truncated at j = 2.
Sample paths consistent with Eq. 18 can be generated using the following set of SDEs
where the w_{ k }(t) are independent Gaussian white noise processes. These equations are often referred to as the chemical Langevin equations. For Eqs. 1 – 3, the explicit form of the SDEs are
BioNetS generates numerical solutions to the SDEs given by Eq. 22 using either an explicit or semiimplicit Euler method. The form of these methods is
where ε = 0 for the explicit method and ε = 1 for the semiimplicit method and the Z_{ k }(t) are independent standard normal random variables. The advantage of using the chemical Langevin equations is that in the appropriate parameter regime, numerical solutions to the set of SDEs given by Eq. 22 can be generated much more efficiently than using the Gillespie algorithm. We expand upon this point in the Examples section. Higher order numerical algorithms for SDEs are available [36], but the noise structure of the chemical Langevin equations makes these schemes very cumbersome to implement. In the Examples, we verify that the Euler methods given by Eq. 26 are sufficient to produce reliable results. We note that the Δ matrix is generally sparse, and BioNetS takes advantage of this sparseness to optimize the efficiency of the two Euler methods (see Numerical Methods, below).
Hybrid schemes
It is often desirable to allow some of the chemical species to be treated as continuous random variables and some to be treated discretely. This is particularly true for the case of transcriptional regulation by transcription factors. In this situation there can be as few as one DNA/transcription factor binding site and mRNA abundances can be as small as 10 or fewer. In contrast, protein abundances can be in the thousands. The technical difficulty with implementing hybrid schemes that include both discrete and continuous random variables is that the Gillespie method requires constant transition rates between reactions. This may not be the case, if some of the chemical species are evolving continuously in time. BioNetS overcomes this problem in one of two ways.
Let N_{ d }<N be the number of discrete chemical species and M_{ d }≤ M the number of reactions that produce a change in one of the N_{ d }chemical species. The overall reaction rate at time t_{ j }for the discrete set of chemical species is
If the time step Δt for the SDEs is small enough such that
then p_{ t }is approximately the probability of a transition in Δt. In the above equation ε is a user specified tolerance. The probability of two discrete transitions in Δt is proportional to (Δt)^{2}. Choosing ε < 0.1, which means the probability of two reactions in Δt is less than 0.01, generally produces good results. However, this should be verified on a case by case basis. At each time step, BioNetS checks to verify that Ineq. 28 is satisfied for the specified ε. If so, a uniform random number R is generated and compared against p_{ t }. If R < p_{ t }, then a transition occurred and the conditional probability R/p_{ t }is used to determine which of the discrete transitions occurred. If p_{ t }> ε, then the discrete reactions determine the fastest time scale in the system. In this case the Gillespie algorithm is used to update the discrete reactions, and the random time stepΔt_{ j }is used to update the SDEs.
A pseudocode description of the above algorithm is displayed in Table 3:
Numerical methods
BioNetS generates code that is tailored to efficiently simulate biochemical reactions. The optimization techniques used by BioNetS allows the software to simulate large systems in reasonable times without requiring highend computational hardware.
Techniques used to optimize the Gillespie method are:

For the discrete variables, the program uses data structures that allow only the chemical species and reaction rates that are affected by the current reaction to be updated.

A bisection search is used to determine which reaction occurred.
The code has both an explicit and a semiimplicit solver, for simulating the chemical Langevin equations. The user specifies at runtime which method to use. By default the semiimplicit solver will be used. The semiimplicit solver uses Newton's method to solve the implicit equations, and for that the program needs to compute the Jacobian and solve a linear system at each iteration. For updating the chemical Langevin equations and hybrid models optimization techniques include:

The sparse nature of the stoichiometric matrix is used to efficiently store and per form matrix operations.

After every reaction, only the species and reaction rates affected by that reaction are updated. This can be seen in the Rates.cpp file, where all the different cases have been worked out and written for optimal execution speed.

The Jacobian is sparse, and the code takes full advantage of this fact. The program solves and factorizes the Jacobian using sparse methods. Before the code generation, BioNetS computes the entries in the Jacobian symbolically and finds a permutation that decreases the number of fillins during the LU factorization. As a result, no zero entries are saved, and the sparse structure is fully exploited. The sparse structure is then used in the LU solve. In the code, no pivots are visible, and no ifstatements are left.
Results and discussion
Dimerization
We begin with a simple system that consists of the following two reactions:
In this system, monomer molecules M are produced at an average rate γ and degraded at an average rate δ_{ m }M(t). Two monomers can then bind to form a dimer molecule D. The average forward and backward rates for the this reaction are k_{1}M(t)(M(t)  1) and k_{2}D(t), respectively. The dimers are degraded at a rate δ_{ d }. We will treat two cases. In the first case the cell volume is assumed to be constant, and in the second case the cell is allowed to grow and divide. To model cell growth, the cell volume V_{ c }is treated as a random variable V_{ c }=αV, where V is a nonnegative discrete random variable and α represents a unit of volume. The random variable V is governed by the reaction
The above reaction causes V to grow exponentially fast with an average rate of k_{3}. Note that logistic growth is produced when the backward reaction in Eq. 32 is included.
Constant volume
We start by considering the simple case in which the volume of the cell remains constant. To use BioNetS follow these steps. Copy BioNetS onto your machine, and double click to launch. Help is included as part of the program, and accessed from the Help menu. The Help document will walk you through all the steps needed to enter reactions and run the simulator.
To run BioNetS as a BioSpice agent, you need to move the source directory onto a OAAsupported system. Once there, open up the MakeOAA file and specify the locations of your oaalib folder. Then just type "make f MakeOAA" (without the quotes) to create the agent.
Cell growth and division
In this section we describe how cell growth and division can be modeled using BioNetS. We will assume that the cell is experiencing exponential growth up until the time it divides. As discussed above, the cell volume V_{ c }= αV is treated as a random variable. In this model cell division occurs when V exceeds a threshold value V_{ max }. Note that the choice of V_{ max }influences the degree of variability observed in the cell division times: cells growing from V = 1 to V_{ max }= 2 will have a large amount of variability in their division times, while those growing from V = 100 to V_{ max }= 200 will have less variable times, and those ranging from V = 1000 to V_{ max }= 2000 will be still less variable. Changing the range of V in this way requires rescaling the relationship of V to the cell volume by adjusting the value of α. When cell division occurs the volume is halved, and the proteins are randomly divided between the two cells using a binomial distribution. Only one of the daughter cells is tracked. Because second order reactions require two molecules to collide, the rate constants for these reactions should scale like k_{1} = k'_{1}/V_{ c }. We also assume that the production rate of monomers scales as γ = γ'V_{ c }. This is a reasonable assumption, because as the cell grows the transcription and translation machinery increases. These assumptions produce the following rate equations for the concentrations
The terms in Eqs. 35 and 36 that involve k_{3} arise because of dilution due to cell growth. We use the same parameter values as in the constant volume case except δ_{ m }= 1 and δ_{ d }= 0. The cell growth rate is k_{3} = 0.02 (assuming a scaling of 1 time unit to one minute, this yields an average cell division time of ln 2/k_{3} ≃ 35 minutes, typical for bacteria), the scale factor for the cell volume is α = 1 (for simplicity), and V_{ max }= 100. With these choices of parameter values, Eqs. 35 and 36 are identical with Eqs. 33 and 34, and we expect the average behavior of this system to be similar to that of the constant volume case.
A chemical oscillator
We next use BioNetS to simulate a two gene system that has been studied in the literature [37]. In this system, the protein A coded for by gene a acts as an activator for gene a and gene r, by binding to the promoter regions, P_{ a }and P_{ r }, of the respective gene. This increases the rate of mRNA_{ a }and mRNA_{ r }production by a factor α_{ a }and α_{ r }, respectively. The protein R, acts as a represser for both genes by binding to A to form the inactive complex A_R. All gene products, mRNA and protein, are actively degraded. However, the heterodimer A_R protects the R subunit from degradation. The system consists of 9 chemical species and the following 14 biochemical reactions:
An engineered promoter system
Using standard techniques in modern molecular biology, it is possible to design novel systems of promotergene pairs, such that virtually any desired regulatory network architecture may be instantiated; such networks are often called "synthetic gene networks." Recent implementations have included direct negative [22] and positive [23] feedback, a bistable switch [12], an oscillator [11], an intercellular communication system [38], and a bimodal selfactivating system [39].
The processes to be captured by the model are: transcription and degradation of mRNA strands; translation of mRNA into protein; degradation of protein; formation of protein multimers (dimers in the case of CI, tetramers in the case of LacI); LacI binding to isopropylβDthiogalactopyranoside (IPTG), a chemical inducer that reduces LacI's binding affinity for O_{ lac }; and proteinDNA binding at the O_{ R }O_{ lac }promoter's operator sites. We define the following chemical species: G, GFP; M_{ g }, mRNA coding for GFP; X, CI monomer; X_{2}, CI dimer; M_{ x }, mRNA coding for CI; D_{ x }, the arabinoseinducible pBAD promoter site producing CI; Y, LacI monomer; Y_{2}, LacI dimer; Y_{4}, LacI tetramer; I_{0}, IPTG (present in massive excess and thus taken to be constant); Y_{ I }, LacI tetramer bound to IPTG; M_{ y }, mRNA coding for LacI; and D_{ y }, the P_{ L }tetO 1 site constitutively producing LacI. In addition to these, we define species D_{0} through D_{8}, representing the various permutations of proteins bound to the three operator sites in the O_{ R }O_{ lac }promoter (see Table for a list). There are twelve combinatorial possibilities, but we eliminate three of them on the basis that CI (X_{2}) binding O_{ R }2 but notO_{ R }1 is unlikely, because of the low binding affinitity of CI for O_{ R }2 compared toO_{ R }1. Table also lists the effect on the basal rate of production when the promoter is in each state. This reflects the regulatory effect of the proteins; for example, CI bound to O_{ R }2 leads to a 10fold increase in transcription rate, while LacI bound to O_{ lac }halts transcription completely (note that we assume in the event of simultaneous binding of activator and repressor, repression "wins" and transcription is halted).
The following irreversible reactions represent the processes of transcription, translation, and degradation:
As in previous reactions, the caligraphic letters represent individual molecules of each species. We scale all times and rates by the cell division time.
Experimental measurements generally provide equilibrium rather than rate constants, and thus when writing reversible reactions we use the following notational convention: a reaction with equilibrium constant K has forward rate constant KR and backward rate constant R, where R is a scaling factor which sets the speed at which the reaction approaches equilibrium (we will consider three values of R: 1, 10, and 100). Using this notation, we represent proteinprotein binding with the following set of reactions
Finally, proteinDNA binding is given by:
In all, the system consists of 21 species, participating in 34 reactions. The reactions are entered into BioNetS using the same method described in the previous examples. We use BioNetS' ability to represent individual species as either discrete or continuous to formulate three versions of the model: fully discrete, fully continuous, and a hybrid version in which the DNA species D_{0} through D_{8} are discrete while all other species are continuous. We vary the value of R, the scaling factor for reversible reactions, and keep all other parameters fixed at the following nondimensionalized values: β_{ g }= 0.1, β_{ y }= 1, β_{ x }= 0.5, β_{ T }= 10, γ_{ mrna }= 3.5, γ_{ prot }= 0.7, K_{ y }= 0.01, K_{y 2}= 0.1, K_{ yI }= 2 × 10^{6}, K_{ x }= 0.05, K_{1} = 0.3, K_{2} = 2K_{1},K_{3} = 0.008, K_{4} = 1.4 × l0^{4}K_{3}, I_{0} = 1 × 10^{6}.
We used simulations 200 cell cycles in length to test the speed at which the three model versions ran. In each case, 200 simulations were run using a consistent set of 200 different random seeds; all runs were started with identical initial conditions. For the fully continuous and hybrid systems, the semiimplicit scheme was numerically stable and yielded consistent histograms for all time step sizes between dt = 0.001 and dt = 0.5, but the latter corresponds to just two time points per cell division cycle (recall that all times are scaled by the cell division time), and we chose instead to sample 20 points per cycle and set dt = 0.05. As shown in Table 2, the fully continuous method was always fastest, with the degree of improvement over the exact, fully discrete method depending strongly on the value of R, the scaling factor for the reversible reaction rates. For R = 1, the fully continuous method was only 1.4fold faster than the fully discrete method, but as R is increased this speed advantage increases to over 4fold at R = 10, then to over 30fold at R = 100. (Note that the speed advantage of the fully continuous over the fully discrete method increases with the abundances of the chemical species. Shifting parameters to generate higher protein numbers can yield cases in which the continuum approximation is hundreds of times faster than the discrete approach; runs not shown here.) Use of a hybrid discrete/continuous method did not, for this particular model system, offer any speed gain over the fully discrete approach; the increased time involved in computing the Jacobian for the semiimplicit method is more timeconsuming than simply simulating the reactions directly. Optimizing efficiency requires testing various potential approaches, and BioNetS makes this a simple process.
Conclusions
We have developed BioNetS to be a reliable tool for studying the stochastic dynamics of large chemical networks. The software allows the user to specify which of the chemical species in the network should be treated as discrete random variables and which can be approximated as continuous random variables. The software is highly optimized for speed and should be be able to simulate networks consisting of hundreds of chemical species. We have verified the accuracy of the numerical methods by considering several test systems (a dimerization reaction, a chemical oscillator, and an engineered promoter), each of which shows excellent agreement between the fully discrete version and the fully or partially continuous versions. Our hope is that BioNetS, by providing a simple, userfriendly interface, will allow biological experimentalists to formulate biochemical reaction models of their systems quickly and easily, ideally increasing the number of systems in which direct comparisons are available between models and experimental results. Clearly, not every possible biological system can be captured in the current version of BioNetS, and its capabilities will continue to grow in the future. We wish to encourage users, or potential users, to contact us regarding which additional features would be most helpful to them.
Availability and requirements

Project name: BlOchemical NETwork Stochastic Simulator (BioNetS)

Project home page: http://x.amath.unc.edu/BioNetS

Operating system:
* User interface: Macintosh OS X, version 10.2 or above.
* Generated source code: Ability to compile portable C++ code. Makefiles included for OS X and Linux.

Programming language: C++.

Other requirements: None.

License: BSD license.

Restrictions on use by nonacademics: None.
List of DNAbinding states in the O_{ R }O_{ lac }engineered promoter example.
Species  O _{ lac }  O _{ R } 2  O_{ R }1  Production Rate 

D _{0}        1 
D _{1}      X _{2}  1 
D _{2}    X _{2}  X _{2}  10 
D _{3}  Y _{4}      0 
D _{4}  Y _{4}    X _{2}  0 
D _{5}  Y _{4}  X _{2}  X _{2}  0 
D _{6}  Y _{I}      0 
D _{7}  Y _{I}    X _{2}  0 
D _{8}  Y _{I}  X _{2}  X _{2}  0 
Execution times for three versions of the O_{ R }O_{ lac }promoter model.
Fully continuous  Fully discrete  Hybrid  

R  Mean (s)  Std (s)  Mean (s)  Std (s)  Mean (s)  Std (s) 
1  0.26 [1]  0.01  0.36 [1.4]  0.01  0.69 [2.7]  0.01 
10  0.27 [1]  0.01  1.15 [4.3]  0.01  6.91 [27]  0.12 
100  0.27 [1]  0.01  9.11 [34]  0.28  69.43 [260]  1.13 
Table 3
• Initialize all species and rate constants 

• Compute all reaction rates 
• Loop: 
* Set μ = sum of rates for the discrete reactions 
* if (p_{ t }= μΔt > ε), use Gillespie algorithm: 
* R = a uniform random number in (0,1) 
* Set timeStep = log(R)/μ 
* Find which reaction occurred, update the species involved 
* else, use small Δt approximation: 
* R = a uniform random number in [0,1] 
* timeStep = continuousTimeStep 
* if (R <p_{ t }= μ × timeStep), discrete transition has occurred: 
• Determine which discrete transition occurred: 
• Find the first value of k for which 
• If , the forward reaction occurred, otherwise the backward reaction occurred 
* else, no discrete transition: 
• No discrete reaction occurs, update is entirely due to continuous reactions (below) 
* end if (small Δt method, determination if discrete transition occurred) 
* end if (selection of Gillespie or small Δt method for discrete reactions) 
* Update the continuous species using the Langevin equation, with step size timeStep (where timeStep is either equal to continuousTimeStep or to the step size found by the Gillespie algorithm), using a semiimplicit numerical method 
* Update any rates that have been changed by the continuous reactions and the single discrete reaction 
* Break when userdefined total simulation time is reached 
• end loop 
Declarations
Acknowledgements
This work was supported by DARPA grant F306020120579. D. Adalsteinsson acknowledges support by the Alfred P. Sloan Foundation.
Authors’ Affiliations
References
 Glass L, Kauffman S: The logical analysis of continuous, nonlinear biochemical control networks. J Theor Biol 1973, 39: 103–129.View ArticlePubMedGoogle Scholar
 Kauffman S: The largescale structure and dynamics of gene control circuits: An ensemble approach. J Theor Biol 1974, 44: 167–190.View ArticlePubMedGoogle Scholar
 Savageau M: Comparison of classical and autogenous systems of regulation in inducible operons. Nature 1974, 252: 546–549.View ArticlePubMedGoogle Scholar
 Glass L: Classification of biological networks by their qualitative dynamics. J Theor Biol 1975, 54: 85–107.View ArticlePubMedGoogle Scholar
 Tyson J, Othmer H: The dynamics of feedback control circuits in biochemical pathways. Progr Theor Biol 1978, 5: 1–60.View ArticleGoogle Scholar
 Ackers G, Johnson A, Shea M: Quantitative model for gene regulation by λ phage repressor. Proc Natl Acad Sci USA 1982, 79: 1129–1133.PubMed CentralView ArticlePubMedGoogle Scholar
 Shea M, Ackers G: The O_{R}control system of bacteriophage lambda: A physicalchemical model for gene regulation. J Mol Biol 1985, 181: 211–230.View ArticlePubMedGoogle Scholar
 Reinitz J, Vaisnys J: Theoretical and experimental analysis of the phage lambda genetic switch implies missing levels of cooperativity. J Theor Biol 1990, 145: 295–318.View ArticlePubMedGoogle Scholar
 Keller A: Model genetic circuits encoding autoregulatory transcription factors. J Theor Biol 1995, 172: 169–185. 10.1006/jtbi.1995.0014View ArticlePubMedGoogle Scholar
 Mestl T, Lemay C, Glass L: Chaos in highdimensional neural and gene networks. Physica D 1996, 98: 33–52. 10.1016/01672789(96)000863View ArticleGoogle Scholar
 Elowitz M, Leibler S: A synthetic oscillatory network of transcriptonal regulators. Nature 2000, 403: 335–338. 10.1038/35002125View ArticlePubMedGoogle Scholar
 Gardner T, Cantor C, Collins J: Construction of a genetic toggle switch in Escherichia coli . Nature 2000, 403: 339–342. 10.1038/35002131View ArticlePubMedGoogle Scholar
 Endy D, You L, Yin J, Molineux I: Computation, prediction, and experimental tests of fitness for bacteriophage T7 mutants with permuted genomes. Proc Natl Acad Sci USA 2000, 97: 5375–5380. 10.1073/pnas.090101397PubMed CentralView ArticlePubMedGoogle Scholar
 Hasty J, Isaacs F, Dolnik M, McMillen D, Collins J: Designer gene networks: Towards fundamental cellular control. Chaos 2001, 11: 207–220. 10.1063/1.1345702View ArticlePubMedGoogle Scholar
 Santillán M, Mackey M: Dynamic regulation of the tryptophan operon: A modeling study and comparison with experimental data. Proc Natl Acad Sci USA 2001, 98: 1364–1369. 10.1073/pnas.98.4.1364PubMed CentralView ArticlePubMedGoogle Scholar
 McMillen D, Kopell N, Hasty J, Collins J: Synchronizing genetic relaxation oscillators by intercell signaling. Proc Natl Acad Sci USA 2002, 99: 679–684. 10.1073/pnas.022642299PubMed CentralView ArticlePubMedGoogle Scholar
 Spudich J, Koshland DEJ: Nongenetic individuality: Chance in the single cell. Nature 1976, 262: 467–471.View ArticlePubMedGoogle Scholar
 Ross I, Browne C, Hume D: Transcription of individual genes in eukaryotic cells occurs randomly and infrequently. Immunol Cell Biol 1994, 72: 177–185.View ArticlePubMedGoogle Scholar
 Hume D: Probability in transcriptional regulation and its implications for leukocyte differentiation and inducible gene expression. Blood 2000, 96: 2323–2328.PubMedGoogle Scholar
 Arkin A, Ross J, McAdams H: Stochastic kinetic analysis of developmental pathway bifurcation in phage λ infected Escherichia coli cells. Genetics 1998, 149: 1633–1648.PubMed CentralPubMedGoogle Scholar
 Barkai N, Leibler S: Biological rhythms: Circadian clocks limited by noise. Nature 2000, 403: 267.PubMedGoogle Scholar
 Becskei A, Serrano L: Engineering stability in gene networks by autoregulation. Nature 2000, 405: 590–593. 10.1038/35014651View ArticlePubMedGoogle Scholar
 Becskei A, Séraphin B, Serrano L: Positive feedback in eukaryotic gene networks: cell differentiation by graded to binary response conversion. EMBO J 2001, 20: 2528–2535. 10.1093/emboj/20.10.2528PubMed CentralView ArticlePubMedGoogle Scholar
 Bialek W: Stability and noise in biochemical switches. In Advances in Neural Information Processing Systems Cambridge: The MIT Press 2001, 13: 103–109.Google Scholar
 Elowitz M, Levine A, Siggia E, Swain P: Stochastic gene expression in a single cell. Science 2002, 297: 1183–1186. 10.1126/science.1070919View ArticlePubMedGoogle Scholar
 Hasty J, Pradines J, Dolnik M, Collins J: Noisebased switches and amplifiers for gene expression. Proc Natl Acad Sci USA 2000, 97: 2075–2080. 10.1073/pnas.040411297PubMed CentralView ArticlePubMedGoogle Scholar
 Ko M: A stochastic model for gene induction. J Theor Biol 1991, 153: 181–194.View ArticlePubMedGoogle Scholar
 Kepler T, Elston T: Stochasticity in transcriptional regulation: origins, consequences, and mathematical representations. Biophys J 2001, 81: 3116–3136.PubMed CentralView ArticlePubMedGoogle Scholar
 McAdams H, Arkin A: Stochastic mechanisms in gene expression. Proc Natl Acad Sci USA 1997, 94: 814–819. 10.1073/pnas.94.3.814PubMed CentralView ArticlePubMedGoogle Scholar
 Ozbudak E, Thattai M, Kurtser I, Grossman A, van Oudenaarden A: Regulation of noise in the expression of a single gene. Nature Genet 2002, 31: 69–73. 10.1038/ng869View ArticlePubMedGoogle Scholar
 Thattai M, van Oudenaarden A: Intrinsic noise in gene regulatory networks. Proc Natl Acad Sci USA 2001, 98: 8614–8619. 10.1073/pnas.151588598PubMed CentralView ArticlePubMedGoogle Scholar
 Gillespie D: Exact stochastic simulation of coupled chemical reactions. J Phys Chem 1977, 81: 2340–2361.View ArticleGoogle Scholar
 Gibson M, Bruck J: Efficient exact stochastic simulation of chemical systems with many species and many channels. J Phys Chem 2000, 104: 1876–1889. 10.1021/jp993732qView ArticleGoogle Scholar
 Gillespie D: Approximate accelerated stochastic simulation of chemically reacting systems. J Phys Chem 2001, 115: 1716–1733. 10.1063/1.1378322View ArticleGoogle Scholar
 Gardiner C: Handbook of stochastic methods: For physics, chemistry and the natural sciences Berlin: Springer Verlag 1996.Google Scholar
 Kloeden P, Platen E: Numerical solution of stochastic differential equations Berlin: Springer Verlag 1992.View ArticleGoogle Scholar
 Vilar J, Kueh H, Barkai N, Leibler S: Mechanisms of noiseresistance in genetic oscillators. Proc Natl Acad Sci USA 2002, 99: 5988–5992. 10.1073/pnas.092133899PubMed CentralView ArticlePubMedGoogle Scholar
 Weiss R, Knight T: Engineered communications for microbial robotics. In DNA6: Sixth International Meeting on DNA Based Computers: June 13–17 2000; Leiden, Netherlands (Edited by: Condon A, Rozenberg G). Berlin: Springer Verlag 2001.Google Scholar
 Isaacs F, Hasty J, Cantor C, Collins J: Prediction and measurement of an autoregulatory genetic module. Proc Natl Acad Sci USA 2003, 100: 7714–7719. 10.1073/pnas.1332628100PubMed CentralView ArticlePubMedGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article: verbatim copying and redistribution of this article are permitted in all media for any purpose, provided this notice is preserved along with the article's original URL.