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.

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.

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 (A₀/K_A 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 δ⁻¹. (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, K_A·V_cell = 25 molecules, and the fully activated burst size is b = 4 (black), 9 (dark gray), and 14 (light gray).

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. and ε = 0.1, Δ_bA = 6 × 10⁻² (cross in A). (See Section III-C of SI.) (C) A plot of the noise-to-signal ratio η_T = 〈(〈T〉 − T)²〉^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.