Analytic theory of stochastic oscillations in single-cell gene expression
AAnalytic Theory of Stochastic Oscillations in Single-Cell Gene Expression
Chen Jia , ∗ Michael Q. Zhang , † and Hong Qian ‡ Department of Mathematics, Wayne State University, Detroit, MI 48202, U.S.A. Department of Biological Sciences, University of Texas at Dallas, Richardson, TX 75080, U.S.A. Department of Applied Mathematics, University of Washington, Seattle, WA 98195, U.S.A.
Single-cell stochastic gene expression, with gene state switching, transcription, translation, andnegative feedback, can exhibit oscillatory kinetics that is statistically characterized in terms ofa non-monotonic power spectrum. Using a solvable model, we illustrate the oscillation as astochastic circulation along a hysteresis loop. A triphasic bifurcation upon the increasing strength ofnegative feedback is observed, which reveals how random bursts evolve into stochastic oscillations.Translational bursting is found to enhance the efficiency and the regime of stochastic oscillations.Time-lapse data of the p53 protein from MCF7 single cells validate our theory; the general conclusionsare further supported by numerical computations for more realistic models. These results provide aresolution to R. Thomas’ two conjectures for the single-cell stochastic gene expression kinetics.
The engineering concept of feedback loops has beenextremely useful in understanding nonlinear biologicaldynamics [1] and cellular regulations [2]. In connection togene regulatory networks, the biologist Ren´e Thomas [3]proposed two conjectures in 1981: (a) The existence of apositive feedback loop is a necessary condition for multiplestable states; (b) The existence of a negative feedbackloop is a necessary condition for sustained oscillations.Many efforts since then have been devoted to theinvestigation and proof of these two conjectures. However,almost all previous studies are based on deterministicrepresentations of continuous or discrete dynamicalsystems such as ordinary differential equations (ODEs) orBoolean networks. In particular, the first conjecture hasbeen studied extensively under both the continuous [4–8]and discrete frameworks [9–13]. Studies of the secondconjecture are relatively limited; nice results have beenobtained in the discrete case [11, 14] and only partialresults are available in the continuous case [5, 6].Over the past two decades, large amounts of single-cellexperiments [15–17], some with single-molecule sensitivity,have shown that gene expression in an individual cell is aninherently stochastic process due to small copy numbersof biochemical molecules and their stochastic kinetics [18].To explain “noisy” massive experimental data, significantprogress has been made in the kinetic theory of single-cellstochastic gene expression [18–24]. Yet, there is still alack of an analytic theory of stochastic oscillations at thesingle-cell level; neither is there a resolution of Thomas’ssecond conjecture for stochastic gene expression kineticswith feedback controls.Gene regulatory networks can be tremendously complex,involving numerous feedback loops and signaling steps.However, the situation becomes much simpler if wefocus on a particular gene and the feedback loop thatregulates it [25]. In general, there are three types of grossfeedback topologies for the gene of interest: no feedback,positive feedback, and negative feedback (Fig. 1(a)).Based on the central dogma of molecular biology, geneexpression in an individual cell has a canonical three-stage mRNA gene
EAD B ... and so onprotein C none positive negative active d u v s inactive a b a ( x ) b ( x ) mRNA proteingene FIG. 1. Schematic diagrams of stochastic gene expression inliving cells. (a) Three types of gross feedback topologies. Generegulatory networks in a living cell can be extremely complex,involving numerous feedback loops and signaling steps (graybox). If we focus on a particular gene of interest (red), thenthere are three types of fundamental regulatory relations: nofeedback (none), positive feedback, and negative feedback.The dotted line denotes that there is no link between adjacentnodes. (b) Three-stage model of stochastic gene expressioninvolving gene state switching, transcription, and translation.The gene of interest can switch between an active (red) andan inactive (blue) epigenetic states. Both the mRNA (yellow)and protein (green) can be degraded. representation involving gene state switching betweenactivation and inactivation, transcription, and translation[18] (Fig. 1(b)). Here s is the transcription rate, u is thetranslation rate, and v and d are the degradation rates ofthe mRNA and protein, respectively. The evolution of theprotein concentration x is usually modeled as a hybridswitching ODE, a special class of the so-called piecewisedeterministic Markov process [26–31]:˙ x = sh − dx b ( x ) (cid:15) (active state) , ˙ x = − dx a ( x ) (cid:79) (inactive state) . where h = u/v is the mean burst size of the protein, thatis, the average number of protein copies produced by asingle transcript. In the presence of feedback regulation,the protein abundance x will directly or indirectly affectthe switching rates a ( x ) and b ( x ) of the gene betweenthe active and inactive states. In the simplest case, the a r X i v : . [ q - b i o . M N ] S e p product of the gene may regulate its own expression toform an autoregulatory gene network [21, 23]. If the geneis unregulated, a ( x ) and b ( x ) are constants independentof x . In a positive (negative) feedback network, a ( x ) isan increasing (decreasing) function of x due to epigeneticcontrols such as association (dissociation) of activators,while b ( x ) is an decreasing (increasing) function of x dueto epigenetic controls such as dissociation (association) ofrepressors. To establish an analytic theory of stochasticoscillations, we assume that a ( x ) and b ( x ) are linear withrespect to x [21, 23]: a ( x ) = a − ux, b ( x ) = b + ux, (1)where | u | characterizes the strength of feedback regulationwith u = 0 corresponding to unregulated genes, u > u < a ( x ) < x , we shall automatically set a ( x ) = 0. Thefunctional forms of a ( x ) and b ( x ) introduced aboveare essentially consistent with a popular model for thebiological oscillator of NF- κ B nuclear translocation [32].Let p i ( x, t ) denote the probability density of the proteinconcentration at time t when the gene is in state i , where i = 1 and i = 0 correspond to the active and inactivestates of the gene, respectively. Then the evolution of theMarkovian model is governed by the Kolmogorov forwardequation (cid:40) ∂ t p = d∂ x ( xp ) + b ( x ) p − a ( x ) p ,∂ t p = d∂ x ( xp ) − sh∂ x p + a ( x ) p − b ( x ) p . (2)Let p active ( t ) = (cid:82) ∞−∞ p ( x, t ) dx denote the probability ofthe gene being active at time t and let m ( t ) = (cid:104) x ( t ) (cid:105) denote the mean protein concentration at time t . Byusing (2), the active probability p active ( t ) and proteinmean m ( t ) satisfy the following system of ODEs: ddt (cid:18) p active ( t ) m ( t ) (cid:19) = − T (cid:18) p active ( t ) m ( t ) (cid:19) + (cid:18) a (cid:19) , (3)where T is a matrix defined by T = (cid:18) a + b u − sh d (cid:19) . The two eigenvalues λ and λ of the matrix T are thesolutions to the quadratic equation λ − ( a + b + d ) λ + ( a + b ) d + shu = 0 , which has a discriminant ∆ = ( a + b − d ) − shu . Weassume that the gene network can reach a steady state,which implies that both λ and λ have positive realparts. If the gene is unregulated or positively regulated,then u ≤ >
0, which shows that the twoeigenvalues are both positive real numbers. The oscillations in a stochastic system are usuallycharacterized by the autocorrelation function and powerspectrum. The former C ( t ) = Cov ss ( x (0) , x ( t )) is definedas the steady-state covariance of x (0) and x ( t ) and thelatter G ( ω ) = (cid:82) ∞−∞ C ( | t | ) e − iωt dt is defined as the Fouriertransform of C ( | t | ). In general, oscillations cannot beobserved if G ( ω ) is monotonically decreasing over [0 , ∞ ),while a non-monotonic G ( ω ) over [0 , ∞ ) serves as acharacteristic signal of robust stochastic oscillations withthe maximum point being the dominant frequency. Infact, the autocorrelation function can be represented as C ( t ) = (cid:88) i =0 (cid:90) ∞−∞ xp i ( x, ∞ )[ m i ( x, t ) − m i ( x, ∞ )] dx, where m i ( x, t ) is the conditional mean of the proteinconcentration at time t given that the initial gene stateis i and the initial protein concentration is x . Solving (3)with appropriate initial conditions, we can obtain m ( x, t )and m ( x, t ). In this way, the autocorrelation functioncan also be calculated analytically.If the gene is unregulated, the two eigenvalues of thematrix T are given by λ = a + b and λ = d . In thiscase, both C ( t ) and G ( ω ) can be calculated explicitly as C ( t ) = s h abd ( a + b ) ( d + a + b ) (cid:20) de − ( a + b ) t − ( a + b ) e − dt d − a − b (cid:21) ,G ( ω ) = 2 s h ab ( a + b )( ω + d )( ω + ( a + b ) ) . (4)where C ( t ) is a linear combination of two exponentialfunctions. This result is in full agreement with a recentsingle-cell experiment on nuclear localization of Crz1protein in response to extracellular calcium, where theauthors found that at calcium concentrations greaterthan 100 mM, the autocorrelation function of localizationtrajectories is better fit by a sum of two exponentials [33].In addition, our result also shows that the coefficientsbefore the two exponentials have different signs, whichsuggests that the gene network is in a nonequilibriumsteady state. This is because if a Markovian system is inequilibrium, the autocorrelation function must be a linearcombination of exponential functions with nonpositivecoefficients [34]. From (4), it is easy to see that both C ( t )and G ( ω ) are monotonically decreasing. In a positivefeedback network, similar computations show that C ( t ) = K (cid:20) λ e − λ t − λ e − λ t λ − λ (cid:21) ,G ( ω ) = 2 Kλ λ ( λ + λ )( ω + λ )( ω + λ ) , (5)where K > λ and λ are positivereal numbers, both C ( t ) and G ( ω ) are still monotonicallydecreasing. Thus, no robust oscillations could be observedif the gene is unregulated or positively regulated.More interesting is the situation of negative feedbacknetworks. In the presence of a negative feedback loop,the discriminant ∆ may become negative. In particular,the negative feedback strength u has a critical value u s = ( a + b − d ) sh . If u ≤ u s , then ∆ ≥ λ and λ are positivereal numbers. In this case, both C ( t ) and G ( ω ) have theform of (5) and are monotonically decreasing. If u > u s ,however, then ∆ < λ and λ become conjugatedcomplex numbers. For convenience, we represent themas λ = α − βi and λ = α + βi , where α, β >
0. Inthis case, the autocorrelation function exhibits a dampedoscillation and thus must be non-monotonic: C ( t ) = Ke − αt cos( βt − φ ) , where K > φ satisfiestan φ = α/β . Moreover, the power spectrum is given by G ( ω ) = 2 K ( α cos φ + β sin φ )( α + β )[ α + ( ω − β ) ][ α + ( ω + β ) ] . (6)It is easy to show that G ( ω ) is non-monotonic if and onlyif β > α . Straightforward computations show that β − α = 12 [2 shu − ( a + b ) − d ] . Therefore, there is another critical value u c = ( a + b ) + d sh > u s = ( a + b − d ) sh for the feedback strength u . If u s < u ≤ u c , then β ≤ α and G ( ω ) must be monotonically decreasing. In thiscase, although the autocorrelation displays a dampedoscillation, the power spectrum is still monotonic and thusoscillations cannot be observed. If u > u c , then β > α and thus G ( ω ) becomes non-monotonic. In particular, anegative feedback network with inversely correlated geneswitching rates can exhibit robust oscillations when thefeedback strength is sufficiently large.We have now provided a complete characterization ofthe oscillatory behavior of stochastic gene expression inthree types of gene networks. If the gene is unregulated orpositively regulated, both C ( t ) and G ( ω ) are monotonic,and thus no oscillations could be observed. However, astochastic bifurcation will occur if the gene is negativelyregulated. To gain an intuitive picture of this, we simulatethe trajectories of the Markovian model [35] under threesets of biologically relevant parameters and then usethem to estimate C ( t ) and G ( ω ) (Fig. 2(a)-(c)). Thelatter, as an estimate, is computed by performing fastFourier transforms on the stochastic trajectories and thenapplying the Winner-Kinchin theorem. In fact, the twocritical values u s and u c separate the parameter region time t r a j e c t o r y time au t o c o rr e l a t i onpo w e r s pe c t r u m frequency timefrequency abc time u ≤ u s time u s < u ≤ u c u > u c frequencytime FIG. 2. Stochastic trajectories, autocorrelation functions, andpower spectrums for stochastic gene expression in negativefeedback networks. (a) Trajectories in three phases. (b)Autocorrelation functions in three phases. (c) Power spectrumsin three phases. The model parameters in (a)-(c) are chosenas a = 3 , b = 0 , s = 1 , s = 0 , d = 1. into three phases. In the non-oscillatory phase of u ≤ u s ,both C ( t ) and G ( ω ) are monotonic. In the transitionalphase of u s < u ≤ u c , the former becomes non-monotonicwhile the latter is still monotonic. In the oscillatory phaseof u > u c , both quantities become non-monotonic.Let G peak denote the peak value of the power spectrum.To evaluate the performance of stochastic oscillations, wedefine a quantity called oscillation efficiency: η = G peak − G (0) G peak = 1 − G (0) G peak , which is a number between 0 and 1. If G ( ω ) is monotonic,then G peak = G (0) and thus η vanishes. Therefore, theefficiency η serves as an effective indicator that describesthe robustness of oscillations. When u > u c , the efficiencycan be computed explicitly as η = 1 − (cid:20) αβα + β (cid:21) = (cid:20) sh ( u − u c )( d + a + b ) + 2 sh ( u − u c ) (cid:21) . As the feedback strength u increases, the efficiency η becomes larger and thus the oscillation becomes moreapparent (Fig. 2(a)).Experimentally, the single-cell tracks of gene expressionalways fluctuate stochastically around the mean. For suchtime-lapse data, it is sometimes difficult to distinguishwhether the underlying dynamics belongs to stochasticbursts [36] or stochastic oscillations. This issue has beenpointed out in a recent single-cell experiment on nuclearlocalization of Crz1 protein, where the authors foundthat at very high calcium levels, the stochastic burstsof Crz1 nuclear localization trajectories look similar tosustained oscillations [33]. In fact, our theory leads toa clear distinction between stochastic bursts ( u ≤ u c )and stochastic oscillations ( u > u c ). Typically, stochasticbursts do not have an intrinsic frequency and thus giverise to a monotonic power spectrum, whereas stochasticoscillations have an intrinsic period and are featured by anon-monotonic power spectrum. Therefore, for time-lapsegene expression data, a simple way of distinguishing theunderlying dynamics is to estimate the power spectrumand identify its monotonicity.Thus far, our analytic theory is developed under theassumption that the gene switching rates have the linearexpressions of (1). However, our main results are actuallyinsensitive to the specific functional forms of a ( x ) and b ( x )except their monotonicity. More generally, we may assumethat a ( x ) = a − vx and b ( x ) = b + ux , where u and v aretwo constants jointly characterizing the feedback strength.In this case, it is almost impossible to obtain the analyticsolutions of C ( t ) and G ( ω ). However, according to ournumerical simulations, a negative feedback network alsoundergoes a stochastic bifurcation as u and v increaseswhile keeping u/v as a constant. Furthermore, oscillationsbecome even more robust when a ( x ) and b ( x ) are chosento be nonlinear Michaelis-Menten or Hill functions.Our Markovian model can also be applied to gain adeeper understanding on Thomas’s first conjecture. Bysolving (2), we can obtain the steady-state probabilitydensity of the protein concentration, which turns out tobe the beta distribution p ss ( x ) = Γ( β ) w − β Γ( α )Γ( β − α ) x α − ( w − x ) β − α − , x < w, (7)where w = sh/d is the maximum protein concentration,and α and β are two constants given by α = ad , β = a + bd + shud . In fact, the steady-state protein distribution can be eitherunimodal or bimodal. From (7), bistability occurs if andonly if a < d and u < d ( d − b ) sh . (8)In particular, for unregulated genes, we have u = 0 andthus gene expression shows bistability if and only if a < d and b < d . This suggests that Thomas’s first conjecture isnot always true in living organisms: Bistability can occurin unregulated or even negatively regulated networks whenthe gene switches slowly between the active and inactivestates [31]. Moreover, when b ≥ d , the right side of (8)is nonpositive and thus bistability can only take place inpositive feedback networks, which reinforces Thomas’sfirst conjecture.Thomas’s two conjectures can be understood intuitivelyusing the schematic diagrams depicted in Fig. 3(a),(b).Since the gene switching rates a ( x ) and b ( x ) have distinctmonotonicity in positive and negative feedback networks,the stable and unstable regions of the active and inactivestates are also different in the two types of networks, each forming a hysteresis loop. Here stable and unstableregions are defined as regions with slow and fast geneswitching rates, respectively. In the presence of positivefeedback, the stable attractors of the two gene states liein the stable regions, forming bistability. In the presenceof negative feedback, however, the stable attractors liein the unstable regions. Before the protein level couldapproach a stable attractor, it has already entered theunstable region and thus will switch to the other genestate. With time, the protein level will fluctuate aroundthe hysteresis loop, forming sustained oscillations. Thisunderstanding not only provides a theoretical complementto previous computational studies on relaxation oscillators,where similar hysteretic switch was found [37], but also issupported by previous experimental studies, which showedthat the response of many cell cycle regulatory elementsin Xenopus are hysteretic [38]. unstable stable unstablestable a ( x ) b ( x ) activestateinactivestate unstablestableunstable stable b ( x ) a ( x ) hysteresis loopbistability oscillation a b protein concentration protein concentration FIG. 3. Schematic diagrams of the Markovian dynamics in (a)positive and (b) negative feedback networks. The upper andlower layers represent the active and inactive states of the gene.The gray curves illustrate the graphs of the gene switchingrates a ( x ) and b ( x ). The blue balls show the locations ofthe stable attractors of the two gene states. The stable andunstable regions in the two gene states are depicted as greenand red bars, respectively. The yellow arrows indicate thedirection of the Markovian dynamics. Interestingly, we emphasize that our main results arealso robust with respect to the specific gene expressionmodel used. In some previous works [39–42], the evolutionof the protein concentration x is modeled as a hybridswitching stochastic different equation (SDE):˙ x = ˙ ξ t − dx b ( x ) (cid:15) (active state) , ˙ x = − dx a ( x ) (cid:79) (inactive state) . where ξ t is a compound Poisson process capturing randombursts of the protein. The jump rate of ξ t is exactlythe transcription rate s and the jump distribution of ξ t describes the distribution of the burst size and thus hasthe mean h = u/v . Recently, it has been shown [42]that the switching ODE and SDE models can be viewedas different macroscopic limits of the discrete chemicalmaster equation model of stochastic gene expression [18]as the size of the system tends to infinity. The formerperforms better in the regime of large burst frequencies,while the latter performs better in the regime of largeburst sizes.The the evolution of the switching SDE model isgoverned by the Kolmogorov forward equation (cid:40) ∂ t p = d∂ x ( xp ) + b ( x ) p − a ( x ) p ,∂ t p = d∂ x ( xp ) + sµ ∗ p + a ( x ) p − [ b ( x ) + s ] p . where µ ( dx ) is the probability distribution of the burst sizeand µ ∗ p ( x ) = (cid:82) x p ( x − y ) µ ( dy ) is the convolution of µ and p . If the gene is always active and the burst size hasthe exponential distribution µ ( dx ) = e − x/h /h , then theabove equation reduces to the classical Freidman-Cai-Xierandom bursting model [39]. If the gene is unregulated,both C ( t ) and G ( ω ) can be calculated explicitly as C ( t ) = (cid:98) C ( t ) + (cid:101) C ( t ) , G ( ω ) = (cid:98) G ( ω ) + (cid:101) G ( ω ) . where (cid:98) C ( t ) and (cid:98) G ( ω ) are defined as in (4) and (cid:101) C ( t ) = sσp ( ∞ )2 d e − dt , (cid:101) G ( ω ) = sσp ( ∞ ) d + ω . Here σ = (cid:82) ∞−∞ x µ ( dx ) is the second moment of the burstsize distribution. In this case, both C ( t ) and G ( ω ) arethe sums of two monotonically decreasing functions. Thisleads to the same conclusion that oscillations cannot beobserved if the gene is unregulated. u c u c u c u c u c time t r a j e c t o r y po w e r s pe c t r u m frequency negative feedback strength ac time switching SDEswitching ODE switching SDE switching ODE e ff i c i en cy db t r a j e c t o r y FIG. 4. Stochastic oscillations in negative feedback networksunder the switching ODE and SDE models. (a) Trajectory ofthe switching ODE model when u = u c . (b) Trajectory of theswitching SDE model when u = u c . (c) Power spectrums ofthe switching ODE (blue) and SDE (red) models when u = u c .The power spectrums are normalized so that G (0) = 1. (d)Efficiencies of the switching ODE (blue) and SDE (red) modelsversus the negative feedback strength u . The model parametersare chosen as s = d = 1, b = 0, a = 0 . a = 1in (d). In the switching SDE model, the burst size has theexponential distribution µ ( dx ) = e − x/h /h with h = 1. In the presence of feedback loops, it is almost impossibleto find the analytic expressions of C ( t ) and G ( ω ). Tocompare the oscillations in the switching ODE and SDEmodels, we simulate the trajectories of the two modelsso that they share the same mean field dynamics andthen use them to estimate their power spectrums and oscillation efficiencies (Fig. 4). Based on our simulations,stochastic bifurcations take place in both models. To oursurprise, we find that the critical value ˜ u c for oscillationsin the switching SDE model is always smaller than u c ,suggesting that the switching SDE model needs a smallerfeedback strength u to generate robust oscillations. To seethis, we illustrate the trajectories and power spectrums ofthe two models at the critical value u c in Fig. 4(a)-(c). Itis clear that the switching ODE model has a monotonicpower spectrum and oscillations are not observed, whilethe switching SDE model has a non-monotonic powerspectrum and displays apparent oscillations. This revealsthat random translational bursts are expected to boost theoccurrence and border the region of stochastic oscillations,although the bursting kinetics gives rise to larger geneexpression noise [18]. A reasonable explanation for thiscounterintuitive result is that a single burst will drive theprotein abundance to jump from a low to a much highervalue. Once the burst occurs, the gene will switch rapidlyfrom the active to the inactive state and then the proteinabundance will decay to a lower value again to finish acycle. Therefore, random bursts play an important role inprolonging the cycle times and enhancing the robustnessof stochastic oscillations (Fig. 4(b)). To reinforce theabove results, we depict the efficiencies of the two modelsin Fig. 4(d). Under the model parameters chosen, thecritical value for the switching SDE model is only 0 . u c and for a fixed feedback strength u , the switching SDEmodel possesses a much higher efficiency. Stochasticityincreasing the regime of oscillations has been discussedearlier in [43].To validate our analytic theory, we apply it to one of thebest-studied biological oscillators, the p53-Mdm2 feedbackloop in mammalian cells (Fig. 5(a)). Experiments haveshown that transient DNA damage of double strandbreaks, induced by γ -radiation or radiomimetic drugs,could trigger the oscillatory response of the tumorsuppressor p53 and its negative regulator Mdm2 [44].Specifically, the stress signal of cellular DNA damagefacilitates the phosphorylation of p53 and the degradationof Mdm2. The phosphorylated p53 transcriptionallyactivates Mdm2, which in turn targets p53 for degradation,forming a negative feedback loop. The negative feedbackstrength can be tuned by treating cells with varying dosesof the radiomimetic drug Neocarzinostatin (NCS), whichinduces DNA lesion and elicits p53 oscillations [45]. Thehuman breast cancer cell line MCF7, which containsa stably integrated fluorescent reporter Venus fused top53, was used and the single-cell p53 dynamics underfive different concentrations of NCS was captured usingtime-lapse microscopy [45, 46].In the experiment, fluorescence signals from hundredsof cells were captured every 10 minutes over a timeperiod of 20 hours, where robust stochastic oscillationswere observed under strong DNA damage (Fig. 5(b)).For each dose of NCS, we estimated the autocorrelation time c e ll s po w e r s pe c t r u m frequency bd DNA damagep53 Mdm2mdm2promoter a time (hour) au t o c o rr e l a t i on c
200 ng/mL NCS100 ng/mL NCS0 ng/mL NCS25 ng/mL NCS50 ng/mL NCS 200 ng/mL NCS100 ng/mL NCS25 ng/mL NCS50 ng/mL NCS0 ng/mL NCS
FIG. 5. Experimental validation of the analytic theory. (a)The p53-Mdm2 negative feedback loop in single cells. (b) Heatmap of the single-cell trajectories of p53 fluorescence aftertreatment with 200 ng/mL NCS. Each row represents thetime-lapse measurements for a single cell. The colorbar showshow many standard deviations a data value is away from themean. (c) Normalized autocorrelation functions under fivedifferent doses of NCS: 200 ng/mL (black), 100 ng/mL (gray),50 ng/mL (blue), 25 ng/mL (red), and 0 ng/mL (green). (d)Normalized power spectrums under five doses of NCS. function and power spectrum (Fig. 5(c)-(d)). Accordingto our data analysis, the autocorrelation does not oscillatewhen there is no NCS added and displays an apparentdamped oscillation when the NCS dose is above 25 ng/mL,whereas the power spectrum is monotonically decreasingwhen the NCS dose is below 25 ng/mL and yields anonzero peak when the NCS dose is above 50 ng/mL. Withincreased doses of NCS, the system undergoes a stochasticbifurcation from the non-oscillatory to the transitionalphase at a critical value u s between 0 ng/mL and 25ng/mL and undergoes another stochastic bifurcation fromthe transitional to the oscillatory phase at a higher criticalvalue u c between 25 ng/mL and 50 ng/mL. All theseresults are in full agreement with our theory.There are growing observations showing thatgene expression in individual cells can displaystochastic oscillations [47–51], which are distinct fromrandom fluctuations and should be quantified by anon-monotonic power spectrum [52, 53] with an oscillatoryautocorrelation function [54, 55], or be understood as acircular random walk or stochastic flow [43, 56]. Althoughnegative feedback is widely suspected to be necessarybased on deterministic dynamics [3], more mechanisticnature of this phenomenon is still largely unknown. Ouranalytic theory shows that a negative feedback loop withinversely correlated gene switching rates can give rise to astochastic bifurcation evolving from stochastic bursts tostochastic oscillations for single-cell gene expression. Inthe oscillatory regime, the essence of stochastic oscillationsis revealed to be a circular motion along a stochastichysteresis loop. Even though random translational bursts yield large gene expression noise, it can enhancethe efficiency and the regime of stochastic oscillations.Further elucidations of the biochemical steps for negativefeedback loops in specific cells are expected.We are grateful to X.S. Xie and X.-J. Zhang for helpfuldiscussions and sharing unpublished research results. M.Q.Zhang was supported by NIH Grants MH102616 andMH109665 and also by NSFC Grants 31671384 and91329000. ∗ [email protected] † [email protected] ‡ [email protected][1] J. D. Murray, Mathematical Biology: I. An Introduction ,3rd ed. (Springer, New York, 2007).[2] U. Alon,
An Introduction to Systems Biology: DesignPrinciples of Biological Circuits (Chapman and Hall, NewYork, 2006).[3] R. Thomas, in
Numerical Methods in the Study of CriticalPhenomena (Springer, 1981) pp. 180–193.[4] E. Plahte, T. Mestl, and S. W. Omholt, J. Biol. Syst. ,409 (1995).[5] J.-L. Gouz´e, J. Biol. Syst. , 11 (1998).[6] E. H. Snoussi, J. Biol. Syst. , 3 (1998).[7] O. Cinquin and J. Demongeot, J. Theor. Biol. , 229(2002).[8] C. Soul´e, ComPlexUs , 123 (2003).[9] J. Aracena, J. Demongeot, and E. Goles, IEEETransactions on Neural Networks , 77 (2004).[10] A. Richard and J.-P. Comet, Discrete Appl. Math. ,2403 (2007).[11] ´E. Remy, P. Ruet, and D. Thieffry, Adv. Appl. Math. , 335 (2008).[12] J. Aracena, Bull. Math. Biol. , 1398 (2008).[13] A. Richard, Discrete Appl. Math. , 3281 (2009).[14] A. Richard, Adv. Appl. Math. , 378 (2010).[15] M. B. Elowitz, A. J. Levine, E. D. Siggia, and P. S. Swain,Science , 1183 (2002).[16] L. Cai, N. Friedman, and X. S. Xie, Nature , 358(2006).[17] L. Bintu, J. Yong, Y. E. Antebi, K. McCue, Y. Kazuki,N. Uno, M. Oshimura, and M. B. Elowitz, Science ,720 (2016).[18] J. Paulsson, Phys. Life Rev. , 157 (2005).[19] J. Paulsson and M. Ehrenberg, Phys. Rev. Lett. , 5447(2000).[20] E. Sontag, A. Kiyatkin, and B. N. Kholodenko,Bioinformatics , 1877 (2004).[21] J. Hornos, D. Schultz, G. Innocentini, J. Wang,A. Walczak, J. Onuchic, and P. Wolynes, Phys. Rev.E , 051907 (2005).[22] V. Shahrezaei and P. S. Swain, Proc. Natl. Acad. Sci.USA , 17256 (2008).[23] N. Kumar, T. Platini, and R. V. Kulkarni, Phys. Rev.Lett. , 268105 (2014).[24] Z. Cao and R. Grima, Nat. Commun. , 3305 (2018).[25] I. Lestas, G. Vinnicombe, and J. Paulsson, Nature ,174 (2010).[26] H. Ge, H. Qian, and X. S. Xie, Phys. Rev. Lett. , , 185001 (2015).[28] Y. T. Lin and C. R. Doering, Phys. Rev. E , 022409(2016).[29] P. C. Bressloff, J. Phys. A: Math. Theor. , 133001(2017).[30] Y. T. Lin and N. E. Buchler, J. R. Soc. Interface ,20170804 (2018).[31] H. Ge, P. Wu, H. Qian, and X. S. Xie, PLoS Comput.Biol. , e1006051 (2018).[32] L. Ashall, C. A. Horton, D. E. Nelson, P. Paszek, C. V.Harper, K. Sillitoe, S. Ryan, D. G. Spiller, J. F. Unitt,D. S. Broomhead, et al. , Science , 242 (2009).[33] L. Cai, C. K. Dalal, and M. B. Elowitz, Nature , 485(2008).[34] C. Jia and Y. Chen, J. Phys. A: Math. Theor. , 205001(2015).[35] G. Yin and C. Zhu, Hybrid switching diffusions: propertiesand applications , Vol. 63 (Springer, New York, 2010).[36] D. M. Suter, N. Molina, D. Gatfield, K. Schneider,U. Schibler, and F. Naef, Science , 472 (2011).[37] T. Y.-C. Tsai, Y. S. Choi, W. Ma, J. R. Pomerening,C. Tang, and J. E. Ferrell, Science , 126 (2008).[38] W. Sha, J. Moore, K. Chen, A. D. Lassaletta, C.-S. Yi,J. J. Tyson, and J. C. Sible, Proc. Natl. Acad. Sci. USA , 975 (2003).[39] N. Friedman, L. Cai, and X. S. Xie, Phys. Rev. Lett. ,168302 (2006).[40] M. C. Mackey, M. Tyran-Kaminska, and R. Yvinec,SIAM J. Appl. Math. , 1830 (2013).[41] J. Jedrak and A. Ochab-Marcinek, Phys. Rev. E ,032401 (2016). [42] C. Jia, M. Q. Zhang, and H. Qian, Phys. Rev. E ,040402(R) (2017).[43] M. Vellela and H. Qian, Proceedings of the Royal SocietyA: Mathematical, Physical and Engineering Sciences ,771 (2010).[44] G. Lahav, N. Rosenfeld, A. Sigal, N. Geva-Zatorsky, A. J.Levine, M. B. Elowitz, and U. Alon, Nat. Genet. , 147(2004).[45] E. Batchelor, A. Loewer, C. Mock, and G. Lahav, Mol.Syst. Biol. , 488 (2011).[46] R. Moore, H. K. Ooi, T. Kang, L. Bleris, and L. Ma,PLoS Comput. Biol. , e1004653 (2015).[47] M. B. Elowitz and S. Leibler, Nature , 335 (2000).[48] A. Hoffmann, A. Levchenko, M. L. Scott, andD. Baltimore, Science , 1241 (2002).[49] D. Nelson, A. Ihekwaba, M. Elliott, J. Johnson, C. Gibney,B. Foreman, G. Nelson, V. See, C. Horton, D. Spiller, et al. , Science , 704 (2004).[50] N. Geva-Zatorsky, N. Rosenfeld, S. Itzkovitz, R. Milo,A. Sigal, E. Dekel, T. Yarnitzky, Y. Liron, P. Polak,G. Lahav, et al. , Mol. Syst. Biol. (2006).[51] E. Batchelor, C. S. Mock, I. Bhan, A. Loewer, andG. Lahav, Mol. Cell , 277 (2008).[52] H. Qian and M. Qian, Phys. Rev. Lett. , 2271 (2000).[53] D. A. Potoyan and P. G. Wolynes, Proc. Natl. Acad. Sci.USA , 2391 (2014).[54] H. Qian, S. Saffarian, and E. L. Elson, Proc. Natl. Acad.Sci. USA , 10376 (2002).[55] Y. Cao, H. Wang, Q. Ouyang, and Y. Tu, Nat. Phys. ,772 (2015).[56] J. Wang, L. Xu, and E. Wang, Proc. Natl. Acad. Sci.USA105