Role of DNA binding sites and slow unbinding kinetics in titration-based oscillators
RRole of DNA binding sites and slow unbinding kineticsin titration-based oscillators
Sargis Karapetyan
1, 2 and Nicolas E. Buchler
1, 2, 3 Department of Physics, Duke University, Durham, NC 27708 Center for Genomic & Computational Biology, Durham, NC 27710 Department of Biology, Duke University, Durham, NC 27708 (Dated: October 15, 2018)Genetic oscillators, such as circadian clocks, are constantly perturbed by molecular noise arisingfrom the small number of molecules involved in gene regulation. One of the strongest sources ofstochasticity is the binary noise that arises from the binding of a regulatory protein to a promoter inthe chromosomal DNA. In this study, we focus on two minimal oscillators based on activator titrationand repressor titration to understand the key parameters that are important for oscillations and forovercoming binary noise. We show that the rate of unbinding from the DNA, despite traditionallybeing considered a fast parameter, needs to be slow to broaden the space of oscillatory solutions.The addition of multiple, independent DNA binding sites further expands the oscillatory parameterspace for the repressor-titration oscillator and lengthens the period of both oscillators. This effect isa combination of increased effective delay of the unbinding kinetics due to multiple binding sites andincreased promoter ultrasensitivity that is specific for repression. We then use stochastic simulationto show that multiple binding sites increase the coherence of oscillations by mitigating the binarynoise. Slow values of DNA unbinding rate are also effective in alleviating molecular noise due tothe increased distance from the bifurcation point. Our work demonstrates how the number of DNAbinding sites and slow unbinding kinetics, which are often omitted in biophysical models of genecircuits, can have a significant impact on the temporal and stochastic dynamics of genetic oscillators.
I. INTRODUCTION
Genetic oscillatory networks are ubiquitous in natureand perform important functions. For example, thecell cycle oscillator regulates cell growth and division,whereas the circadian clock regulates the behavior of or-ganisms with respect to daily changes in light. Thesegenetic oscillators are used by living systems to reliablycoordinate various periodic internal processes with eachother as well as with their rhythmic environment. How-ever, this presents a challenge at the cellular level becauseoscillators have to maintain proper timing (temporal co-herence of oscillation) in the presence of stochastic noisethat arises from the small number of regulatory moleculesin cells [1].A simple mechanism to mitigate the effect of molecu-
A I R I (a) (b)
FIG. 1. (a) In the Activator-Titration Circuit (ATC), the ac-tivator is constitutively produced at a constant rate and acti-vates the expression of the inhibitor, which, in turn, titratesthe activator into inactive complex. (b) In the Repressor-Titration Circuit (RTC) the constitutively expressed inhibitortitrates the self-repressing repressor. lar noise would be to increase the number of molecules ofeach species [2–4]. While the number of RNAs and pro-teins made per gene can be large, most cells are funda-mentally constrained to 1-2 gene copies and are subject tobinary noise in the first step of gene regulation (i.e., tran-scription factor binding to DNA) [5, 6]. This binary generegulation noise manifests itself as a stochastic temporalpattern of all-or-none gene activity depending on whetherthe promoter is bound by the regulatory protein or not.Recent work shows that slow DNA binding/unbindingkinetics (also called the non-adiabatic limit) can exacer-bate the binary noise and have profound consequenceson gene expression [7], epigenetic switching [8], and os-cillation [3, 4, 9, 10]. Faster kinetic rates and complexgene promoter architectures have been proposed as a wayto suppress the effect of this binary noise. For example,increasing the DNA binding/unbinding rate can increasetemporal coherence of oscillations via more precise sam-pling of the concentration of transcription factors [3, 9–11] or by increasing the distance from a bifurcation point[4, 12]. However, transcription factors often have slowDNA unbinding rates [13–17], which suggests that thesemechanisms are not generally applicable. The coopera-tive binding of a transcription factor to multiple bindingsites has also been shown to increase temporal coherenceof oscillations [18]. However, multiple binding sites donot always lead to cooperativity and transcription fac-tor binding to a single DNA site may often be enough toeffectively activate or repress transcription.To better understand the potential mechanisms thatsuppress the binary gene regulation noise, in particularthe influence of slow DNA unbinding rates and multi-ple binding sites, we study an Activator-Titration Cir- a r X i v : . [ phy s i c s . b i o - ph ] J a n cuit (ATC) that has been theoretically shown to oscil-late [19]. The ATC consists of a constitutively-expressedactivator that promotes the expression of the inhibitor,which then titrates the activator into an inactive het-erodimer complex (Fig. 1). Studying the ATC has twoadvantages. First, it lies at the core of animal circa-dian clocks [20] and oscillatory NF- κ B signaling [21, 22]and has served as a model of natural genetic oscillators[10, 19, 23–26]. Second, the ATC generates the necessarynonlinearities through protein titration [27] and does notrequire cooperative binding of activator to the inhibitorpromoter. Thus, by studying a titration-based oscilla-tor, we can better explore the kinetic effects of multiplebinding sites on coherence independently of the effectsthat might arise from cooperativity. To obtain generalinsights that are not specific to activation, we also studya Repressor-Titration Circuit (RTC), which consists ofa self-repressor and a constitutively-expressed inhibitor(Fig. 1). This novel titration-based oscillator is analo-gous to the ATC but uses repression instead of activationfor the transcriptional regulation.We first characterize these oscillators and how theydepend on several key parameters in Section II. We de-liberately constrain ourselves to physiological parame-ters found in a simple eukaryote
S. cerevisiae , commonlyknown as budding yeast. We show that, in additionto slow mRNA degradation, slow DNA unbinding ratesof transcription factors are important for providing thenecessary delay in the negative feedback loop for oscilla-tory solutions. Thus, both the DNA unbinding rate andmRNA degradation rate can set the period of oscillation.We then demonstrate that the addition of multiple, inde-pendent binding sites has nontrivial effects on the ATCand the RTC. While multiple binding sites lengthen theperiod of both oscillators due to an effective increase inthe delay of negative feedback, they dramatically increasethe oscillatory solution space of the RTC only. This isbecause multiple, independent binding sites generate ul-trasensitivity (i.e. strong nonlinear response to changesin regulatory protein concentration) in repression-basedpromoters only, and thus only RTC can benefit from thiseffect. In section III, we use stochastic Gillespie simula-tions to understand the extent to which DNA unbindingrates and numbers of binding sites suppress the molec-ular noise in ATC and RTC oscillators. We show thatmultiple binding sites increase the temporal coherence ofoscillations by alleviating the binary noise resulting fromdiscrete gene states. We also show that slower values ofDNA unbinding rates are best for coherent oscillations insimple titration-based oscillations. Last, we compare andcontrast our results on temporal coherence with those ofprevious models of genetic oscillators in Section IV.
II. BIOPHYSICAL MODEL OF ATC AND RTC
Oscillators require negative feedback with nonlinearityand time delays [28]. Mechanistically, negative feedback on gene expression can occur transcriptionally via repres-sors [29–31] or post-transcriptionally via phosphorylation[32–34], degradation [32, 34, 35], or titration of activa-tors into inactive complexes by inhibitors [10, 19, 23–26].The ATC is a minimal two-gene circuit that can oscillatethrough the use of protein titration both as a source ofnonlinearity and indirect negative feedback. In the firstphase of oscillation, high levels of free activator homo-dimerize, bind the promoter, and overproduce inhibitoruntil all free activator has been titrated into inactive het-erodimer. In the second phase of oscillation, the remain-ing activator will unbind from the inhibitor promoter andbe sequestered by inhibitor, thus causing the promoter toreturn to low levels of expression of inhibitor. The levelsof inhibitor will eventually decline to a point where freeactivator can re-accumulate and restart the cycle.In the RTC, protein titration is used exclusively asa source of nonlinearity and the negative feedback isdirectly achieved through auto-repression. In the firstphase of oscillation, high levels of free repressor willhomo-dimerize, bind directly to its own promoter, and re- (a)(b)
Activator
A AA A AI I + δ p ρ βδ m γ ε δ p Inhibitor
I I IA A A + A δ p ρ f , ρ b β δ m γε α θδ p Repressor R R R + R δ p ρ f , ρ b β δ m γε α θδ p R RI I + γ ε δ p Inhibitor δ p ρ βδ m I I IR RR
FIG. 2. A biophysical model for ATC (a) and RTC (b)with explicit transcription, translation, protein-protein andprotein-DNA interactions. Each arrow corresponds to a reac-tion rate in Eqs. 1-3. Neither of these titration-based oscilla-tors have been built or studied by synthetic biologists. press its transcription. The free repressor will be titratedaway by the constitutively expressed inhibitor. In thesecond phase of oscillation when free repressor levels arelow, the remaining repressor will unbind from the pro-moter, returning to high levels of transcription and therapid over-production of free repressor. As we will showbelow, the indirect versus direct nature of negative feed-back in ATC and RTC is responsible for many of thedifferences between these two titration-based oscillators.
A. ATC and RTC oscillators with a single DNAbinding site
Even simple genetic circuits such as the ATC and RTCinclude many reactions and parameters (Fig. 2). An ex-haustive search over all the parameter space was not fea-sible, and we decided to constrain our parameter spaceby studying synthetic gene circuits that could be builtin budding yeast. Synthetic genetic oscillators have beenuseful tools to understand the properties of natural oscil-lators. For example, a synthetic oscillator built in bacte-ria [36] was useful in understanding entrainment capabili- ties of genetic oscillators, as well as elucidating sources ofstochasticity that affected entrainment. Surprisingly, allsynthetic genetic oscillators built to date have neglectedprotein titration, a common mechanism in natural os-cillators. To this end, we built a mathematical modelof ATC and RTC oscillators using a basic leucine zip-per (bZIP) transcription factor that dimerizes and bindsDNA, and a rationally-designed inhibitor that binds freebZIP into an inactive heterodimer. These synthetic com-ponents have been successfully used in yeast [37] and,importantly, many of the protein-protein and protein-DNA binding affinities of this bZIP and inhibitor pairare known [38, 39]; see Table I. We fixed these param-eters and scanned through a range of other biophysicalparameters to understand which ones affect oscillation.Our results should help guide future experimental im-plementation of synthetic ATC and RTC oscillators inyeast.The biophysical model of our ATC and RTC circuitsis based on chemical mass-action kinetics where the dy-namic variables are the mean concentrations of all molec-ular species. The ODEs that correspond to the reactionsin Fig. 2 are the following: d [ G ] dt = − α [ G ][ X ] + θ [ G ] (1a) d [ G ] dt = α [ G ][ X ] − θ [ G ] (1b) d [ I ] dt = β [ r I ] − γ [ X ][ I ] + (cid:15) [ XI ] − δ p [ I ] (1c) d [ X ] dt = β [ r X ] − γ [ X ][ I ] + (cid:15) [ XI ] − γ [ X ] + 2 (cid:15) [ X ] − δ p [ X ] (1d) d [ XI ] dt = γ [ X ][ I ] − (cid:15) [ XI ] − δ p [ XI ] (1e) d [ X ] dt = γ [ X ] − (cid:15) [ X ] − δ p [ X ] − α [ G ][ X ] + θ [ G ] (1f)With [ r X ] and [ r I ] described by: d [ r X ] dt = ρ [ G T ] − δ m [ r X ] (2a) d [ r I ] dt = ρ f [ G ] + ρ b ([ G T ] − [ G ]) − δ m [ r I ] (2b)for the ATC, where X=A (activator), and d [ r X ] dt = ρ f [ G ] + ρ b ([ G T ] − [ G ]) − δ m [ r X ] (3a) d [ r I ] dt = ρ [ G T ] − δ m [ r I ] (3b)for the RTC, where X=R (repressor). The first twoequations represent the dynamics of promoter DNA where [ G ] and [ G ] are the mean concentrations of freeand bound promoter, respectively. The molar concentra-tion of total DNA [ G T ] = [ G ]+[ G ] = 1 / ( N A · V ) = 1 / N A is the Avogadro constant and V is theyeast cell volume; see Table 1. Here, we consider onlya single DNA binding site, but we will later expand ouranalysis to include multiple binding sites. At any in-stant, the promoter is either free or bound. The proba-bility of free or bound promoters is equal to the ratio ofconcentrations [ G ] / [ G T ] or [ G ] / [ G T ], respectively. Theother equations describe the mean concentration dynam-ics of the respective molecular species such as mRNA( r I , r X ), monomeric protein ( I or X ), and dimeric pro-teins ( X , XI ), where X stands for the activator A orrepressor R , respectively. The regulatory homodimer X associates with G at a rate α to form G , which disso- FIG. 3. (Color online) Parameter space of oscillatory solutions on a logarithmic scale for RTC (top) and ATC (bottom) withincreasing DNA binding sites. The colormap shows the number of ρ f values that exhibited oscillations for each combinationof f , δ m , and θ . (a) and (b), single DNA binding site for RTC and ATC, (c) and (d) three independent binding sites for RTCand ATC, (e) and (f), three synergistic binding sites for RTC and ATC, where activation/repression strength is f when morethan one activator/repressor is bound. ciates at the rate θ . The r I and r X are the inhibitor andactivator/repressor mRNAs. For the ATC, the activa-tor mRNA [ r X ] is transcribed constitutively at the rate ρ , where as the inhibitor mRNA is transcribed at rates ρ f and ρ b from free and bound DNA, respectively (Eqs.2a,b). In contrast, for the RTC, the repressor mRNA istranscribed at rates ρ f and ρ b (Eqs. 3a,b) while inhibitormRNA r I is constitutively transcribed at the rate ρ . Weassume that all mRNA species are degraded at the samerate δ m and translated into proteins with the same rate β . The activator/repressor X protein dimerizes into ac-tive homodimer X and forms inactive heterodimer XI with the inhibitor protein I at the same rate γ . The ho-modimer and heterodimer dissociation rates are (cid:15) and (cid:15) , respectively. We assume that all protein species arestable and diluted by cell growth at rate δ p . B. DNA unbinding kinetics influence oscillation
Our parameters were restricted to physiological val-ues from yeast (see Table I). Most parameters were keptfixed, but we varied four key parameters. The first pa-rameter was the mRNA production rate ( ρ f ) of free, un-bound promoter because a desired expression level caneasily be selected from existing promoter libraries [40].Second, we varied the activation/repression strength ( f ),which is the ratio of the larger ρ divided by the smaller ρ . Thus, f = ρ b /ρ f for the ATC and f = ρ f /ρ b for the RTC. The ratio f can be tuned by appropriate choice ofactivation or repression domains fused to our bZIP tran-scription factor [41–43]. The third parameter was themRNA degradation rate ( δ m ), which is known to set thetime scale of the ATC oscillator [19]. Last, we varied theDNA unbinding rate ( θ ) because it is our point of focusand this parameter can vary between different transcrip-tion factors. The DNA dissociation constant ( K d ) fixesthe DNA binding rate α = θK d ; see Appendix for details.We divided the physiological range of each variable pa-rameter into 30 values (on a logarithmic scale) and evalu-ated the long-term dynamics of a total of (30) parametersets per circuit. We solved the ODEs over time for eachset of ( ρ f , f , δ m , θ ). A solution was classified as oscil-latory if the trough of activator or repressor homodimerconcentration was below the K d of DNA-binding and ifthe peak was above 2 K d ; see Appendix for justification.We noticed that ρ f had the smallest effect on the numberof oscillatory solutions and, thus, we plot the marginalfrequency distribution of oscillatory solutions over f , δ m ,and θ in Fig. 3. We see that strong activators (large f forthe ATC), stable mRNAs (small δ m ), and slow DNA un-binding rates (small θ ) generally favor oscillation. Thelast two parameters dictate the timescale of the delayin the negative feedback loop. Increased delay supportsoscillation and, thus, the largest number of oscillatory so-lutions occur at the smallest θ and δ m for both RTC (Fig.3a) and ATC (Fig. 3b). The parameter space of stableoscillations is larger in ATC relative to RTC for a single FIG. 4. The DNA unbinding rate θ sets the period of theoscillations for RTC (a) and ATC (b) at slow unbinding rates.The mean period of oscillatory solutions for a given θ is shown(solid black line) with the shaded area representing the rangeof periods. binding site because of the additional step and delay inthe negative feedback loop: negative feedback throughthe activator in the ATC is indirect (i.e. activator regu-lates the expression of inhibitor, which then inhibits itsactivity), where as the self-repressor in the RTC is direct(i.e. repressor regulates its own expression).The period of oscillation τ should be set by thetimescale of the slowest parameters in the delay. Thenegative feedback in our circuits is dominated by DNAunbinding rate θ and mRNA degradation rate δ m [19].This relationship can be seen in Figure 4 where the DNAunbinding rate sets the oscillation period at low θ . Anincrease in θ leads to mRNA degradation rate ( δ m ) be-coming the slower timescale at which point τ becomesflatter and less dependent on θ . Eventually, a bifurca-tion occurs at a critical, maximum value of θ max whichleads to loss of the stable limit cycle. A similar relation-ship exists for the mRNA degradation rate δ m ; see FigureS2. C. Multiple DNA binding sites affect ATC andRTC oscillators differently
This role of DNA unbinding rate in generating de-lays led us to hypothesize that multiple DNA bindingsites should increase the parameter space of oscillationsand lengthen the period of the oscillators. We reasoned X X X ρ b G G G G Same Transcription Rate X X ρ b X X ρ b X X ρ b X ρ b X ρ b X ρ b ρ f FIG. 5. Transitions between promoter states for multipleDNA binding sites (n=3). G i denotes the set of promoterstates with i out of total n binding sites occupied by activa-tor (X=A) or repressor (X=R) dimers. There are n − i ways ofswitching from G i to [ G i +1 ] via the binding of X , and thereare i ways of switching from G i to G i − by the unbinding of X . Our model conservatively assumes that the binding of X does not affect the binding or unbinding of the next tran-scription factor to an adjacent site (no cooperativity). Wealso assume that the transcription rate is equal to ρ b for G , G , . . . G n promoter states. that if the occupancy of any binding site by a transcrip-tion factor activates or represses transcription, then theeffective unbinding rate ( θ n ) from a state of saturatedDNA binding to the unbound DNA state ( G , where thetranscription rate changes) should decrease with the in-creasing number of binding sites ( n ). We can show that θ n = θ/H n , where H n is the n-th harmonic number (seeSupplement [44]).The addition of multiple DNA binding sites to ourmodel will modify Eqs. (1a-b, 1f) by increasing the num-ber of promoter states that can be bound by X ; see Fig-ure 5 and Supplement [44]. For three binding sites (n=3),our new Eqs (1a-b) are: d [ G ] dt = − α · [ G ][ X ] + θ [ G ] d [ G ] dt = 3 α · [ G ][ X ] − ( θ + 2 α · [ X ])[ G ] + 2 θ · [ G ] d [ G ] dt = 2 α · [ G ][ X ] − (2 θ + α · [ X ])[ G ] + 3 θ · [ G ] d [ G ] dt = α · [ G ][ X ] − θ · [ G ] (4)where the total concentration of DNA [ G T ] = [ G ] +[ G ] + . . . + [ G n ] is fixed to 1 / ( N A · V ) = 1 /
24 nM. Forthree binding sites (n=3), the term − α [ G ][ X ] + θ [ G ]in Eq. (1f) is replaced with: − X i =0 (3 − i ) α · [ G i ][ X ] + X i =0 i · θ [ G i ] (5) G i denotes the promoter states with i out of total n binding sites occupied by activator (X=A) or repressor(X=R) dimers. The factors in front of each term repre-sents the amount of degeneracy of each state, i.e [ G i ] has i bound sites, thus i ways of switching to [ G i − ]. There-fore, we have the term i · θ [ G i ]. At the same time, [ G i ] has n − i vacant sites, so it has n − i ways of switching to thestate [ G i +1 ]. Thus, we have the term ( n − i ) α · [ G i ][ X ].The addition of two more independent DNA bindingsites dramatically increased the oscillatory space of theRTC (Fig. 3c), while slightly decreasing the oscillatoryspace of the ATC (Fig. 3d). These opposite results arisefrom a compound effect. First, two extra binding sitesdecreased the effective unbinding rate for the promoterto be fully vacated by half ( θ/H ≈ θ/ θ increased the delay and resulted in some im-provement in oscillations for both RTC and ATC. Thiseffect is best observed in the increased period of bothoscillators (Fig. 4). The second, more dominant ef-fect is the fundamental difference in how the promotersensitivity changes with multiple, independent bindingsites. It is well established that nonlinear promoter re-sponses facilitate oscillation [28]. We use the logarithmicsensitivity ( S ) to quantify the nonlinearity in the pro-moter response, where S = dlog P/ dlog[ X ] [45]. P isthe synthesis rate of the promoter and [ X ] is the acti-vator/repressor homodimer concentration that regulatesthe promoter. As shown previously [45, 46], an increasein the number of independent repressive binding sites willincrease the magnitude of S and create an ultrasensi-tive promoter response, (i.e. | S | >
1, see Supplement[44]). However, increasing the number of independentactivating binding sites cannot generate an ultrasensi-tive promoter response ( | S | ≤ G and G ) are acti-vated/repressed f -fold instead of f -fold because theycan interact with RNA polymerase at several interfaces[45]. Although this synergy increased the activation orrepression strength, it did not significantly change theoscillatory parameter space (Fig. 3e,f). III. STOCHASTIC SIMULATIONS
Deterministic simulations were useful for understand-ing how DNA unbinding rate and the number of bind-ing sites affect the phase space and period of oscillation.However, they cannot provide insights into the loss oftemporal coherence that arises from stochastic gene ex-pression. To this end, we used the Gillespie algorithm
Time (min) r R ( n M ) (a) − − Delay (min) C o rr e l a t i on (b) Time (min) r I ( n M ) (c) − − Delay (min) C o rr e l a t i on (d) FIG. 6. Sample stochastic trajectories for one and three bind-ing sites for RTC (a) and ATC (c), and their autocorrelationfunctions (b,d). The variable parameter values ( θ , δ m , ρ f , f )were fixed to (0.02 min − , 0.0159 min − , 0.8928 min − , 3.63)for the RTC and (0.02 min − , 0.0186 min − , 0.1781 min − ,30) for the ATC. We chose parameters that produced oscilla-tion over the largest range of DNA unbinding rates. The restof the parameters are given in Table 1. [47] to simulate stochastic chemical reaction kinetics andinvestigate how DNA binding/ unbinding dynamics andthe addition of binding sites affect the temporal coher-ence of ATC and RTC oscillators.For each ATC and RTC, we quantified temporal co-herence by calculating the autocorrelation function ofmRNA transcripts levels (repressor mRNA for the RTCand inhibitor mRNA for the ATC); see Fig. 6. In the ab-sence of noise, an undamped oscillatory signal will havean undamped, periodic autocorrelation function. Thepresence of noise will stochastically perturb period andphase, such that the autocorrelation now exhibits damp-ening or loss of temporal coherence. We quantified theloss of coherence by measuring the rate of exponentialdecay ( e − t/τ ) of the envelope of a periodic (cos(2 πt/τ ))autocorrelation function (see Appendix for details). Sim-ilar to other studies [4, 18], our metric for temporal co-herence is the normalized autocorrelation function decayrate, which is the ratio of timescales τ /τ . A larger ra-tio indicates better temporal coherence. We varied theDNA unbinding rate ( θ ) and number of binding sites ( n )to understand the role of each feature in resisting molec-ular noise. A. DNA Unbinding Rate
Our results show that ATC and RTC oscillators withsmaller DNA unbinding rates exhibit better temporal co-
FIG. 7. The normalized autocorrelation function decay ratefor the RTC (a) and ATC (b) for varying θ and number ofbinding sites. All parameters, except θ , are the same as inFig. 6. Boxed, outlined regions are parameters past thebifurcation point ( θ max ) where deterministic oscillations areunsustainable and damped, yet exhibit stochastic excitableoscillations. herence (Fig. 7). Lower θ increases the temporal co-herence of the oscillations because of the increased dis-tance of the dynamical system from the bifurcation point( θ max ); see Figure 8. Eventually there is another bifur-cation at small θ min , but these unbinding rates are un-physiological and do not affect our conclusions regard-ing biophysical ATC and RTC oscillators. Strikingly,some θ > θ max , which do not show deterministic oscil-lation, exhibit oscillation in the presence of noise. Thisphenomenon is consistent with coherence resonance [48]which has been observed in other excitable, genetic cir-cuits near oscillatory bifurcation points [12, 24]. B. Multiple DNA Binding Sites
Increasing the number of binding sites ( n ) also in-creased the temporal coherence of ATC and RTC oscil-lators over all DNA unbinding rates (Fig. 7). To bet-ter understand this result, we must consider the effectof stochastic binding and unbinding of regulators on thevariance of gene expression. In the phase of changing ac-tivator/repressor concentrations, the binding sites startbeing occupied or vacated. Each additional binding siteintroduces an additional DNA binding state. Becausewe treat the expression level of all bound DNA states asequivalent ( ρ b ), the spontaneous binding and unbindingevents that occur between states that have at least onebinding site occupied have no effect on transcription; seeFig. 5. These “protected” states act as a buffering mech-anism to mitigate the effects of binary noise on temporalcoherence. IV. DISCUSSION
We analyzed the properties of two titration based ge-netic oscillators, the activator-titration circuit (ATC)and the repressor-titration circuit (RTC). The focus ofour study was to understand how the number of DNA − − − − (a) (cid:101) (min − ) r R ( n M ) − − − − (b) (cid:101) (min − ) r I ( n M ) FIG. 8. Bifurcation diagram of the RTC (a) and ATC (b)oscillators as a function of DNA unbinding rate ( θ ). All pa-rameters, except θ , are the same as in Fig. 6. There are twobifurcation points ( θ max , θ min ) and the amplitude of mRNAoscillation is shown by the upper and lower branches. Phys-iological values of θ are to the right of the dashed verticalline. binding sites and slow unbinding kinetics in promotersmitigate or exacerbate the binary gene regulation noise.First, we showed that multiple DNA binding sites andslow unbinding kinetics were important for providing thenecessary delay in the negative feedback loop for oscilla-tory solutions. The role of slow DNA binding/unbindingin providing delay for oscillations is consistent with priorwork on a small negative feedback oscillator [9]. Sec-ond, we used stochastic simulation to show that slowerDNA unbinding rates exhibited better temporal coher-ence, a result which appears at odds with previous workon circadian clocks and NF- κ B oscillators [3, 10, 12] andwhich is more in line with the results obtained for a smallnegative feedback oscillator model [9]. This previouswork showed that slower DNA unbinding kinetics neg-atively affected the temporal coherence for two reasons.First, slow DNA unbinding increased the stochasticity ofgene expression due to imprecise concentration sampling,which decreased the temporal coherence of oscillations[3, 10]. Second, slower DNA unbinding ( θ ) pushed thedynamical system towards θ min bifurcation point, whichmade it less robust to noise [12]. These results are differ-ent from ours because the delays in the circadian clockand NF- κ B models rely on slow intermediate steps (e.g.phosphorylation and/or nuclear transport) in the nega-tive feedback loop. Unlike our titration-based oscillatorsin Fig. 8, these models do not have θ max and still oscil-late at infinitely fast unbinding rates where the promoterdynamics are in steady-state.We observed the opposite effect for our titration-basedoscillators because physiological θ overlaps the θ max bi-furcation point for ATC/RTC. Thus, lowering θ alwaysincreases the robustness in ATC/RTC because the dy-namical system is moving away from θ max and deeper intooscillatory parameter space. This phenomenon likely ex-plains the similar results presented in [9]. The influenceof DNA unbinding rate on temporal coherence dependson the structure of the underlying bifurcation diagram ofeach oscillator as a function of θ . Changes in topology,mechanism, and parameters can change the bifurcationdiagram and, thus, the influence of DNA unbinding rateon temporal coherence of oscillation may also change.Last, we demonstrate that multiple independent bind-ing sites consistently increased the temporal coherenceof oscillations by alleviating the binary noise result-ing from binary gene states. Our results agree withprevious work, which showed that multiple, coopera-tive DNA binding sites increased the coherence of cir-cadian clocks [18]. However, in contrast to our results,the temporal coherence of circadian clocks peaked at 3binding sites and then decreased with additional sites.The difference likely arises from our slower DNA bind-ing/unbinding rates, where the ATC and RTC oscilla-tors spend significant time in protected states that arebuffered against molecular noise. In contrast, the circa-dian clock model spends very little time in the interme-diate protected states between fully free or fully boundpromoters and, therefore, the increased coherence is onlydue to cooperativity [18]. The idea of buffering to re-duce noise in gene circuits has been discussed in the con-text of decoy binding sites [49]. However, this requiresfast DNA binding/unbinding, where as buffering throughpromoter states requires slow DNA binding/unbinding.We note that increased temporal coherence due to pro-tected states is a stochastic effect because the additionof binding sites consistently increased the coherence ofATC oscillators, despite occasionally pushing it past thebifurcation point at θ max (Fig. 6b). Appendix A: Parameter Values
To constrain the physiological parameters of our mod-els, we used data from large-scale studies of the yeasttranscriptome and proteome; see Table I. These data pro-vide typical ranges and values for our parameters. First,we converted numbers of molecules into nanomolar (nM)concentrations using the cell volume V = 40 fL for hap-loid yeast. For the ATC, we assumed that the basalmRNA transcription would be low. Thus, ρ f for the ATCwas constrained to values from the bottom 5th percentileto the median of all mRNA synthesis rates [50]. Simi-larly, ρ f for the RTC was constrained to values from themedian to top 95th percentile. In the case of the ATC,the constraint ρ f < ρ < ρ b ensured that the inhibitorcan completely titrate the constitutively-expressed acti-vator when the inhibitor is maximally produced at ρ b , butnot when it is expressed at the basal rate ρ f . Similarly, for the RTC, the constraint ρ b < ρ < ρ f ensures thatconstitutively-expressed inhibitor can completely titratethe repressor when the repressor is produced at the re-pressed rate ρ b , but not ρ f . We set ρ = √ ρ f ρ b to satisfyboth conditions. mRNA degradation rate ranged fromthe bottom 5th percentile to top 95th percentile valuesfor all genes [50]. To obtain a rough approximation forthe translation rates, we assumed a constitutive gene ex-pression model for all genes: d [ r ] dt = ρ − δ m [ r ] (A1) d [ P ] dt = β [ r ] − δ p [ P ] (A2)At steady state, β = [ P ] δ p δ m ρ . Protein concentrationsand degradation rates were taken from [51, 52]. We cal-culated β for all genes and used the median value in ourmodel. We also assumed that our activators, repressors,and inhibitors are not actively degraded and are dilutedby growth. Thus, δ p = ln (2) /T , where T = 90 mins isthe duration of the yeast cell cycle. The proteins in ourmodels were based on a mammalian transcription fac-tor basic leucine zipper (bZIP) protein C/EBP α and itsdominant-negative inhibitor (3HF) [37]. We used previ-ously measured rates for protein-protein interaction ki-netics [37]. Since we did not know the DNA unbind-ing rate for C/EBP α , we considered the range for theknown DNA unbinding rates for other bZIP proteins[13–16]. The thermodynamic dissociation constant ( K d )of C/EBP α to its specific DNA binding site is known[39]. We set the DNA association rate to α = θ/K d .Finally, we varied the activation/repression strength f from 1 to 30, to consider both strong and weak activa-tion/repression. Appendix B: Methods
We scanned the parameter space for oscillations byrunning simulations on MATLAB (Mathworks) usingode15s for 2000 min and recording the minima and max-ima of the activator/repressor homodimer during the last1000 min. We imposed the following restrictions: 1) thelast minimum should be below K d so that DNA-bindingis not saturated, 2) the last maximum should be greaterthan 2 K d so that the change in transcription is noticeablyaltered. While this restriction slightly underestimates thenumber oscillatory solutions, it ensures that a syntheticversion of these circuits would produce detectable oscilla-tions. We verified that our definition gave similar resultsto a less stringent criterion for oscillation.We used the direct Gillespie method to perform thestochastic simulations [53]. We ran the simulations for10 min and recorded the concentration of the regu-lated mRNA (inhibitor for the ATC and the repres-sor for the RTC). We then normalized the concentra-tion such that the average would be zero and evaluated TABLE I. Parameter valuesParameter Min Max Reference θ (min − ) 0.0188 34.5 [13–16]ATC ρ f (min − ) 0.0509 0.1781 [50]RTC ρ f (min − ) 0.1781 0.8928 [50] δ m (min − ) 0.0159 0.1516 [50] f ρ b (min − ) f · ρ f RTC ρ b (min − ) ρ f /fα (nM − min − ) θ/ .
344 nM [39] ρ (min − ) √ ρ f ρ b β (min − ) 14.1 [50–52] δ p (min − ) 0.0077 γ (nM min − ) 0.6 [37] (cid:15) (min − ) 6 [37] (cid:15) (min − ) 0.024 [37] V (fL) 40 [27] the autocorrelation function. We then fit the function C ( t ) = e − t/τ · cos(2 πt/τ ) to the first 1500 min of theautocorrelation function to measure the decay constant τ and period τ . The ratio τ /τ describes how rapidlythe envelope of autocorrelation function decays per oscil-lation period. ACKNOWLEDGMENTS
We thank Joshua Socolar for helpful comments. Thiswork was funded by an NIH Director’s New InnovatorAward (DP2 OD008654-01) and Burroughs WellcomeFund CASI Award (BWF 1005769.01). An SBML ver-sion of the ATC and RTC oscillators can be downloadedfrom the BioModels Database (MODEL1512100000 andMODEL1512100001, respectively). [1] E. Schr¨odinger,
What is life? The physical aspect of theliving cell (Cambridge University Press, 1944).[2] D. Gonze, J. Halloy, and A. Goldbeter, J. Biol. Phys. , 637 (2002).[3] D. B. Forger and C. S. Peskin, Proc Natl Acad Sci USA , 321 (2005).[4] D. Gonze and A. Goldbeter, Chaos , 026110 (2006).[5] I. Golding, J. Paulsson, S. M. Zawilski, and E. C. Cox,Cell , 1025 (2005).[6] L. Cai, N. Friedman, and X. S. Xie, Nature , 358(2006).[7] J. E. M. Hornos, D. Schultz, G. C. P. Innocen-tini, J. Wang, A. M. Walczak, J. N. Onuchic, andP. G. Wolynes, Phys. Rev. E , 051907 (2005).[8] A. M. Walczak, J. N. Onuchic, and P. G. Wolynes, ProcNatl Acad Sci USA , 18926 (2005).[9] H. Feng, B. Han, and J. Wang, Biophys. J , 1001(2012).[10] D. A. Potoyan and P. G. Wolynes, Proc Natl Acad SciUSA , 2391 (2014).[11] D. Labavi´c, H. Nagel, W. Janke, and H. Meyer-Ortmanns, Phys. Rev. E , 062706 (2013).[12] D. Gonze, J. Halloy, and A. Goldbeter, Physica A ,221 (2004).[13] H. Kwon, S. Park, S. Lee, D.-K. Lee, and C.-H. Yang,Eur. J. Biochem. , 565 (2001).[14] M. Kyo, T. Yamamoto, H. Motohashi, T. Kamiya,T. Kuroita, T. Tanaka, J. D. Engel, B. Kawakami, andM. Yamamoto, Genes Cells , 153 (2004).[15] Y. Okahata, K. Niikura, Y. Sugiura, M. Sawada, andT. Morii, Biochemistry , 5666 (1998).[16] M. Geertz, D. Shore, and S. J. Maerkl, Proc Natl AcadSci USA , 16540 (2012).[17] P. Hammar, M. Walld´en, D. Fange, F. Persson, ¨O. Bal- tekin, G. Ullman, P. Leroy, and J. Elf, Nat Genet ,405 (2014).[18] D. Gonze, J. Halloy, and A. Goldbeter, Proc Natl AcadSci USA , 673 (2002).[19] P. Fran¸cois and V. Hakim, Phys. Rev. E. , 031908(2005).[20] J. Menet, K. Abruzzi, J. Desrochers, J. Rodriguez, andM. Rosbash, Genes Dev , 358 (2010).[21] A. Hoffmann, A. Levchenko, M. L. Scott, and D. Balti-more, Science , 1241 (2002).[22] D. E. Nelson, A. E. C. Ihekwaba, M. Elliott, J. R. John-son, C. A. Gibney, B. E. Foreman, G. Nelson, V. See,C. A. Horton, D. G. Spiller, S. W. Edwards, H. P. Mc-Dowell, J. F. Unitt, E. Sullivan, R. Grimley, N. Benson,D. Broomhead, D. B. Kell, and M. R. H. White, Science , 704 (2004).[23] N. Barkai and S. Leibler, Nature , 267 (2000).[24] J. M. G. Vilar, H. Y. Kueh, N. Barkai, and S. Leibler,Proc Natl Acad Sci USA , 5988 (2002).[25] S. Krishna, M. H. Jensen, and K. Sneppen, Proc NatlAcad Sci USA , 10840 (2006).[26] J. K. Kim and D. B. Forger, Mol Syst Biol (2012).[27] N. E. Buchler and M. Louis, J Mol Biol , 1106 (2008).[28] B. Novak and J. J. Tyson, Nat Rev Mol Cell Biol , 981(2008).[29] M. B. Elowitz and S. Leibler, Nature , 335 (2000).[30] M. R. Atkinson, M. A. Savageau, J. T. Myers, and A. J.Ninfa, Cell , 597 (2003).[31] J. Stricker, S. Cookson, M. R. Bennett, W. H. Mather,L. S. Tsimring, and J. Hasty, Nature , 516 (2008).[32] B. Novak and J. J. Tyson, J Theor Biol , 101 (1993).[33] M. Nakajima, H. Ito, and T. Kondo, FEBS Letters ,898 (2010).[34] Q. Yang and J. E. Ferrell, Nature , 519 (2013).[35] W. W. Wong, T. Y. Tsai, and J. C. Liao, Mol Syst Biol , 130 (2007).[36] O. Mondrag´on-Palomino, T. Danino, J. Selimkhanov,L. Tsimring, and J. Hasty, Science , 1315 (2011).[37] N. E. Buchler and F. R. Cross, Mol Syst Biol , 272(2009).[38] D. Krylov, M. Olive, and C. Vinson, EMBO J , 5329(1995).[39] Z. Cao, R. M. Umek, and S. L. McKnight, Genes Dev. , 1538 (1991).[40] L. Keren, O. Zackay, M. Lotan-Pompan, U. Barenholz,E. Dekel, V. Sasson, G. Aidelberg, A. Bren, D. Zeevi,A. Weinberger, U. Alon, R. Milo, and E. Segal, MolSyst Biol (2013).[41] U. Baron, M. Gossen, and H. Bujard, Nucleic Acids Res , 2723 (1997).[42] G. Bell´ı, E. Gar´ı, L. Piedrafita, M. Aldea, and E. Her-rero, Nucleic Acids Res , 942 (1998).[43] P. Perez-Pinera, D. G. Ousterout, J. M. Brunger, A. M.Farin, K. A. Glass, F. Guilak, G. E. Crawford, A. J.Hartemink, and C. A. Gersbach, Nat Methods , 239(2013).[44] See Supplemental Material at URL for the derivation ofpromoter sensitivity, unbinding from multiple sites and additional figures.[45] L. Bintu, N. E. Buchler, H. G. Garcia, U. Gerland,T. Hwa, J. Kondev, T. Kuhlman, and R. Phillips, CurrOpin Genet Dev , 125 (2005).[46] I. M. Lengyel, D. Soroldoni, A. C. Oates, and L. G.Morelli, Papers in Physics (2014).[47] D. Gillespie, J. Phys. Chem. , 2340 (1977).[48] A. S. Pikovsky and J. Kurths, Phys. Rev. Lett. , 775(1997).[49] A. Burger, A. M. Walczak, and P. G. Wolynes, Phys.Rev. E , 041920 (2012).[50] C. Miller, B. Schwalb, K. Maier, D. Schulz, S. D¨umcke,B. Zacher, A. Mayer, J. Sydow, L. Marcinowski,L. D¨olken, D. E. Martin, A. Tresch, and P. Cramer,Mol Syst Biol , 458 (2011).[51] S. Ghaemmaghami, W.-K. Huh, K. Bower, R. W. How-son, A. Belle, N. Dephoure, E. K. O’Shea, and J. S.Weissman, Nature , 737 (2003).[52] A. Belle, A. Tanay, L. Bitincka, R. Shamir, and E. K.O’Shea, Proc Natl Acad Sci USA , 13004 (2006).[53] S. Hoops, S. Sahle, R. Gauges, C. Lee, J. Pahle, N. Simus,M. Singhal, L. Xu, P. Mendes, and U. Kummer, Bioin-formatics , 3067 (2006). upplementary Material Sargis Karapetyan
1, 2 and Nicolas E. Buchler
1, 2, 31
Department of Physics, Duke University, Durham, NC 27708 Center for Genomic & Computational Biology, Durham, NC 27710 Department of Biology, Duke University, Durham, NC 27708 a r X i v : . [ phy s i c s . b i o - ph ] J a n . PROMOTER WITH MULTIPLE BINDING SITES We first derive the chemical kinetic equations for 2 binding sites and then generalize to n binding sites [1]. Similar to previous work [2], we calculate the promoter sensitivity ofactivators and repressors in the limit of thermodynamic equilibrium. A. Promoter with 2 equal binding sites
We first begin with an explicit derivation that includes each binding state combination: d [ G ] dt = − α · [ G ][ X ] + θ · [ G ] − α · [ G ][ X ] + θ · [ G ] (1) d [ G ] dt = α · [ G ][ X ] − θ · [ G ] − α · [ G ][ X ] + θ · [ G ] (2) d [ G ] dt = α · [ G ][ X ] − θ · [ G ] − α · [ G ][ X ] + θ · [ G ] (3) d [ G ] dt = α · [ G ][ X ] − θ · [ G ] + α · [ G ][ X ] − θ · [ G ] (4)where [ G T ] = [ G ] + [ G ] + [ G ] + [ G ] and [ X ] is the homodimer transcription factorconcentration. We rewrite these equations in a compact form: d [ G ] dt = − α · [ G ][ X ] + θ [ G ] (5) d [ G ] dt = 2 α · [ G ][ X ] − θ [ G ] − α [ G ][ X ] + 2 θ [ G ] (6) d [ G ] dt = α · [ G ][ X ] − θ [ G ] (7)where [ G ] = [ G ], [ G ] = [ G ] + [ G ], and [ G ] = [ G ]. This simplification is validwhen all binding sites have equivalent kinetics and the regulator X binds independently toeach site.
1. Promoter synthesis rate at equilibrium
At equilibrium, Eqs. (5-7) are equal to 0. Thus,2[ G ][ X ] = K d · [ G ] (8)( K d + [ X ])[ G ] = 2[ G ][ X ] + 2 K d · [ G ] (9)[ G ] · [ X ] = 2 K d · [ G ] (10)where K d = θ/α . Solving these equations with constraint [ G T ] = [ G ] + [ G ] + [ G ]:[ G ] = [ G T ] · [ X ] K d ) (11)[ G ] = [ G T ] · · ( [ X ] K d )(1 + [ X ] K d ) (12)[ G ] = [ G T ] · ( [ X ] K d ) (1 + [ X ] K d ) (13)The promoter synthesis rate P where any bound state has the same gene expression ( ρ b )and the unbound state has gene expression ( ρ f ) is given by: P ( X ) = [ G T ] · (cid:20) ρ b [2 · X + ( X ) ] + ρ f (1 + X ) (cid:21) = [ G T ] · (cid:20) ρ b + ( ρ f − ρ b )(1 + X ) (cid:21) (14)where X = [ X ] K d . B. Promoter with n equal binding sites
By generalizing our previous derivation for an arbitrary number of binding sites, we get n + 1 differential equations, where n is the total number of binding sites and the index i spans 1 ≤ i ≤ n − d [ G ] dt = − nα · [ G ][ X ] + θ [ G ] (15)... d [ G i ] dt = ( n − i + 1) α · [ G i − ][ X ] − ( i · θ + ( n − i ) α · [ X ])[ G i ] + ( i + 1) θ · [ G i +1 ] (16)... d [ G n ] dt = α · [ G n − ][ X ] − nθ · [ G n ] (17)3 i denotes the promoter states with i out of total n binding sites occupied by activator(X=A) or repressor (X=R) dimers. The factors in front of each term represents the amountof degeneracy of each state, i.e [ G i ] has i bound sites, thus i ways of switching to [ G i − ].Therefore, we have the term i · θ [ G i ]. At the same time, [ G i ] has n − i vacant sites, so it has n − i ways of switching to the state [ G i +1 ]. Thus, we have the term ( n − i ) α · [ X ][ G i ]. Notethat our model assumes that the binding of a transcription factor does not affect the bindingor unbinding of the next transcription factor to an adjacent site (no cooperativity). We alsoassume that transcription is fully activated/repressed for any bound state of the promoter,i.e. transcription is equal to ρ b for G , G , . . . , G n . Although not shown, the ODE for therate of change in dimer concentration Eq. (1f) in the main text is also modified to reflectbinding and unbinding from each state G i .
1. Promoter synthesis rate at equilibrium
Setting Equations (15-17) equal to zero and solving with the constraint [ G T ] = [ G ] +[ G ] + . . . + [ G n ]: [ G ] = [ G T ] · [ X ] K d ) N (18)...[ G i ] = [ G T ] · N ! i !( N − i )! · ( [ X ] K d ) i (1 + [ X ] K d ) N (19)...[ G N ] = [ G T ] · ( [ X ] K d ) N (1 + [ X ] K d ) N (20)The promoter synthesis rate P n is thus given by: P n ( X ) = [ G T ] · (cid:20) ρ b + ( ρ f − ρ b )(1 + X ) n (cid:21) (21)where X = [ X ] K d . 4 . Logarithmic sensitivity analysis We compute the logarithmic sensitivity of the promoter: S = dlog P dlog X = X P ( X ) · dPdX : S ( X ) = n ( ρ b − ρ f ) X (1 + X ) · ( ρ f + ρ b [(1 + X ) n −
1. Activation
We consider the simplified case of ideal activation, where ρ f = 0. Thus, S ( X ) = nX (1 + X ) · [(1 + X ) n −
1] (23) S is a monotonically decreasing function of X with boundaries S (0) = 1 and S ( ∞ ) =0. Thus, activation is never ultrasensitive (i.e. | S | ≤ X , n , or ρ b ).Additionally, | S | is also a monotonically decreasing function of n for all X >
0; see FigureS1. The boundaries are at | S (0) | = X (1+ X ) · ln (1+ X ) and | S ( ∞ ) | = 0. This explains why theATC oscillator phase space decreases with the addition of multiple DNA binding sites.
2. Repression
We consider the case of ideal repression where ρ b = 0. Thus, S ( X ) = − nX (1 + X ) (24) S is a monotonically decreasing function of R with boundaries S (0) = 0 and S ( ∞ ) = − n .Promoter repression is ultrasensitive (i.e. | S | >
1) for X > n − and n >
1; see Figure S1. | S | is a monotonically increasing function of n for all X >
0. The boundaries are at | S (0) | = 0 and | S ( ∞ ) | = ∞ . Thus, the RTC oscillator phase space increases with theaddition of multiple DNA binding site. II. KINETICS OF COMPLETE UNBINDING FROM MULTIPLE SITES
The probability of unbinding for a single binding site is given by a Poisson distribution:5 ( t ) = θe − θt (25)where θ is the unbinding rate. The average unbinding time for a single site is given by < τ > = R ∞ t · p ( t ) dt and is equal to θ . Consider the case where all n binding sites are filledand where the concentration of regulator abruptly changes to 0 (such that no rebindingoccurs). We calculate the probability distribution of time for complete unbinding from all n binding sites. The probability that the slowest unbinding event is τ and the other n − τ is: P ( τ | n, θ ) = (cid:20)Z τ p ( t ) dt (cid:21) n − · p ( τ ) (26) P ( τ | n, θ ) = (cid:2) − e − θτ (cid:3) n − · θe − θτ (27)When we normalize the distribution ( R ∞ P ( τ | n, θ ) dt = 1): P ( τ | n, θ ) = n (cid:2) − e − θτ (cid:3) n − · θe − θτ (28)The average < τ > can be calculated by R ∞ τ · P ( τ ) dτ : < τ > = H n θ = 1 θ n (29)where H n = 1 + + + . . . + n (Harmonic number) and θ n is the effective DNA unbindingrate. 6 II. SUPPORTING FIGURES | S | UltrasensitiveSubsensitiveActivationRepression
FIG. S1. Multiple binding sites affect activation- and repression-based promoters differently. Theabsolute value of logarithmic sensitivity S is plotted as a function of number of binding sites N ,with X = 2. The points of 1 and 3 binding sites considered in this study are marked with circles.FIG. S2. The mRNA degradation rate δ m sets the period of the oscillations for RTC (a) and ATC(b) at slow mRNA degradation rates. The mean period of oscillatory solutions for a given δ m isshown (solid black line) with the shaded area representing the range of periods.
1] I. M. Lengyel, D. Soroldoni, A. C. Oates, and L. G. Morelli, Papers in Physics (2014).[2] L. Bintu, N. E. Buchler, H. G. Garcia, U. Gerland, T. Hwa, J. Kondev, T. Kuhlman, andR. Phillips, Curr Opin Genet Dev , 125 (2005)., 125 (2005).