Biochemical Network Stochastic Simulator (BioNetS): software for stochastic modeling of biochemical networks

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 efficiently and 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 solves the 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.

Figure 1

A) A screen shot of BioNetS that illustrates how the dimerization example with constant volume is entered. The tabs at the bottom of the screen allow the user to enter various options. B) A screen shot that illustrates how the dimerization example with cell growth and division is entered into BioNetS. In each screen shot, the area to the right of the reaction entry interface is a a slide-out testing panel that allows the user to enter rate constants and initial conditions, then view the results of a run from within BioNetS.

Figure 2

A) A single realization of M(t) for the discrete process. B) A single realization of M(t) produced by the chemical Langevin equations. The solid line in both panels is the result produced from Eqs. 33 and 34. The parameter values used to generate these figures were γ = 50, δ_m = 1.02, δ_d = 0.02, k₁ = 0.01, and k₂ = 0.1. The initial conditions used were M(0) = 0 and D(0) = 0 and a time step of 0.01 was used for the continuous case.

Figure 3

PDFs for the dimer concentration at various times. The staircase plots are the results for the discrete case and the continuous lines are the results for the continuous case. Panel A is for the case of constant volume and Panel B is the varying volume case. The parameter values and initial conditions are the same as in Figs. 2 and 4, for the constant volume case and varying volume case, respectively. Each PDF consists of 10, 000 realizations of the process.

Figure 4

A) Time series for the volume and monomer number for the case of cell growth and division. The parameter values used to generate these figures were γ = 50, δ_m = 1.0, δ_d = 0.0, k₁ = 0.01, k₂ = 0.1, k₃ = 0.02, and V_max = 100. The initial conditions used were M(0) = 0, D(0) = 0 and V(0) = 50. B) The monomer concentration M(t)/V(t) as a function of time. The solid line is the result from Eqs. 35–36.

Figure 5

A screen shot of the BioNetS user interface for the chemical oscillator example. The parameter values shown in the figure are the ones used to produce the results shown in Figs. 6, 7, 8. The area to the right of the reaction entry interface is a a testing panel provided by BioNetS to allow the user to enter rate constants and initial conditions, then examine the results of a run.

Figure 6

Sample paths for the repressor protein number. Panel A is the fully discrete case and Panel B is the hybrid model.

Figure 7

2-D histograms of repressor mRNA number versus repressor protein number. Red indicates regions of large frequency (i.e., where the system spends a lot of time) and blue indicates regions of low frequency. Panel A is for the discrete system and Panel B is the hybrid model. Roughly 10, 000 oscillations were used to produce both plots.

Figure 8

Power spectra for the repressor protein number. Panel A is the discrete case and B is the hybrid model.

Figure 9

Schematic of the O_RO_lac engineered promoter. The promoter fuses three operator sites, one of which (O_lac) yields repression when the LacI tetramer is bound to it, while binding of the CI dimer to O_R2 yields approximately ten-fold activation (the O_R1 site assists activation by cooperatively enhancing CI binding to O_R2 when CI is bound to O_R1). In the experimental system, CI and LacI are produced in the cell by external promoters, not shown here. The Green Fluorescent Protein (GFP) product is then used to monitor the output of the promoter by flow cytometry.

Figure 10

A) Probability densities of GFP molecules, for three versions of the O_RO_lac promoter model. Densities were generated by accumulating statistics in runs 250000 cell cycles in duration. (Solid line) Fully discrete model. (Dashed line) Fully continuous model. (Dash-dotted line) Hybrid model: the DNA binding states are represented by discrete variables, while all other species are continuous. All three methods produce virtually identical probability distributions. B) Probability density for species D₈, for the continuous model. The hybrid and fully discrete methods produce identical distributions in which D₈ spends 99.996% of the time at zero. This histogram shows that the continuous model artificially allows negative numbers (see inset time series).