Tuning Genetic Clocks Employing DNA Binding Sites

Periodic oscillations play a key role in cell physiology from the cell cycle to circadian clocks. The interplay of positive and negative feedback loops among genes and proteins is ubiquitous in these networks. Often, delays in a negative feedback loop and/or degradation rates are a crucial mechanism to obtain sustained oscillations. How does nature control delays and kinetic rates in feedback networks? Known mechanisms include proper selection of the number of steps composing a feedback loop and alteration of protease activity, respectively. Here, we show that a remarkably simple means to control both delays and effective kinetic rates is the employment of DNA binding sites. We illustrate this design principle on a widely studied activator-repressor clock motif, which is ubiquitous in natural systems. By suitably employing DNA target sites for the activator and/or the repressor, one can switch the clock “on” and “off” and precisely tune its period to a desired value. Our study reveals a design principle to engineer dynamic behavior in biomolecular networks, which may be largely exploited by natural systems and employed for the rational design of synthetic circuits.


Introduction
Periodic oscillations are essential for biological phenomena such as cell cycle regulation and circadian rhythms [1,2]. Several studies attribute these oscillations to bio-molecular clocks composed of genes arranged in feedback networks [3,4]. Of the several arrangements that may produce oscillations, activator-repressor motifs are recurrent in several natural systems [3,5]. These motifs comprise an activator module that is self activated and that activates a repressor module. The repressor module, in turn, represses the activator (Figure 1a). This motif has been shown to be remarkably robust to biological noise [5], leading to synthetic implementations as model systems to study natural clocks [6][7][8][9].
Independently of the specific topology of the network, the presence of delays in feedback loops has long been recognized as a key mechanism to obtain periodic behavior and to tune the clock period (see the review by [1] and the study by [10]). Similarly, a key (related) parameter controlling periodic behavior is the relative value among protein decay rates [11,12]. For the activatorrepressor motif, for example, analytical studies have demonstrated that a crucial mechanism for sustained oscillations is the time-scale difference between the activator and the repressor dynamics, that is, the repressor dynamics should be sufficiently slower than the activator dynamics [13,14]. This is, to some extent, qualitatively similar to having a delay in the negative feedback from the repressor to the activator. How does nature realize and tune delays and kinetic rates in feedback motifs? Known ways to increase a delay in a feedback or to make the feedback slower include either decreasing the decay rates of species involved in the negative feedback and/or increasing the number of steps in the feedback loop (see, for example, [10,14,15]).
Recent studies of modularity in biomolecular circuits have revealed that excess of DNA targets to a protein can slow down the protein's dynamics [16][17][18]. This effect, called retroactivity, is a consequence of changes in the dynamics of the system due to the sequestration of the protein from the network of interactions composing the system. Basically, the protein is ''busy'' in binding the targets and hence takes longer to perform its function in the system to which it belongs. In the context of modularly designing circuits in synthetic biology, this is an undesired effect (similar to impedance in electrical circuits) that occurs when two modules are interconnected by a transcription factor of one module binding to DNA target sites in the other module. From the perspective of a natural system, however, this loading effect may provide a simple method to tune delays and change the effective kinetic rates without changing the ''hardware'' of the network.
In this work, we demonstrate that indeed DNA target sites can be employed as a powerful design parameter to finely tune and control the dynamic behavior of a biomolecular circuit, the activator-repressor clock of Figure 1a in particular. Specifically, we illustrate how one can change the dynamics of an activatorrepressor clock utilizing DNA binding sites (load) with affinity to each of the species. Initially, a mechanism to switch an oscillator ''on'' or ''off'' is shown depending on which node (the activator or repressor) the load is being added to. Robustness of this behavior to intrinsic noise is verified by employing stochastic simulation of a mechanistic model of the clock. Finally, a method to tune the period of the clock by employing a carefully chosen amount of load to both nodes is demonstrated.

Results
We consider a general model for a two-component clock incorporating both positive and negative feedback loops based on the activator-repressor configuration of [6] and illustrated in Figure 1a. Oscillations for activator-repressor clocks often arise from Hopf bifurcation, wherein a stable equilibrium point bifurcates into an unstable equilibrium and a stable periodic orbit when a key parameter is changed [9,13,14,19,20]. In the models surveyed in the literature, the fundamental mechanism responsible for this oscillatory behavior is well captured by a reduced twodimensional model that describes the rate of change of the activator and repressor concentrations. This model is obtained by taking into account that the period of oscillations occurs in a timescale slower than the dynamics of multimerization, binding and dissociation interactions, so that quasi-steady state approximations can be made [6,9,20]. Additionally, it has been shown that transcription and translation can be lumped into a one-step expression model with no impact to the dynamics of interest [5,14]. Following these prior works, we also focus on a reduced two-dimensional model.
In the system of Figure 1a, activator protein A promotes its own expression as well as the expression of repressor protein R. Protein R, in turn, represses expression of protein A. Let K m1 be the apparent dissociation constant between the activator protein and its DNA binding site and K m2 be the apparent dissociation constant between repressor protein and its DNA binding site [21] (see SI for details). For any species X, we denote in italics X its concentration. Consider the concentration of A and R given in units of their respective dissociation constants a :~A=K m1 and r :~R=K m2 . Considering a one-step model for protein expression, the dynamics for this system can be represented by in which d A and d R model protein decay (due to either dilution or degradation) and functions f 1 and f 2 model expression rates and take the form of the standard Hill functions [2] f 1 (a,r)~b 1 a m zb 2 1za m zr n and f 2 (a)~b 3 a m zb 4 1za m , in which b 1 and b 3 are the maximal expression rates, b 2 and b 4 represent the basal expression, and m and n are the Hill coefficients of the affinity between the proteins A and R and their respective binding sites. The mathematical derivation of this reduced nondimensional model is given in the SI. In the sequel, we refer to system (1) as the isolated system. We assume that the values of the parameters are such that system (1) has a unique equilibrium point. We give conditions for which this assumption holds when either m~1 or m~2 in the SI. In particular, it is shown that when m~1, the system always presents a unique and stable equilibrium and, therefore, no oscillatory behavior can be observed. When m~2 the uniqueness of the equilibrium is guaranteed under the following conditions: (i) the value of b 2 must be sufficiently smaller than the maximal expression rate of the activator, which is proportional to b 1 ; (ii) b 2 must be non-zero; (iii) the maximal expression rate of the repressor must be larger than the maximal expression rate of the activator; (iv) the smaller b 2 becomes, the smaller b 4 must be. In the general case (mw2), results related to existence and uniqueness of equilibria require a case by case analysis, which is out of the scope of this work. The results from this paper, do not explicitly impose conditions on the Hill coefficients m and n and only assume the uniqueness of the equilibrium (a Ã ,r Ã ) for system (1).
Since system (1) is a two-dimensional system, Poincaré-Bendixson theorem [22] can be employed to obtain conditions for the existence of a periodic orbit. Specifically, one must show that the trajectories of the system are bounded in a compact set and that the equilibrium point is unstable and not locally a saddle.
The following proposition shows that the trajectories of system (1) are bounded in a compact set. Proposition 1. There exists a constant D[R z such that the set K~f(a,r)[R 2 z Da 2 zr 2 ƒD 2 g is a positively invariant set under the vector field defined by system (1) and its equilibrium (a Ã ,r Ã )[K.
Proof. Note that f 1 (a,r) and f 2 (a) are positive bounded functions in the domain R 2 z . Let M 1~s up a,r ð Þ[R 2 z ff 2 (a,r)g and M 2~s up a[Rz ff 2 (a)g. First, notice that for a~0, _ a aw0 according to (1). Similarly, for r~0, _ r rw0. The quadrant R 2 z is, therefore, a positively invariant set. Define d Ã :~minfd A ,d R g and M :~maxfM 1 ,M 2 g. Consider the positive definite function v(a,r)~a 2 =2zr 2 =2. Using the chain rule, we obtain dv(a,r) dt~{ d A a 2 {d R r 2 zaf 1 (a,r)zrf 2 (a) ƒ{d Ã a 2 {dr 2 zaMzrM From the above, it is clear that _ v v(a,r)v0 on the exterior of a circle with center (M=2d Ã ,M=2d Ã ) and radius ffiffi ffi 2 p M=2d Ã . Therefore, for any Dw maxf ffiffi ffi 2 p M=d Ã ,a Ã ,r Ã g, _ v v(a,r)v0 along the arc defined by the boundary of K. Hence, K is a positively invariant set that contains the equilibrium (a Ã ,r Ã ).
To show that the equilibrium point is unstable and not locally a saddle, consider the Jacobian matrix of system (1) calculated at the equilibrium: and denote by (J 0 ) and det (J 0 ) the trace and the determinant of J 0 , respectively. The eigenvalues of the Jacobian are given by hence the equilibrium point is unstable and not locally a saddle if tr(J 0 )w0 and det (J 0 )w0. Given the specific expression of the Jacobian in (3), the equilibrium (a Ã ,r Ã ) of system (1) is unstable and not locally a saddle if the following conditions are fulfilled: System (1) satisfying conditions (i) and (ii) presents periodic orbits and will be referred to as Functional Clock.
Condition (ii) highlights a crucial design principle for the activator-repressor clock. In fact, assume that which is satisfied if the self activation is sufficiently strong. Then, This, in turn, implies that the timescale of the activator dynamics are sufficiently faster than that of the repressor dynamics. Hence, a central mechanism for the appearance of a limit cycle is a fast activator dynamics compared to the repressor dynamics. Retroactivity on a species due to downstream binding sites has been shown to slow down the species dynamics [16,17]. It follows that downstream binding sites can be employed to vary the relative speeds between the activator and the repressor dynamics.
Hence, we will also consider the non-oscillating version of system (1) that does not satisfy condition (ii), referred to as Non-Functional Clock. The non-functional clock is given by system (1) in which, in addition to condition (i), the following condition is satisfied: In this work, we study how the addition of binding sites to the repressor or activator can switch system (1) between the Functional Clock and the Non-Functional Clock behavior, with no change to the parameters of the original system (1).

Switching the Clock Off by Loading the Activator
In this section, we show the effect of additional DNA binding sites for the activator in a Functional Clock. Specifically, consider system (1) satisfying conditions (i) and (ii). The addition of DNA binding sites q A with affinity to the activator A, which binds as homomers, illustrated in Figure 1b, is modeled by the following chemical reaction in which D 1 represents the complex formed by A and q. In order to model the addition of DNA binding sites that are identical copies of the ones in the operator, we assume that the affinity between the DNA site q and the activator protein A is given by the apparent dissociation constant K m1~ffi ffiffiffiffiffiffiffiffiffiffiffiffiffiffi k b1 =k a1 m p , identical to the affinity of A to the promoters in the isolated clock. The impact in the dynamics from retroactivity can be obtained by employing binding sites with different affinities as long as the quantity of binding sites is adjusted accordingly [16]. Additionally, we assume the total concentration of binding sites q q A~( q A zD 1 )=K m1 to be constant. Let the complex concentration D 1 be given in units of K m1 using the nondimensional variable d 1~D1 =K m1 . The dynamics of the system after nondimensionalization are given by _ a a~{d A azf 1 (a,r)zmG 1 in which G 1~kb1 =d A models the timescale separation between the dissociation rate and the protein degradation. A mathematical derivation for this model is found in the SI. Since binding and unbinding reactions can occur in the order of milliseconds, they are in a timescale significantly faster than expression and degradation of proteins, which occur in the order of minutes [2]. As a result, parameter G 1 is very large. This fact allows to employ a singular perturbation argument [23,24] to facilitate the analysis of this system. To this end, define the small parameter E :~1=G 1 and re-write system (5) as In order to reduce this system to standard singular perturbation form, we perform the change of variables y~md 1 za, so that system (6) becomes _ y y~{d A (y{md 1 )zf 1 (y{md 1 ,r) ð7Þ which is in standard singular perturbation form. Setting E~0 one obtains from (9) the solution d 1~ q q A a m a m z1 :~w 1 (a). This equation defines the slow manifold, which can be shown to be locally exponentially stable (see SI). Hence, system (7) is well approximated by the reduced system obtained by replacing d 1 by its expression on the slow manifold w 1 (a). Specifically, we have that from which we obtain that _ a a~1 1zm dw 1 (a) da ({d A azf 1 (a,r)): 1zm q q A a m{1 (1za m ) {2 , the reduced system in the original coordinates is given by Since S A (a, q q A )=0, the equilibria of (10) are the same as the ones of (1). Therefore, if (1) has a unique equilibrium (a Ã ,r Ã ), this will also be a unique equilibrium of (10). Also, we have that 0vS A (a, q q A )ƒ1 and that S A (a, q q A ) is a strictly monotonically decreasing function of the amounts of DNA binding sites q q A . Hence, in system (10), the dynamics of the activator have been slowed down compared to the original isolated system (1). That is, the effective kinetic rate of the activator dynamics is now decreased by a factor equal to S A (a, q q A ). Note additionally that lim q q A ??
S A (a, q q A )~0 and S A (a,0)~1: The Jacobian of system (10) calculated at the equilibrium is given by in which we use the shorthand notation S Ã A ( q q A ) :~S A ( q q A ,a Ã ). We have det (J A ( q q A ))~S Ã A ( q q A ) det (J 0 )w0 from condition (i) and that Hence, while the addition of load does not change the sign of the determinant of the Jacobian, it can change the sign of the trace. For large enough load, because of (11), the trace becomes negative and the equilibrium point becomes stable. Hence, the periodic orbit disappears (see the SI for details). Figure 3 a shows the effect of load on system (5).
For the value of q q A for which (J A ( q q A ))~0, the eigenvalues of the Jacobian are imaginary, hence the system goes through a Hopf bifurcation. A continuation study shows that the Hopf bifurcation is present also in the full three-state system (5). In particular, the amounts of load needed to switch the clock off is about four times the amplitude of the activator oscillations. For the specific choice of parameters in this example, the amount of load required to stop this clock is of the same order of the dissociation constant K m1 , which usually amounts to a low concentration. For example, for the NRI activator used in the oscillator in [25], K m1 &10pM [6] which amounts to approximately 10 copies of the binding site per cell in E. coli.

Switching the Clock on by Loading the Repressor
We now consider a Non-Functional Clock and show how it can be turned into a Functional Clock by adding load to the repressor. Specifically, consider system (1) satisfying conditions (i) and (ii)'. Following the idea in the previous system, we model here the addition of DNA binding sites q R with affinity to the repressor R, identical to the binding sites found in the original clock. This interaction, illustrated in Figure 1c, is modeled by the following chemical reaction in which D 2 represents the complex formed by the R and q R . Let the affinity between the repressor and the binding sites is given by the apparent dissociation constant K m2~ffi ffiffiffiffiffiffiffiffiffiffiffiffiffiffi k b2 =k a2 n p . Let d 2 :~D 2 =K m2 be the nondimensional concentration of complexes and q q R~( q R zD 2 )=K m2 be the total nondimensional concentration of binding sites. The nondimensionalized dynamics of the system are given by in which G 2 :~k b2 =d R models timescale separation between the dissociation rate of the complex D 2 and the repressor decay rate. It is possible to reduce the order of system (14) by a similar technique employed in the previous section. To this end, define E :~G {1 2 . Define also the variable y :~rznd 2 , system (14) can be taken to the standard singular perturbation form _ a a~{d A azf 1 (a,y{nd 2 ) _ y y~{d R (y{nd 2 )zf 2 By setting E~0, one obtains the reduced system in the original coordinates, which, since the slow manifold is locally exponentially stable (see the SI), is a good approximation of system (14). This reduced system is given by _ a a~{d A azf 1 (a,r) _ r r~S R (r, q q R )({d R rzf 2 (a)) ð16Þ in which S R (r, q q R )~1 1zn q q R r n{1 (1zr n ) {2 : Since S R (r, q q R )=0, the equilibrium points of (16) are the same as the ones of the isolated system (1). Therefore the unique equilibrium point (a Ã ,r Ã ) of (1) is also the unique equilibrium point of (16). We employ the shorthand notation S Ã R ( q q R ) :~S R (r, q q R ). It is easy to verify that 0vS Ã R ( q q R )ƒ1 and that S Ã R ( q q R ) is a strictly monotonically decreasing function of q q R . Furthermore, we have that lim q q R ??
Hence, the addition of the load to the repressor makes the dynamics of the repressor slower compared to that of the isolated system (1). That is, the repressor effective kinetic rates are now smaller by a factor equal to S Ã R ( q q R ), which can be arbitrarily decreased by increasing the amounts of sites q q R . The Jacobian of system (16) calculated at the equilibrium (a Ã ,r Ã ) is given by Thus, the addition of load to the repressor does not change the sign of the determinant of the Jacobian as det (J R ( q q R ))~S Ã R ( q q R ) det (J 0 )w0. However, it can change the sign of the trace from negative to positive as condition (ii)' is satisfied and condition (17) holds. Hence, the equilibrium point can become unstable with sufficient addition of the load and the system begins oscillating (see the formal derivations in the SI). Figure 4a shows the effect of load on system (14). Note that the parameters were chosen so that the system satisfies conditions (i) and (ii)'. When tr(J R ( q q R ))~0, a Hopf bifurcation occurs since both eigenvalues are imaginary. A continuation analysis can be used to show that this Hopf bifurcation is also present in the full system (14). Figure 4b illustrates that the amount of load required for the Hopf bifurcation is given by q q R~1 :32 in units of K m2 . Hence, the amounts of load needed to switch the clock on is on the same order of the amounts of repressor at the equilibrium. For the LacI repressor employed in [6], K m2 &1pM [26], which amounts to few copies per cell of the load. Figure 4c shows that the addition of load increases the period of oscillation. This suggests the possibility that the load can be employed not only for switching an oscillator ''on'' and ''off'' but for also tuning the period. However, the increase in period is accompanied by an increase in the amplitude of the oscillation (Figure 4b), which may be undesired. We discuss how the period can be changed while maintaining the amplitude through simultaneous addition of activator and repressor loads in Section ''Tuning the Clock period''.

Stochastic Simulations of the Switching Behavior
In order to understand how robust the switching behavior is to intrinsic noise, we employ stochastic simulations of the system. An implementation of the Gillespie algorithm [27] was employed to produce realizations of trajectories of an activator repressor clock in which both activator and repressor bind to DNA as dimers (m~n~2).
In these simulations, we assumed the presence of 5 copies of each activator and repressor gene to emulate the situation in which the circuit is present in a low copy number plasmid. Expression rates and degradation rates were chosen based on the values used in the deterministic models to obtain a functional and a nonfunctional oscillator. The association and dissociation rates between proteins and dimers were chosen so that the apparent dissociation constants K m1~Km2~1 , which consider a bacterial transcription factor with apparent dissociation constant on the order of picomolars. A detailed description of this model is given in the SI. Figure 5a shows that addition of binding sites with affinity to the activator can eliminate oscillations from a functional clock. Figure 5b shows how the addition of binding sites with affinity to the repressor can generate sustained more robust oscillations in a non-functional clock. In both situations, the amount of loads employed to switch the clock is on the order of 10 2 {10 3 copies of binding sites per cell, which can be achieved by inserting small arrays in high copy number plasmids.

Tuning the Clock Period
As noticed in Figure 4c, addition of binding sites to the repressor increases the period of the limit cycles of the system. However, this may cause an increase in the amplitude of the cycle (Figure 4b), which may be undesirable. In this section, we illustrate how the simultaneous addition of load to both the activator and repressor can be employed to vary the period as desired with little impact on the cycle amplitude.
Consider the nondimensional model for the system with DNA binding sites for both the activator and the repressor as shown in Figure 1d: Here, k b1 ,k u1 model the association and dissociation rates between the activator protein and its corresponding DNA binding site q A , k b2 ,k u2 model the association and dissociation rates between the repressor protein and its corresponding DNA binding site q R , g 1 (A,R), g 2 (A) represent the dimensional version of the Hill functions (see SI), and q A,T , q R,T represent the total concentration of activator and repressor DNA sites.
This system can be nondimensionalized, by setting the nondimensional states a~A=K m1 , r~R=K m2 , d 1~D1 =K m1 and k 2~D2 =K m2 , as shown in the SI, to obtain system _ a a~{d A azf 1 (A,R)zmG 1 (14) with two different amounts of load. The parameters in the simulation were b 1~b3~1 00, b 2~: 04, b 4~: 004, d A~1 , d R~1 :5, G 2~1 00, m~2 and n~4. The amount of DNA binding sites in the system with no load is q q R~0 whereas in the system with repressor load is q q R~2 0. (b) Hopf Bifurcation with q q R as a parameter. A continuation of the equilibrium as a function of the load parameter q q R shows that, for this set of parameters, the amount of load required to activate the clock is in the same order of magnitude as that of the the affinity coefficient K m2 , with bifurcation occurring at q q R~1 :32. This plot was obtained via continuation of system (14) with the same parameters as before. Solid lines indicate a stable trajectory (limit cycle to the right of the Hopf bifurcation and the equilibrium to its right). The dotted line indicates an unstable equilibrium point. (c) Period increases as a function of the repressor load q q R . doi:10.1371/journal.pone.0041019.g004 Figure 5. Effect of the load on clock holds under intrinsic noise. The plots above are stochastic realizations of an activator-repressor clock with m~n~2 and containing 5 copies of activator and repressor genes. (a) Functional clock stops with load to the activator. We show that, with the chosen parameters, it is possible to stop the clock with an amount of load that is roughly 100 times higher than the copy number of the circuit. (b) Load to the repressor leads to robust oscillation. We show that, the it is possible to generate robust oscillation with roughly 400 times the number of circuit genes with the choice of parameters above. doi:10.1371/journal.pone.0041019.g005 as given in expressions (2), q q A~qA,T =K m1 and q q R~qR,T =K m2 , and G 1 and G 2 are as defined before. In order to employ a singular perturbation argument similar to what was done in the previous sections, define E~1=G 1 , n~G 2 =G 1 to model the explicit timescale separation present in this system. Define also the following change of variables y 1 :~azmd 1 and y 2~r znd 2 . Substituting these in (20), one obtains the system in standard singular perturbation form: By setting E~0, one obtains the slow manifold , q q R r n r n z1 :~w 1 (a),w 2 (r) ð Þ : Since the slow manifold is locally exponentially stable (see SI), the reduced system is a good approximation of system (21). Since _ y y 1~_ a azm dw 1 (a) da _ a a and _ y y 2~_ r rzn dw 2 (r) dr _ r r, this reduced system, in the original variables, takes the form in which 1z q q A m 2 a m{1 (1za m ) {2 and S R (r, q q R ) :~1 1zn dw 2 (r) dr~1 1z q q R n 2 r n{1 (1za n ) {2 : Let the activator and repressor loads be added at a fixed ratio r~ q q A = q q R and define F (a,r, q q R ) :~S R (r, q q A =r) S A (a, q q A ) . System (22) can be re-written as Since S A (a, q q A )w0, this system is orbitally equivalent [19] to the system _ a a~{d A azf 1 (a,r) ð Þ _ r r~{d R rzf 2 (a) ð Þ F (a,r, q q A ): ð24Þ Hence, if system (23) has a periodic orbit, system (24) will have a corresponding periodic orbit with identical trajectories. The corresponding periodic signals, however, will have different periods whose values depend on function S A (a, q q A ). Thus, if the value of F(a,r, q q A ) does not appreciably change when q q A changes, the addition of the load will affect the period of oscillations without impacting their amplitudes. Since we have that for large values of q q A , LF (a,r, q q A ) L q q A &0. Under these conditions, since the function S A (a, q q A ) is a monotonically decreasing function of q q A , the periodic orbits of system (23) will display decreasing periods as q q A increases, while maintaining the same amplitude, due to orbital equivalence between system (24) and system (23) (see the SI for a formal proof). Figure 6a illustrates this result. The addition of repressor load to a functioning clock increases the period but also leads to a higher amplitude. This effect in the amplitude is not observed when both activator and repressor loads are added. Figure 6b shows this behavior for increasing amount of load. When only repressor load is added, there is an increase in the period of the limit cycles along with an increase in the amplitude, as it was seen in the previous section (Figure 4(b) and (c)). However, if a sufficient amount of activator load is simultaneously added along with the repressor load, the increase of the period occurs with very little impact on the amplitude of oscillations.

Discussion
Effective kinetic rates are crucial parameters for the dynamic behavior of biomolecular networks. In particular, delays in negative feedback loops have been shown to be a fundamental mechanism for periodic oscillations both in electronic circuits [28] and in biomolecular networks [1,10,15]. Research has shown that in natural systems these delays are realized by the number of steps, such as transcription, translation, and post-translational modifications, involved in the implementation of the feedback loop. More steps lead to larger delays. Hence, adding a delay involves engineering the structure and length of a pathway. In this paper, we have revealed that a different mechanism exists for adding and carefully tuning delays and effective kinetic rates: the addition of DNA targets. In natural systems, transcription factors can have large numbers of DNA binding sites, several of which do not even have regulatory functions (see [29] and [30], for example). Our study suggests that a role of these DNA binding sites is to carefully tune effective kinetic rates to realize the desired dynamics in genetic networks.
As an example, consider the regulation network of cellular resources such as ribosomes or RNA polymerase (RNAP). Since both molecules need themselves to be assembled, there is a self activating loop. Additionally, it has been shown [31,32] that RNAP and ribosomes are negatively regulated through transcriptional repression. Hence, the regulation motif of these species has the form of Figure 1a, in which we can view A as the resource (RNAP or ribosome) and R as a repressor system. This motif, as we have shown, can present sustained oscillations, which would be undesired for RNAP or ribosomes. However, due to the large demand by the cellular environment, RNAP and ribosomes are  (20) increases when we add exclusively DNA binding sites with affinity to the repressor ( q q A~0 , q q R~1 0). However, if we simultaneously add DNA binding sites with affinity to the activator, the amplitude is not affected as much ( q q A~ q q R~1 0). (b) The period of system (20) can be changed with no effect on the amplitude when DNA binding sites with affinity to both the repressor and the activator are added simultaneously. The upper plot shows that a similar increase of period observed via the addition of repressor load can be obtained via the simultaneous addition of activator and repressor load. This second method has the advantage of not generating an increase in the amplitude, as shown in the lower plot. In this simulation we assumed the ratio r~ q q A = q q R~1 . Parameters of the activator repressor system used in the simulation were b 1~b3~1 00, b 2~: 04, b 3~: 004, d A~1 , d R~0 :5, G 1~G2~1 00 and m~2, n~4. In the traces showing only repressor load r~0, while the traces showing simultaneous repressor and activator load, r~1. doi:10.1371/journal.pone.0041019.g006