SpringerOpen Newsletter

Receive periodic news and updates relating to SpringerOpen.

Open Access Research

Modeling stochasticity and variability in gene regulatory networks

David Murrugarra12*, Alan Veliz-Cuba3, Boris Aguilar4, Seda Arat12 and Reinhard Laubenbacher12

Author Affiliations

1 Department of Mathematics, Virginia Tech, Blacksburg, VA 24061-0123, USA

2 Virginia Bioinformatics Institute, Virginia Tech, Blacksburg, VA 24061-0477, USA

3 Department of Mathematics, University of Nebraska-Lincoln, Lincoln, NE 68588, USA

4 Department of Computer Science, Virginia Tech, Blacksburg, VA 24061-0123, USA

For all author emails, please log on.

EURASIP Journal on Bioinformatics and Systems Biology 2012, 2012:5  doi:10.1186/1687-4153-2012-5


The electronic version of this article is the complete one and can be found online at: http://bsb.eurasipjournals.com/content/2012/1/5


Received:14 November 2011
Accepted:6 June 2012
Published:6 June 2012

© 2012 Murrugarra et al; licensee Springer.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Modeling stochasticity in gene regulatory networks is an important and complex problem in molecular systems biology. To elucidate intrinsic noise, several modeling strategies such as the Gillespie algorithm have been used successfully. This article contributes an approach as an alternative to these classical settings. Within the discrete paradigm, where genes, proteins, and other molecular components of gene regulatory networks are modeled as discrete variables and are assigned as logical rules describing their regulation through interactions with other components. Stochasticity is modeled at the biological function level under the assumption that even if the expression levels of the input nodes of an update rule guarantee activation or degradation there is a probability that the process will not occur due to stochastic effects. This approach allows a finer analysis of discrete models and provides a natural setup for cell population simulations to study cell-to-cell variability. We applied our methods to two of the most studied regulatory networks, the outcome of lambda phage infection of bacteria and the p53-mdm2 complex.

1 Introduction

Variability at the molecular level, defined as the phenotypic differences within a genetically identical population of cells exposed to the same environmental conditions, has been observed experimentally [1-4]. Understanding mechanisms that drive variability in molecular networks is an important goal of molecular systems biology, for which mathematical modeling can be very helpful. Different modeling strategies have been used for this purpose and, depending on the level of abstraction of the mathematical models, there are several ways to introduce stochasticity. Dynamic mathematical models can be broadly divided into two classes: continuous, such as systems of differential equations (and their stochastic variants) and discrete, such as Boolean networks and their generalizations (and their stochastic variants). This article will focus on stochasticity and discrete models.

Discrete models do not require detailed information about kinetic rate constants and they tend to be more intuitive. In turn, they only provide qualitative information about the system. The most general setting is as follows. Network nodes represent genes, proteins, and other molecular components of gene regulation, while network edges describe biological interactions among network nodes that are given as logical rules representing their interactions. Time in this framework is implicit and progresses in discrete steps. More formally, let x1, ..., xn be variables, which can take values in finite sets X1, ..., Xn, respectively. Let X = X1 × ⋯ × Xn be the Cartesian product. A discrete dynamical system (DDS) in the variables x1, ..., xn is a function

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M1">View MathML</a>

where each coordinate function fi: X Xi is a function in a subset of {x1, ..., xn}. Dynamics is generated by iteration of f, and different update schemes can be used for this purpose. As an example, if Xi = {0, 1} for all i, then each fi is a Boolean rule and f is a Boolean network where all the variables are updated simultaneously. We will assume that each Xi comes with a natural total ordering of its elements (corresponding to the concentration levels of the associated molecular species). Examples of this type of dynamical system representation are Boolean networks, logical models and Petri nets [5-7].

To account for stochasticity in this setting several methods have been considered. Probabilistic Boolean networks (PBNs) [8,9] introduce stochasticity in the update functions, allowing a different update function to be used at each iteration, chosen from a probability space of such functions for each network node. For other approaches, see [10-12]. These models will be discussed in more detail in the next section. In this article we present a model type related to PBNs, with additional features. We show that this model type is natural and a useful way to simulate gene regulation as a stochastic process, and is very useful to simulate experiments with cell populations.

1.1 Modeling stochasticity in gene regulatory networks

Gene regulation processes are inherently stochastic. Accurately modeling this stochasticity is a complex and important goal in molecular system biology. Depending on the level of knowledge of the biological system and the availability of data for it one could follow different approaches. For instance, viewing a gene regulatory network as a biochemical reaction network, the Gillespie algorithm can be applied to simulate each biochemical reaction separately generating a random walk corresponding to a solution of the chemical master equation of the system [13,14]. At an even more detailed level one could introduce time delays into the Gillespie simulations to account for realistic time delays in activation or degradation such as in circadian rhythms [15-17]. At a higher level of abstraction, stochastic differential equations [18] contain a deterministic approximation of the system and an additional random white noise term. However, all these schemes require that all the kinetic rate constants to be known which could represent a strong constraint due to the difficulty of measuring kinetic parameters, limiting these approaches to small systems.

As mentioned in the introduction, discrete models are an alternative to continuous models, which do not depend on rate constants. In this setting, several approaches to introduce stochasticity have been proposed. Specially for Boolean networks, stochasticity has been introduced by flipping node states from 0 to 1 or vice versa with some flip probability [12,19-21]. However, it has been argued that this way of introducing stochasticity into the system usually leads to over-representation of noise [11]. The main criticism of this approach is that it does not take into consideration the correlation between the expression values of input nodes and the probability of flipping the expression of a node due to noise. In fact, this approach models the stochasticity at a node regardless of the susceptibility to noise of the underlying biological function [11].

Probabilistic Boolean networks [8,9,22] is another stochastic method proposed within the discrete strategy. PBNs model the choice among alternate biological functions during the iteration process, rather than modeling the stochasticity of the function failure itself. We have adopted a special case of this setting, in which every node has associated to it two functions: the function that governs its evolution over time and the identity function. If the first is chosen, then the node is updated based on its logical rule. When the identity function is chosen, then the state of the node is not updated. The key difference to a PBN is the assignment of probabilities that govern which update is chosen. In our setting, each function gets assigned two probabilities. Precisely, let xi be a variable. We assign to it a probability <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M2">View MathML</a>, which determines the likelihood that xi will be updated based on its logical rule, if this update leads to an increase/activation of the variable. Likewise, a probability <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M3">View MathML</a> determines this probability in case the variable is decreased/inhibited. The necessity for considering two different probabilities is that activation and degradation represent different biochemical processes and even if these two are encoded by the same function, their propensities in general are different. This is very similar to what is considered in differential equations modeling, where, for instance, the kinetic rate parameters for activation and for degradation/decay are, in principle, different.

Note that all these approaches only take account of intrinsic noise which is generated from small fluctuations in concentration levels, small number of reactant molecules, and fast and slow reactions. Another source of stochasticity is related to extrinsic noise such as a noisy cellular environment and temperature. For more about intrinsic vs extrinsic noise see [3,23].

2 Method

Our aim is to model stochasticity at the biological function level under the main assumption that even if the expression levels of the input nodes of an update function guarantee activation or degradation there is a probability that the process will not occur due to stochasticity, for instance, if some of the chemical reactions encoded by the update function may fail to occur. This is similar to models based on the chemical master equation. This model type introduces activation and degradation propensities. More formally, let x1, ..., xn be variables which can take values in finite sets X1, ..., Xn, respectively. Let X = X1 × ⋯ × Xn be the Cartesian product. Thus, the formal definition of a stochastic discrete dynamical system (SDDS) in the variables x1, ..., xn is a collection of n triplets

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M4">View MathML</a>

where

fi: X Xi is the update function for xi, for all i = 1, ..., n.

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M2">View MathML</a> is the activation propensity.

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M3">View MathML</a> is the degradation propensity.

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M2">View MathML</a>, <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M5">View MathML</a>.

We now proceed to study the dynamics of such systems and two specific models as illustration.

2.1 Dynamics of SDDS

Let <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M6">View MathML</a> be a SDDS and consider x X. For all i we define πi, x(xi fi(x)) and πi, x(xi xi) by

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M7">View MathML</a>

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M8">View MathML</a>

That is, if the possible future value of the i-th coordinate is larger (smaller, resp.) than the current value, then the activation (degradation) propensity determines the probability that the i-th coordinate will increase (decrease) its current value. If the i-th coordinate and its possible future value are the same, then the i-th coordinate of the system will maintain its current value with probability 1. Notice that πi, x(xi yi) = 0 for all yi ∉ {xi, fi(x)}.

The dynamics of F is given by the weighted graph X which has an edge from x X to y X if and only if yi ∈ {xi, fi(x)} for all i. The weight of an edge x y is equal to the product

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M9">View MathML</a>

By convention we omit edges with weight zero. See Additional file 1 for pseudocodes of algorithms to compute dynamics of SDDS. Software to test examples is available at http://dvd.vbi.vt.edu/adam.html webcite[24] as a web tool (choose SDDS in the model type).

Additional file 1. Additional File 1contains Supporting Material.

Format: PDF Size: 314KB Download file

This file can be viewed with: Adobe Acrobat ReaderOpen Data

Given <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M6">View MathML</a> a SDDS, it is straightforward to verify that F has the same steady states (fixed points) as the deterministic system <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M10">View MathML</a> (see Additional file 1). It is also important to note that the dynamics of F includes the different trajectories that can be generated from G using other common update mechanisms such as the synchronous and asynchronous schemes (see Additional file 1).

2.1.1 Example

Let n = 2, X = {0, 1} × {0, 1}, F = (f1, f2): X X, where Table 4 represents the regulatory rules for x1 and x2 and Table 5 represents their propensity parameters

and

Pr(01 → 10) = (.1)(.9) = .09, Pr(01 → 00) = (1 - .1)(.9) = .81

Pr(01 → 01) = (1 - .1)(1 - .9) = .09, Pr(01 → 11) = (.1)(1 - .9) = .01

Pr(10 → 10) = (1 - .2)(1 - .5) = .4, Pr(10 → 01) = (.2)(.5) = .1

Pr(10 → 00) = (.2)(1 - .5) = .1, Pr(10 → 11) = (1 - .2)(.5) = .4

Pr(11 → 11) = (1)(1 - .9) = .1, Pr(11 → 10) = (1)(.9) = .9

Pr(00 → 00) = (1)(1) = 1.

Figure 1 shows that there is a 9% chance that the system will transition from 01 to 10. Similarly, there is an 81% chance that the system will transition from 01 to 00. The latter was expected because there is a high degradation propensity for f2. Note that 00 is a fixed point, i.e., there is 100% chance of staying at this state.

thumbnailFigure 1. State Space Diagram. This diagram depicts all trajectories to follow from any given initial state of the network. The numbers next to the edges specify the transition probabilities. Note that dynamics here is not deterministic, most of the states have different trajectories to follow.

3 Applications

We illustrate the advantages of this model type by applying it to two widely studied biological systems, the regulation of the p53-mdm2 network and the control of the outcome of phage lambda infection of bacteria. These regulatory networks were selected because stochasticity plays a key role in their dynamics.

3.1 Regulation in the p53-Mdm2 network

The p53-Mdm2 network is one of the most widely studied gene regulatory networks. Abou-Jaude et al. [25] proposed a logical four-variable model to describe the dynamics of the tumor suppressor protein p53 and its negative regulator Mdm2 when DNA damage occurs. The wiring diagram of this model is represented in Figure 2, where P denotes cytoplasmic p53, nucleic p53, and the gene p53. Mc and Mn stand for cytoplasmic Mdm2 and nuclear Mdm2, respectively. DNA damage caused by ionic irradiation decreases the level of nucleic Mdm2 which enables p53 to accumulate and to remain active, playing a key role in reducing the effect of the damage. There is a negative feedback loop involving three components: p53 increases the level of cytoplasmic Mdm2 which, in turn, increases the level of nuclear Mdm2. Nucleic Mdm2 reduces p53 activity. This model also contains a positive feedback loop involving two components where p53 inhibits its negative regulator nucleic Mdm2. Note the dual role of P, as it positively regulates nucleic Mdm2 through cytoplasmic Mdm2. On the other hand, P negatively regulates nucleic Mdm2 by inhibiting Mdm2 nuclear translocation [25]. For more about the p53-Mdm2 system (see [4,25,26]).

thumbnailFigure 2. Four-variable model for the p53-Mdm2 regulatory network. P, Mc, and Mn stand for protein p53, cytoplasmic Mdm2, and nuclear Mdm2, respectively.

The dynamic behavior of the system is represented in a network of transitions called its state space (see Figure 3). This specifies the different paths to follow and the probabilities of following a specific trajectory from a given state. Dynamics here is not deterministic, i.e., most of the state vectors have different trajectories they can follow. The propensity parameters in Table 1 determine the likelihood of following certain paths. The state 0010 is a steady state, which is differentiated from the others by its oval shape.

thumbnailFigure 3. State space diagram for parameters described in Table 1. The numbers next to the edges encode the transition probabilities. The order of the variables in each vector state is P, Mc, Mn, DNA damage. Self-loops are not depicted. States with darker background comprise the cycle with DNA damage. A second cycle with a lighter shaded background corresponds to the cycle with no DNA damage. The oval shaped state is a steady state.

Table 1. Propensity probabilities for the p53-Mdm2 regulatory network

The state space for this model is specified by [0, 2] × [0, 1] × [0, 1] × [0, 1], that is, except for the first variable P which has three levels {0, 1, 2}, all other variables are Boolean. The update functions for this model are provided in Additional file 1 and also in the model repository of our web tool at http://dvd.vbi.vt.edu/adam.html webcite.

Individual cell simulations render plots similar to the ones shown in Figure 4. Each subfigure shows oscillations as long as the damage is present with a variability in the timing of damage repair. On the other hand, cell population simulations, Figure 5, exhibit damped oscillations of the expression level of p53 as the degradation propensities of the damage increases. This is correlated with the fact that, if the intensity of the damage is increased, more cells exhibit oscillations in the level of p53 which was experimentally observed in [4]. The initial state for all simulations was 0011 which represents the state when DNA damage is introduced (0010 is the steady state without perturbation).

thumbnailFigure 4. Individual cell simulations for parameters described in Table 1. Each subfigure shows oscillations as long as the damage is present. This figure shows variability in the timing of damage repair and in the period of the oscilations. Each frame was generated from a single simulation with sixty time steps. The x-axis represents discrete time steps and the y-axis the expression level. The initial state for all simulations is 0011.

thumbnailFigure 5. Cell population simulations. Each subfigure was generated from 100 simulations, each representing a single cell with sixty time steps. Starting from the top left frame to the right bottom frame the degradation propensity for DNA damage was increased by 5%, i.e. <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M11">View MathML</a> (top left), <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M12">View MathML</a> (top right), <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M13">View MathML</a> (bottom left), and <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M14">View MathML</a> (bottom right). The x-axis represents discrete time steps and the y-axis the average expression level. The initial state for all simulations was 0011. This figure shows that, if the intensity of the damage is increased more cells exhibit oscillations in the level of p53, in agreement with experimental observations [4].

Table 2. Propensity parameters for Figure 7 (top frame)

Table 3. Propensity parameters for Figure 7 (bottom frame)

To highlight the features of our approach we compare our model with the one presented in [25] in which variability has been analyzed. The main difference between these two models is in the way the simulations are performed. In [25], the transition from one state to the next is determined by parameters called "on" and "off" time delays. For instance, to transition from 2001 to 2101 it is required that <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M15">View MathML</a> which means that the "on" delay for Mc (time for activating) is less than the "off" delay (time for degrading) of the damage. Otherwise, if <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M16">View MathML</a> the system will transition from 2001 to 2000. In this article, transitions from one state to others are given as probabilities which are determined from the propensity probabilities. Therefore, the complexity of the model presented here is at the level of the wiring diagram (i.e. the number of variables) while the complexity of the model in [25] is at the level of the state space (i.e. number of possible states) which is exponential in the number of variables. Another key difference is the way DNA damage repair is modeled. In [25], a delay parameter <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M17">View MathML</a> is associated with the disappearance of the damage, and this is decreased by a certain amount τ at each iteration so that <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M18">View MathML</a> where n is the number of iterations. In order to simulate DNA damage with this approach it is required to estimate τ, n, and <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/5/mathml/M19">View MathML</a>. Within our model framework a single parameter, the degradation propensity, is used to model the damage repair which is a more natural setup.

3.2 Phage lambda infection of bacteria

Control of the outcome of phage lambda infection is one of the best understood regulatory systems [3,27,28]. Figure 6 depicts its core regulatory network that was first modeled by Thieffry and Thomas [28] using a logical approach. This model encompasses the roles of the regulatory genes CI, CRO, CII, and N. From experimental reports [3,28-30] it is known that, if the gene CI is fully expressed, all other genes are off. In the absence of CRO protein, CI is fully expressed (even in the absence of N and CII). CI is fully repressed provided that CRO is active and CII is absent.

thumbnailFigure 6. Wiring diagram for phage lambda infection model.

The dynamics of this network is a bistable switch between lysis and lysogeny, Figure 7. Lysis is the state where the phage will be replicated, killing the host. Otherwise, the network will transition to a state called lysogeny where the phage will incorporate its DNA into the bacterium and become dormant. It has been suggested [28,31] that these cell fate differences are due to spontaneous changes in the timing of individual biochemical reaction events.

thumbnailFigure 7. State space for phage lambda model. The order of variables in each vector state is CI, CRO, CII, N. The steady state 2000 represents lysogeny where CI is fully expressed while other genes are off. The cycle between 0200 and 0300 represents lysis where CRO is active and other genes are repressed. Self-loops are not depicted.

The state space for this model is specified by [0, 2] × [0, 3] × [0, 1] × [0, 1], that is, the first variable, CI, has three levels 0, 1, 2, the second variable, CRO, has four levels {0, 1, 2, 3}, and the third and fourth variables, CII and N, are Boolean. Update functions for this model are available in our supporting material, Additional file 1. This model has a steady state, 2000, and a 2-cycle involving 0200 and 0300. The steady state 2000 represents lysogeny where CI is fully expressed while the other genes are off. The cycle between 0200 and 0300 represents lysis where CRO is active and other genes are repressed.

Cell population simulations were performed to measure the cell-to-cell variability. Figure 8 was generated using the probabilities given in Tables 2 (top frame) and 3 (bottom frame). The x-axis in both subfigures represents discrete time steps while the y-axis captures the average expression level. The initial state for all simulations was 0000 which represents the state of the bacterium at the moment of phage infection. Figure 8 shows variability in developmental outcome, some of the networks transition to lysis while others transition to lysogeny. To measure how sensitive the dynamics of the network is to changes in the propensity probabilities, we have plotted the outcome of lysis-lysogeny percentages for different choices of these parameters. Figure 9 shows the variation in developmental outcome as a function of the propensity parameters of CI and CRO. Star points indicate the percentage of networks that transition to lysogeny and circle shaped points indicate the percentage of networks that end up in lysis. The bottom x-axis contains activation propensities for CI and degradation propensities for CRO while the top x-axis contains activation propensities for CRO and degradation propensities for CI. The activation and degradation propensities for CII and N were all set equal to .9. Although the probability distributions for CI and CRO are very symmetric in Figure 9, it gives a good idea of how the variability in developmental outcome will change as the propensity parameters change.

thumbnailFigure 8. Cell population simulations. Both figures were generated from 100 simulations, each representing a single cell iteration of ten time steps. Top frame for parameters in Table 2 shows 93% lysis and 7% lysogeny while bottom frame for parameters in Table 3 shows 4% lysis and 96% lysogeny. The x-axis represents discrete time steps while the y-axis shows the average expression level. The initial state for all the simulations is 0000. Solid (circle) points correspond to the average of CI (CRO), and dashed lines represent standard deviations.

thumbnailFigure 9. Variation in developmental outcome as a function of the propensity parameters. Star points indicate the percentage of networks that transition to lysogeny and circle shaped points indicate the percentage of networks that end up in lysis. Bottom axis represents the activation (and degradation) propensities for CI (CRO) in increasing order. Likewise, the top axis represents the activation (and degradation) propensities for CRO (CI) in decreasing order.

4 Conclusions

Using a discrete modeling strategy, this article introduces a framework to simulate stochasticity in gene regulatory networks at the function level, based on the general concept of PBNs. It accounts for intrinsic noise due to spontaneous differences in timing, small fluctuations in concentration levels, small numbers of reactant molecules, and fast and slow reactions. This framework was tested using two widely studied regulatory networks, the regulation of the p53-Mdm2 network and the control of phage lambda infection of bacteria. It is shown that in both of these examples the use of propensity probabilities for activation and degradation of network nodes provides a natural setup for cell population simulations to study cell-to-cell variability. The new features of this framework are the introduction of activation and degradation propensities that determine how fast or slow the discrete variables are being updated. This provides the ability to generate more realistic simulations of both single cell and cell population dynamics. In the example of the p53-Mdm2 system, one can see that individual simulations show sustained oscillations when DNA damage is present, while at the cell population level these individual oscillations average to a damped oscillation. This agrees with experimental observations [4]. In the second example, λ-phage infection of bacteria, it is observed that differences in developmental outcome due to intrinsic noise can be captured with this framework. Due to the lack of experimental data we are unable to calibrate the model so that it reproduces the correct difference in percentages due to intrinsic noise. So instead we present a plot of the difference in developmental outcome as a function of the propensity parameters.

It is worth noting that this article addresses only intrinsic noise generated from small fluctuations in concentration levels, small numbers of reactant molecules, and fast and slow reactions. Extrinsic noise is another source of stochasticity in gene regulation [3,23], and it would be interesting to see if this framework or a similar setup can be adapted to account for extrinsic stochasticity under the discrete approach. This framework also lends itself to the study of intrinsic noise and it is useful for the study of developmental robustness. For instance, one could ask what the effect of this type of noise is on the dynamics of networks controlled by biologically inspired functions.

Relating the propensity parameters to biologically meaningful information or having a systematic way for estimating them is very important. A preliminary analysis shows that it is possible to relate the propensity parameters in this framework with the propensity functions in the Gillespie algorithm under some conditions (see Additional file 1 where for a simple degradation model, the degradation propensity is correlated by a linear equation with the decay rate of the species being degraded). More precisely, in the Gillespie algorithm [13,14], if one discretizes the number of molecules of a chemical species into discrete expression levels such that within these levels the propensity functions for this species do not change significantly, then one obtains the setup of the framework presented here as a discrete model. That is, simulation within the framework presented here can be viewed as a further discretization of the Gillespie algorithm, in a setting that does not require exact knowledge of model parameters. For a similar approach see [10].

Competing interests

The authors declare that they have no competing interests.

Acknowledgements

DM and RL were partially supported by NSF grant CMMI-0908201. RL and DM thank Ilya Shmulevich for helpful suggestions. The authors thank the anonymous reviewers for many suggestions that improved the article.

References

  1. E Avigdor, M Elowitz, Functional roles for noise in genetic circuits. Nature 467, 167–173 (2010). PubMed Abstract | Publisher Full Text OpenURL

  2. M Acar, J Mettetal, A van Oudenaarden, Stochastic switching as a survival strategy in fluctuating environments. Nat Gen 40, 471–475 (2008). Publisher Full Text OpenURL

  3. F St-Pierre, D Endy, Determination of cell fate selection during phage lambda infection. PNAS 105, 20705–20710 (2008). PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  4. N Geva-Zatorsky, N Rosenfeld, S Itzkovitz, R Milo, A Sigal, E Dekel, T Yarnitzky, Y Liron, P Polak, G Lahav, U Alon, Oscillations and variability in the p53 system. Mol Syst Biol 2 (2006) (2006, 2006), . 0033, doi: 10.1038/msb4100068 OpenURL

  5. D Irons, Logical analysis of the budding yeast cell cycle. J Theor Biol 257, 543–559 (2009). PubMed Abstract | Publisher Full Text OpenURL

  6. R Thomas, R D'Ari, Biological Feedback (CRC Press, Boca Raton, 1990)

  7. C Chaouiya, E Remy, B Mossé, D Thiery, Qualitative Analysis of Regulatory Graphs: A Computational Tool Based on a Discrete Formal Framework (Lecture Notes in Control and Information Sciences, 2003) 294, pp. 830–832

  8. I Shmulevich, E Dougherty, S Kim, W Zhang, Probabilistic Boolean networks: a rule based uncertainty model for gene regulatory networks. Bioinformatics 18(2), 261–274 (2002). PubMed Abstract | Publisher Full Text OpenURL

  9. I Shmulevich, E Dougherty, Probabilistic Boolean Networks: The Modeling and Control of Gene Regulatory Networks (SIAM, Philadelphia, 2010)

  10. S Teraguchi, Y Kumagai, A Vandenbon, S Akira, D Standley, Stochastic binary modeling of cells in continuous time as an alternative to biochemical reaction equations. Phys Rev E 84(4) (2011) 062903 OpenURL

  11. A Garg, K Mohanram, A Di Cara, G De Micheli, I Xenarios, Modeling stochasticity and robustness in gene regulatory networks. Bioinformatics 15;25(12), i101–i109 (2010). PubMed Abstract | Publisher Full Text OpenURL

  12. AS Ribeiro, SA Kauffman, Noisy attractors and ergodic sets in models of gene regulatory networks. J Theor Biol 247, 743–755 (2007). PubMed Abstract | Publisher Full Text OpenURL

  13. D Gillespie, Exact stochastic simulation of coupled chemical reactions. J Phys Chem 81(25), 2340–2361 (1977). Publisher Full Text OpenURL

  14. D Gillespie, Stochastic simulation of chemical kinetics. Annu Rev Phys Chem 58, 35–55 (2007). PubMed Abstract | Publisher Full Text OpenURL

  15. D Bratsun, D Volfson, LS Tsimring, J Hasty, delay-induced stochastic oscillations gene regulation. PNAS 102(41), 14593–14598 (2005). PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  16. AS Ribeiro, Stochastic and delayed stochastic models of gene expression and regulation. Math Biosci 223(1), 1–11 (2010). PubMed Abstract | Publisher Full Text OpenURL

  17. AS Ribeiro, R Zhu, SA Kauffman, A general modeling strategy for gene regulatory networks with stochastic dynamics. J Comput Biol 13(9), 1630–1639 (2006). PubMed Abstract | Publisher Full Text OpenURL

  18. T Toulouse, P Ao, I Shmulevich, S Kauffman, Noise in a small genetic circuit that undergoes bifurcation. Complexity 11(1), 45–51 (2005). PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  19. ER Álvarez-Buylla, A Chaos, M Aldana, M Benítez, Y Cortes-Poza, C Espinosa-Soto, DA Hartasánchez, RB Lotto, D Malkin, GJ Escalera Santos, P Padilla-Longoria, Floral morphogenesis: stochastic explorations of a gene network epigenetic landscape. PLoS ONE 3(11), e3626 (doi:10, 2008), . 1371/journal.pone.0003626 PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  20. MI Davidich, S Bornholdt, Boolean network model predicts cell cycle sequence of fission yeast. PLoS ONE 3(2), e1672 (doi:10, 2008), . 1371/journal.pone.0001672 PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  21. K Willadsen, J Wiles, Robustness and state-space structure of Boolean gene regulator models. J Theor Biol 249(4), 749–765 (2008)

  22. R Layek, A Datta, R Pal, ER Dougherty, Adaptive intervention in probabilistic Boolean networks. Bioinformatics 25(16), 2042–2048 (2009). PubMed Abstract | Publisher Full Text OpenURL

  23. SS Peter, BE Michael, DS Eric, Intrinsic and extrinsic contributions to stochasticity in gene expression. PNAS 99(20), 12795–12800 (2002). PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  24. F Hinkelmann, M Brandon, B Guang, R McNeill, G Blekherman, A Veliz-Cuba, R Laubenbacher, ADAM: analysis of the dynamics of algebraic models of biological systems using computer algebra. BMC Bioinf 12, 295 (2011). BioMed Central Full Text OpenURL

  25. W Abou-Jaoudé, D Ouattara, M Kaufman, From structure to dynamics: frequency tuning in the p53-mdm2 network: I. logical approach. J Theor Biol 258(4), 561–577 (2009). PubMed Abstract | Publisher Full Text OpenURL

  26. E Batchelor, A Loewer, G Lahav, The ups and downs of p53: understanding protein dynamics in single cells. Nat Rev Cancer 9, 371–377 (2009). PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  27. M Ptashne, A Genetic Switch: Phage λ and Higher Organisms (Cell Press and Blackwell Scientific Publications, Cambridge, 1992)

  28. D Thiery, R Thomas, Dynamical behaviour of biological regulatory networks-II. Immunity control in bacteriophage lambda. Bull Math Biol 57, 277–295 (1995). PubMed Abstract OpenURL

  29. L Reichardt, D Kaiser, Control of λ repressor synthesis. PNAS 68, 2185–2189 (1971). PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  30. P Kourilsky, Lysogenization by bacteriophage lambda I. Multiple infection and the lysogenic response. Mol Gen Gen 122, 183–195 (1973). Publisher Full Text OpenURL

  31. A Arkin, J Ross, H McAdams, Stochastic kinetic analysis of developmental pathway bifurcation in Phage-infected Escherichia coli cells. Genetics 149, 1633–1648 (1998). PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL