pmc.ncbi.nlm.nih.gov

Deterministic characterization of stochastic genetic circuits

Abstract

For cellular biochemical reaction systems where the numbers of molecules is small, significant noise is associated with chemical reaction events. This molecular noise can give rise to behavior that is very different from the predictions of deterministic rate equation models. Unfortunately, there are few analytic methods for examining the qualitative behavior of stochastic systems. Here we describe such a method that extends deterministic analysis to include leading-order corrections due to the molecular noise. The method allows the steady-state behavior of the stochastic model to be easily computed, facilitates the mapping of stability phase diagrams that include stochastic effects, and reveals how model parameters affect noise susceptibility in a manner not accessible to numerical simulation. By way of illustration we consider two genetic circuits: a bistable positive-feedback loop and a negative-feedback oscillator. We find in the positive feedback circuit that translational activation leads to a far more stable system than transcriptional control. Conversely, in a negative-feedback loop triggered by a positive-feedback switch, the stochasticity of transcriptional control is harnessed to generate reproducible oscillations.

Keywords: intrinsic noise, phase diagram, synthetic biology, oscillations

Parallel advances in the conceptual understanding of gene regulation along with technological advances in molecular biology have given rise to the possibility of system-level quantitative kinetic measurements of living organisms (1) and synthetic genetic circuit designs (2, 3). Interpretation of time-series data from complex networks and reliable forward-design of gene circuits depend on detailed quantitative mathematical models (2, 4). These models generally take one of two largely exclusive forms: either deterministic formulations with reactant concentration varying continuously in time and governed by a system of rate equations or stochastic formulations that explicitly include the discrete and probabilistic change in reactant molecule numbers as each subsequent reaction occurs (5). Both approaches have benefits and associated limitations.

The great practical advantage of rate equation models is the ease with which the qualitative behavior of the system can be extracted. By focusing on the long-term behavior, the model dynamics are simplified, and one is able to gain insight into the expected response of the system (6). Rate equation models, however, neglect the fact that chemical reaction networks are composed of species that evolve on discrete space, jumping from some number of molecules to another as each reaction occurs (7). The resulting deviation from the deterministic formulation is called the “intrinsic noise” in the system (because the fluctuations arise from the reaction dynamics themselves and not from some external source) (8, 9). In cellular systems with small numbers of reactant molecules, the relative magnitude of the intrinsic noise can be large, and it can give rise to qualitatively different behavior than what rate equation models would predict. A system that has several possible stable states, for example, may be induced to spontaneous transitions between them as a result of intrinsic noise (10, 11), leading to a stochastic switching of states. In an excitable system, noise may cause oscillations to occur in a model that is otherwise stable (12–14). With a given set of physical parameters, it is possible to simulate explicitly the individual chemical reaction events, including the effect of intrinsic noise (5). Nevertheless, the design of synthetic circuits, or therapeutics aimed at altering an existing network, require knowledge of the “phase diagram,” which involves a systematic mapping of the parameter space. There, stochastic simulation becomes prohibitively time-consuming even for reasonably simple genetic circuits involving two or three genes (see below), and analytical methods are needed.

A number of analytic studies have been done recently to model intrinsic noise in genetic circuits, much of them built on the linear noise approximation (15) and focused on the noise property itself, e.g., “noise propagation” through genetic networks (16, 17), the equilibrium distribution of fluctuations about multiple steady-states (18), and constructive effects of noise in signal processing (13, 19). There has been comparatively little work, however, aimed at providing tools to study the effect of intrinsic noise on the stability of systems where stochastic models exhibit qualitatively different behavior from their deterministic counterparts (20). Under these conditions, the linear noise approximation alone cannot predict qualitative changes in the observable dynamics of the system, as, for example, in the case of noise-induced oscillations (21). Here we present an analytic method, which we call the “effective stability approximation” (ESA), that extends the applicability of existing deterministic methods to include stochastic effects. The method is an extension of the linear noise approximation, including correction of stochasticity to the deterministic equations to the order 1/N (where N is the number of molecules in the system). It conveniently connects deterministic and stochastic descriptions, allowing for systematic exploration of parameter space while at the same time including the essential effect of intrinsic fluctuations. For the two model systems examined here, we find the ESA to capture reliably the essential features of those systems, correctly estimating the effect of intrinsic noise on the phase diagrams of systems dominated by as little as a few dozen molecules.

ESA can be applied to generic models of genetic circuits, and a brief tutorial is presented in Methods with the hope that the approach can be used by other investigators to include stochastic effects in deterministic models. The full mathematical details are presented in supporting information (SI). We illustrate the power of the method below by considering two examples: an autoregulator with positive feedback (an “autoactivator”) (22) and an excitable genetic oscillator linking positive and negative feedback loops (12, 23). The behavior of both circuits is conveniently visualized by means of a phase diagram that cannot be practically constructed by using numerical simulations if stochastic effects are to be included. Furthermore, the analysis reveals that the system behavior is completely governed by a few dimensionless combinations of model parameters: combinations that would be very difficult to infer from simulation data alone. We hope that our presentation of the ESA method will make it accessible to modelers, bioengineers, and synthetic circuit designers for the analysis of various molecular circuits, while our description of the behaviors of the two model systems will provide quantitative-minded biologists with a concrete sense of the effect of stochasticity as well as a succinct means of characterization (e.g., a phase diagram with reduced variables).

Results and Discussion

Autoactivator.

Perhaps the simplest circuit motif able to exhibit multiple stable states is the autoactivating positive-feedback loop (Fig. 1A) (24). The circuit consists of a single gene encoding an activator. Several autoactivator circuits have been experimentally characterized, including the autoactivation of CI protein by the PRM promoter of phage λ studied by Isaacs et al. (22) and the autoactivation of NtrC by the glnAp promoter of Escherichia coli studied by Atkinson et al. (23). The autoactivator circuit is expected to exhibit either a HIGH state, characterized by an elevated level of protein synthesis, or a LOW state, characterized by a low basal level of production. We simplify the model by assuming that the activator binding and mRNA turnover are fast compared with the lifetime of the protein activator. The effect of the activator is quantified by the “activation function” g(A/KA, f), where A is the activator concentration, KA is the equilibrium dissociation constant of the activator and its cognate binding site, and f is the maximum fold-activation in the circuit. As a particular example, we assume a Hill-form for the activation function g(A/KA, f),

graphic file with name zpq01807-6053-m01.jpg

with cooperative activation (n = 2) (4). The resulting model is a single kinetic equation governing the activator concentration A(t) (22, 25) (Fig. 1A),

graphic file with name zpq01807-6053-m02.jpg

where γ is the fully activated rate of protein synthesis and δ is the protein degradation rate (which in prokaryotes is often estimated from the growth rate due to growth-mediated dilution).

Fig. 1.

Fig. 1.

Two example circuit motifs. (A) A positive-feedback loop capable of maintaining two stable states (22). (B) An excitable oscillator that exhibits noise-induced oscillations (12, 23). The autoactivator triggers the production of a repressor R that provides negative feedback control. Dashed arrows, lumped transcription and translation; bold filled arrows, activation; blunt arrow, repression; wavy arrows, degradation.

In the deterministic limit, when the number of reactant molecules is very large, we expect Eq. 2 to adequately describe the system behavior. Once initial transients have died out, the system will approach a steady state, and A reaches its steady-state value As where the rate of synthesis and degradation balance, i.e., γ·g(As) = δ·As. The stability of the steady state is determined by the response of the system to a small perturbation Ap, found by linearizing Eq. 2 about As,

graphic file with name zpq01807-6053-m03.jpg

The expression in the square brackets λ ≡ [γ·g′(As) − δ] is a constant that depends on the model parameters. If λ is positive, the small perturbations will grow in time (As is an “unstable state”), whereas if λ is negative, the small perturbation will decay (As is a “stable state”). In the stable case, the long-term state of the system can be thought of as a point located at the bottom of a valley (or basin of attraction): the more negative the constant λ, the steeper the valley. As the model parameters are varied, the valley may become more flat (λ ≈ 0) or even develop into a mountain (λ > 0), resulting in a loss of stability. The parameter space is divided into regions of different qualitative behavior (as in Fig. 2A, black curve); the threshold between these domains indicates where λ has changed sign and is called the “phase boundary.” Although the model seems to depend on a large family of parameters (γ, δ, KA, etc.), the stability of the deterministic model is actually described by two dimensionless combinations of these parameters: the ratio of the protein concentration with fully activated promoter (A0 = γ/δ) to the dissociation constant, A0/KA, and the fold-activation, f.

Fig. 2.

Fig. 2.

Stability phase plot for the autoactivator (Fig. 1A), including the effect of intrinsic noise. (A) The black dashed curve is the phase boundary of the deterministic model with transcriptional activation (A0/KA is the fully activated protein concentration scaled by the activator/DNA dissociation constant). Increasing the level of intrinsic noise by increasing the discreteness parameter Δb (i.e., increasing the “burstiness” of translation or decreasing the number of molecules) diminishes the parameter regime of reliable bistability (Re[λ′] < 0). Here, Δb = 0.1 (black solid line), 0.2 (dark gray), and 0.3 (light gray). (B) The average escape time from the stable state is an indicator of the permanence of the bistability. Here, the dark gray curve from A corresponds to an escape time of approximately τ = 6, where time has been scaled relative to the protein lifetime δ−1. (C) As in A, but now with translational activation. The range of bistability is considerably widened as transitions from the LOW to the HIGH state are suppressed. Here, KA·Vcell = 25 molecules, and the fully activated burst size is b = 4 (black), 9 (dark gray), and 14 (light gray).

The ESA we propose is an approximation that allows the average effect of intrinsic noise to be expressed as a positive correction to λ,

graphic file with name zpq01807-6053-m04.jpg

(see Eq. 13 below). The correction reflects an effective “flattening” of the local landscape by stochastic fluctuations, making it easier for the system to escape from the basin of attraction. Adopting this perspective allows the analysis used to study the deterministic model to be extended to the stochastic model with only minor modification. With λ′ corrected to include the effect of the intrinsic noise, the new phase boundaries are drawn to coincide with points in parameter space where λ′ = 0.

A major source of intrinsic noise in gene regulatory networks is so-called translational bursting (7, 26), where each mRNA transcript is translated into several peptides before the message is degraded, leading to a “burst” of protein synthesis. Typical values of the “burst size” b can vary from close to zero for poorly translated genes (27) up to several dozen (28, 29) depending on the rate of translation and the lifetime of the transcript. When intrinsic noise is included in the autoactivator model, and the procedure described in detail in Section IIIA of SI is applied, we find the correction to λ is λcorr ∝ Δb/λ2, where

graphic file with name zpq01807-6053-m05.jpg

is a third dimensionless quantity we call the “discreteness parameter.” This parameter captures the average change in protein number when a synthesis or degradation event occurs, scaled relative to the protein number required to initiate activation NA = KA × Vcell, where KA is the activator dissociation constant and Vcell is the cell volume. Increasing the discreteness parameter Δb increases the magnitude of the discrete change in activator numbers and, therefore, increases the relative magnitude of the perturbation to the system caused by the intrinsic noise. One would expect the circuit to switch more readily from stable state to stable state as the magnitude of the intrinsic noise is increased, thereby reducing the average stability of the circuit. On the other hand, as the number of activator molecules increases (NA → ∞), the discreteness parameter vanishes and the behavior of the system is fully described by the deterministic model. Thus, the discreteness parameter Δb represents a distillation of the complicated effect of intrinsic noise on the model behavior, captured in a compact expression that would be difficult to extract from numerical simulation data.

As shown in Fig. 2A, for the autoactivator the parameter space is divided into regions of bistability (two stable states) and monostability (one stable state). The bistability is most easily lost near the phase boundary separating the bistable and monostable states. The circuit parameters of Isaacs et al. (22) lie close to the left-hand tip of the black triangle in Fig. 2A (f ≈ 10), and as they observed in their experiments, the noise overwhelms bistability in such a system (cf. figure 2A of ref. 22), leading to rapid transitions between the stable states. A much greater fold-activation is required to maintain two distinct stable states (as likewise noted by the authors).

Actually, once noise is allowed in the autoactivator model, one no longer has stability in the strictest sense because there is always a chance that a perturbation will switch the system from one steady state to the other. With noise, it is not a question of stability, but rather the average escape time from the steady state (10, 11). The longer the escape time (compared with other time scales in the problem), the more “stable” the system. To emphasize the effect of the intrinsic noise on the stability phase plot, we consider a system with a small number of activator proteins (KA·Vcell = 25 molecules). Using the parameters γ = 2 protein·min−1, δ−1 = 30 min (a half-life of ≈20 min), KA = 25 nM, and a burst size of b = 10, the discreteness parameter in E. coli (Vcell ≈ 1 μm3) is Δb ≈ 0.2. From Fig. 2A (dark gray curve), a maximum fold-activation of f ≥ 40 is necessary to ensure long-lived bistable states (shown as a cross on the plot). It is possible to explicitly compute the average escape time from the stable states for this simple model (see refs. 30 and 31 and Section II of SI). Fig. 2B compares the average escape time as a function of A0/KA and f for the case above, with Δb = 0.2 (dark gray curve in Fig. 2A). Along the dark gray curve, the escape time is τ = (3 ± 0.5) h, which is approximately six times longer than the protein lifetime (which sets the basic time scale of the system's “memory”).

The escape time is an indirect measure of the system's stability. We have developed a more direct method that measures the effective rate of divergence of an ensemble of stochastic trajectories. This method is of general applicability and allows a direct evaluation of the accuracy of the ESA. The details of that calculation are reserved for SI (see Section III-A.2). Comparing λ′ with the effective rate of divergence in the stochastic simulations of the autoactivator, the ESA is found to be accurate for systems with Δb ≲ 0.25.

The burst size b can be reduced by decreasing the rate of translation, and indeed Ozbudak et al. (27) suggest that many poorly translated genes in E. coli could be the result of evolutionary selection against burst noise. Alternatively, the method of control in the circuit can be shifted from transcriptional to translational activation. Although the simple deterministic model remains unchanged for either choice of transcriptional or translational control, the resulting stochastic model exhibits improved stability for translational activation.

Fig. 2C shows the result of putting the translation rate under control of the activator. Decreasing the translation rate in the LOW state has the effect of shifting the upper branch of the phase boundary, indicating a decrease in transitions from the LOW to the HIGH state. The translational autoactivator can tolerate a larger range of transcription rates (i.e., higher γ) and a lower maximum fold-activation (f ≥ 20), even for large burst size. As above, with γ = 2 protein·min−1, δ−1 = 30 min, KA = 25 nM, Vcell = 1 μm3 and a fully activated burst size b = 10, a fold-activation of f ≥ 25 is required to sustain the bistability (shown as a cross on the plot), almost half that required in the transcriptional autoactivator above.

To generate the phase plot for a given stochastic model requires division of the parameter space of the model (Eq. 2) into a fine grid, with stochastic simulation performed at each point. Even after several such simulations are generated, it is unlikely that the discreteness parameter Δb will suggest itself as a key measure of the magnitude of intrinsic noise. The ESA method provides not only a rapid overview of the parameter space, but also compact expressions characterizing the effect of intrinsic noise on the observable dynamics. In the next section, we apply ESA to the analysis of a more elaborate circuit model.

Genetic Oscillator.

Oscillating systems underlie many physiological processes in the cell, from circadian rhythms (32) to the cell cycle itself (33). In addition to the natural systems, several synthetic genetic oscillator designs have been studied, including the mutually repressing ring-oscillator (Repressilator) of Elowitz and Leibler (34) and the activator–repressor design of Atkinson and coworkers (23) (which has a great deal in common with the model discussed below). A recurring motif in experimentally characterized networks is a negative-feedback loop serving as a system reset (32, 35). Without some time delay or intervening mechanism to prevent reversibility, the system will rapidly approach an intermediate equilibrium, and it is found both theoretically (36) and experimentally (33) that a negative feedback loop alone is not sufficient to maintain reliable oscillations. If, however, the feedback repressor is controlled by a bistable autoactivator, the oscillations become more robust and coherent because the bistable switch acts as a ratchet that “locks” into the HIGH state generating a large amount of repressor to feedback and reset the system to the LOW state where the system remains until the activator accumulates over a critical threshold to initiate another cycle (37). This motif is highly represented in natural gene networks (35), and we will use the ESA to ascertain the contribution of intrinsic noise to the performance of such an oscillator.

We consider the generic model proposed by Vilar et al. (12) to describe circadian rhythms in eukaryotes, with a transcriptional autoactivator driving expression of a repressor that provides negative control by sequestering activator proteins through dimerization (13, 38). The repressor and activator form an inert complex until the activator degrades, recycling repressor back into the system. In their model, the degradation rate of the activator, δA, is the same irrespective of whether it is bound in the inert complex or free in solution. We simplify their original model somewhat, and as in the previous section, we assume fast activator/DNA binding and rapid mRNA turnover, leading to a reduced set of rate equations governing the concentration of activator A, repressor R, and the inert dimer C,

graphic file with name zpq01807-6053-m06.jpg

We further assume no cooperativity in activator binding (n = 1 in the activation function g) and the nominal parameter set used in ref. 12. For this more complicated system, there is a larger number of dimensionless combinations of parameters that characterize the system dynamics. The scaled repressor degradation rate ϵ = δR/δA is a key control parameter in the model because oscillations occur in the deterministic system only for an intermediate range of this parameter. For the nominal parameter set used in ref. 12, the deterministic model exhibits oscillations over the range 0.12 < ϵ < 40 (Fig. 3A, black region). We will focus on the parameter regime near to the phase boundary at ϵ ≈ 0.12 and examine the role intrinsic noise plays in generating regular oscillations from a deterministically stable system.

Fig. 3.

Fig. 3.

Noise induced oscillations in the excitable oscillator. (A) Stability phase plot as a function of the scaled repressor degradation rate ε = δR/δA for the circuit shown in Fig. 1B. The discreteness in the activator synthesis, ΔbA, characterizes the average discrete change in activator concentration during each reaction, and consequently the magnitude of the intrinsic noise. The intrinsic noise expands the region of instability (gray), extending the parameter range over which oscillations are expected to occur. The deterministic phase boundary is located at ε ≈ 0.12 (dashed line separating the black and gray regions). The solid line is the phase boundary predicted from the roots of Eq. 12, and filled circles denote the phase boundary found by stochastic simulation (see text). The model and parameters are as in Vilar et al. (12). (B) The circuit exhibits noise-induced oscillations (dotted line) with interspike time T. The parameters used in the simulation correspond to a deterministically stable system (solid line). Numerical simulation data were generated by using Gillespie's direct method (5), with parameters as used in ref. 12 and ε = 0.1, ΔbA = 6 × 10−2 (cross in A). (See Section III-C of SI.) (C) A plot of the noise-to-signal ratio ηT = 〈(〈T〉 − T)2〉1/2/〈T〉 as a function of ε. The oscillations are regular when ηT is small (the region of noise-induced oscillations predicted by the ESA is gray), and ηT was calculated by using at least 200 spikes for each point.

Applying the ESA to the oscillator model, the parameter ΔbA = (bA + 1)/(2·KA·Vcell) emerges as an important measure quantifying the discreteness in activator synthesis (see Eq. 36 in SI). Here again, bA is the burst size in the activator synthesis, KA is the activator/DNA dissociation constant, and Vcell is the cell volume. (Here, Vcell = 100 μm3 as is appropriate for eukaryotic cells.)

Using the nominal parameter set of Vilar et al. (12) in our reduced model leads to a burstiness in activator synthesis of bA = 5 (giving ΔbA = 6 × 10−2) and a burstiness in repressor synthesis of bR = 10. The phase boundary predicted by the ESA is shown as a solid line in Fig. 3A, bounding a region of parameter space between the deterministic phase boundary where qualitatively different behavior is expected from the stochastic model. We examine the system behavior in this region by running a stochastic simulation using the parameter choice ε = 0.1 and ΔbA = 6 × 10−2 (denoted by a cross in Fig. 3A). With this choice, the deterministic model is stable (Fig. 3B, black line). Nevertheless, a stochastic simulation of the same model, including protein bursting and stochastic dimerization, clearly shows oscillations (Fig. 3B, dotted line).

The time between successive peaks in the stochastic simulation of Fig. 3B is denoted by T. As is clear from Fig. 3B, T is itself a random variable. Each simulation run generates a collection of interspike times from which the mean 〈T〉 and the variance 〈(〈T〉 − T)2〉1/2 can be calculated. Following Steuer et al. (13), the quality of the noise-induced oscillations is measured by using the noise-to-signal ratio ηT = 〈(〈T〉 − T)2〉1/2/〈T〉, and the system is said to exhibit regular oscillations where ηT is small (13, 38). The dependence of ηT on the repressor degradation rate ε is shown in Fig. 3C, with the discreteness parameter ΔbA = 6 × 10−2 (as in Fig. 3B), using at least 200 spikes to calculate ηT. At low repressor degradation rate, the noise-to-signal ratio is high, indicating large variance in the interspike time T and corresponding to a stable (i.e., nonoscillatory) system. As the repressor degradation rate is increased, the variance in the interspike time T decreases with a consequent decrease in the noise-to-signal ratio ηT, indicative of a more regularly oscillating system. Physically, the intrinsic noise in this parameter range is sufficient to drive the system away from the deterministically stable steady state, yet the noise is not so strong that the return trajectory through phase space is much affected.

As in the autoactivator model, it is useful to compare the phase boundary predicted by the ESA to some independent measure of stability, in this case ηT. In Fig. 3C, the ESA phase boundary (for ΔbA = 6 × 10−2) is denoted by the interface between the white and gray regions, corresponding to a value of ηT ≈ 0.2. Using data such as that shown in Fig. 3C, the points in the phase plot with ηT = 0.2 can be found for a range of discreteness parameter ΔbA (Fig. 3A, filled circles). These points correspond very well to the phase boundary calculated by using the ESA (Fig. 3A, solid line). The results are as one would expect: near the deterministic phase boundary, very little molecular noise is required to sustain oscillations, and reasonable periodicity persists even for small values of the discreteness parameter (ΔbA → 0, bR ≠ 0). As the repressor degradation rate ε is decreased to a region favoring stability, more noise is required to overcome the deterministic stability of the system and initiate the autoactivator trigger. It is illustrative to remark that each data point in Fig. 3A, obtained from stochastic simulation (5), took ≈1 day to generate on a dual-processor desktop computer because at low repressor degradation rate, a large separation of timescales is introduced, necessitating long stochastic simulation runs to capture the slowly varying dynamics of the system. By contrast, the solid line generated from the roots of Eq. 12 took <1 h to produce on the same machine. Thus, even for a two-gene circuit with several degrees of freedom, the ESA affords a compact and convenient means to survey the phase space, drawing attention to those regions of particular interest that may be probed in more detail by more realistic (although also more computationally costly) stochastic simulation methods.

Methods

The ESA can be applied to generic models of genetic circuits in a straightforward way. Here, a brief outline of the method is provided. A self-contained tutorial on stochastic modeling and the ESA is provided in SI.

A useful abstraction of genetic regulatory networks is as a system of ordinary differential equations (39, 40). (Here, and throughout, we shall assume a spatially homogeneous environment.) We denote the concentrations of the reactants of interest by the state vector x, where the xi values correspond to the concentration of mRNA, transcription factors, protein products, etc. The kinetic equation governing the evolution of the system takes the form dx/dt = f(x), where f is a vector of nonlinear functions of the state variables. We can estimate the long-time, or steady-state, behavior of the model by first computing the equilibrium points xs that satisfy the algebraic constraint f(xs) = 0. We then Taylor expand the reaction rate vector about the equilibrium point by making the substitution x = xs + xp (where xp is an infinitesimal perturbation away from xs), and retain only linear terms in xp. The resulting dynamics of xp are given by (d/dt)xp = J·xp, where J is the Jacobian or response matrix: Jij = ∂fi/∂xj. The eigenvalues of J are the matrix analogue of the parameter λ introduced in Eq. 3, and in a similar fashion if the eigenvalues all have negative real-part, then xs is a stable steady state. [There are, of course, limitations to how far one can trust the linearization (6), but for our purposes it is sufficient as a first approximation.]

To include stochastic effects in the mathematical model, chemical reaction rates must be rewritten in terms of the reaction “propensity” and “stoichiometry” (5). For example, in the positive autoactivator example above, with the individual synthesis and degradation stoichiometries written explicitly, the deterministic model equations (Eq. 2) read

graphic file with name zpq01807-6053-m07.jpg

We encode this information concisely as the “propensity vector” ν = [ν1, ν2] = [γ·g(A)/b, δ·A] and the stoichiometry matrix S = [b, −1]. The discrete change in molecule numbers after the completion of a chemical reaction causes a deviation from the deterministic solution because the deterministic model assumes an infinitesimally small and continuous change in the state. (Consequently, the deterministic model only applies to systems with large numbers of molecules.) We denote the deviation of the stochastic model from the deterministic model by the fluctuating quantity ω·α(t), where ω=1/Vcell and α(t) describes the stochastic deviation in each species x. The Vcell scaling arises from the observation that the relative magnitude of the intrinsic noise scales roughly as the inverse square root of the number of molecules (15). Elf and Ehrenberg (21) developed an algorithmic expression for the statistics of α using the linear noise approximation of vanKampen (15). In that formulation, the mean and covariance of the fluctuations about the deterministic state are written compactly in terms of the propensity vector ν and the stoichiometry matrix S; here, we shall apply their method to characterize the fluctuations about the stable state. The first step in the calculation of the moments of the fluctuations ωα(t) is to construct the auxiliary matrices Γ and D, evaluated at the stable state xs,

graphic file with name zpq01807-6053-m08.jpg

The drift matrix Γ = J is the response matrix (or Jacobian) described above and reflects the local stability of the deterministic system to small perturbations (41). The diffusion matrix D captures the strength of the fluctuations and is related to the magnitude of the reaction step size (21, 42). It is straightforward to show that to leading order in ω the mean of the fluctuations is zero (〈α〉 = 0) and the variance, denoted by the symmetric matrix Ξ = 〈α·αT〉, is determined by the solution of the system of algebraic equations (15), Γ· Ξ + Ξ·ΓT + D = 0. Because the fluctuations about the stable state are stationary, the time autocorrelation function depends on the time difference only and is given by the matrix exponential

graphic file with name zpq01807-6053-m09.jpg

The effect of the fluctuations on the deterministic steady-state is calculated by including an additional term in the deterministic linearization above: x = xs + xp + ωα. Linearizing J in ω, we have a stochastic differential equation governing the decay of the perturbation modes xp

graphic file with name zpq01807-6053-m10.jpg

The fluctuations affect the decay of the infinitesimal disturbance xp as well as the dynamics of the average 〈xp〉, which (provided ωJ(1)(t) ≪ J(0)) is approximately governed by the convolution equation (43, 44)

graphic file with name zpq01807-6053-m11.jpg

where Jc(t − τ) = 〈J(1)(t)eJ(0)(t − τ)J(1)〉 is made up of linear combinations of the cross-correlations 〈αi(t)αj(τ)〉 given by the ith row and the jth column of the right-hand side of Eq. 9. In the noiseless case, the stability of the perturbation xp is determined by the eigenvalues of J(0): diag{λi} = P−1·J(0)·P where the matrix P is made of the eigenvectors of J(0). The analogues of the eigenvalues for the convolution equation above are found from the poles of the Laplace transform, denoted λ′, which solve the resolvent equation (45)

graphic file with name zpq01807-6053-m12.jpg

Here ω2 has been replaced by Vcell−1 and Ĵc(s) = ∫0∞Jc(t)e−stdt is the Laplace transform of Jc(t). If the deterministic eigenvalues are distinct, we can further approximate the effective eigenvalue λ′i by

graphic file with name zpq01807-6053-m13.jpg

where [·]ii denotes the ith diagonal entry of the matrix. Physically, we interpret the leading-order noise correction as the power in the fluctuations at eigenfrequency λi projected in the eigendirection of λi. Because the correction term is quadratic, it is always positive and thus destabilizes the eigenmode on which it is projected. (Hence, in Eq. 4 we write λcorr > 0.)

It often happens that out of the term 1/Vcell[P−1·Ĵc(λi)·P]ii there appears a small parameter that quantifies the effect of the intrinsic noise. (For the two examples above, the small parameters are Δb and ΔbA, each characterizing the discreteness of the protein change.) In the limit that this parameter goes to zero, the effect of the intrinsic noise becomes negligible, at least in that particular eigenmode.

Finally, the ESA can be easily implemented in a symbolic computational environment, without attending to the mathematical details (see Section IV of SI). A version of the ESA coded in Mathematica is freely available from the authors upon request.

Supplementary Material

Supporting Information

Acknowledgments

We thank Jian Liu, Francis Poulin, and Stefan Klumpp for critical reading and constructive comments on the manuscript. M.S. was supported by a postdoctoral fellowship from Canada's Natural Sciences and Engineering Research Council (NSERC). B.I. is supported by an NSERC Discovery grant. This work was supported in part by National Science Foundation Grant MCB0417721 (to T.H.) and by the Physics Frontiers Center-sponsored Center for Theoretical Biological Physics Grants PHY-0216576 and PHY-0225630.

Abbreviation

ESA

effective stability approximation.

Footnotes

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

References

  • 1.Kobiler O, Rokney A, Friedman N, Court DL, Stavans J, Oppenheim AB. Proc Natl Acad Sci USA. 2005;102:4470–4475. doi: 10.1073/pnas.0500670102. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 2.Hasty J, McMillen D, Collins JJ. Nature. 2002;420:224–230. doi: 10.1038/nature01257. [DOI] [PubMed] [Google Scholar]
  • 3.Kaern M, Blake WJ, Collins JJ. Ann Rev Biomed Eng. 2003;5:179–206. doi: 10.1146/annurev.bioeng.5.040202.121553. [DOI] [PubMed] [Google Scholar]
  • 4.Bintu L, Buchler NE, Garcia HG, Gerland U, Hwa T, Kondev J, Phillips R. Curr Opin Gen Dev. 2005;15:116–124. doi: 10.1016/j.gde.2005.02.007. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 5.Gillespie DT. J Phys Chem. 1977;81:2340–2361. [Google Scholar]
  • 6.Strogatz SH. Nonlinear Dynamics and Chaos. New York: Westview–Perseus; 1994. [Google Scholar]
  • 7.Kaern M, Elston TC, Blake WJ, Collins JJ. Nat Rev Genet. 2005;6:451–464. doi: 10.1038/nrg1615. [DOI] [PubMed] [Google Scholar]
  • 8.Swain PS, Elowitz MB, Siggia ED. Proc Natl Acad Sci USA. 2002;99:12795–12800. doi: 10.1073/pnas.162041399. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 9.vanKampen NG. Stochastic Process in Physics and Chemistry. Amsterdam: North–Holland–Elsevier; 1992. [Google Scholar]
  • 10.Aurell E, Sneppen K. Phys Rev Lett. 2002;88:048101. doi: 10.1103/PhysRevLett.88.048101. [DOI] [PubMed] [Google Scholar]
  • 11.Walczak AM, Onuchic JN, Wolynes PG. Proc Natl Acad Sci USA. 2005;102:18926–18931. doi: 10.1073/pnas.0509547102. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 12.Vilar JMG, Kueh HY, Barkai N, Leibler S. Proc Natl Acad Sci USA. 2002;99:5988–5992. doi: 10.1073/pnas.092133899. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 13.Steuer R, Zhou C, Kurths J. BioSystems. 2003;72:241–251. doi: 10.1016/j.biosystems.2003.07.001. [DOI] [PubMed] [Google Scholar]
  • 14.Suel GM, Garcia-Ojalvo J, Liberman LM, Elowitz M. Nature. 2006;440:545–550. doi: 10.1038/nature04588. [DOI] [PubMed] [Google Scholar]
  • 15.vanKampen NG. Adv Chem Phys. 1976;34:245–308. [Google Scholar]
  • 16.Tanase-Nicola S, Warren PB, tenWolde PR. Phys Rev Lett. 2006;97:068102. doi: 10.1103/PhysRevLett.97.068102. [DOI] [PubMed] [Google Scholar]
  • 17.Pedraza JM, vanOudenaarden A. Science. 2005;307:1965–1969. doi: 10.1126/science.1109090. [DOI] [PubMed] [Google Scholar]
  • 18.Tomioka R, Kimura H, Kobayashi TJ, Aihara K. J Theor Biol. 2004;229:501–521. doi: 10.1016/j.jtbi.2004.04.034. [DOI] [PubMed] [Google Scholar]
  • 19.Paulsson J, Berg OG, Ehrenberg M. Proc Natl Acad Sci USA. 2000;97:7148–7153. doi: 10.1073/pnas.110057697. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 20.DeVille REL, Muratov CB, Vanden-Eijnden E. J Chem Phys. 2006;124:231102. doi: 10.1063/1.2217013. [DOI] [PubMed] [Google Scholar]
  • 21.Elf J, Ehrenberg M. Genome Res. 2003;13:2475–2484. doi: 10.1101/gr.1196503. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 22.Isaacs FJ, Hasty J, Cantor CR, Collins JJ. Proc Natl Acad Sci USA. 2003;100:7714–7719. doi: 10.1073/pnas.1332628100. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 23.Atkinson MR, Savageau MA, Myers JT, Ninfa AJ. Cell. 2003;113:597–607. doi: 10.1016/s0092-8674(03)00346-5. [DOI] [PubMed] [Google Scholar]
  • 24.Ferrell JE., Jr Curr Opin Cell Biol. 2002;14:140–148. doi: 10.1016/s0955-0674(02)00314-9. [DOI] [PubMed] [Google Scholar]
  • 25.Keller AD. J Theor Biol. 1995;172:169–185. doi: 10.1006/jtbi.1995.0014. [DOI] [PubMed] [Google Scholar]
  • 26.Thattai M, vanOudenaarden A. Proc Natl Acad Sci USA. 2001;98:8614–8619. doi: 10.1073/pnas.151588598. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 27.Ozbudak EM, Thattai M, Kurtser I, Grossman AD, vanOudenaarden A. Nat Genet. 2002;31:69–73. doi: 10.1038/ng869. [DOI] [PubMed] [Google Scholar]
  • 28.Kennell D, Riezman H. J Mol Biol. 1977;114:1–21. doi: 10.1016/0022-2836(77)90279-0. [DOI] [PubMed] [Google Scholar]
  • 29.Cai L, Friedman N, Xie XS. Nature. 2006;440:358–362. doi: 10.1038/nature04599. [DOI] [PubMed] [Google Scholar]
  • 30.Gardiner CW. Handbook of Stochastic Methods. 3rd Ed. New York: Springer; 2004. [Google Scholar]
  • 31.Kepler TB, Elston TC. Biophys J. 2001;81:3116–3136. doi: 10.1016/S0006-3495(01)75949-8. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 32.Goldbeter A. Biochemical Oscillations and Cellular Rhythms. Cambridge, UK: Cambridge Univ Press; 1997. [Google Scholar]
  • 33.Pomerening JR, Kim SY, Ferrell JE., Jr Cell. 2005;122:565–578. doi: 10.1016/j.cell.2005.06.016. [DOI] [PubMed] [Google Scholar]
  • 34.Elowitz MB, Leibler S. Nature. 2000;403(6767):335–338. doi: 10.1038/35002125. [DOI] [PubMed] [Google Scholar]
  • 35.Dunlap JC. Cell. 1999;96:271–290. doi: 10.1016/s0092-8674(00)80566-8. [DOI] [PubMed] [Google Scholar]
  • 36.Pomerening JR, Sontag ED, Ferrell JE., Jr Nat Cell Biol. 2003;5:346–351. doi: 10.1038/ncb954. [DOI] [PubMed] [Google Scholar]
  • 37.Cross FR, Siggia E. Dev Cell. 2005;9:309–310. doi: 10.1016/j.devcel.2005.08.006. [DOI] [PubMed] [Google Scholar]
  • 38.Guantes R, Poyatos JF. PLoS Comp Biol. 2006;2:188. doi: 10.1371/journal.pcbi.0020030. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 39.Conrad ED, Tyson JJ. In: System Modeling in Cellular Biology. Szallasi Z, Stelling J, Periwal V, editors. Cambridge, MA: MIT Press; 2006. pp. 97–124. Chap 6. [Google Scholar]
  • 40.Kaern M, Weiss R. In: System Modeling in Cellular Biology. Szallasi Z, Stelling J, Periwal V, editors. Cambridge, MA: MIT Press; 2006. pp. 269–296. Chap 13. [Google Scholar]
  • 41.Ali F, Menzinger M. Chaos. 1999;9:348–356. doi: 10.1063/1.166412. [DOI] [PubMed] [Google Scholar]
  • 42.Scott M, Ingalls B, Kaern M. Chaos. 2006;16:026107. doi: 10.1063/1.2211787. [DOI] [PubMed] [Google Scholar]
  • 43.Bourret RC. Can J Phys. 1965;43:619–639. [Google Scholar]
  • 44.vanKampen NG. Phys Rep. 1976;24:171–228. [Google Scholar]
  • 45.Grossman SI, Miller RK. J Diff Eq. 1973;13:551–566. [Google Scholar]

Associated Data

This section collects any data citations, data availability statements, or supplementary materials included in this article.

Supplementary Materials

Supporting Information