Bimodality in gene expression without feedback: From Gaussian white noise to log-normal coloured noise
BBimodality without feedback:Noise-induced transitions in the repressed gene with extrinsic log-normal noise
Gerardo Aquino and Andrea Rocco ∗ Department of Microbial and Cellular Sciences,Faculty of Health and Medical Sciences, University of Surrey,GU2 7XH, Guildford, UK
Extrinsic noise-induced transitions have been largely investigated in a variety of chemical andphysical systems. The key properties that make these systems mathematically tractable is that thenoise appears linearly in the dynamical equations and it is assumed Gaussian and white. Here weconsider the simplest unit of gene regulation, the repressed gene, which is instead characterized bynonlinear dependencies on external parameters, leading to corresponding nonlinear extrinsic noise.We adopt a general methodology to address this nonlinearity based on assuming that the fluctuationsexhibit a finite and large correlation time. We also relax the Gaussian assumption, and consider themore biologically relevant log-normal noise. We find that in contrast to Gaussian noise, log-normalnoise leads to noise-induced transitions, with the system becoming bimodal. We emphasize thatno feedback loops are present in the system, and therefore our findings identify a novel class ofnoise-induced transitions, based on both features of noise nonlinearity and boundedness.
PACS numbers:
Bistability is a feature dramatically relevant for thefunctioning of living systems [1]. Deterministic bistabil-ity is realized in gene regulatory networks by the inclusionof specific topological features, such as feedback loops,which provide for the simultaneous existence of two, ormore, stable steady states. In fact, feedback loops areknown to underlie a number of processes involved in cel-lular decision making, ranging from the LAC operon [2],to quorum sensing [3], to cell differentiation [4].However, in physical and chemical systems topologicalfeatures are not the only known determinants of bistablebehaviour. Extrinsic noise has been investigated now fora while, and it has emerged as an active dynamical player,which can produce so-called noise-induced transitions tobimodal dynamics in systems otherwise deterministicallymonostable [5]. In contrast, no noise-induced transitionshave been identified in biological systems not includingfeedack loops, despite the effect of extrinsic noise havingbeen investigated extensively [6–8].In most sytems, the search and the possible emergenceof such transitions relies on the multiplicative nature ofthe corresponding stochastic dynamics, which is mathe-matically well characterized when the noise appears lin-early in the system’s dynamical equations, and it is gaus-sian and white.These assumptions are not free of criticism. First,regulatory processes of gene expression are usually de-scribed in terms of strongly nonlinear dynamics, such asHill functions. If the parameter affected by noise appearsnonlinearly in the dynamical equations, the mathemat-ical treatement of the corresponding stochastic dynam-ics for white noise becomes very hard, if not impossi-ble [5]. Second, the gaussian approximation is of limitedapplicability for fluctuations affecting parameters whichare strictly positive. This problem has been recognized in some recent literature dealing with so-called boundednoise [9–11], and in the case of linear noise, noise-inducedtransitions have been identified in different systems [12–14]. And third, the white character of extrinsic noise islargely questionable for gene regulatory networks, wherethe typical correlation times of extrinsic noise can extendup to and well above typical cell cycle times [15].In this Letter we take up these issues, and show thatnonlinear log-normal noise provokes noise-induced tran-sitions in a gene regulatory system in absence of feed-back loops. We do so by extending to nonlinear noise awell established formalism, first developed in [16] and ap-plied later in a variety of chemical and physical systems[17]. The formalism is based on relaxing the assumptionof white noise, and on focusing instead on fluctuationscharacterized by a correlation time much longer than allother relaxational timescales in the system. This is in-deed the typical situation when modelling extrinsic noisein gene regulation [15].We start by considering the repressed gene, a genewhich is constitutively expressed at its maximal expres-sion rate, but it is regulated by a transcription factor,the repressor, that represses protein synthesis by actingas an external control parameter (Fig. 1).
FIG. 1: The repressed gene. The expression level of gene G x is downregulated by the binding of the repressor R to the theDNA binding site BS. We make the assumption that whenactive, gene G x synthesizes the protein x in one single step ofcombined transcription and translation. a r X i v : . [ phy s i c s . b i o - ph ] S e p We describe this system by decomposing the rate equa-tion for the protein x into a production term, describedphenomenologically in terms of a Hill function, and adegradation term, straightforwardly written according tothe Law of Mass Action: dxdt = f ( x, R ) = g ρR − kx. (1)Here ρ is the association constant of R to DNA, which de-scribes the strength of the repressor binding to its bindingsite, g is the maximal expression rate for gene G x , and k is the degradation rate of the protein x . These dynamicsdescribe the down regulation of the gene G x by the re-pressor R , which acts as an external control parameter.Due to the much shorter half-life of mRNA with respectto that of proteins [18], we here assume that transcrip-tion and translation are lumped together in one singlestep, so that when the gene is in the active state it syn-thesises directly the protein x . The dynamics specifiedby the function f ( x, R ) is characerized by a timescale τ x ,which in the specific case considered here is simply givenby the protein inverse degradation rate, τ x = 1 /k .Let us now consider fluctuations acting on R . The na-ture of these fluctuations is multiple, and their modellingis entirely phenomenological, as we are not assumingor describing any specific molecular mechanism underly-ing their occurrence. While the Central Limit Theoremwould suggest to model them as a Gaussian process, thischoice is questionable for fluctuations affecting a strictlypositive parameter. In fact, extrinsic fluctuations are welldescribed by log-normal distributions [15, 19], which arealso regarded as phenomenological possible candidates tointerpret universal features in bacteria and yeast [20, 21].Furthermore, irrespective of the chosen distribution, thenonlinear dependency of Eq. (1) on R makes the cor-responding white noise limit mathematically ill-defined[5].In order to tackle these issues, first we properly definethe nonlinear noise terms in (1) by relaxing the whiteproperties of the noise, and considering coloured noisesinstead, characterized by a finite correlation time τ (cid:29) τ x .Secondly, we describe log-normal fluctuations, and ana-lyze their effect on the system, compared to the Gaussiannoise case.Thus we supplement Eq. (1) with a complementaryequation that describes the fluctuations affecting the pa-rameter R : dRdt = 1 τ µ ( R ) + (cid:114) Dτ ν ( R ) ξ. (2)Here the functions µ ( R ) and ν ( R ) are kept generic forthe moment, but can be chosen so as to reproduce dif-ferent types of fluctuations, including Gaussian and nonGaussian noises. The rescaling by the factors 1 /τ is intro-duced so as to slow-down the fluctuations and allows usto carry out a systematic expansion in powers of 1 / √ τ , for τ (cid:29) τ x . Finally, the variable ξ = ξ ( t ) is gaussian,zero average, white noise, with correlator given by (cid:104) ξ ( t ) ξ ( t (cid:48) ) (cid:105) = 2 δ ( t − t (cid:48) ) . (3)The Master Equation of the system specified by Eqs.(1) and (2) results in ∂w t ( x, R ) ∂t = − ∂∂x [ f ( x, R ) w t ( x, R )] − τ ∂∂R [( µ ( R ) + Dν ( R ) ν (cid:48) ( R )) w t ( x, R )]+ Dτ ∂ ∂R (cid:2) ν ( R ) w t ( x, R ) (cid:3) , (4)where w t ( x, R ) is the joint probability density for the x and R variables. The term proportional to ν ( R ) ν (cid:48) ( R ) isthe so-called Stratonovich drift, corresponding to havingassumed the Stratonovich interpretation for the multi-plicative noise equation (2).As it stands, Eq. (4) is difficult to solve, but it canbe solved in stationary conditions when the stochasticfluctuations are slow. By following [16, 17], we apply theso-called “switching-curve approximation”, and expandthe stationary solution w S ( x, R ) of (4) in inverse powersof the correlation time τ of the R process: w S ( x, R ) = w ( x, R ) + 1 τ / w ( x, R )+ 1 τ w ( x, R ) + O (cid:18) τ / (cid:19) . (5)We then obtain to the zeroth order in 1 /τw ( x, R ) = w ( R ) δ ( x − u ( R )) , (6)where u ( R ) is defined so that f ( u ( R ) , R ) = 0. By assum-ing that w ( x, R ) is normalized, w ( R ) in (6) can be identi-fied by solving to order τ − the stationary Fokker-PlanckEquation associated to (2) and supplemented with theStratonovich interpretation:0 = D ∂ ∂R (cid:2) ν ( R ) w ( R ) (cid:3) − ∂∂R [( µ ( R ) + Dν ( R ) ν (cid:48) ( R )) w ( R )] . (7)The marginalized probability density p ( x ) for the x pro-cess can then be obtained as p ( x ) = (cid:90) w ( x, R ) dR = (cid:90) w ( R ) δ ( x − u ( R )) dR. (8)By using δ ( ϕ ( x )) = δ ( x − x ) / | ϕ (cid:48) ( x ) | , where x is theroot of ϕ ( x ), ϕ ( x ) = 0, we readily obtain: p ( x ) = w ( u − ( x )) (cid:12)(cid:12)(cid:12)(cid:12) du − ( x ) dx (cid:12)(cid:12)(cid:12)(cid:12) . (9)The probability density w ( R ) can be computed from (7),and we obtain the central methodological result of thisLetter, namely that p ( x ) = N ν ( u − ( x )) (cid:12)(cid:12)(cid:12)(cid:12) du − ( x ) dx (cid:12)(cid:12)(cid:12)(cid:12) exp (cid:40) D (cid:90) u − ( x ) µ ( y ) dyν ( y ) (cid:41) , (10)where N is a normalization constant.The mode(s) of the probability density p ( x ) representthe natural extension of the stable steady state of the cor-responding deterministic system, and can be computedexplicitly by differentiating p ( x ). This leads readily to µ − Dνν (cid:48) − Dν (cid:18) du − ( x ) dx (cid:19) u (cid:48)(cid:48) = 0 , (11)where µ = µ ( u − ( x )), ν = ν ( u − ( x )), and u (cid:48)(cid:48) = u (cid:48)(cid:48) ( u − ( x )).Eq. (11) is an extension of the standard equationidentifying the modes of a multiplicative stochastic pro-cess supplemented with the Stratonovich prescription.The extra term − Dν (cid:0) du − ( x ) /dx (cid:1) u (cid:48)(cid:48) in (11) accountsfor dynamical modifications of the deterministic systemwhich are induced by the slow fluctuations. This term isidentically zero only in the case when the the relation-ship between the variable R and the variable x is linear(so that u (cid:48)(cid:48) = 0), while in all other cases nontrivial mod-ifications of the dynamics can be expected. However,whether these modifications will correspond to a mereshift of parameter values, a change of stability proper-ties of the deterministic attractor states, or a change intheir number (namely a pure noise-induced transition)will need to be assessed case by case, once the dynamicsof the x and R processes are specified. We are going toshow here two cases in which a noise-induced transitionoccurs, and does not, respectively.Let us consider first the repressed gene model given byEq. (1), with log-normal noise on R . In this case thefunction u ( R ) is defined by g ρR − ku ( R ) = 0 ⇒ u − ( x ) = gkρx − ρ , (12)and we define R → R ( t ) = ¯ Re η ( t ) e − D/ , (13)where η ( t ) is the standard Ornstein-Uhlenbeck noise,with the Langevin representation dηdt = − ητ + (cid:114) Dτ ξ ( t ) . (14)Here ξ ( t ) is zero average Gaussian white noise. The defi-nition (13) guarantees that (cid:104) R (cid:105) = ¯ R , since (cid:104) e η ( t ) (cid:105) = e D/ .The fluctuations on R can then be described by dRdt = − Rτ (cid:20) ln R − ln ¯ R + D (cid:21) + (cid:114) Dτ Rξ ( t ) (15) FIG. 2: Log-normal distributions for R and correspondingprobability distributions for x as from Eqs. (16) and (17)respectively for different values of noise intensity D . Fromtop to bottom, D = 1, D = 3, and D = 5. Other parameters: g = 0 .
01 s − , ρ = 10 nM − , k = 10 − s − , ¯ R = 0 . τ = 10 s. For D = 5, modes are located at x = 1 .
14 and x = 99 . and are characterized by the stationary log-normal dis-tribution: w ( R ) = 1 √ πD R exp (cid:34) − D (cid:18) ln R − ln ¯ R + D (cid:19) (cid:35) . (16)By direct integration of Eq. (10), we readily obtain p ( x ) = g √ πD x ( g − kx )exp (cid:40) − D (cid:20) ln (cid:18) gkρx − ρ (cid:19) − ln ¯ R + D (cid:21) (cid:41) (17)and, for x (cid:54) = 0 and x (cid:54) = g/k , from Eq. (11) g ln (cid:18) gkρx − ρ (cid:19) − g ln ¯ R − gD Dkx = 0 . (18)As an illustration of these dynamics, we consider a pro-totypical bacterial gene with the following biologicallygrounded parameter values. We estimate the effectivemaximal gene expression parameter as g ≈ g P g M /k M =10 − s − , where g P and g M are transcriptional and trans-lational rates respectively, and k M is the degradation FIG. 3: Stationary-state response curves for different valuesof the noise intensity, as determined by Eq. (18). The case D = 0 corresponds to the deterministic system. Parametervalues are as in Fig. 2. rate of mRNA. This choice is compatible with a moder-ate transcriptional and translational activity and typicalmRNA degradation rates in bacteria [22]. Also, we set k = 10 − s − , which corresponds to a protein half-life inthe range of 2 hours, and a typical dissociation constantof trancription factors to DNA ρ = 10 nM − [23]. Wealso fix ¯ R = 0 . ρ gives on average a moderate repression activity of R ongene G x . The timescale of the fluctuations τ is crucialfor the matching of our theoretical predictions with thenumerical results. Given k = 10 − s − , we set τ = 10 s,to capture the dynamics of slow fluctuations of R . Noiseintensities D in the range 1 −
10 are chosen to reproducefairly large fluctuations, as observed in single cells [24].Simulations show an excellent agreement with the the-oretical predictions, as shown in Fig. 2. Further to theexcellent matching of the full distribution for R , the pre-dictions for the modes as given by Eq. (18) is also ver-ified. In Fig. 3 we show the extrema of the stationaryprobability (17) as function of ¯ R for different noise in-tensities, obtained by the numerical solution of Eq. (18)with parameter values as in Fig. 2. The curve associ-ated to D = 0 (deterministic case) is single valued forall values of ¯ R , with further increase of D provoking theappearance of a second stable mode for intermediate ¯ R values.For comparison, we also consider Gaussian fluctuationson the parameter R , by assuming R → ¯ R + η ( t ), where η ( t ) is the standard Ornstein-Uhlenbeck process, whosedynamics is given by Eq. (14). Again by direct integra-tion of Eq. (10), we readily obtain p ( x ) = 1 √ πD (cid:18) gkρx (cid:19) exp (cid:40) − D (cid:18) gkρx − ρ − ¯ R (cid:19) (cid:41) . (19)The resulting p ( x ) is non trivial in this case as well, asthe original Gaussian noise is transformed nonlinearlyinto the non-Gaussian probability density (19). Simula- FIG. 4: Left panel: Ornstein-Uhlenbeck Gaussian distribu-tion for R and corresponding probability distribution for x asfrom Eqs. (19) for D = 2. Right panel: Corresponding bifur-cation diagram. The deterministic branch emanates from thedeterministic solution at D = 0. The stochastic branch onlyexists for D >
0. Parameter values as in Fig. 2. tions agree also in this case very well with the analyticalpredictions. However, in this case only one positive modeis present, which extends the deterministic solution. Infact, for x (cid:54) = 0, Eq. (11) becomes2 Dk ρ x + g (1 + ρ ¯ R ) kx − g = 0 (20)which for D (cid:54) = 0 accounts for the two real solutions: x , = − g Dkρ (1 + ¯ Rρ ) ± (cid:115)(cid:20) g Dkρ (1 + ¯ Rρ ) (cid:21) + g Dk ρ (21)The solution of (20) for D = 0 coincides trivially withthe solution of the deterministic case, with x being thestochastic continuation of it. The x branch emerges in-stead for D >
0, signalling that no transition is takingplace in the system at finite D (the only critical pointbeing the trivial one, D = D c = 0). Despite the appear-ance of a second mode, this is to be biologically discardedbecause at negative x values.In conclusion, we have introduced an efficient method-ology to deal with nonlinear extrinsic noise when the cor-relation time of the fluctuations is large. Our resultsgeneralize the concept of noise-induced transitions ex-clusively due to multiplicative noise terms, and relatedto the so-called Stratonovich drift in the Gaussian andwhite limit. In fact, in the two examples analysed here,the Stratonovich drift does not provoke any qualitativechange in the dynamics of the system, as the drivingnoise is unimodal in all cases. The extra correction termsidentified are instead responsible for the possible onsetof noise-induced transitions that may or may not occurdepending on the general features of the driving noisedistribution. While nonlinear Gaussian noise does notprovoke any biologically relevant transition to bimodal-ity, nonlinear log-normal noise appears to produce moreinteresting effects.Our results show that particular care should be takenwhen reconstructing gene regulatory networks underthe evidence of bimodal distributions of gene expressionlevels. The usual implication that these require a wiringdiagram including feedback loops may be false, asthe same bimodality may be produced by nonlinearextrinsic noise. The novel mechanism for noise-inducedtransitions here described is promising and worthy offurther analysis and experimental validation. It willcontribute to our fundamental understanding of therelevance of noise in biological systems.This work was supported by the UK Biotechnologyand Biological Sciences Research Council (BBSRC), un-der grant BB/L007789/1. ∗ Electronic address: [email protected][1] J.R. Pomerening, Curr. Opin. Biotechnol. , 381 (2008).[2] M. Santillan, M.C. Mackey, E.S. Zeron, Biophysical J. , 3830 (2007).[3] J.W. Williams, X. Cui, A. Levchenko, A.M. Stevens, Mol.Sys. Biol. , 234 (2008).[4] J.E. Ferrell, Curr. Biol. , R458 (2012).[5] W. Horsthemke and R. Lefever, Noise-Induced Transi-tions , Springer (2006).[6] M. Samoilov, S. Plyasunov, A.P. Arkin, PNAS , 2310(2005). [7] E. Pujadas and A.P. Feinberg, Cell , 1123 (2012).[8] M. Weber and J. Buceta, PLOS One , e73487 (2013).[9] S. de Franciscis, G. Caravagna, A. d’Onofrio, Nat. Com-put. , 297 (2014).[10] L. Borland, Phys. Lett. A , 67 (1998).[11] R.V. Bobryk, A. Chrzeszczyk, Phys A , 263 (2005).[12] H.S. Wio, R. Toral, Physica D , 161 (2004).[13] A. d’Onofrio, Phys. Rev. E , 021923 (2010).[14] A. d’Onofrio, A. Gandolfi, Phys. Rev. E , 061901(2010).[15] N. Rosenfeld, J.W. Young, U. Alon, P.S. Swain, M.B.Elowitz, Science , 1962 (2005).[16] L. Arnold, W. Horsthemke, R. Lefever, Z. Physics B ,367 (1978).[17] L.E. Reichl and W.C. Schieve, Instabilities, bifurca-tion, and Fluctuations in Chenical Systems , Universityof Texas Press, Austin (1982).[18] N. Rosenfeld, M.B. Elowitz, U. Alon, J. Mol. Biol. ,785 (2002).[19] V. Shahrezaei, J.F. Ollivier, P.S. Swain, Mol. Sys. Bio. , 196 (2008).[20] H. Salman et al., Phys. Rev. Lett. , 238105 (2012)[21] N. Brenner et al., Phys. Rev. E , 042713 (2015).[22] M Thattai and A. van Oudenaarden, PNAS , 8614(2001).[23] A. Loinger, A. Lipshtat, N.Q. Balaban, O. Biham, PhysRev. , 021904 (2007).[24] L. Cai, N. Friedman, X.S. Xie, Nature440