Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Open Access

Peer-reviewed

Research Article

Are Quasi-Steady-State Approximated Models Suitable for Quantifying Intrinsic Noise Accurately?

Abstract

Large gene regulatory networks (GRN) are often modeled with quasi-steady-state approximation (QSSA) to reduce the huge computational time required for intrinsic noise quantification using Gillespie stochastic simulation algorithm (SSA). However, the question still remains whether the stochastic QSSA model measures the intrinsic noise as accurately as the SSA performed for a detailed mechanistic model or not? To address this issue, we have constructed mechanistic and QSSA models for few frequently observed GRNs exhibiting switching behavior and performed stochastic simulations with them. Our results strongly suggest that the performance of a stochastic QSSA model in comparison to SSA performed for a mechanistic model critically relies on the absolute values of the mRNA and protein half-lives involved in the corresponding GRN. The extent of accuracy level achieved by the stochastic QSSA model calculations will depend on the level of bursting frequency generated due to the absolute value of the half-life of either mRNA or protein or for both the species. For the GRNs considered, the stochastic QSSA quantifies the intrinsic noise at the protein level with greater accuracy and for larger combinations of half-life values of mRNA and protein, whereas in case of mRNA the satisfactory accuracy level can only be reached for limited combinations of absolute values of half-lives. Further, we have clearly demonstrated that the abundance levels of mRNA and protein hardly matter for such comparison between QSSA and mechanistic models. Based on our findings, we conclude that QSSA model can be a good choice for evaluating intrinsic noise for other GRNs as well, provided we make a rational choice based on experimental half-life values available in literature.

Citation: Sengupta D, Kar S (2015) Are Quasi-Steady-State Approximated Models Suitable for Quantifying Intrinsic Noise Accurately? PLoS ONE 10(9): e0136668. https://doi.org/10.1371/journal.pone.0136668

Editor: Attila Csikász-Nagy, Fondazione Edmund Mach, Research and Innovation Centre, ITALY

Received: April 5, 2015; Accepted: August 6, 2015; Published: September 1, 2015

Copyright: © 2015 Sengupta, Kar. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited

Data Availability: The data underlying the results described in our manuscript are available in the body of the manuscript as well as in all the related supporting information uploaded with the manuscript.

Funding: This work was supported by 13IRCCSG008, http://www.ircc.iitb.ac.in, Industrial research and consultancy centre, IIT Bombay, SK; and EMR/2014/000500, http://www.dst.gov.in, Department of science and technology, India, SK.

Competing interests: The authors have declared that no competing interests exist.

Introduction

Most of the events in and around the biological cell are inherently noisy. The sources of noise in the biological systems are diverse and can be mainly classified into two different classes [1–5], namely intrinsic and extrinsic noise. Recently, there were a number of experimental findings [6–11], which revealed that intrinsic and extrinsic noise levels in a particular cell can dictate its ultimate cell fate to a great extent. The intrinsic noise is mainly classified as the molecular noise that is present in any cell type because of low copy numbers of mRNA molecules corresponding to any particular protein [1,3,12,13]. Further, it is well established [14,15] that the coefficient of variation (CV) of the distribution of a particular protein not only depends on the average number of mRNA and protein molecules but also relies on the absolute values of the half-lives of the corresponding mRNA and protein molecules. Gillespie’s stochastic simulation algorithm (SSA) [16] provides fundamentally the most accurate and correct description of the intrinsic noise for a biological network if all the terms in the corresponding model are based on mass action kinetics. This means as the biological network grows in size and complexities, it will be computationally highly time consuming to simulate every reaction event using SSA. To overcome this challenge, different kinds of efforts had been made to avoid continuously simulating the reactions that happen in a faster time scale in case of SSA [15,17–19]. Quasi-steady-state approximation (QSSA) is one of such methods where one reduces the detailed mechanistic model of a complex biological network depending on the time scales of the reactions in the concerned system and as a result effectively cut short the number of molecular species and reactions to simulate stochastically [15,19]. In this context, one should keep in mind that well-separated timescales in a deterministic model do not guarantee a reduced QSSA model to be valid and there are instances where additional conditions needed to be fulfilled [20]. In last few decades, quasi-steady-state approximated models of different biological systems were proved to be quite successful to furnish all the deterministic dynamical features of the corresponding biological systems under different biological conditions and even in some cases the stochastic behavior as well [15,19,20]. There were evidences that some of these QSSA models were faithfully reproducing the intrinsic noise for simple biological systems [15,19,20]. In this context, we want to investigate what happens if we use the QSSA approach to quantify the intrinsic noise for relatively complex biological networks? Under what conditions the QSSA approach will quantify the intrinsic noise as accurately as the SSA performed using the mechanistic version of the model and when the result will differ significantly?

The organization of the paper is as follows. In the first section, we develop a mechanistic and the corresponding QSSA model for a simple positive feedback module. We employ Gillespie’s stochastic algorithm [16] to simulate the QSSA model as well as the mechanistic model to identify the conditions where both the models will give identical measure of intrinsic noise for the corresponding motif and where the results will start to differ from each other significantly. First, we compare the stochastic results obtained from the mechanistic and the corresponding QSSA models to investigate the role played by the mRNA and protein half-lives while keeping the protein and mRNA abundance levels fixed. In this section we changed the half-life values following two different approaches. In approach 1, we change the ratio of the mRNA and protein half-lives by keeping the mRNA half-life fixed [15] and in approach 2, we vary the absolute values of mRNA and protein half-lives by keeping the half-life ratio fixed. Next, we have done similar comparison between the two models as a function of number of molecules of protein (mRNA number of molecules remained fixed) and number of molecules of mRNA (protein number of molecules remained fixed) by following approach 1 and approach 2. We further perform similar studies with extended GRNs to show that the conclusions made from a simple positive feedback motif holds good even for relatively bigger networks as well. Finally, in the discussion section, we summarize our results and discuss the measures to be taken before using a QSSA model to faithfully quantify intrinsic noise for relatively complex biological networks.

Creating mechanistic and QSSA models for a simple positive feedback module

We have constructed a toy model of a simple positive feedback motif [21–26] that is frequently observed in many biological regulations (Fig 1A, detailed network shown in the right panel) [25,26] in terms of mass action kinetics (mechanistic model, right panel, Table 1). We further made some appropriate steady state approximations on the variables except the total protein (X) and mRNA (MP) and obtained the corresponding QSSA model (left panel, Table 1) for the given GRN.

thumbnail
Download:
Fig 1. Positive feedback motif under consideration and the corresponding bifurcation diagrams.

(A) QSSA positive feedback motif (left panel) and the detailed mechanistic scheme (right panel) are given where solid and dotted lines represent biochemical reactions and catalytic effects respectively. G and Ga denote inactive and active states of gene, Mp denotes mRNA, P denotes protein, P2 is the dimer of protein and X is the total protein. The protein (P) molecules form dimers P2. P2 binds to the promoter region of the inactive gene (G) and activates G to Ga. km and kp are the degradation rates of MP and P respectively. (B) Bifurcation diagrams (left panel, QSSA model and right panel, mechanistic model) of total protein (X) are plotted as a function of J0. The solid lines represent stable steady states and dotted lines represent unstable steady states; J0 is the basal rate of mRNA (MP) synthesis and acts as the bifurcation parameter. In both the models deterministic mean of X (molecules) = 457.9 at J0 = 3 min-1. (C) Bifurcation diagrams (left panel, QSSA model and right panel, mechanistic model) of mRNA (MP) are plotted as a function of J0. In both the models deterministic mean of MP (molecules) = 185.2 at J0 = 3 min-1.

https://doi.org/10.1371/journal.pone.0136668.g001

thumbnail
Download:
Table 1. Equations of the QSSA and Mechanistic models.

https://doi.org/10.1371/journal.pone.0136668.t001

The parameter values for the model concerned are given in Table 2. Here we have kept km>>J3, as it was shown by Thomas et al. [20] that in order for the reduced QSSA model to be accurate one needs to fulfill this additional condition to be satisfied on top of the appropriate time scale separation. The details of the steady state approximations are provided in the S1 Text. The identical bifurcation diagrams shown in Fig 1B (left panel, QSSA model and right panel, mechanistic model) for the total protein (X) as a function of J0 (the basal rate of the MP synthesis), signify that the two models provided in Table 1 are deterministically similar (the corresponding bifurcation diagrams of MP as a function of J0 are given in Fig 1C). This clearly shows that both the deterministic models will lead to identical steady states deterministically under different values of J0 and sets the stage to investigate whether the QSSA model shown in Table 1 (left panel) can quantify the intrinsic noise as efficiently and accurately as the detailed mechanistic model (right panel, Table 1) or not.

thumbnail
Download:
Table 2. Parameters of the QSSA and Mechanistic models.

https://doi.org/10.1371/journal.pone.0136668.t002

Results

Comparison of the stochastic distributions at different values of J0

Once we created the deterministically similar models, we employed the Gillespie’s simulation algorithm [16] at J0 = 3 min-1 (deterministic means of X = 457.9 molecules and MP = 185.2 molecules) to quantify the intrinsic noise from both the models (The reaction scheme followed for Gillespie simulation for both the models are provided in the S1 Text). We have performed the stochastic simulation at J0 = 3 min-1 where we do have only one stable steady state and measured the steady state distribution of the total protein (X) to compare the QSSA (left panel, Fig 2A) and mechanistic (right panel, Fig 2A) models (The corresponding time courses for the two situations are given in S1A Fig). The result looks extremely promising from the perspective of the QSSA model as the statistical quantities (mean = 457.2 and standard deviation = 42.98 in number of molecules, Coefficient of variation (CV) = 9.4%) obtained from the total protein distribution resembles quite well with the statistical values (mean = 462.3 and standard deviation = 43.54 in number of molecules, CV = 9.42%) obtained from the corresponding total protein distribution of the mechanistic model.

thumbnail
Download:
Fig 2. Stochastic simulation based on Gillespie Algorithm.

(A) Steady state distributions of X (total protein) are plotted at J0 = 3 min-1 for the QSSA (left panel, mean = 457.2, standard deviation = 42.98, CV = 9.4%) and the mechanistic models (right panel, mean = 462.3, standard deviation = 43.54, CV = 9.42%). (B) Steady state distributions of MP (mRNA) are plotted at J0 = 3 min-1 for the QSSA (left panel, mean = 185.2, standard deviation = 15.05, CV = 8.13%) and the mechanistic models (right panel, mean = 185.6, standard deviation = 27.26, CV = 14.68%). All the values (except CVs) are in number of molecules.

https://doi.org/10.1371/journal.pone.0136668.g002

This evidently shows that under the chosen parameter regime the QSSA model does a good job to accurately quantify the molecular noise at the protein level in comparison to the mechanistic model. This kind of Gillespie runs is also performed with different sets of random number sequences (one of them shown in S1C Fig) and we got similar results. Additionally, we make sure that the parameter set used for the steady state approximation to obtain the QSSA model is a sensible choice. Since, if one performs Gillespie simulation with a wrong choice of parameters (for example, ka = 8E-03 min-1 and kd = 5E-03 min-1 instead of ka = 8 min-1 and kd = 5 min-1 as given in Table 2, Fig 2A), the stochastic results (S1 Table) for the mechanistic model (S2 Fig, right panel) will be totally different than that of the stochastic simulation result (S1 Table) of the corresponding QSSA model (S2 Fig, left panel). Stochastic results for the QSSA model are not at all affected by the inaccuracies made in the steady state approximations; where as mechanistic model can sense such changes quite adequately.

Although we got similar total protein distributions for the two models (Fig 2A), the corresponding mRNA (MP) steady state distributions that are given in Fig 2B (left panel, QSSA model and right panel, mechanistic model) appear to be quite different in nature (The corresponding time courses for the two situations are given in S1B Fig). We observe that the statistical quantities such as standard deviation and the CV at the mRNA level (mean = 185.2, standard deviation = 15.05, CV = 8.13%, for QSSA model and mean = 185.6, standard deviation = 27.26, CV = 14.68%, for mechanistic model) are quite different for the two models under the same parametric domain even though the average abundance level of MP is same for both the cases. It clearly indicates that the protein (P) in the given biological network is not sensing the molecular noise present at the mRNA level adequately in case of the mechanistic model. As a consequence the intrinsic noise calculated in the form of CV of the total protein distribution for QSSA and mechanistic models seems to be in agreement. For other values of J0 (for example J0 = 1.5 min-1, S3 Fig) corresponding to the bi-stable domain sometimes we have found subtle differences (S3A Fig) and sometimes the results match quite well (S3B–S3D Fig) between the two models depending on whether we start the simulations from lower or upper steady states for a particular random number sequence. Since in the bi-stable domain it is hard to quantify and compare the noise statistics in a systematic manner, at this juncture we concentrate to understand why in case of the mechanistic model, the total protein level is not sensing the molecular noise at the mRNA level adequately at J0 = 3 min-1 for the current parameter set given in Table 2? We believe, this will eventually lead to the answer that under what conditions the QSSA model for the simple positive feedback motif will give similar result stochastically and will efficiently quantify the intrinsic noise like its detailed mechanistic version.

Comparison of intrinsic noise as a function of the half-lives of the mRNA (MP) and total protein (X)

In literature [3,4,14,15], it is well known that the half-lives of the mRNA and protein play a crucial role to determine the noise at the protein level. Pedraza and Paulsson [14] had shown that for simple transcription and translational events the intrinsic noise at the protein level is given by, (1) where, CVP is the coefficient of variation of the protein distribution, which not only depends on the average number of protein (<NP>) and mRNA (<NM>) molecules but it is also highly dependent on the half-lives τP and τM of the corresponding protein and mRNA molecules. Keeping these observations in mind, we wanted to verify whether the dependence on half-lives still exist in a gene regulatory network with feedback regulation or not and whether the stochastic calculation with QSSA model can capture such effects effectively or not like a detailed mechanistic model. At this point, we have taken two different approaches to vary the half-lives of mRNA and protein molecules to compare the two models.

Approach 1: Varying the ratio of the half-lives (KC) by keeping the absolute value of the mRNA half-life fixed.

We performed a systematic study where we have calculated the CV for the total protein and mRNA distributions as a function of different ratio () of mRNA and protein half-lives keeping the abundance levels of the protein (~ 458 molecules) and mRNA (~ 185 molecules) fixed (deterministically) for all the ratios of half-lives used (Fig 3 and S2 Table).

thumbnail
Download:
Fig 3. Plot of CV of protein and mRNA versus ratio of half-lives (KC) of MP and P at J0 = 3 min-1 and the corresponding sensitivity analysis.

(A) Coefficient of variation of total protein (X) abundance versus KC. (B) Coefficient of variation of MP abundance versus KC. km = 1E-01 min-1 is considered in all the cases. When KC = 0.01, then kp = 1E-03 min-1, J3 = 9.22E-04 min-1. When KC = 0.02, then kp = 2E-03 min-1, J3 = 1.85E-03 min-1. When KC = 0.025, then kp = 2.5E-03 min-1, J3 = 2.31E-03 min-1. When KC = 0.033, then kp = 3.34E-03 min-1, J3 = 3.08E-03 min-1. When KC = 0.046, then kp = 4.67E-03 min-1, J3 = 4.3E-03 min-1. When KC = 0.1, then kp = 1E-02 min-1, J3 = 9.22E-03 min-1. When KC = 0.2, then kp = 2E-02 min-1, J3 = 1.85E-02 min-1. When KC = 0.33, then kp = 3.34E-02 min-1, J3 = 3.08E-02 min-1. When KC = 1.0, then kp = 1E-01 min-1, J3 = 9.22E-02 min-1. When KC = 7.0, then kp = 7E-01 min-1, J3 = 6.5E-01 min-1. All other parameters are same as Table 2. (C) Sensitivity analysis of the parameters involved at KC = 1E-02 and J0 = 3 min-1. Sensitivity is measured on the basis of CV of stochastic mean of X as a function of rate constants involved in QSSA model and mechanistic model. Here CVX refers to the coefficient of variation in stochastic mean of X resulting due to parameter variation. It is divided by the CV of our standard stochastic model (Fig 2A) using our model-parameter set (Table 2). All the parameters are increased individually at an amount of 10% of the model parameters (Table 2) keeping all other parameters constant.

https://doi.org/10.1371/journal.pone.0136668.g003

This approach is similar to some earlier study performed by Shahrezaei and Swain [15] for simple transcription and translational regulations of a protein (without any feedback regulation) where they have shown that if γ > 1 (where according to our definition) then the analytical expression obtained by them for the noise at the protein level after QSSA on the mRNA dynamics, nicely fits the data for budding yeast experimental observations [13,27,28]. To obtain Fig 3A and 3B, we fixed the value of the mRNA half-life at τM = 7 min [29] (in Fig 2 we have used τM = 7 min, as km = 1E-01 min-1) and varied protein half-life τP (from 700 min to 1 min) [30] that eventually varies the ratio KC. From Fig 3A, it is quite evident that for KC<1, i.e. when the half-life of protein is much greater than the mRNA half-life, the stochastic calculation obtained from QSSA model resembles the stochastic result generated from mechanistic model if we compare the CV of total protein (X). CVX in case of stochastic calculation from QSSA model only starts to deviate significantly from that of the mechanistic version when the protein half-life is taken 35 mins (KC = 0.2) (Fig 3A and S2 Table). This indicates that in our case as the ratio KC becomes greater than 0.1, the mechanistic model starts to sense the molecular fluctuations due to mRNA more efficiently and hence the stochastic quantifications done from the two models start to differ. As the ratio KC is increased the CVX is also rising in case of the stochastic calculations done from the QSSA model and it shows similar trend as that of the mechanistic model but the rise is not as steep as it is in case of the mechanistic model simulations (Fig 3A). This is expected as we can see in Fig 3B (S2 Table) that the stochastic simulations from QSSA model cannot capture the molecular noise at the mRNA level even for the lowest magnitude of the ratio KC in comparison to the simulations performed with mechanistic model. We start to observe the effect of the differences in the mRNA CV levels found in Fig 3B at the CVX for the two models (Fig 3A), only for higher values of KC. This shows that the stochastic runs with QSSA model for the considered network are good to understand the trend of the intrinsic noise variation for the total protein level but it appears from Fig 3A that they are not appropriate for accurate quantification of the intrinsic noise if the ratio of the half-lives of the corresponding mRNA and protein are such that we have KC>0.1 but for KC≤0.1 the stochastic results from QSSA model will be quite reliable. At this juncture we did perform sensitivity analysis [31,32] of the rate constants involved in the QSSA and mechanistic models by taking the corresponding associated CVs for the total protein (X) (Fig 3C) at J0 = 3 min-1. We found that for KC = 0.01 situation, both the QSSA and mechanistic model parameters show quite identical sensitivities (Fig 3C) for all the rate constants related to the system, so no wonder that the stochastic calculation executed with the QSSA model gives similar results like the mechanistic model for KC = 0.01 (Fig 3A). The sensitivity analysis performed with KC = 0.2 (S4A Fig) and KC = 1 (S4B Fig) start to differ for QSSA and mechanistic models, which kind of supports why the CVs of total protein quantified from the two models start to vary as we increase the value of KC further in Fig 3A.

Approach 2: Varying the absolute values of the mRNA and protein half-lives keeping the ratio KC fixed.

Till now we have discussed if we change the ratio KC (by keeping τM = 7 min and changing only τP) how good a QSSA model performs in comparison to the mechanistic model to quantify the intrinsic noise. This kind of calculation makes sense for the case of budding yeast where for more than 80% of genes the value of KC<1 [15,27,28]. However, it is well-known that the absolute values of the half-lives (τM and τP) for mRNA and protein are equally inportant to the level of gene expression noise [3,4,14].

This implies that one can have fixed value of KC for different values of τM and τP, which can eventually change the burst size and burst frequencies and impact the overall intrinsic noise significantly. The important question is whether stochastic calculation from QSSA model can capture that feature accurately or not? To address this issue, we keep on changing the absolute values of the half-lives (τM and τP) for mRNA and protein by keeping the ratio KC fixed at 0.01 (Fig 4A and 4D), 1 (Fig 4B and 4E) and 100 (Fig 4C and 4F) (the average number of the total protein (~ 458 molecules) and mRNA (~ 185 molecules) are also kept fixed for all the cases deterministically, parameters are given in S2 Table). The way we varied the absolute values of the half-lives τM and τP (with KC fixed at 0.01, 1, 100) span almost all the possibilities of half-life combinations that can exist in reality for mammalian cells found in recent experiment [30]. The stochastic simulation results obtained for both the models performed equally well to quantify intrinsic noise at the total protein (X) level for all the KC values used in Fig 4A–4C for a certain range of absolute values of τM and τP. For KC = 0.01 (Fig 4A), the intirnsic noise quantified from QSSA model at the protein (X) level starts to differ from mechanistic model when τM = 0.1 min and τP = 10 min. Similarly for KC = 1 (Fig 4B), the deviation starts at τM = 7 min and τP = 7 min and for KC = 100 (Fig 4C), the deviation starts at τM = 10 min and τP = 0.1 min. Interestingly, we found that in Fig 4D–4F for certain combinations of absolute values of τM and τP (at different fixed values of KC), the stochastic QSSA model can also quantify the intrinsic noise at the mRNA (MP) level with reasonable accuracy in comparison to the mechanistic model. In all the cases (Fig 4A–4F), the CVs obtained from stochastic calculation of the QSSA model hardly changes for a fixed value of KC when the absolute values of the τM and τP are changed whereas the CVs calculated from the mechanistic model starts to deviate from QSSA model as soon as the value of either τM or τP or both τM and τP are such that they can significantly change the burst frequency and burst size related to the intrinsic noise of the concerned network. This evidently shows an important fact that even for ratio KC>0.1 the QSSA model can quantify the intrinsic noise (for protein as well as for mRNA) as accurately as the mechanistic model and we can have situations where for ratio KC<0.1 the QSSA model can fail to quantify the intrinsic noise (even at the protein level) as precisely as the mechanistic model. This result is quite different than what has been shown by Shahrezaei and Swain [15] earlier as they did not concentrate on changing the absolute values of τP and τM while performing their analysis.

thumbnail
Download:
Fig 4. Impact of absolute values of τP and τM on CV of X (Total protein) and MP (mRNA) as well as on ZX and (where Zi = % deviation of the measured intrinsic noise between the QSSA and the mechanistic models for i = X or MP).

Plot of CV of X versus τP and τM for (A) KC = 1E-02 (B) KC = 1.0 (C) KC = 100.0. Plot of CV of MP versus τP and τM for (D) KC = 1E-02 (E) KC = 1.0 (F) KC = 100.0. (G) Impact of absolute values of τP and τM on ZX. (H) Impact of absolute values of τP and τM on . τP and τM are the half-lives of P and MP respectively and are given in minutes.

https://doi.org/10.1371/journal.pone.0136668.g004

To understand this issue in more comprehensive manner, we further plotted the % deviation of the measured intrinsic noise (Zi) (Where Zi = [(CVi (mechanistic model)–CVi (QSSA model)) / CVi (mechanistic model)] x 100, i = X or MP) at the protein (X) (Fig 4G) and mRNA (MP) (Fig 4H) level in 2-dimentional heatmaps as a function of absolute values of τP and τM. In Fig 4G and 4H, we have defined 4 different regions by considering the stabilities of the mRNA and the corresponding protein. We have made such classification by keeping the experimental data available for budding yeast [27,28] as well as the data for Mammalian cell [30] in mind. It is quite evident from Fig 4G that in region I (unstable mRNA/unstable protein) the stochastic results obtained from QSSA model deviates (>40–50% deviation) the most from the mechanistic model and in region IV (stable mRNA/stable protein), stochastic results generated from QSSA model are as good as (within 5–10% deviation) the mechanistic model with orders of magnitude reduction in computational time. The stochastic results obtained for the case of region II (stable protein/unstable mRNA) and region III (unstable protein/stable mRNA) by using the QSSA model are moderately satisfactory (10–30% deviation) in comparison to stochastic results from mechanistic model. Fig 4H also evidently shows that the intrinsic fluctuations at the mRNA level can only be captured faithfully by stochastic QSSA model for combinations of τP and τM values falling within region IV (within 10–20% deviation) and for a part of region (III) where the mRNA is highly stable. This clearly shows that even if there is enough separation of timescale between the mRNA and protein dynamics (as in case of region II and III), if either of them has a short absolute value of half-life then bursting frequency generated because of that will be hard to capture by stochastic simulation using QSSA model. On the contrary, The QSSA model will do a fine job to quantify the intrinsic noise for both mRNA and protein even if there is no separation of time scale between mRNA and protein dynamics (as in case of region IV) as long as the absolute values of τP and τM are such that the bursting frequencies generated due to them are negligible. As soon as the absolute values of τP and τM are quite small (region I) the effective increase in bursting frequency will change the intrinsic noise significantly and QSSA model will fail to capture that effect.

Comparison of intrinsic noise as a function of number of protein molecules

In this section we wanted to understand how the QSSA model competes with the mechanistic model to capture the intrinsic noise as a function of number of protein molecule for a fixed number of mRNA molecules and fixed value of KC (Fig 5A–5D, approach 1). We will further investigate what happens if we change the absolute values of the τP and τM by keeping the KC fixed (Fig 5E and 5F, Approach 2).

thumbnail
Download:
Fig 5. Plot of CV of total protein (X) versus total protein (X) abundance following approach 1 and approach 2 by keeping deterministic mean of MP = 185.2 molecules for all the cases.

(A) KC = 1E-02, τP = 700 min (approach 1), (B) KC = 1E-01, τP = 70 min (approach 1), (C) KC = 1.0, τP = 7 min (approach 1), (D) KC = 100.0, τP = 7.0E-02 min (approach 1). In all the cases we kept τM = 7 min. (E) KC = 1E-02, τP = 10 min, τM = 1E-01 min (approach 2). (F) KC = 100.0, τP = 7 min, τM = 700 min (approach 2).

https://doi.org/10.1371/journal.pone.0136668.g005

The changes made in the parameter values (than given in the Table 2) to keep the deterministic mean of mRNA fixed ~ 185 molecules in all the cases for the Fig 5, are given in S3 Table. In Fig 5, we have shown how the CV of the total protein level will vary as a function of total protein abundance level keeping the deterministic mean of mRNA fixed at ~185 molecule at J0 = 3.0 min-1 for four different KC values. For KC = 0.01 (Fig 5A), the CV at the total protein level from the stochastic calculation obtained by using QSSA model seems to be identical as that obtained from mechanistic model as a function of total protein number. In both the cases the intrinsic noise at the total protein level decreases as the average number of protein molecule is increasing in a similar quantitative fashion. As the value of KC is increased to 0.1 (Fig 5B), 1.0 (Fig 5C) and 100.0 (Fig 5D), the statistical results quantified from stochastic calculation of QSSA model started to differ from the stochastic results obtained from the mechanistic model quantitatively. The qualitative feature of the stochastic calculation from QSSA model still showed similar trend as that of the stochastic result of the mechanistic model. We have performed similar studies to compare the intrinsic noise at the total protein level as a function of average values of total protein by fixing the number of mRNA molecules at lower and higher levels (S5 Fig, related changes in the parameter values given in S4 Table) and it strengthens the conclusions made from Fig 5. From this study it is evident that if the half-lives ratio of mRNA and protein remains as KC≤0.1, then the stochastic results using QSSA model will definitely quantify the intrinsic noise at the total protein level as accurately as stochastic calculation done with the mechanistic model even if the protein number starts to vary significantly in the corresponding biological network. It does not mean that the intrinsic noise calculated from QSSA model for KC = 100.0 (i.e., KC>0.1case) will always be inaccurate in comparison to mechanistic model and for KC<0.1 will always match. As we observed in case of Fig 4G, here too, if we measure the intrinsic noise as a function of total protein number using QSSA and mechanistic models by changing the absolute values of the half-lives τP and τM keeping the KC fixed, we found that even for KC<0.1 (Fig 5E, KC = 0.01, τM = 0.1 min and τP = 10 min) the results can differ significantly and for KC>0.1(Fig 5F, KC = 100.0, τM = 700 min and τP = 7 min) the results can match quite well. This clearly demonstrates the fact that the absolute values of the protein and mRNA half-lives are the important dictating factors to decide whether the QSSA model is good enough to capture the intrinsic noise for a gene regulatory network in comparison to the corresponding mechanistic model or not, number of molecules of total protein hardly matters.

Comparison of intrinsic noise as a function of number of mRNA molecules

In most of the biological systems, the abundance level of mRNAs is found to be much less than the corresponding protein molecules [1,3,12,30]. Keeping this in mind, next we wanted to investigate, how far the stochastic calculation from QSSA model can capture the effect due to the change in the mRNA numbers in comparison to the stochastic results obtained from mechanistic model at a fixed number of protein molecules and for different values of τP and τM (following both approach 1 and 2).

Following approach 1, we fixed the values of KC and kept the total protein level fixed at ~ 458 molecules in all the cases for the Fig 6A–6D (parameter values are given in the S5 Table). From Fig 6A it is evident that the CV calculated for the total protein level are comparable between the stochastic calculation done with QSSA and mechanistic models as a function of mRNA population when KC (0.01) is small. As the KC value is increased systematically by keeping the mean value of total protein fixed at ~ 458 molecules, the differences between the CV values at the total protein level (as a function of mRNA abundance) increased between the stochastic calculations performed with the QSSA and mechanistic models (Fig 6B with KC = 0.1, Fig 6C with KC = 1.0 and Fig 6D with KC = 100.0). The qualitative nature of the variations in CV values remain identical in both the cases but quantitatively they again start to differ as a function of average number of mRNA molecules for KC>0.1. In this case too, it is quite clear from Fig 6E and 6F that if we change the absolute values of τP and τM (by keeping KC fixed), the stochastic results calculated from QSSA model can differ for KC = 0.01 (Fig 6E, τM = 0.1 min and τP = 10 min) and can be similar for KC = 100.0 (Fig 6F, τM = 700 min and τP = 7 min) as a function of mRNA number.

thumbnail
Download:
Fig 6. Plot of CV of total protein (X) versus mRNA (MP) abundance following approach 1 and approach 2 by keeping deterministic mean of X = 457.9 molecules for all the cases.

(A) KC = 1E-02, τP = 700 min (approach 1), (B) KC = 1E-01, τP = 70 min (approach 1), (C) KC = 1.0, τP = 7 min (approach 1), (D) KC = 100.0, τP = 7.0E-02 min (approach 1). In all the cases we kept τM = 7 min. (E) KC = 1E-02, τP = 10 min, τM = 1E-01 min (approach 2). (F) KC = 100.0, τP = 7 min, τM = 700 min (approach 2).

https://doi.org/10.1371/journal.pone.0136668.g006

This clearly shows that before using any QSSA model for quantifying the molecular noise faithfully, one must carefully check the absolute values of the half-lives of the proteins and mRNAs involved in that network and the number of molecules of mRNA and protein will be less influencing factors for this purpose.

Considering the effect of additional feedback loops in the simple positive feedback module

To generalize the observations made in the previous sections, we have further systematically constructed mechanistic and QSSA models of two more modules (Fig 7). We have assumed that the protein (X) involved in the positive feedback module can further activate a protein (YP) which inturn can deactivate (Fig 7A) or activate (Fig 7B) the production of protein (X) by giving a feedback at the protein (X) level (Fig 7A, negative feedback and Fig 7B, positive feedback).

thumbnail
Download:
Fig 7. Various gene regulatory networks under consideration.

(A) Positive feedback with additional negative feedback motif, where X (total protein) activates the synthesis of another protein YP which negatively regulates X population level. (B) Double positive feedback motifs, where X (total protein) activates the synthesis of another protein YP which positively regulates the synthesis of X. In all the cases X autoregulates its own synthesis positively. Here MP denotes mRNA of P, YP denotes another protein and YM is the mRNA of protein YP. Detailed mechanistic models are given in (S6 Fig).

https://doi.org/10.1371/journal.pone.0136668.g007

We often encounter this kind of biological networks in different biological systems [33–35] and it is imperative to understand the reliability of the stochastic calculations using QSSA models to accurately evaluate the intrinsic noise regulations for such systems. The detailed version of the modules shown in Fig 7 are provided in S6 Fig and the corresponding QSSA models and the mechanisitc models along with the parameter values are given in the supplementary information (S6 and S8 Tables). We are interested in understanding how the molecular noise at the total protein (X) and mRNA (MP) levels will be influenced by the absolute values of (half-life of YM), (half-life of YP) and their ratio KS (where()) by following the two approaches discussed in the earlier sections. Further, we would like to investigate how far the stochastic simulations based on the QSSA models constructed for the two modules shown in Fig 7 can reproduce the stochastic simulation results obtained from the corresponding mechanistic models under these situations.

Comparison of intrinsic noise as a function of the half-lives of proteins and mRNAs (for Module 1).

With the small positive feedback module studied earlier, we have observed that the intrinsic noise at the total protein level (X) as well as at the mRNA (MP) level can be evaluated nicely by the QSSA model in comparison to the mechanistic model under certain restricted domain of absolute values of the corresponding protein and mRNA half-lives (approach 2) involved in the network. To further investigate in this direction, we took a step forward and considered the existence of an additional negative feedback loop of YP on X in presence of positive feedback of X on its own transcription (for example, in mammalian cell cycle E2F autocatalyzes its own transcription and also activates CycA which together with Cdk2 activates the degradation process of E2F [33]). Before doing the stochastic simulations, we first did the bifurcation analysis (S7A and S7B Fig (left panel QSSA model and right panel mechanistic model, parameters are given in S6 Table)) of the corresponding deterministic models to show that both the models are deterministically similar even when the additional negative feedback loop at the total protein (X) level is operative in the system. We have calculated the molecular noise (following approach 1) at the total protein (X) level by using the stochastic version of both the models (S1 Text) given in S7 Table for module 1 (Fig 7A) as a function of KS at two different fixed values of KC = 0.01(τM = 7 min and τP = 700 min, Fig 8A) and KC = 1.0 (τM = 7 min and τp = 7 min, S8A Fig).

thumbnail
Download:
Fig 8. Plot of CV of the total protein (X) and YP versus ratio of half-lives (KS) of YM and YP as well as the absolute values of the half-lives of YM and YP at KC = 1E-02 in the module 1 with additional negative feedback on X.

(A) τP = 700 min, τM = 7 min and = 7 min (B) τP = 10 min, τM = 1E-01 min and = 7 min. (C) KS = 1E-02, τp = 700 min and τM = 7 min. (D) KS = 1.0, τP = 700 min and τM = 7 min. Values of and are given in minutes.

https://doi.org/10.1371/journal.pone.0136668.g008

In Fig 8A, the CV for the total protein (X) calculated from QSSA model resembles the mechanistic model calculation but the result for the CV of YP protein differs as we go on changing the KS (following approach 1, fixed at 7 min) at a fixed value of the ratio KC (parameters used to keep the protein and mRNA numbers fixed are given in S7 Table). We have performed the similar comparison for KC = 1 (τM = 7 min and τP = 7 min) as well and found that for this value of KC (S8A Fig), QSSA model can not quantify the intrinsic noise precisely for both X and YP. This result is in agreement with what we have discussed earlier for the simple positive feedback motif and consistent with the findings of Shahrezaei and Swain [15].

At this point, we changed the absolute values of the half-lives (approach 2) of the protein (X) and mRNA (MP) by keeping the ratio KC same in Fig 8B (τM = 0.1 min and τP = 10 min for KC = 0.01) as well as in S8B Fig (τM = 700 min and τP = 700 min for KC = 1) and again performed same comparison between QSSA and mechanistic models. We can clearly see that now the stochastic results from QSSA model fails to capture the intrinsic noise accurately for both X and YP protein in comparison to mechanistic model in Fig 8B. On the contrary, S8B Fig shows that intrinsic noise at the level of X are quite comparable even for KC = 1 but in case of YP there is still disagreement. The reason behind this is quite clear from Fig 4G. In case of Fig 8A, the values of the τM (7 min) and τP (700 min) used fall in the top part of the region (II) where stochastic result from QSSA model are reasonably in good agreement with mechanistic model calculation whereas the values of the τM (0.1 min) and τP (10 min) used for Fig 8B fall in the region (I) where the disagreement between QSSA and mechanistic models is evident. In case of S8A Fig, the values used for τM and τP are both 7 mins respectively, which corresponds to region (I) and consequently there was disagreement and as soon as the values of τM and τP are both changed to 700 mins in S8B Fig, (coresponding to region IV) there was agreement between the stochastic results performed with QSSA and mechanistic models.

Further we wanted to investigate if we change the absolute values of the half-lives of the protein YP and mRNA YM at fixed values of KS and KC, how stochastic calculation with QSSA model performs in comparison to mechanistic model. To do this we considered two different cases KS = 0.01, KC = 0.01 (Fig 8C) and KS = 1, KC = 0.01 (Fig 8D). Fig 8C (left panel) and Fig 8D (left panel), vividly show that whatever may be the absolute values of and (for KS either fixed at 0.01 or at 1), the intrinsic noise at the level of protein X, can always be quantitatively reproduced by the QSSA model in comparison to mechanistic model even if the two pairs of absolute values of and fall in region (I) defined in Fig 4G. Under the same situation the QSSA model fails to quantify the intrinsic noise accurately at the level of protein YP (Fig 8C and 8D, right panel) whenever the pair of absolute values of and falls in region (I) corresponding to Fig 4G. Fig 8C (right panel) and Fig 8D (right panel) also clearly show that the QSSA model will quantify the intrinsic noise appropriately at the level of protein YP in comparison to mechanistic model if the absolute values of and fall in region (II) and (IV) defined in Fig 4G. To show that the stochastic QSSA model can even capture the mRNA fluctuations for module 1, we performed the comparison between stochastic QSSA and SSA calculation from mechanistic model for KS = 1 and KC = 1 (S8C Fig) by keeping the absolute values of τP = 700 min and τM = 700 min (to be in the region (IV) of Fig 4H). It is evident that the intrinsic fluctuations for the mRNA (MP) (S8C Fig, top right panel) can be captured beautifully by stochastic QSSA model for all the combinations of absolute values of and in different domains. Moreover, The stochastic QSSA model can satisfactorily quantify the intrinsic noise for mRNA (YM) (S8C Fig, bottom right panel) if the combinations of half-life values ( and ) correspond to region (IV) as mentioned in Fig 4H and the quantification starts to differ as soon as the combination of half-life values tend to fall in region (I) or (II).

We have performed similar studies with the Module 2 (with an additional positive feedback motif) and the detail results are provided in S2 Text and in the supplementary figures (S9–S11 Figs, S9 Table) referred there on. The stochastic simulation studies performed with Module 1 and Module 2, reiterate the fact that absolute values of the protein and mRNA half-lives will be the major factor to decide whether the stochastic simulation performed by a QSSA model can quantify the intrinsic noise as efficiently as mechanistic model or not.

Discussion

In important biological processes such as cell cycle regulation [6,9,10], maintainance of stemness [36,37], differentiation regulation [36,37] etc., intrinsic noise has been shown to play crucial role to determine the ultimate cell fate. Not only that, intrinsic noise can even drive a system to a completely different dynamical regime which is not permitted if we analyze the corresponding deterministic model under same conditions [38,39]. One of the standard practice to tackle intrinsic noise theoretically is to construct a proper mechanistic model in terms of mass action kinetic terms for the corresponding gene regulatory network under consideration and employ the Gillespie stochastic simulation algorithm (SSA) [16] to accurately quantify the molecular noise present in the system. Unfortunately, in most of the cases the Gillespie simulation for a detailed mechanistic model is highly expensive in terms of computational time. In literature, efforts had been made [15,17–19] to make appropriate approximations to simulate only the slow reactions involved in a large biological network by SSA. Among them the QSSA method found to be quite successful in reducing the computational time for SSA without losing much the accuracy of the result for simple biological systems [15,19]. In this regard the pertinent question that we raised is that whether it is appropriate to use these QSSA models to quantify the intrinsic noise for gene regulatory networks with feedback motifs? How much reliable they are to quantify the molecular noise in comparison to the detailed mechanistic models? This is an important question because if the QSSA models do a decent job to quantify the intrinsic noise in comparison to mechanistic models, then computationally SSA calculations will become much faster and easy to execute for larger networks.

By following approach 1, We have demonstrated that for the simple positive feedback motif under consideration, the stochastic calculations of QSSA models can quantify the intrinsic noise for a particular protein in a network with reasonable accuracy as a function of the ratio (KC) of the half-lives of mRNA and protein, provided both the protein and mRNA are reasonably stable (which means less bursting frequency) and under the condition KC<0.1. We observed that after a certain specific value of the ratio KC (KC>0.1), the stochastic results obtained from the QSSA model simulation for the protein (X) started to differ from the stochastic results of the mechanistic model simulation (Fig 3) which is quite in accordance with the observations made by Shahrezaei and Swain [15]. Although our approach 1 makes sense in the context of budding yeast where for more than 80% genes the protein half-life [13] is greater than the corresponding mRNA half-life [27,28], we have to keep in mind that even if the lifetimes of mRNA and protein are well separated in time scale (for the values KC<0.1) it might not guarantee that the reduced QSSA model will always be able to capture the intrinsic noise efficiently. We have further shown by implementing approach 2 that we can have disagreement for KC<0.1 and agreement for KC>0.1 between the two models and these things highly depend on the absolute values of the mRNA and protein half-lives. This result is quite unique and Shahrezaei and Swain [15] did not perform such analysis.

In this context, we have convincingly demonstrated using 2D heat maps (Fig 4G and 4H) that absolute values of the mRNA and protein half-lives are the crucial determining factors when one tries to compare the quality of the intrinsic noise quantified from a stochastic QSSA model in comparison to SSA performed using mechanistic model. We can have perfect agreement between the two calculations both at the protein and mRNA levels if the absolute values of the half-lives of the correposnding protein and the mRNA are located within region (IV) (as defined in Fig 4G or 4H). There will be a significant level of disagreement between the stochastic calculations performed from QSSA and mechanistic models if the absolute values of the half-lives correspond to region (I) i.e., when both the protein and mRNA are highly unstable causing a high degree of bursting frequency, which stochastic calculation performed from QSSA model will not be able to capture efficiently. As the half-lives of the proteins and mRNAs related to the gene regulatory network start to approach region (II) or (III), the success of the stochastic calculations using a QSSA model in comparison to mechanistic model will vary from case to case as shown in Fig 4G and 4H for both protein and mRNA respectively. This is a very important result in the context of recent experiments performed with mammalian cell [30] where they have shown that half-life combinations of protein and mRNA rarely turns out to be in the region (I) (as defined in Fig 4G) but most of the time the half-life combinations will fall in region (II), (III) and specially in region (IV) and beyond. This means that more often and not we can employ QSSA simplification for the corresponding mechanistic models of any gene regulatory network related to mammalian cell and have a good idea about the intrinsic noise if we take care of the additonal requirements as suggested by Thomas et al. [20].

We have further shown that the abundance level of protein (Fig 5) or mRNA (Fig 6) hardly influences the difference that we observe between the stochastic results obtained from QSSA and mechanistic models. To generalize our observations, we have further considered two different modules [33,34] (Fig 7) with additional negative and positive feedback loops and executed similar kind of studies as performed in case of simple positive feedback module shown in Fig 1A. These studies performed with module 1 and module 2 further corroborate the fact that absolute values of the half-lives of the proteins and mRNAs related to the gene regulatory network critically control the success of QSSA model in comparison to mechanistic model in the context of accuracy level achieved for intrinsic noise quantification. If the absolute values of all the half-lives (related to all the proteins and mRNAs involved in the extended GRNs) are such that they correspond to the region (IV) defined in Fig 4G and 4H, then we can definitely use the QSSA model to faithfully calculate the intrinsic noise for both protein and mRNA for a gene regulatory network. If they fall in region (II) or (III) then depending on other additional conditions the QSSA model might quantify the intrinsic noise with reasonable accuracy in comparison to its mechanistic counterpart but it will definitely fail to do so when the half-life values correspond to region (I).

In conclusion, our analysis with few frequently observed GRNs suggests that care must be taken before using any QSSA model for stochastic calculations although they are computationally less expensive. We must have an idea about the experimental half-lives of the proteins and mRNAs involved in the corresponding biological network before using a QSSA model to reliably quantify intrinsic noise for a GRN. We strongly believe that our finding about the importance of the absolute values of the half-lives while considering the efficiency of the stochastic QSSA model in comparison to SSA performed with mechanistic model is quite generic and will be applicable for even bigger GRNs. We hope that this work of us will shine some light to systematically use QSSA method for accurate measurement of intrinsic noise for large networks.

Materials and Methods

Deterministic simulations

The complete gene auto-regulatory network (Fig 1A) was expressed in terms of 5 ordinary differential equations in case of mechanistic model and 2 ordinary differential equations in case of QSSA model (Table 1). For simplicity we neglected the contribution of GP2 to total protein. As G+Ga+GP2 = 1, so GP2 level was very low. Thus we can safely neglect the contribution of GP2 in total protein. The deterministic models were encoded as .ode files and deterministic simulations were done by the freely available software XPP-AUT. The bifurcation diagrams were drawn using AUTO facility available in XPP-AUT. The bifurcation diagrams shown in Fig 1 and S5 Fig were drawn in MATLAB using the data points generated by XPP-AUT.

In a similar fashion each of the other modules (Fig 7) was expressed in terms of 9 ordinary differential equations in case of mechanistic models and 4 ordinary differential equations in case of QSSA models (S6 and S8 Tables). Based on these ordinary differential equations we have constructed the .ode files. For simplicity contribution of GP2 and GP2s to total protein population were neglected in all the cases. The bifurcation diagrams shown in S7 and S9 Figs were generated using the procedure discussed above.

Stochastic simulations

We employed stochastic simulation based on Gillespie’s algorithm [16]. In the traditional deterministic analysis reaction constants are considered as reaction rates but in stochastic analysis the reaction constants are considered as ‘probabilities per unit time’ [16]. The mechanistic model corresponding to Fig 1 consisted of 11 reactions and the QSSA model consisted of 5 reactions. The mechanistic models of the modules mentioned in Fig 7 were expressed in terms of 21 reactions and the corresponding QSSA models were expressed in terms of 11 reactions in each case. Corresponding reactions and propensities of the reactions were given in S1 Text. We have stochastically simulated the models using different random number sequences. The plots shown in the figures were drawn in MATLAB using the data points generated by stochastic simulation.

Supporting Information

S1 Fig. Stochastic simulation based on Gillespie Algorithm.

(A) In the left and right panels, we plot variation in population of X (total protein) with time in the QSSA model and mechanistic model respectively at J0 = 3 min-1. (B) In the left and right panels, we plot variation of MP with time in the QSSA model and mechanistic model respectively at J0 = 3 min-1. (C) Stochastic simulation based on Gillespie Algorithm using another random number sequence. In the left and right panels, we plot steady state distribution of X (total protein) in the QSSA model at J0 = 3 min-1, mean = 457.3, standard deviation = 42.69, CV = 9.34% and steady state distribution of X (total protein) in the mechanistic model, mean = 461.4, standard deviation = 43.46, CV = 9.42% respectively. All the values (except CVs) are in number of molecules.

https://doi.org/10.1371/journal.pone.0136668.s001

(TIFF)

S2 Fig. Steady state analysis of QSSA and mechanistic models.

In the left panel, we plot steady state distribution of X (total protein) in the QSSA model at J0 = 3 min-1, when ka and kd are 8E-03 min-1 and 5E-03 min-1, mean = 457.2, standard deviation = 42.98, CV = 9.4%. In the right panel, we plot steady state distribution of X (total protein) in the mechanistic model at J0 = 3 min-1, when ka and kd are 8E-03 min-1 and 5E-03 min-1, mean = 523.7, standard deviation = 197.6, CV = 37.7%. All the values (except CVs) are in number of molecules.

https://doi.org/10.1371/journal.pone.0136668.s002

(TIFF)

S3 Fig. Steady state distribution of X at J0 = 1.5 min-1.

(A) In the left and right panels, we plot the steady state distributions of X (total protein) in the QSSA model and mechanistic model respectively starting from upper steady state. (B) In the left and right panels, we plot the steady state distributions of X (total protein) in the QSSA model and mechanistic model respectively starting from lower steady state. Using another random number sequence in the left and right panels, we plot the steady state distribution of X (total protein) in the QSSA model and mechanistic model respectively starting from (C) upper steady state (D) lower steady state.

https://doi.org/10.1371/journal.pone.0136668.s003

(TIFF)

S4 Fig. Sensitivity analysis of the parameters involved at J0 = 3 min-1.

(A) KC = 0.2 (τP = 700 min, τM = 35 min) and (B) KC = 1.0 (τP = 700 min, τM = 7 min). For both cases sensitivity is measured on the basis of CVs of stochastic mean of X as a function of rate constants involved in QSSA model and mechanistic model. Here CVX refers to the coefficient of variation in stochastic mean of X resulting due to parameter variation. It is divided by the CV of our standard stochastic model (Fig 3A) using our model-parameter set at KC = 0.2. In all the cases parameters are increased individually at an amount of 10% of the model parameters keeping all other parameters constant.

https://doi.org/10.1371/journal.pone.0136668.s004

(TIFF)

S5 Fig. Noise (CV) level in X with increase in X population at various MP population levels.

(A) Deterministically similar QSSA and mechanistic models, deterministic mean of X (molecules) = 457.9 and MP (molecules) = 46.3 at J0 = 7.5E-01 min-1. The solid lines represent stable steady states and dotted lines represent unstable steady states; J0 is the basal rate of MP synthesis and acts as the bifurcation parameter. (B) Coefficient of variation of total protein (X) versus X abundance keeping deterministic mean of MP = 46.3 molecules at J0 = 7.5E-01 min-1 with various KC. CVs of X are plotted at KC = 1E-02, KC = 1E-01, KC = 2E-01 and KC = 1.0 respectively. (C) Deterministically similar QSSA and mechanistic models, deterministic mean of X (molecules) = 457.9 and MP (molecules) = 370.4 at J0 = 6 min-1. (D) Coefficient of variation of X versus X abundance keeping deterministic mean of MP = 370.4 molecules at J0 = 6 min-1 with various KC. CVs of X are plotted at KC = 1E-02, KC = 1E-01, KC = 2E-01 and KC = 1.0 respectively. (E) Deterministically similar QSSA and mechanistic models with deterministic mean of X (molecules) = 457.9 and MP (molecules) = 740.9 at J0 = 12 min-1. (F) Coefficient of variation of X versus X abundance keeping deterministic mean of MP = 740.9 molecules at J0 = 12 min-1 with various KC. CVs of X are plotted at KC = 1E-02, KC = 1E-01, KC = 2E-01 and KC = 1.0 respectively. In all the cases τM~7 min, km = 1E-01 min-1. Deviation between QSSA model and mechanistic model is more pronounced at higher KC.

https://doi.org/10.1371/journal.pone.0136668.s005

(TIFF)

S6 Fig. Mechanistic frameworks of different gene regulatory networks under consideration.

(A) Gene regulatory network with positive and additional negative feedback motifs under consideration where G and Ga denote inactive and active states of gene, MP denotes mRNA, P denotes protein, P2 is the dimer of protein. The protein (P) molecules form dimers P2. P2 binds to the promoter region of the inactive gene and activates G to Ga as well as P2 binds to the promoter region of another inactive gene Gs and activates Gs to Gas. mRNAs (YM) are synthesized at a basal rate J4. km and kp are the degradation rates of mRNA of P (MP) and P and kym and kyp are the degradation rates of mRNA of YP (YM) and YP. YP activates the degradation of protein P at Kn rate. (B) Gene regulatory network with double positive feedback motifs under consideration where G and Ga denote inactive and active states of gene, Mp denotes mRNA, P denotes protein, P2 is the dimer of protein. The protein P form dimers P2. P2 binds to the promoter region of the inactive gene and activates G to Ga as well as P2 binds to the promoter region of another inactive gene Gs and activates Gs to Gas. mRNAs (YM) are synthesized at a basal rate J4. km and kp are the degradation rates of MP and P and kym and kyp are the degradation rates of YM and YP. YP activates the synthesis of protein P at Kn rate.

https://doi.org/10.1371/journal.pone.0136668.s006

(TIFF)

S7 Fig. Deterministically similar QSSA and mechanistic models of gene regulatory networks with additional positive and negative feedback motifs.

The solid lines represent stable steady states and dotted lines represent unstable steady states; J0 is the basal rate of MP synthesis and acts as the bifurcation parameter. (A) In the left and right panels, we plot steady state population of X in QSSA and mechanistic models respectively. (B) In the left and right panels, we plot steady state population of YP in QSSA and mechanistic models respectively at KC = 1E-02, KS = 1E-02, τP = 700 min, τM = 7 min [S6 Table]. Values of Kn, given in the plots, are in molecule-1min-1. Deterministic mean of X = 390.9 molecules, YP = 149.6 molecules, MP = 174.2 molecules and YM = 162.3 molecules at J0 = 3 min-1 when Kn = 2.5E-07 molecule-1min-1. Deterministic mean of X = 457.9 molecules, YP = 161.1 molecules, MP = 185.2 molecules and YM = 174.7 molecules at J0 = 3 min-1 when Kn = 0 molecule-1min-1.

https://doi.org/10.1371/journal.pone.0136668.s007

(TIFF)

S8 Fig. Plot of CV of proteins and mRNAs versus ratio of half-lives (KS) of YM and YP as well as the absolute values of the half-lives of YM and YP at KC = 1.0 in the module 1 with additional negative feedback on X.

(A) τP = 7 min, τM = 7 min and = 7 min (B) τP = 700 min, τM = 700 min and = 7 min. (C) Plots of CV of X, MP, YP and YM abundance in QSSA and mechanistic models versus absolute values of the half-lives of YM and YP at KS = 1.0, keeping τP = 700 min, τM = 700 min. Values of and are given in minutes.

https://doi.org/10.1371/journal.pone.0136668.s008

(TIFF)

S9 Fig. Deterministically similar QSSA and mechanistic models of gene regulatory network with double positive feedback motifs (additional positive feedback is working on X).

The solid lines represent stable steady states and dotted lines represent unstable steady states; J0 is the basal rate of MP synthesis and acts as the bifurcation parameter. (A) In the left and right panels, we plot steady state population of X in QSSA and mechanistic models respectively. (B) In the left and right panels, we plot steady state population of YP in QSSA and mechanistic models respectively at KC = 1E-02, KS = 1E-02, τP = 700 min, τM = 7 min [S8 Table]. Values of Kn, given in the plots, are in min-1. Deterministic mean of X = 543.6 molecules, YP = 172.6 molecules, MP = 196.3 molecules and YM = 187.1 molecules at J0 = 3 min-1 when Kn = 5E-05 min-1. Deterministic mean of X = 457.9 molecules, YP = 161.1 molecules, MP = 185.2 molecules and YM = 174.7 molecules at J0 = 3 min-1 when Kn = 0 min-1.

https://doi.org/10.1371/journal.pone.0136668.s009

(TIFF)

S10 Fig. Plot of CV of the total protein (X) and YP versus ratio of half-lives (KS) of YM and YP as well as the absolute values of the half-lives of YM and YP at KC = 1E-02 in the module 2 with additional positive feedback on X.

(A) τP = 700 min, τM = 7 min and = 7 min (B) τP = 10 min, τM = 1E-01 min and = 7 min. (C) KS = 1E-02, τP = 700 min and τM = 7 min. (D) KS = 1.0, τP = 700 min and τM = 7 min. Values of and are given in minutes.

https://doi.org/10.1371/journal.pone.0136668.s010

(TIFF)

S11 Fig. Plot of CV of proteins and mRNAs versus ratio of half-lives (KS) of YM and YP as well as the absolute values of the half-lives of YM and YP at KC = 1.0 in the module 2 with additional positive feedback on X.

(A) τP = 7 min, τM = 7 min and = 7 min (B) τP = 700 min, τM = 700 min and = 7 min. (C) Plots of CV of X, MP, YP and YM abundance in QSSA and mechanistic models versus absolute values of the half-lives of YM and YP at KS = 1.0, keeping τP = 700 min, τM = 700 min. Values of and are given in minutes.

https://doi.org/10.1371/journal.pone.0136668.s011

(TIFF)

S1 Table. Steady state analysis of the QSSA and mechanistic models.

[S2 Fig].

https://doi.org/10.1371/journal.pone.0136668.s012

(DOCX)

S2 Table. (A) Stochastic results of Total protein (X) and mRNA (MP) at different KC in QSSA model and mechanistic model. [Fig 3] (B) Parameters used in Fig 4.

https://doi.org/10.1371/journal.pone.0136668.s013

(DOCX)

S3 Table. Parameters used in Fig 5.

https://doi.org/10.1371/journal.pone.0136668.s014

(DOCX)

S4 Table. Parameters used in S5 Fig.

https://doi.org/10.1371/journal.pone.0136668.s015

(DOCX)

S5 Table. Parameters used in Fig 6.

https://doi.org/10.1371/journal.pone.0136668.s016

(DOCX)

S6 Table. Equations and parameters of the positive feedback with additional negative feedback on Protein level (Module 1).

[S7A and S7B Fig].

https://doi.org/10.1371/journal.pone.0136668.s017

(DOCX)

S7 Table. Parameters used in Fig 8 and S8 Fig.

https://doi.org/10.1371/journal.pone.0136668.s018

(DOCX)

S8 Table. Equations and parameters of the additional positive feedback on Protein level (Module 2).

[S9A and S9B Fig].

https://doi.org/10.1371/journal.pone.0136668.s019

(DOCX)

S9 Table. Parameters used in S10–S11 Figs.

https://doi.org/10.1371/journal.pone.0136668.s020

(DOCX)

S1 Text. Mathematical modeling and stochastic version of different modules.

https://doi.org/10.1371/journal.pone.0136668.s021

(DOCX)

S2 Text. Comparison of intrinsic noise as a function of the half-lives of proteins and mRNAs (for Module 2).

https://doi.org/10.1371/journal.pone.0136668.s022

(DOCX)

Acknowledgments

Thanks are due to IRCC, IIT Bombay (13IRTAPSG005) for a fellowship to (DS). This work is supported by the funding agencies IRCC, IIT Bombay (13IRCCSG008) and DST grant (EMR/2014/000500).

Author Contributions

Conceived and designed the experiments: DS SK. Performed the experiments: DS. Analyzed the data: DS SK. Contributed reagents/materials/analysis tools: DS SK. Wrote the paper: DS SK.

References

  1. 1. Kar S, Baumann WT, Paul MR, Tyson JJ. Exploring the roles of noise in the eukaryotic cell cycle. Proc Natl Acad Sci. 2009;106: 6471–6476. pmid:19246388
  2. 2. Meyer R, D’Alessandro LA, Kar S, Kramer B, She B, Kaschek D, et al. Heterogeneous kinetics of AKT signaling in individual cells are accounted for by variable protein concentration. Front Physiol. 2012;3 NOV: 1–14.
  3. 3. Swain PS, Elowitz MB, Siggia ED. Intrinsic and extrinsic contributions to stochasticity in gene expression. Proc Natl Acad Sci. 2002;99: 12795–12800. pmid:12237400
  4. 4. Blake WJ, KÆrn M, Cantor CR, Collins JJ. Noise in eukaryotic gene expression. Nature. 2003;422: 633–637. pmid:12687005
  5. 5. Kaern M, Elston TC, Blake WJ, Collins JJ. Stochasticity in gene expression: from theories to phenotypes. Nat Rev Genet. 2005;6: 451–464. pmid:15883588
  6. 6. Yao G, Lee TJ, Mori S, Nevins JR, You L. A bistable Rb-E2F switch underlies the restriction point. Nat Cell Biol. 2008;10: 476–482. pmid:18364697
  7. 7. Yao G, Tan C, West M, Nevins JR, You L. Origin of bistability underlying mammalian cell cycle entry. Mol Syst Biol. Nature Publishing Group; 2011;7: 485.
  8. 8. Di Talia S, Skotheim JM, Bean JM, Siggia ED, Cross FR. The effects of molecular noise and size control on variability in the budding yeast cell cycle. Nature. 2007;448: 947–951. pmid:17713537
  9. 9. Spencer SL, Cappell SD, Tsai FC, Overton KW, Wang CL, Meyer T. The proliferation-quiescence decision is controlled by a bifurcation in CDK2 activity at mitotic exit. Cell. Elsevier Inc.; 2013;155: 369–383.
  10. 10. Overton KW, Spencer SL, Noderer WL, Meyer T, Wang CL. Basal p21 controls population heterogeneity in cycling and quiescent cell cycle states. Proc Natl Acad Sci. 2014;111: E4386–E4393. pmid:25267623
  11. 11. Niepel M, Spencer SL, Sorger PK. Non-genetic cell-to-cell variability and the consequences for pharmacology. Curr Opin Chem Biol. 2009;13: 556–561. pmid:19833543
  12. 12. Zenklusen D, Larson DR, Singer RH. Single-RNA counting reveals alternative modes of gene expression in yeast. Nat Struct Mol Biol. 2008;15: 1263–1271. pmid:19011635
  13. 13. Belle A, Tanay A, Bitincka L, Shamir R, O’Shea EK. Quantification of protein half-lives in the budding yeast proteome. Proc Natl Acad Sci. 2006;103: 13004–13009. pmid:16916930
  14. 14. Pedraza JM, Paulsson J. Effects of molecular memory and bursting on fluctuations in gene expression. Science. 2008;319: 339–343. pmid:18202292
  15. 15. Shahrezaei V, Swain PS. Analytical distributions for stochastic gene expression. Proc Natl Acad Sci. 2009;106: 17256–17261.
  16. 16. Gillespie DT. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J Comput Phys. 1976;22: 403–434.
  17. 17. Cao Y, Gillespie DT, Petzold LR. Accelerated stochastic simulation of the stiff enzyme-substrate reaction. J Chem Phys. 2005; 123: 144917. pmid:16238434
  18. 18. Cao Y, Gillespie DT, Petzold LR. The slow-scale stochastic simulation algorithm. J Chem Phys. 2005; 122: 014116.
  19. 19. Rao C V. Stochastic chemical kinetics and the quasi-steady-state assumption: Application to the Gillespie algorithm. J Chem Phys. 2003;118: 4999–5010.
  20. 20. Thomas P, Straube A V, Grima R. The slow-scale linear noise approximation: an accurate, reduced stochastic description of biochemical networks under timescale separation conditions. BMC Syst Biol. BioMed Central Ltd; 2012; 6: 39.
  21. 21. Karmakar R, Bose I. Positive feedback, stochasticity and genetic competence. Phys Biol. 2007;4: 29–37. pmid:17406083
  22. 22. Smits WK, Kuipers OP, Veening J-W. Phenotypic variation in bacteria: the role of feedback regulation. Nat Rev Microbiol. 2006;4: 259–271. pmid:16541134
  23. 23. Isaacs FJ, Hasty J, Cantor CR, Collins JJ. Prediction and measurement of an autoregulatory genetic module. Proc Natl Acad Sci. 2003;100: 7714–7719. pmid:12808135
  24. 24. Ferrell JE. Self-perpetuating states in signal transduction: Positive feedback, double-negative feedback and bistability. Curr Opin Cell Biol. 2002;14: 140–148. pmid:11891111
  25. 25. Maamar H, Dubnau D. Bistability in the Bacillus subtilis K-state (competence) system requires a positive feedback loop. Mol Microbiol. 2005;56: 615–624. pmid:15819619
  26. 26. Smits WK, Eschevins CC, Susanna KA, Bron S, Kuipers OP, Hamoen LW. Stripping Bacillus: ComK auto-stimulation is responsible for the bistable response in competence development. Mol Microbiol. 2005;56: 604–614. pmid:15819618
  27. 27. Grigull J, Mnaimneh S, Pootoolal J, Mark D, Hughes TR, Robinson MD. Genome-Wide Analysis of mRNA Stability Using Transcription Inhibitors and Microarrays Reveals Posttranscriptional Control of Ribosome Biogenesis Factors Genome-Wide Analysis of mRNA Stability Using Transcription Inhibitors and Microarrays Reveals Posttran. Mol Cell Biol. 2004;24: 5534–5547. pmid:15169913
  28. 28. Wang Y, Liu CL, Storey JD, Tibshirani RJ, Herschlag D, Brown PO. Precision and functional specificity in mRNA decay. Proc Natl Acad Sci. 2002;99: 5860–5865. pmid:11972065
  29. 29. Holstege FC, Jennings EG, Wyrick JJ, Lee TI, Hengartner CJ, Green MR, et al. Dissecting the regulatory circuitry of a eukaryotic genome. Cell. 1998;95: 717–728. pmid:9845373
  30. 30. Schwanhäusser B, Busse D, Li N, Dittmar G, Schuchhardt J, Wolf J, et al. Global quantification of mammalian gene expression control. Nature. 2011;473: 337–342. pmid:21593866
  31. 31. Sarkar AX, Sobie EA. Regression analysis for constraining free parameters in electrophysiological models of cardiac cells. Plos Comput Biol. 2010;6: e1000914. pmid:20824123
  32. 32. Perumal TM, Gunawan R. Understanding dynamics using sensitivity analysis: caveat and solution. BMC Syst Biol. BioMed Central Ltd; 2011;5: 41.
  33. 33. Gérard C, Goldbeter A. Temporal self-organization of the cyclin/Cdk network driving the mammalian cell cycle. Proc Natl Acad Sci. 2009;106: 21643–21648. pmid:20007375
  34. 34. Aguda BD, Kim Y, Piper-Hunter MG, Friedman A, Marsh CB. MicroRNA regulation of a cancer network: consequences of the feedback loops involving miR-17-92, E2F, and Myc. Proc Natl Acad Sci. 2008;105: 19678–19683. pmid:19066217
  35. 35. Novák B, Tyson JJ. Design principles of biochemical oscillators. Nat Rev Mol Cell Biol. 2008;9: 981–991. pmid:18971947
  36. 36. Kalmar T, Lim C, Hayward P, Muñoz-Descalzo S, Nichols J, Garcia-Ojalvo J, et al. Regulated fluctuations in Nanog expression mediate cell fate decisions in embryonic stem cells. Plos Biol. 2009;7: 33–36.
  37. 37. Garcia-ojalvo J, Elowitz MB. Tunability and noise dependence in Differentiation Dynamics. Science. 2007;315: 1716–1719. pmid:17379809
  38. 38. McKane AJ, Newman TJ. Predator-prey cycles from resonant amplification of demographic stochasticity. Phys Rev Lett. 2005;94: 1–4.
  39. 39. Cantini L, Cianci C, Fanelli D, Massi E, Barletti L. Linear noise approximation for stochastic oscillations of intracellular calcium. J Theor Biol. Elsevier; 2014;349: 92–99.