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

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.


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][2][3][4][5], namely intrinsic and extrinsic noise. Recently, there were a number of experimental findings [6][7][8][9][10][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][18][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][22][23][24][25][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 (M P ) and obtained the corresponding QSSA model (left panel, Table 1) for the given GRN.
The parameter values for the model concerned are given in Table 2. Here we have kept k m >>J 3 , 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 J 0 (the basal rate of the M P synthesis), signify that the two models provided in Table 1 are deterministically similar (the corresponding bifurcation diagrams of M P as a function of J 0 are given in Fig 1C). This clearly shows that both the deterministic models will lead to identical steady states deterministically and catalytic effects respectively. G and G a denote inactive and active states of gene, M p denotes mRNA, P denotes protein, P 2 is the dimer of protein and X is the total protein. The protein (P) molecules form dimers P 2 . P 2 binds to the promoter region of the inactive gene (G) and activates G to G a . k m and k p are the degradation rates of M P 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 J 0 . The solid lines represent stable steady states and dotted lines represent unstable steady states; J 0 is the basal rate of mRNA (M P ) synthesis and acts as the bifurcation parameter. In both the models deterministic mean of X (molecules) = 457.9 at J 0 = 3 min -1 . (C) Bifurcation diagrams (left panel, QSSA model and right panel, mechanistic model) of mRNA (M P ) are plotted as a function of J 0 . In both the models deterministic mean of M P (molecules) = 185.2 at J 0 = 3 min -1 .

Comparison of the stochastic distributions at different values of J 0
Once we created the deterministically similar models, we employed the Gillespie's simulation algorithm [16] at J 0 = 3 min -1 (deterministic means of X = 457.9 molecules and M P = 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 J 0 = 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. 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, k a = 8E-03 min -1 and k d = 5E-03 min -1 instead of k a = 8 min -1 and k d = 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 (M P ) 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 M P 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 J 0 (for example J 0 = 1. 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 J 0 = 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 (M P ) 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, where, CV P is the coefficient of variation of the protein distribution, which not only depends on the average number of protein (<N P >) and mRNA (<N M >) 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 (K C ) 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 (K C ¼ t M = t P ) 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). Coefficient of variation of total protein (X) abundance versus K C . (B) Coefficient of variation of M P abundance versus K C . k m = 1E-01 min -1 is considered in all the cases. When K C = 0.01, then k p = 1E-03 min -1 , J 3 = 9.22E-04 min -1 . When K C = 0.02, then k p = 2E-03 min -1 , J 3 = 1.85E-03 min -1 . When K C = 0.025, then k p = 2.5E-03 min -1 , J 3 = 2.31E-03 min -1 . When K C = 0.033, then k p = 3.34E-03 min -1 , J 3 = 3.08E-03 min -1 . When K C = 0.046, then k p = 4.67E-03 min -1 , J 3 = 4.3E-03 min -1 . When K C = 0.1, then k p = 1E-02 min -1 , J 3 = 9.22E-03 min -1 . When K C = 0.2, then k p = 2E-02 min -1 , J 3 = 1.85E-02 min -1 . When K C = 0.33, then k p = 3.34E-02 min -1 , J 3 = 3.08E-02 min -1 . When K C = 1.0, then k p = 1E-01 min -1 , J 3 = 9.22E-02 min -1 . When K C = 7.0, then k p = 7E-01 min -1 , J 3 = 6.5E-01 min -1 . All other parameters are same as Table 2. (C) Sensitivity analysis of the parameters involved at K C = 1E-02 and J 0 = 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 CV X 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 modelparameter set ( Table 2). All the parameters are increased individually at an amount of 10% of the model parameters ( 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 g ¼ 1 = K C 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 k m = 1E-01 min -1 ) and varied protein half-life τ P (from 700 min to 1 min) [30] that eventually varies the ratio K C . From Fig 3A, it is quite evident that for K C <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). CV X 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 (K C = 0.2) (Fig 3A and S2 Table). This indicates that in our case as the ratio K C 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 K C is increased the CV X 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 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 K C 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 CV X for the two models ( Fig 3A), only for higher values of K C . 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 K C >0.1 but for K C 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 J 0 = 3 min -1 . We found that for K C = 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 K C = 0.01 ( Fig 3A). The sensitivity analysis performed with K C = 0.2 (S4A Fig) and K C = 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 K C further in Fig 3A. Approach 2: Varying the absolute values of the mRNA and protein half-lives keeping the ratio K C fixed. Till now we have discussed if we change the ratio K C (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 K C <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 K C 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 K C 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 K C 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 K C values used in Fig 4A-4C for a certain range of absolute values of τ M and τ P . For K C = 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 K C = 1 (Fig 4B), the deviation starts at τ M = 7 min and τ P = 7 min and for K C = 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 K C ), the stochastic QSSA model can also quantify the intrinsic noise at the mRNA (M P ) 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 K C 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 K C >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 K C <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.
To understand this issue in more comprehensive manner, we further plotted the % deviation of the measured intrinsic noise (Z i ) (Where Z i = [(CV i (mechanistic model)-CV i (QSSA model)) / CV i (mechanistic model)] x 100, i = X or M P ) at the protein (X) ( Fig 4G) and mRNA (M P ) ( 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 K C (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 K C fixed (Fig 5E and 5F, Approach 2).
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 J 0 = 3.0 min -1 for four different K C values. For K C = 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 K C 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 K C 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 K C = 100.0 (i.e., K C >0.1case) will always be inaccurate in comparison to mechanistic model and for K C <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 K C fixed, we found that even for K C <0.1 (Fig 5E, K

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 K C and kept the total protein level fixed at4 58 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 K C (0.01) is small. As the K C 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 K  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 (Y P ) 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).
We often encounter this kind of biological networks in different biological systems [33][34][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  information (S6 and S8 Tables). We are interested in understanding how the molecular noise at the total protein (X) and mRNA (M P ) levels will be influenced by the absolute values of t Y M (half-life of Y M ), t Y P (half-life of Y P ) and their ratio K S (where(K S ¼ t Y M = t Y P )) 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 (M P ) 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 Y P 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 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 K S at two different fixed values of K C = 0.01(τ M = 7 min and τ P = 700 min, Fig 8A) and K C = 1.0 (τ M = 7 min and τ p = 7 min, S8A Fig).
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 Y P protein differs as we go on changing the K S (following approach 1, t Y M fixed at 7 min) at a fixed value of the ratio K C  (parameters used to keep the protein and mRNA numbers fixed are given in S7 Table). We have performed the similar comparison for K C = 1 (τ M = 7 min and τ P = 7 min) as well and found that for this value of K C (S8A Fig), QSSA model can not quantify the intrinsic noise precisely for both X and Y P . 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 (M P ) by keeping the ratio K C same in Fig 8B (τ M = 0.1 min and τ P = 10 min for K C = 0.01) as well as in S8B Fig (τ M = 700 min and τ P = 700 min for K C = 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 Y P 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 K C = 1 but in case of Y P there is still disagreement. The reason behind this is quite clear from Fig 4G. In case of Fig 8A,  Further we wanted to investigate if we change the absolute values of the half-lives of the protein Y P and mRNA Y M at fixed values of K S and K C , how stochastic calculation with QSSA model performs in comparison to mechanistic model. To do this we considered two different cases K S = 0.01, K C = 0.01 ( Fig 8C) and K S = 1, K C = 0.01 (Fig 8D). Fig 8C (left panel) and Fig  8D (left panel), vividly show that whatever may be the absolute values of t Y M and t Y P (for K S 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 t Y M and t Y P 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 Y P (Fig 8C and 8D 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 Applicability of QSSA Models to Quantify Intrinsic Noise 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][18][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 (K C ) 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 K C <0.1. We observed that after a certain specific value of the ratio K C (K C >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 K C <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 K C <0.1 and agreement for K C >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.

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 GP 2 to total protein. As G+G a +GP 2 = 1, so GP 2 level was very low. Thus we can safely neglect the contribution of GP 2 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 GP 2 and GP 2s 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. . 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 CV X 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)  (A) Gene regulatory network with positive and additional negative feedback motifs under consideration where G and G a denote inactive and active states of gene, M P denotes mRNA, P denotes protein, P 2 is the dimer of protein. The protein (P) molecules form dimers P 2 . P 2 binds to the promoter region of the inactive gene and activates G to G a as well as P 2 binds to the promoter region of another inactive gene G s and activates G s to G as . mRNAs (Y M ) are synthesized at a basal rate J 4. k m and k p are the degradation rates of mRNA of P (M P ) and P and k ym and k yp are the degradation rates of mRNA of Y P (Y M ) and Y P . Y P activates the degradation of protein P at K n rate. (B) Gene regulatory network with double positive feedback motifs under consideration where G and G a denote inactive and active states of gene, M p denotes mRNA, P denotes protein, P 2 is the dimer of protein. The protein P form dimers P 2 . P 2 binds to the promoter region of the inactive gene and activates G to G a as well as P 2 binds to the promoter region of another inactive gene G s and activates G s to G as . mRNAs (Y M ) are synthesized at a basal rate J 4. k m and k p are the degradation rates of M P and P and k ym and k yp are the degradation rates of Y M and Y P . Y P activates the synthesis of protein P at K n rate. 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 Y P in QSSA and mechanistic models respectively at K C = 1E-02, K S = 1E-02, τ P = 700 min, τ M = 7 min [S6 Table]. Values of K n , given in the plots, are in molecule -1 min -1 . Deterministic mean of X = 390.9 molecules, Y P = 149. 6