Nonlinear Biochemical Signal Processing via Noise Propagation

# Nonlinear Biochemical Signal Processing via Noise Propagation

Kyung Hyuk Kim, Hong Qian, and Herbert M. Sauro
Department of Bioengineering, University of Washington, William H. Foege Building, Box 355061, Seattle, WA 98195, U.S.A.
Department of Applied Mathematics, University of Washington, Lewis Hall, Box 353925, Seattle, WA 98195, U.S.A.
August 23, 2019
###### Abstract

Single-cell studies often show significant phenotypic variability due to the stochastic nature of intra-cellular biochemical reactions. When the numbers of molecules, e.g., transcription factors and regulatory enzymes, are in low abundance, fluctuations in biochemical activities become significant and such “noise” can propagate through regulatory cascades in terms of biochemical reaction networks. Here we develop an intuitive, yet fully quantitative method for analyzing how noise affects cellular phenotypes based on identifying a system’s nonlinearities and noise propagations. We observe that such noise can simultaneously enhance sensitivities in one behavioral region while reducing sensitivities in another. Employing this novel phenomenon we designed three biochemical signal processing modules: () A gene regulatory network that acts as a concentration detector with both enhanced amplitude and sensitivity. () A non-cooperative positive feedback system, with a graded dose-response in the deterministic case, that serves as a bistable switch due to noise-induced bimodality. () A noise-induced linear amplifier for gene regulation that requires no feedback. The methods developed in the present work allow one to understand and engineer nonlinear biochemical signal processors based on fluctuation-induced phenotypes.

## I Introduction

Quantifying biochemical processes at the cellular level is becoming increasingly central to modern molecular biology Elowitz et al. (2002); Pedraza and van Oudenaarden (2005). Much of this attention can be attributed to the development of new methodologies for measuring biochemical events. In particular, the use of green fluorescent protein and related fluorophores, high sensitivity light microscopy and flow cytometry have provided researchers with greater details of the dynamics of cellular processes at the level of single cells Elowitz et al. (2002); Pedraza and van Oudenaarden (2005); Xie et al. (2008) with single-molecule sensitivity Yu et al. (2006).

Quantitative laboratory measurements demand theoretical models that are able to describe and interpret the data. The traditional approach to modeling cellular processes has been to use deterministic kinetic equations with continuous dynamic variables: It is assumed that concentrations of metabolites, proteins, etc. vary deterministically in a continuous manner. This may be a reasonable assumption when one deals with molecules that occur in relatively large numbers. In other cases, the number of molecules can be much lower; for example, the number of LacI tetrameric repressor proteins in a single E. coli bacterium has been estimated to be of the order of 10 to 50 molecules Oehler et al. (1994). Such small numbers suggest that a continuous-variable model is no longer an appropriate description of the reality. In fact, recent experimental measurements have clearly highlighted the stochastic nature of biochemical processes in individual cells Elowitz et al. (2002); Choi et al. (2008); Qian (2012).

Stochastic reaction processes can be analyzed in terms of probability distributions from an isogenic population or from stochastic trajectories of individual cells. The former has been theoretically investigated in terms of the chemical master equation (CME) Gillespie (1992); McQuarrie (1967). For general reaction systems, analytically solving this equation is extremely challenging because the rates of biochemical reactions are often nonlinear functions of concentrations. In addition, numerical solutions are equally impractical because of the significant increase in the number of states; progress on the computational approach to the CME has been slow Liang and Qian (2010). Alternatively, stochastic trajectories can be computed via Monte Carlo procedures. Methods such as Gillespie algorithm Gillespie (1977) can also be highly intensive computationally for reaction systems involving several molecular species. Stochastic dynamics, therefore, are often studied with approximations Van Kampen (2001); Gómez-Uribe and Verghese (2007); Scott, Hwa, and Ingalls (2007); Munsky and Khammash (2006).

Gómez-Uribe et al. Gómez-Uribe and Verghese (2007) provided an approximate method called mass fluctuation kinetics. Recall that in a biological reaction system, the rate of a reaction is usually a nonlinear function of reactant concentrations, where is a vector of concentrations. In terms of cellular biological functions, is usually referred to as a ‘signal’, and we shall adopt this terminology. When the signal is noisy, the expected value of the nonlinear function, , is different from the function of the expectation of the concentrations, , as illustrated in Fig. 1 Kim and Price (2010). (Note that for a uni-molecular reaction, with being a constant. Then the mean rate is the same as the rate of the mean concentration: .) Gómez-Uribe et al. Gómez-Uribe and Verghese (2007) tried to quantify this difference in terms of both the the noise strength of a signal and the degree of reaction nonlinearity with a first-order approximation, assuming the noise strength is sufficiently small. We revisit this nonlinear signal processing of noise and investigate this problem more thoroughly.

Noise in biochemical signals following biological reaction pathways is dependent not only on the ‘local’ reaction rate functions Fell (1992); Kacser and Burns (1995); Fell (1996); Paulsson (2004); Sauro (2011) but also on the ‘global’ structure of the pathway that the signal ‘propagates’ through Paulsson (2004); Pedraza and van Oudenaarden (2005); Hooshangi, Thiberge, and Weiss (2005). The propagation of noise has been studied mainly by identifying the mathematical analytical structure Paulsson (2004); Tănase-Nicola, Warren, and ten Wolde (2006); Bruggeman, Blüthgen, and Westerhoff (2009). As the system size increases, such a structure quickly becomes intractable. Here, we provide a more manageable approach for understanding noise propagation and propose a simple principle we call stochastic focusing compensation. This principle states that a sensitivity increase (stochastic focusing) Paulsson, Berg, and Ehrenberg (2000) in one parameter region is typically achieved by reducing sensitivity in another region due to a geometric constraint. To show its wide applicability, we design three different synthetic biochemical networks that exhibit noise-induced phenotypes: (1) noise-enhanced concentration detection, (2) noise-induced bistable ON-OFF switching, and (3) noise-induced linear amplification. These examples illustrate the proposed method as a powerful tool for designing, analyzing, and understanding noise-induced phenotypes.

## Ii Results

### ii.1 Model systems.

When molecules undergo chemical reactions particularly at very low concentrations, changes in amounts are more conveniently described using discrete numbers to represent individual molecules Van Kampen (2001). The random nature of reaction events also dictates the changes being stochastic. Thus, the time evolution of a system is described by discrete-state continuous-time stochastic processes. Under appropriate conditions (spatially homogeneous and statistically independent reaction events), the stochastic processes are fully described by the CME Gillespie (1992), or more precisely the Delbrück-Gillespie process Qian (2011). The CME describes the time evolution of a probability distribution function, which represents the probability that one finds the number of molecules for each species at a given time.

### ii.2 Elasticity.

Our analysis is based on the concept of elasticity (), which quantifies the change in a reaction rate due to a change in a reactant, product or effector concentration. Mathematically, the sensitivity is defined as the ratio of the fold change in the rate to the fold change in the concentration :

 εd=sv∂v∂s=∂lnv∂lns.

has been widely used in metabolic control analysis Fell (1992); Kacser and Burns (1995); Fell (1996) and here we expand its application to the stochastic regime.

### ii.3 Elasticity vs. noise propagation.

The mean values of reaction rates can be affected by concentration fluctuations. This effect was formulated in terms of concentration variances and the curvatures of reaction rate functions Van Kampen (2001); Gómez-Uribe and Verghese (2007); Kim and Price (2010). Consider a reaction step shown in Fig. 1a. Although the concentration of , denoted by , is assumed to fluctuate ‘symmetrically’ with respect to its mean value, the reaction rate fluctuates ‘asymmetrically’ Ochab-Marcinek and Tabaka (2010). Let us assume that an increase in by the amount causes to increase by an amount . A decrease in by the same amount , however, will not cause to decrease by the same amount . Instead it is larger than , due to the down-curved shape of the reaction rate function (Fig. 1b). As the curvature of the reaction rate function increases, the fluctuations in become more asymmetric and the mean value of the rate function gets smaller, resulting in more deviation from the deterministic prediction . Here, the angle bracket denotes ensemble average. As the distribution function of gets narrower, the distribution function of also becomes narrower and the difference between and smaller. The first order correction to the deterministic rate was expressed proportional to both concentration variance (covariance in general cases) and the curvature of a rate function, and the estimate of the true mean rate was expressed Van Kampen (2001); Gómez-Uribe and Verghese (2007):

 w(s∗)=v(s∗)+12∂2v∂s2∣∣∣s=s∗Σ∗, (1)

where is the estimate of the true mean concentration and that of concentration variance (Methods). The variance in is dependent on how is synthesized and degraded as well as other sources of random fluctuations from the connected networks. Thus, the interplay between the nonlinearity in the rate function and the noise in the concentration causes a change in the mean reaction rate, resulting in changes in elasticities. Since elasticities together with stoichiometry are the determinants of dynamic behavior, changes in the elasticities due to noise will result in changes in system dynamics. This result allows us to relate deterministic properties of a network to the nonlinear transmission of noise. We will define elasticities that depend on noise propagation as:

 εs=∂ln⟨v⟩∂ln⟨s⟩.

and its first-order approximation as .

This approximation can be used to derive analytical solutions of mean values and covariances of concentrations, as a next level of the linear noise approximation Van Kampen (2001), but needs to be verified on a case by case basis; for example, the level of in Fig. 3a was assumed to be constant but it can be a function of other species concentrations (e.g., Fig. S2c), which could result in additional noise being propagated into the synthesis rate.

### ii.4 Stochastic Focusing Compensation.

Consider the curvature of in Fig. 1c and 2a which represent an inhibitory regulation. The curvature is positive for all except when or . Thus, from Eq. (1) the mean reaction rate is shown to be larger than the deterministic rate . When the curvature approaches zero and when the variance of vanishes. Thus, in these limits becomes equal to . Inspection of Fig. 1c and 2a shows that the slope of the mean rate ( or ) is smaller than that of when is small, implying . This elasticity decrease due to noise will be called stochastic de-focusing (SD) (Fig. 2b). However, as increases the slope of the mean rate becomes larger than , implying that . This sensitivity increase will be called stochastic focusing (SF) (Fig. 2b). A minor difference from the original definition of SF Paulsson, Berg, and Ehrenberg (2000) is discussed in the Appendix.

The previous discussion shows that elasticities that increase due to noise in one parameter region will typically come with elasticity decreases in another region. Consider the case that satisfies the Hill equation with Hill coefficient 4 (Fig. 2c). The curvature of changes from positive to negative as increases from zero, and the curvature vanishes as or . Thus, defocusing appears in between two focused regions (Fig. 2d). In general, focusing does not, however, always come with defocusing. For example, when is a Michaelis-Menten-type rate function, only focusing appears as shown in Fig. 2e and f.

When satisfies the Poisson distribution, focusing and defocusing can occur strongly. With a different type of distribution such as the Gamma distribution, these can appear even more strongly (Fig. 2). The Gamma distribution has been frequently found in the case when translation events occur in a burst manner Taniguchi et al. (2010). One of the important differences between the Poisson and Gamma distributions is that in the Gamma distribution one can control the noise (coefficient of variation; CV) and mean levels independently, while one cannot do so in the Poisson distribution Kim and Sauro (2012). Furthermore, the Poisson distribution is the asymptotic case of the Gamma distribution as the burst size in the translation events reduces to zero. These two facts imply that the noise strength (CV) increases from that of the Poisson case by turning on the bursting events, resulting in a larger deviation from the deterministic case based on Eq. (1). This is the reason that stronger focusing and defocusing were observed for the case of the Gamma distribution as shown in Fig. 2.

In this section, we have shown that higher elasticities in one region of parameters can be achieved by reducing elasticities in another region, and have explained this stochastic focusing compensation by using the curvature-variance correction to a mean reaction rate (Eq. (1)). The compensation was intuitively understood as a geometric effect of the reaction rate. In the following sections we will show how this effect can be used to modify the behavior in gene regulatory networks.

### ii.5 Application 1: Noise-improved Concentration Detector.

Let us consider an incoherent feed-forward gene regulatory network that functions as a concentration detector Entus, Aufderheide, and Sauro (2007); Kaplan et al. (2008). We can use noise to enhance the detection amplitude and sensitivity of this system. Consider a reaction process where a protein ( in Fig. 3a) is regulated by two different pathways: either via direct transcriptional activation by or indirect inhibition by . When the concentration of is zero, is not expressed. As increases, increases together (this region of will be named OFFON as shown in Fig. 3b). When becomes larger than a thresh-hold point, it begins, however, to decrease and is eventually dominated by ’s inhibition (this region of will be named ONOFF as shown in Fig. 3b). Thus, one can detect a specific range of the concentration of by monitoring the concentration of . The plot of vs. (Fig. 3b) will be used to characterize the concentration detection by using the slope as sensitivity, its width at the half maximum as bandwidth, and the peak value of as amplitude.

We will apply stochastic focusing compensation to enhance the detection amplitude and sensitivity. In this system, noise in the concentration levels of and is generated from a series of stochastic reaction events of synthesis and degradation of and . Since regulates the synthesis of , noise in the concentration level of introduces additional noise into , i.e., the noise propagates. Based on the fact that the regulation is inhibitory, we can conclude that both strong SF and SD can appear in the ONOFF region (refer to Fig. 2a,b). For simplicity, here we assumed that the concentration of does not fluctuate stochastically and thus the activation by does not have any noise effect. Therefore, we focus on the inhibition pathway and aim to increase the detection sensitivity in the ONOFF region.

For this aim, strong SF is required to appear in the ONOFF region, while minimizing the appearance of SD in the same region. In the hyperbolic inhibition case (Fig. 2a), strong SF was achieved for the value of  nM and for that of between 1 and 10 nM (for the Poisson case in Fig. 2a,b) SI . Thus, we tuned the values of the terms corresponding to , , to less than 0.1 nM and the mean value of , , between 1 and 10 nM. These two conditions provided the possible range of : , where the function MAX chooses the larger value of its parameter arguments.

The above region was placed to the right side of the detection peak but not too far away from the peak to maximize the appearance of strong SF and to minimize the effect of strong SD in the region. Then, concentration detection was enhanced 5 times in the ONOFF sensitivity and 7 times in the amplitude (Fig. 3b) due to the strong SF SI . However, the bandwidth stayed the same (arrows in Fig. 3b), which is related to the stochastic focusing compensation; SD appears as decreases away from the region of SF, resulting in reduced detection sensitivity and eventually preventing the detection bandwidth from decreasing.

To enhance the detection amplitude and sensitivity further, we used the Gamma distribution to allow larger noise levels in . As the burst size – number of translation events from a single mRNA during its average lifetime – increases, the amplitude of the detection is significantly enhanced: 20 times for , 83 times for compared to the deterministic case. However, when we allow to fluctuate stochastically, its signal noise turned out to slightly decrease the amplitude due to the positive curvature contribution in the activation by as well as the noise correlation between and SI . This further analysis shows that our curvature-variance-based analysis can successfully predict the change in the system properties.

### ii.6 Application 2: Noise-induced Bistable ON-OFF switch.

Bistability is a common feature of many biological systems Ozbudak et al. (2004); Brandman et al. (2005). In the stochastic framework, bistability can sometimes be observed as a bimodal distribution because the system can jump from one stable state to the other. In this section we will show how a non-bistable system can be made bistable through stochastic focusing compensation. This example shows the predictive power of our method in system design and furthermore provides a novel mechanism of noise-induced bistability.

We consider the model of mutual inhibition between two genes (g1 and g2 in Fig. 4a), which regulate each other non-cooperatively (Hill coefficient = 1). Under these conditions it is impossible for the system to show bistability with first-order degradation in the deterministic case Sauro (2009). Now we consider including noise in the system and show noise-induced bistability. The dissociation constant of , , was varied while that of , , was fixed at  nM. The reason that this value was chosen is to achieve strong SF and SD in the synthesis rate of , . For  0.1 nM, another strong SF and SD in appears. This doublet of SF-SD can lead to ultra-sensitive positive feedback (Fig. 4c), leading to bimodality (Fig. 4b). As increases, the SF-SD in becomes weaker and the bimodal distribution becomes uni-modal.

To understand how the doublet of SF-SD causes bimodality, the nullclines of and were plotted in Fig. 4d. The nullclines were computed using the Poisson distributions for and as a first level of approximation. In Fig. 4d, the solid line (NC1) corresponds to with  nM and was significantly distorted due to the strong SF-SD (cf. Fig. 2a). With  nM, the other nullcline (NC2) corresponding to was also significantly distorted. These two distorted curves, NC1 and NC2, crossed each other by making a shape of . As increases, NC2 became straightened out due to weaker SF-SD and shifted upward. Thus the bimodality disappeared.

The above analysis is limited to unimodal distributions and thus significant care may need to be taken when the distribution deviates from unimodality and enters into bimodality. Therefore, we checked the numerical accuracy of our analysis on this system. The phase diagram for bimodality and unimodality were plotted for different values of and . The predicted phase diagram (gray and white areas in Fig. 4e) matched well with the results based on exact stochastic simulations (dots in Fig. 4e). This shows that the noise-induced bimodality appears due to the interplay between the system nonlinearity and stochasticity that are described by the CME description and that our method based on stochastic focusing compensation can be used for predicting bimodality in a rational way.

### ii.7 Application 3: Noise-induced Linear Amplifier.

A linear amplifier is a device that transfers an input signal to an output without distorting the signal shape, resulting in reliable information transfer. Such devices have been observed in biological systems, for example, in signal transduction pathways in yeast cells (Saccharomyces cerevisiae) Yu et al. (2008) and mammalian cells Sturm et al. (2010). This linear amplification has been attributed to negative feedback mechanisms Sauro and Ingalls (2007); Sturm et al. (2010). Here we present a novel way to achieve linear amplification by exploiting noise without resorting to any feedback.

The SF-SD makes a sigmoidal curve less steep as illustrated in Fig. 2c. The contribution becomes larger as the strength of noise increases and the SF-SD can appear more strongly (Eq. (1)), i.e., a transfer curve can become less sigmoidal with the increase of noise strength. Based on this prediction, we investigated how gene expression noise can be leveraged to build a linear amplifier without using any feedback.

We considered two different activation reaction systems (Fig. 5a): a Hill equation with Hill coefficient = 12 (Fig. 5b) and a piece-wise linear function (Fig. 5c). An activator is assumed to satisfy the negative binomial distribution – a discrete version of the Gamma distribution:

 Ps=Γ(a+s)Γ(s+1)Γ(a)(b1+b)s(1−b1+b)a (2)

where is a burst frequency (number of transcription events during protein lifetime), ranging between 0.1 and 11 in E. coli with a YFP-fusion tag Taniguchi et al. (2010), and a burst size, ranging between 1 and 900 Taniguchi et al. (2010). Here the mean of the distribution is equal to and the strength of noise (standard deviation divided by mean) is for . We varied the strength of noise by increasing (more precisely for a fixed mean value, increasing corresponds to decreasing , resulting in an increase in the noise strength) and observed the change in . It turned out that became dramatically linearized as increased (Fig. 5b; cf. Fig. 3 in Thattai and van Oudenaarden (2002)). More surprisingly, this linearized response was observed not only for the Hill function-type activation with Hill coefficients ranging from 1 to 20 SI , but also holds for various types of activation patterns SI including piece-wise linear functions (Fig. 5c).

To confirm this noise-induced linear amplification, stochastic simulations were performed for a three-gene cascade (Fig. 5d) by using the Gillespie stochastic simulation algorithm, with transcription processes included. We achieved a linear response for the transfer curve for vs. (Fig. 5e). This result implies that gene expression noise generated upstream can linearize downstream nonlinear transcriptional responses. This leads to the non-intuitive result that signals upstream can be transferred reliably without distortion as a result of noise. This noise-induced linearization may be observed in protein signaling pathways as well if the protein concentration distributions are wide Birtwistle et al. (2012a); Thattai and van Oudenaarden (2002). This could serve as an alternative mechanism for dose-response alignment phenomenon Yu et al. (2008); Sturm et al. (2010) without strongly resorting to negative feedback.

## Iii Discussion

We analyzed how noise propagation, in the steady state of stochastic biochemical reaction systems Paulsson (2004); Pedraza and van Oudenaarden (2005); Hooshangi, Thiberge, and Weiss (2005), affects the population average of reaction rates and thus their elasticities (sensitivities). The change in the elasticity was explained in terms of a curvature-variance contribution, which led us to discover stochastic focusing compensation: A (stochastic) elasticity increase in one parameter region due to noise can lead to an elasticity decrease in another region.

This stochastic focusing compensation was used to enhance the functional properties of a concentration detector having an incoherent feedforward network structure. We designed strong stochastic focusing to a desired detection region and improved the detection amplitude and sensitivity. The detection bandwidth was, however, not enhanced due to stochastic defocusing In addition, we applied the same idea of stochastic focusing compensation in designing a noise-induced bistable ON-OFF switch that has a non-cooperative positive feedback loop (thus showing a single stable state in the absence of the noise). The underlying mechanisms of the emergent bimodality were related to the doublet of stochastic focusing-defocusing that appear in two individual inhibitory interactions, with the two as a whole comprising a positive feedback loop. Lastly, the stochastic focusing compensation principle was also applied to the design of a noise-induced linear amplifier. Surprisingly, the biologically-relevant distribution of proteins (Gamma distribution) were shown to have important properties that linearize a wide range of sigmoidal responses.

In our linear amplification study, it seems paradoxical that greater noise enables more reliable information transfer. However, the fact that the linearization becomes saturated as the noise strength (CV) increases (refer to SI ) implies that there will be an optimal level of noise satisfying optimal signal transfer much like stochastic resonance Gammaitoni et al. (1998) but with a completely different mechanism.

Experimental evidence of the Gamma distribution of intracellular protein copy numbers Paulsson and Ehrenberg (2000); Friedman, Cai, and Xie (2006); Shahrezaei and Swain (2008); Taniguchi et al. (2010); Birtwistle et al. (2012b) implies that when measuring the Hill coefficient of transcriptional regulation, significant care may need to be taken to include this stochastic focusing compensation. Furthermore, in signal transduction networks, cell-to-cell variability in receptor occupancies and intra-cellular enzyme copy numbers Birtwistle et al. (2012a, b) could lead to lowering the observed Hill coefficients of downstream regulation and contributing to dose-response alignment. It would be important to investigate the degree of individual contributions of this noise effect (SF-SD) and negative feedback loops Yu et al. (2008).

## Iv Conclusions

Noise propagation can have significant effects on cellular phenotypes Arkin, Ross, and McAdams (1998); Maamar, Raj, and Dubnau (2007); Wang and Zhang (2011). It is important to be able to explain propagation mechanisms and to use them for designing or adjusting phenotypes Kim and Sauro (2012). We provide a manageable, yet fully quantitative approach for studying noise propagation and its effect, in terms of the degree of a system’s nonlinearity and the strength of noise. We show that an increased sensitivity in one parameter region often comes with a decreased sensitivity in another parameter region. We applied this sensitivity compensation characteristics in designing three gene regulatory networks: a concentration detector, a noise-induced ON-OFF switch, and a noise-induced linear amplifier. These examples illustrate how our sensitivity-based approach can be used to design, analyze, and optimize the functionalities of stochastic reaction networks by exploiting noise.

## acknowledgments

We acknowledge the support from the National Science Foundation (Theoretical Biology 0827592 and Molecular and Cellular Biosciences 1158573; for Preliminary studies, FIBR 0527023).

## Appendix A Model Systems

Random chemical reactions are often modeled by stochastic processes governed by the chemical master equation (CME) Gillespie (1992), under the assumption that the reactions occur uniformly in position space and independently. This equation describes the time evolution of a probability distribution function of reactant numbers for all species at a given time.

We consider species of molecules involved in reactions:

 nl1S1+⋯+nlmSmVl−→ml1S1+⋯+mlmSm,

where the molecule numbers are denoted by with and the rate of reaction by with . We also assume that the rate of reactions , denoting {}, can be controlled by changing a set of system parameters . The parameters are non-fluctuating (bold symbols represent matrices and vectors). The number change in a species by a single event of the above reaction is described by a reduced stoichiometry matrix: . The molecule numbers evolve stochastically and their evolution can be described by the CME:

 ∂P(\boldmathS,t)∂t = n∑j=1[P({Si−NRij},t)Vj({Si−NRij,\boldmathp}) (3) −P(\boldmathS,t)Vj(\boldmathS,% \boldmathp)],

where denotes with , and is a probability distribution function of reactant numbers for all species at a given time.

## Appendix B Time evolution of concentration mean values and covariances

We switch the representation of states from numbers to concentrations of molecules with a system volume, since this concentration representation has a direct correspondence to deterministic macroscopic kinetics.

The mean concentration of a species is given as

 ⟨si⟩t≡[m∏j=1∫∞sj=0dsj]siP(\boldmaths,t),

where is the probability distribution function of molecule concentrations, , at a given time . We will use the angle bracket to represent the ensemble average at a given time from now on. The evolution of the mean concentration of a species is governed by the following equation SI :

 d⟨\boldmaths⟩tdt=\boldmathNR⟨\boldmathv(\boldmaths,\boldmathp)⟩t, (4)

where is a reduced stoichiometry matrix Reder (1988); Sauro and Ingalls (2004). is the number of linearly independent rows of a stoichiometry matrix. represents rate functions: with . We note that is a function of with and .

A concentration covariance between two species ( and ) is defined as

 Σtij≡⟨(si−⟨si⟩t)(sj−⟨sj⟩t)⟩t.

The angle brackets again represents ensemble averages at a given time . The evolution of the covariance matrix is described by the following equation SI :

 d\boldmathΣtdt=⟨(\boldmathNR\boldmathv)(\boldmaths−⟨\boldmaths⟩t)+(\boldmaths−⟨\boldmaths⟩t)T(% \boldmathNR\boldmathv)T+\boldmathDΩ⟩t, (5)

where the diffusion coefficient matrix is defined by with a diagonal matrix .

It is difficult to solve equations (4) and (5) unless the rate function is linear with . Under the assumption that the third and higher order moments of are negligible, the time evolution of concentration mean values and covariances can be described by the following equations Gómez-Uribe and Verghese (2007) SI :

 d⟨\boldmaths⟩tdt = \boldmathNR(\boldmathv(⟨% \boldmaths⟩t,\boldmathp)+∑ij12∂2\boldmathv(⟨\boldmaths⟩t,% \boldmathp)∂⟨si⟩t∂⟨sj⟩tΣtij), (6)
 d\boldmathΣtdt = \boldmathJ\boldmathΣt+\boldmathΣTt\boldmathJT+\boldmathD(⟨%\boldmath$s$⟩t)Ω, (7)

where is an element of the Jacobian matrix defined as

 Jij≡n∑k=1NRik∂vk(⟨% \boldmaths⟩t,\boldmathp)∂⟨sj⟩t.

The validity of the above equations should be checked on the basis of each different model, since noise in adjacent systems can be propagated into the rate function Paulsson (2004); Pedraza and van Oudenaarden (2005); Hooshangi, Thiberge, and Weiss (2005). We note that Eq. (7) is different from the results of Gómez-Uribe, et al. Gómez-Uribe and Verghese (2007): We have neglected all the terms of the order of that have been kept in Eq. 11 of Gómez-Uribe, et al. Gómez-Uribe and Verghese (2007). This is consistent within the approximation of truncation of third and higher order moments.

## Appendix C Mean reaction rates and covariance matrices at the stationary state

From Eq. (6), the mean rate of reaction can be approximated as:

 \boldmathw≡\boldmathv(⟨\boldmaths⟩,\boldmathp)+∑ij12∂2% \boldmathv∂si∂sj∣∣∣\boldmaths=⟨\boldmaths⟩Σtij. (8)

At the stationary state, Eqs. (6) and (7) becomes:

 \boldmathNR\boldmathw=0
 \boldmathJ\boldmathΣ+\boldmathΣT% \boldmathJT+\boldmathD(⟨\boldmaths⟩)Ω=0. (9)

By solving the above two equations one can estimate the mean values of the concentrations and their covariance matrix.

## Appendix D Definition of stochastic focusing

Stochastic focusing (SF) that we have discussed here has minor conceptual differences from the one described in Paulsson et al. Paulsson, Berg, and Ehrenberg (2000). Consider the following reactions:

 cs  −−−−→Sγss−−→\o,  →v(s)=1/(KM+s)Pγpp−−→\o, (10)

with a constant synthesis rate of and a synthesis rate of . The SF that we could discuss for the above system is independent of how fast the concentration fluctuates. In Paulsson, Berg, and Ehrenberg (2000), SF is however the only effect of rapid fluctuations in concentrations. This difference is due to the fact that the focus in Paulsson et al. Paulsson, Berg, and Ehrenberg (2000) was on what is the most probable state of , and we focus on the mean value of , more precisely since the flux balance at the stationary state leads to , resulting in being proportional to . To understand this difference, we need to understand the dynamics of that is correlated with the fluctuation in . When all reaction rates , , , and are in the same order of magnitude, and fluctuate on similar time scales. If hits zero, the inhibition on is removed and thus the number of can rapidly increase to a very large number. When increases to , however, the inhibition on appears and the number of rapidly decreases. Therefore, the time series profile of shows a flat lower bound at zero with many large sharp spikes. This time series profile changes as the synthesis and degradation of become faster. If the parameter values of and increase such that fluctuates much faster than , sees the averaged behavior of and is unlikely to hit zero. Thus, the strong/weak inhibition by is averaged out. Due to this averaging, the sensitivity enhancement, i.e., SF, manifests itself. This is why the SF was claimed to be the effect of a rapidly fluctuating Paulsson, Berg, and Ehrenberg (2000). However if one can observe the time series profile for a sufficiently long time to obtain good statistics of the spike heights, the mean value of becomes actually independent from how fast and are, if the ratio is presumed to be kept constant. This means that the mean rates of and also becomes time-scale independent. This is why the SF defined here becomes time-scale independent.

## References

• Elowitz et al. (2002) M. B. Elowitz, A. J. Levine, E. D. Siggia,  and P. S. Swain, Science 297, 1183 (2002).
• Pedraza and van Oudenaarden (2005) J. M. Pedraza and A. van Oudenaarden, Science 307, 1965 (2005).
• Xie et al. (2008) X. S. Xie, P. J. Choi, G.-W. Li, N. K. Lee,  and G. Lia, Annu. Rev. Biophys. 37, 417 (2008).
• Yu et al. (2006) J. Yu, J. Xiao, X. Ren, K. Lao,  and X. S. Xie, Science 311, 1600 (2006).
• Oehler et al. (1994) S. Oehler, M. Amouyal, P. Kolkhof, B. von Wilcken-Bergmann,  and B. Müller-Hill, EMBO J. 13, 3348 (1994).
• Choi et al. (2008) P. J. Choi, L. Cai, K. Frieda,  and X. S. Xie, Science 322, 442 (2008).
• Qian (2012) H. Qian, Annu. Rev. Biophys. 41, 179 (2012).
• Gillespie (1992) D. T. Gillespie, Physica A 188, 404 (1992).
• McQuarrie (1967) D. A. McQuarrie, J. Appl. Probab. 4, 413 (1967).
• Liang and Qian (2010) J. Liang and H. Qian, J. Comput. Sci. Technol. 25, 154 (2010).
• Gillespie (1977) D. T. Gillespie, J. Phys. Chem. 81, 2340 (1977).
• Van Kampen (2001) N. G. Van Kampen, Stochastic Processes in Physics and Chemistry., 3rd ed. (North Holland, 2001).
• Gómez-Uribe and Verghese (2007) C. A. Gómez-Uribe and G. C. Verghese, J. Chem. Phys. 126, 24109 (2007).
• Scott, Hwa, and Ingalls (2007) M. Scott, T. Hwa,  and B. Ingalls, Proc. Natl. Acad. Sci. U.S.A. 104, 7402 (2007).
• Munsky and Khammash (2006) B. Munsky and M. Khammash, J. Chem. Phys. 124, 44104 (2006).
• Kim and Price (2010) P.-J. Kim and N. D. Price, Phys. Rev. Lett. 104, 9 (2010).
• Fell (1992) D. A. Fell, Biochem. J. 286, 313 (1992).
• Kacser and Burns (1995) H. Kacser and J. A. Burns, Biochem. Soc. Trans. 23, 341 (1995).
• Fell (1996) D. A. Fell, Understanding the Control of Metabolism. (London, Portland Press, 1996).
• Paulsson (2004) J. Paulsson, Nature 427, 415 (2004).
• Sauro (2011) H. M. Sauro, Enzyme Kinetics for Systems Biology, 2nd ed. (Ambrosius Publishing and Future Skill Software, 2012).
• Hooshangi, Thiberge, and Weiss (2005) S. Hooshangi, S. Thiberge,  and R. Weiss, Proc. Natl. Acad. Sci. U.S.A. 102, 3581 (2005).
• Tănase-Nicola, Warren, and ten Wolde (2006) S. Tănase-Nicola, P. B. Warren,  and P. R. ten Wolde, Phys. Rev. Lett. 97, 68102 (2006).
• Bruggeman, Blüthgen, and Westerhoff (2009) F. Bruggeman, N. Blüthgen,  and H. Westerhoff, PLoS Comput. Biol. 5, e1000506 (2009).
• Paulsson, Berg, and Ehrenberg (2000) J. Paulsson, O. G. Berg,  and M. Ehrenberg, Proc. Natl. Acad. Sci. U. S. A. 97, 7148 (2000).
• Qian (2011) H. Qian, Nonlinearity 24, R19 (2011).
• Ochab-Marcinek and Tabaka (2010) A. Ochab-Marcinek and M. Tabaka, Proc. Natl. Acad. Sci. U. S. A. 107, 22096 (2010).
• Taniguchi et al. (2010) Y. Taniguchi, P. J. Choi, G.-W. Li, H. Chen, M. Babu, J. Hearn, A. Emili,  and X. S. Xie, Science 329, 533 (2010).
• Kim and Sauro (2012) K. H. Kim and H. M. Sauro, PLoS Comput. Biol. 8, e1002344 (2012).
• Entus, Aufderheide, and Sauro (2007) R. Entus, B. Aufderheide,  and H. M. Sauro, Syst. Synth. Biol. 1, 119 (2007).
• Kaplan et al. (2008) S. Kaplan, A. Bren, E. Dekel,  and U. Alon, Mol. Syst. Biol. 4, 203 (2008).
• Ozbudak et al. (2004) E. M. Ozbudak, M. Thattai, H. N. Lim, B. I. Shraiman,  and A. Van Oudenaarden, Nature 427, 737 (2004).
• Brandman et al. (2005) O. Brandman, J. E. Ferrell, Jr., R. Li,  and T. Meyer, Science 310, 496 (2005).
• Sauro (2009) H. M. Sauro, in Computational Systems Biology (Humana Press, New York, 2009) Chap. 13, pp. 269–309.
• Yu et al. (2008) R. C. Yu, C. G. Pesce, A. Colman-Lerner, L. Lok, D. Pincus, E. Serra, M. Holl, K. Benjamin, A. Gordon,  and R. Brent, Nature 456, 755 (2008).
• Sturm et al. (2010) O. E. Sturm, R. Orton, J. Grindlay, M. Birtwistle, V. Vyshemirsky, D. Gilbert, M. Calder, A. Pitt, B. Kholodenko,  and W. Kolch, Sci. Signal. 3, ra90 (2010).
• Sauro and Ingalls (2007) H. M. Sauro and B. Ingalls, arXiv:0710.5195[q-bio.MN]  (2007) .
• Thattai and van Oudenaarden (2002) M. Thattai and A. van Oudenaarden, Biophys. J. 82, 2943 (2002).
• Birtwistle et al. (2012a) M. R. Birtwistle, J. Rauch, A. Kiyatkin, E. Aksamitiene, M. Dobrzyński, J. B. Hoek, W. Kolch, B. A. Ogunnaike,  and B. N. Kholodenko, BMC Syst. Biol. 6, 109 (2012a).
• Gammaitoni et al. (1998) L. Gammaitoni, P. Hänggi, P. Jung,  and F. Marchesoni, Rev. Mod. Phys. 70, 223 (1998).
• Paulsson and Ehrenberg (2000) J. Paulsson and M. Ehrenberg, Phys. Rev. Lett. 84, 5447 (2000).
• Friedman, Cai, and Xie (2006) N. Friedman, L. Cai,  and X. Xie, Phys. Rev. Lett. 97, 1 (2006).
• Shahrezaei and Swain (2008) V. Shahrezaei and P. S. Swain, Proc. Natl. Acad. Sci. U. S. A. 105, 17256 (2008).
• Birtwistle et al. (2012b) M. R. Birtwistle, A. von Kriegsheim, M. Dobrzyński, B. N. Kholodenko,  and W. Kolch, Mol. Biosyst. 8, 3068 (2012b).
• Arkin, Ross, and McAdams (1998) A. P. Arkin, J. Ross,  and H. H. McAdams, Genetics 149, 1633 (1998).
• Maamar, Raj, and Dubnau (2007) H. Maamar, A. Raj,  and D. Dubnau, Science 317, 526 (2007).
• Wang and Zhang (2011) Z. Wang and J. Zhang, Proc. Natl. Acad. Sci. U. S. A. 108, E67 (2011).
• Reder (1988) C. Reder, J. Theor. Biol. 135, 175 (1988).
• Sauro and Ingalls (2004) H. M. Sauro and B. Ingalls, Biophys. Chem. 109, 1 (2004).
• Sauro and Fell (2000) H. M. Sauro and D. A. Fell, in Animating the Cellular Map: Proceedings of the 9th International Meeting on BioThermoKinetics (Stellenbosch University Press, 2000) pp. 221–228.
• R Development Core Team (2008) R Development Core Team, R: A language and environment for statistical computing (R Foundation for Statistical Computing, Vienna, Austria, 2008).
• () J. Quinn, J. Beal, S. Bhatia, P. Cai, J. Chen, K. Clancy, N. Hillson, M. Galdzicki, A. Maheshwari, P Umesh, M. Pocock, C. Rodriguez, G.-B. Stan, D. Endy, BioBricks Foundation RFP 93  (2013).
• () See Supplementary Material Document No.xxxxxx for more detailed discussions on the derivation of Eqs. (8) and (9) and other numerical results. For information on Supplementary Material, see http://www.aip.org/pubservs/epaps.html.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters