Theory on the Dynamics of Oscillatory Loops in the Transcription Factor Networks

We develop a detailed theoretical framework for various types of transcription factor gene oscillators. We further demonstrate that one can build genetic-oscillators which are tunable and robust against perturbations in the critical control parameters by coupling two or more independent Goodwin-Griffith oscillators through either -OR- or -AND- type logic. Most of the coupled oscillators constructed in the literature so far seem to be of -OR- type. When there are transient perturbations in one of the -OR- type coupled-oscillators, then the overall period of the system remains constant (period-buffering) whereas in case of -AND- type coupling the overall period of the system moves towards the perturbed oscillator. Though there is a period-buffering, the amplitudes of oscillators coupled through -OR- type logic are more sensitive to perturbations in the parameters associated with the promoter state dynamics than -AND- type. Further analysis shows that the period of -AND- type coupled dual-feedback oscillators can be tuned without conceding on the amplitudes. Using these results we derive the basic design principles governing the robust and tunable synthetic gene oscillators without compromising on their amplitudes.


Introduction
Transcription factors (TFs) regulate the quantitative levels of several proteins inside a living cell [1][2][3][4]. TF networks present across various organisms ranging from prokaryotes to higher eukaryotes and consist of fundamental building blocks such as autoregulatory loops, cascades and single input modules, feedforward and feedback loops, dense overlapping regulons and oscillatory loops [5][6][7]. Feedback loops act as bistable switches and feedforward loops have been shown to act as efficient filters for transient external signals [8], [10][11][12]. Positive self-regulatory loops seem to play important roles in the maintenance of cellular memory [3] and subsequent reprogramming of the cellular states whereas negative auto regulatory loops have been shown [11] to speed up the response times against an external stimulus [8][9][10], [12]. Oscillatory loops drive the developmental as well as mitotic cell-cycle dynamics [13] and circadian-rhythms [14], [15] associated with the intracellular concentration of various types of proteins, metabolites and other cell-signaling molecules. Understanding of the detailed dynamics of oscillatory loops associated with the TF networks is a central topic in biophysics, synthetic and systems biology.
The minimalist TF network model that can generate selfsustained oscillations is the well-known Goodwin-Griffith oscillator which has a single gene that codes for a TF protein that negatively auto-regulates its own transcription [16][17][18]. In this model the TF protein-product undergoes a one-step modification that yields the matured or active end-product and subsequently n numbers of this end-product bind with the cis-regulatory modules (CRMs) of the associated promoter that in turn results in down-regulation. Here n is the Hill coefficient associated with the cooperative type binding. Detailed studies on this minimalist model showed [17] that the inequality condition n.8 is necessary to generate selfsustained oscillations in the levels of mRNA and protein. This result was obtained with the assumptions that the rate constants associated with the synthesis and decay of the protein and endproduct are equal and the binding-unbinding of the end-product with the promoter is much faster than the rate of change in the synthesis and degradation of mRNA, protein and end-product. Further it was assumed that the decay of mRNA and protein product follows a first order type reaction.
It was realized later that the inequality condition n.8 is unlikely [19] under in vivo conditions since the formation of such large multimeric protein complexes via pure three dimensional diffusion (3D) limited collisions ( Figure 1) is almost an improbable event and several other modifications over the Goodwin-Griffith model were proposed to reduce the required value of n. Most of these modifications were mainly associated with the insertion of (a) a temporal delay in the negative auto-regulatory loop either explicitly as a time-delay in between the synthesis and binding of end-product at the promoter [19], [20] or implicitly via inserting additional reaction steps [19] in the formation of end-product that interacts with its own promoter and (b) a non-linear Michaelis-Menten type kinetics in the decay of mRNA and protein products despite of the first order type kinetics and (c) additional TF gene members in the negative feedback loop which again indirectly acts as temporal delay in the overall negative feedback. The delayed negative feedback may also be coherently or incoherently amplified [21][22][23] via the insertion of a positively regulated intermediate. Here the temporal delay is connected with the overall time that is required for the transport of fully transcribed mRNAs from nucleus to cytoplasm, post-translational modifications and subsequent transport of active TF proteins into the nucleus through 3D diffusion. When the decay of mRNA and protein product follows a Michaelis-Menten type kinetics then the Goodwin-Griffith (GG) oscillator seems to produce self-sustained oscillations [19] even for n = 1. The genetic oscillator module with three TF genes connected in a cyclic negative feedback loop is named as repressilator [24]. Though these motifs were shown to be oscillatory through deterministic and stochastic simulations, significant fraction of cells containing the constructs of these motifs were not showing any oscillations under in vivo experimental conditions. It was argued that it could be partially due to the noisy nature of intracellular environment [18], [24]. Here one should note that most of the simulation studies were performed with constant parameter values which may not be true under in vivo conditions. In this context it is essential to investigate how the oscillatory dynamics of these motifs reacts to perturbations in the system parameter values.
Most of the earlier studies on GG and other oscillator models assumed a quasi-equilibrium condition for the binding-unbinding dynamics of the negatively autoregulated TF proteins at their own promoters. This is mainly to reduce the four or higher dimensional Jacobian matrix associated with the non-linear system of differential rate equations into a three dimensional one to ease further analysis since there is an additional rate equation corresponding to the promoter state dynamics apart from the rate equations associated with mRNA, protein and end-product. However this assumption is valid [8], [9] only when the timescales associated with the synthesis and degradation of mRNAs and TF proteins are much slower than the timescales associated with the binding-unbinding of regulatory TFs at the respective promoters. Recent studies [8] on feedforward loops suggested that the binding-unbinding dynamics of TF protein at the promoter can be ignored only when the cellular volume V c ( = volume of nucleus in Figure 1. Goodwin-Griffith genetic oscillator model. The transcription factor (TF) gene A is transcribed and translated into the TF protein product that in turn is converted to the active end-product. The end-product (or its oligomer as in case of lacI repressor negative-feedback-only system that was constructed in reference [25]) is the key molecule that locates the respective cis-regulatory elements associated with the promoter of TF gene A through a combination of one-dimensional (1D) and three-dimensional (3D) routes as that of typical site-specific DNA-protein interactions.
Here either monomers of the end-product directly assemble at the corresponding regulatory elements in a combinatorial manner (I) or the fullyformed complex of n a -mer binds with the respective regulatory sites (II). Assembly of combinatorial TF molecules results in the looping of DNA segment that is present in between the promoter and cis-regulatory elements. ARPC is the assembled repressor-promoter complex that in turn results in down-regulation. Our analysis shows that out of these two competing pathways, the pathway I is the most probable one since it takes shortest time. The corresponding set of differential equations is given in Eqs (4)(5). This system is well characterized by parameters of Group I, II and III. Group I consists of parameters w a ,v a ,e a ð Þwhereas Group II consists of equilibrium parameters l a ,m a ð Þ and Group III consists of ordinary type perturbation parameters s a ,k a ,x a ð Þ . Most of the earlier studies assumed zero values to Group II parameters apart from assuming zero for v a that controls the promoter state dynamics. Blue colored spheres are the dimers of lac repressor. Here a.a denotes amino acids and n.a denotes nucleotides. doi:10.1371/journal.pone.0104328.g001 case of eukaryotes) is comparable with that of the prokaryotes [8] such as E. coli (V c ,10 218 m 3 ) and the influence of the promoter state fluctuations on the overall dynamics of feedforward/feedback loops seems to significantly increase as the nuclear volumes increases as in eukaryotic cells across yeast to human. Further, the Michaelis-Menten type degradation kinetics associated with mRNA and protein is a valid assumption only when the concentrations of these species are much higher than the concentration of the corresponding nucleases and proteases. Nevertheless in most of the in vivo conditions, the intracellular levels of mRNA and protein of a particular TF gene will be much lesser than the corresponding levels of the non-specific nucleases and proteases. When the latter is true then the enzyme mediated decay of mRNA and protein will eventually follow a first order type kinetics. In this article, using a combination of theoretical and simulation tools (a) we develop a generalized theoretical framework of various types of genetic oscillators by explicitly incorporating the promoter state dynamics and other chemical reaction balances in detail. Using our detailed model (b) we identify and classify various critical control parameters and compute their physiological ranges which are required to generate self-sustained oscillations in the intracellular levels of mRNAs and transcription factor proteins and (c) explore various possibilities of coupling independent gene oscillators and fine-tuning the period of such coupled system. We further (d) demonstrate that by coupling two or more independent Goodwin-Griffith oscillators one can design oscillatory network architectures which are tunable and also robust against perturbations in system parameters.

Theoretical framework of transcription factor gene oscillators
The Goodwin-Griffith oscillator consists of a negatively selfregulated gene (we denote it as TF gene A) which codes for a transcription factor protein ( Figure 2A). We denote the cellular concentrations (mol/lit, M) of its mRNA as m a , protein as p a , the transformed end-product as z a and the complex of promoter with the end-product as x a . Here the total cellular concentration of promoter is d za and the overall promoter state occupancy by the end-product will be X a = x a /d za where X a [ (0, 1). Though there is only one copy of the promoter by definition, we use a continuous type probability variable X a to describe promoter state occupancy mainly to account for its partially occupied status [8], [9]. The transcription and translation rates are denoted as k ma (Ms 21 ) and k pa (s 21 ) respectively. The first order decay rate constants corresponding to mRNA and TF protein are c ma (s 21 ) and c pa (s 21 ) respectively. The first order on-and off-rates associated with the transformation of protein into the matured end-product are denoted as l af (s 21 ) and l ar (s 21 ) and the corresponding dimensionless dissociation constant is l a~lar l af . The overall forward and reverse rate constants associated with the binding and unbinding of n a numbers of end-product molecules with the respective cis-regulatory modules (CRMs) of the promoter of TF gene A are k af (M {na s {1 ) and k ar (s 21 ) and the corresponding dissociation constant is defined as K arf~kar k af (M na ). To simplify the analysis further we introduce the following scaling transformations to project the time and concentration variables into the dimensionless space. t~c pa t; P a~pa =p as ; M a~ma =m as ; Z a~za =p as ; X a~xa =d az ð1Þ In these equations t denotes the dimensionless time that is measured as the number of lifetimes of the protein product of TF gene A and P a , M a , Z a and X a are respectively the dimensionless concentrations of protein, mRNA, end-product and promoter complexes. We also should note that (P a , M a , Z a and X a ) [ (0, 1) by definition. Here the steady state values of mRNA and protein in the absence of negative self-regulation can be defined as follows [8], [9].
We further transform the parameter associated with the multimerization of end-product and subsequent binding events as follows. m a~Karf p na as ; K arf~kar k af ð3Þ Using the scaling transformations given by Eqs (1-3) one can write the deterministic rate equations corresponding to the temporal evolution of dimensionless concentration variables (X a , M a , P a , and Z a ) over dimensionless time variable t as follows.
The initial conditions are X a ,Z a ,M a ,P a ð Þ0 at t~0. When (v a = 0, s a = 0 and k a = 0), then this system reduces to the usual GG oscillator model for three concentration variables. Here we have defined the dimensionless ordinary perturbation parameter k a~cza l af where c za (s 21 ) is the first order decay rate constant associated with the protein end-product Z a . Since c za &l af will be true in most of the physiological conditions and k a is an ordinary perturbation parameter one can assume k a &0. When there is a dimerization of z a (we denote the dimer z a -z a as y a ) as in case of Lac repressor system that has been constructed and studied in reference [25] (negative-feedback-only model using lacI gene) then the first and last equations of Eqs (4) will be modified as follows.
Here various parameters associated with the dimerization of end-product of TF protein A and subsequent assembly of this endproduct at the own promoter are defined as follows.
x ya~p na{1 as d az k af l fya ; e ya~cpa l fya p as ; l ya~lrya p as l fya ; s ya~pas l fya l fa Here we have defined Y a~ya =p as and l fya (M 21 s 21 ) is the forward rate constant associated with the dimerization reaction and l rya (s 21 ) is the corresponding reverse rate constant. One should note that for a fully functional Lac repressor system the Hill coefficient will be n a = 4 since an octamer of Lac repressor protein (which is a dimer) is involved in the overall looping of DNA that results in strong repression of lacI. The system of Eqs (4)  Here one should note that the parameters m a ,x a ,v a ð Þ are functions of n a that can be simplified by assuming an in vivo protein level as p as~1 . To simplify the analysis further we can classify these dimensionless control parameters into Group I, II and III. Group I consists of v a ,w a ,e a ð Þwhich are all singular type  perturbation parameters since they multiply the first order derivative terms. One should note that this set of parameters directly controls the dynamics of changes in the cellular concentrations of active promoter, mRNA and end-product respectively. Group II consists of (s a , x a , k a ) those are ordinary type perturbation parameters. In Group II, s a controls the coupled dynamics associated with the concentrations of TF protein A and its end-product whereas x a controls the coupled dynamics of changes in the concentrations of end-product and it's binding with the promoter sequence. The lifetime of end-product is controlled by k a . One should note that almost all the earlier studies assumed that (s a , x a , k a ) = 0. Group III consists of the equilibrium and promoter affinity parameters (l a , m a ). In this l a controls the equilibrium associated with the formation of end-product and m a controls the equilibrium associated with the binding of n a molecules of end-product with CRMs of the promoter of TF gene A.

Biophysical modeling of promoter state dynamics
The total time required to initiate transcription consists of at least two different components viz. the time required (proportional to 1/k af ) for the assembly of n a numbers of TFs at the respective CRMs of promoter and the time required for subsequent looping of DNA and subsequent distal action of TFs on RNAPII-promoter complex. The time component associated with the looping and distal action along with the time required for elongation and termination steps are included in the definition of transcription rate (the total time required for transcribing a full length mRNA will be equal to 1/k ma ).The kinetics of interaction of n a molecules of end-product with the sequentially located CRMs can occur in two different possible ways namely binding of the full-fledged complex of n a molecules of end-product (pathway II) or sequential assembly of the monomers of end-product at the corresponding cis-regulatory DNA-binding sites similar to that of a combinatorial binding of TFs with CRMs (pathway I) as in eukaryotic systems ( Figure 1). Though the pathway I resembles a (n a +1) th order reaction it is still an operable one since the length scale of the genomic DNA is much higher than the combinatorial binding TF proteins. Binding of n a numbers of transcription factors in a sequential manner or n a -mer of end-product leads to looping of the DNA segment that is present in between promoter and CRMs of TF gene A that results in the spatial or distal communication between the end-product present at CRMs and the already formed RNAPII-promoter complex which in turn activates (positive) or deactivates (negative) the initiation of transcription depending on the type of self-regulation [1], [3], [4], [26], [27]. In case of activation or positive regulation, the combinatorial transcription factors bound at CRMs enhance the initiation of transcription by strengthening the RNAPII-promoter interactions through their distal action (positive arrows in Figure 2) whereas in case of negative regulation, the RNAPII-promoter complex will be destabilized by the combinatorial TFs present at CRMs (negative arrows in Figure 2). Here the destabilization of RNAPIIpromoter complex may be through the formation of stem and loop structures. In prokaryotes, these types of up and down regulations generally do not involve recruitment or combinatorial binding of several TFs and the regulator transcription factor directly influences the RNAP-promoter interactions as in case of negatively self-regulated oscillatory motifs constructed with a lacrepressor gene. Here binding of lac-repressor at the Operator sequence directly destabilizes RNAP-promoter complex that in turn lead to the down regulation of transcription [1], [3]. The total time t d,na required for the formation of a full-fledged n a -mer via 3D diffusion-controlled collisions and subsequent binding with the cis-regulatory sites can be calculated as follows.
In this equation n a t s,1 is the time required for the searching and binding of the entire n a -mer at the corresponding CRMs on DNA (for a monomer it will be t s,1 ) via a combination of 1D and 3D diffusion, t d,1 is the time that is required for the formation of a dimer of the end-product through 3D diffusion under in vivo conditions, Y n a ð Þ~d ln C n a ð Þ=dn a where C n a ð Þ is Gamma function and c e~0 :5772157: is the Euler-Mascheroni constant.
Here t d,1~1 0 {3 16pRD T N A C Z (s) is the minimum possible 3D diffusion controlled bimolecular collision time inside the cellular volume where R is average radius of the monomers of endproduct, C Z (mol/lit) is the concentration of end-product inside the cellular volume, N A is the Avogadro's number, D T~kB T=6pQR ( m 2 s 21 ) is the 3D diffusion coefficient associated with the dynamics of monomers of end-product in aqueous medium where k B is the Boltzmann constant, T is the absolute temperature (K), Q is the viscosity of the medium. In the calculation of t d,1 , we have assumed that the reaction radius between two monomers is ,2R. Since the overall maximum radius of the m-mer will be ,mR, we find that t d,m !m= 1zm ð Þ and subsequently the total time that is required to form a n a -mer in a sequential manner via 3D diffusion will be given by the sum t d, 1 P na m~1 m= 1zm ð Þ. We should note that this is the maximum possible search time since we have assumed a maximum possible radius for the n a -mer complex and also we have not considered the possibility of formation of the final n a -mer through non-sequential and random pathways and the steric factor associated with the multimerization reaction. The total search-time that is required by the monomer of end-product to find its cognate site on DNA is defined as t s,1~N t L zt ns ð Þ =L where the overall 1D sliding time is defined as t L~L 2 6x d [28]. In this calculation we have assumed that the end-product searches for its binding sites on DNA via a combination of 1D and 3D diffusion controlled collision routes. Monomers of the end-product undergo at least N/L numbers of cycles of 3D diffusion mediated association that is followed by 1D scanning and dissociation where N is the size of the genomic DNA (base-pairs, bps) and L (bps) is the average 1D sliding-length between non-specific 3D association and dissociation. Here t L is the time that is required by the monomers of end-product to scan L bps of the genomic DNA via 1D diffusion along the DNA chain, x d (bps 2 s 21 ) is the 1D diffusion coefficient associated with the sliding of monomers on DNA (this will be scaled down to x d =n a for n a -mer) and t ns is the time that is required for non-specific binding of end-products with the genomic DNA via 3D collisions under in vivo conditions. When all the n a monomers of end-product search for their binding sites on DNA in a parallel manner, then the total time [26] that is required (t s,na ) for all these n a monomers to assemble at the sequentially located cis-acting elements can be derived from the theory of combinatorial binding of transcription factors [27][28][29] with DNA as follows.
From this equation we find that the 1D scanning time increases with the number of monomers n a in a power law manner as n a a where typical value of the exponent seems to be a2 =5 2=5 [26]. From Eqs (6) and (7) one can compute the following ratio.
From the theory of site-specific DNA-protein interactions we find that t d,1 =t s,1 ð Þ &10 2 for n a = 1 [26][27][28][29] which suggests that h na §10 2 for all values of n a . Eqs (7) and (8) suggest the pathway I is more efficient than pathway II. This means that though the diffusion limited multimerization of the monomers of the endproduct is not a reasonable assumption for large values of n a , direct assembly of the monomers of the end-products of TFs on the sequentially located cis-regulatory DNA binding sites via a combination of 1D and 3D diffusion controlled routes can be still a reasonable assumption even for higher values of n a . In this context we can replace the combinatorial binding of n a numbers of TFs with CRMs of template DNA with a single step (n a +1) th order reaction as given by the first equation in Eqs (4). One should note that unlike the prokaryotic systems, most of the eukaryotic promoters are activated through a combinatorial binding of several TFs at the corresponding CRMs [27]. This observation suggests that GG oscillator can be still a feasible model that can be used to generate limit-cycle oscillations in the cellular levels of negatively self-regulated TF proteins especially in eukaryotic systems.

Steady-state analysis of Goodwin oscillator
The system of Eqs (4) has a fixed point P a~ga which is a real solution of the following polynomial equation of the order (n a +1).
The steady state values of other concentration variables (Z a , M a and X a ) can be calculated using the fixed point g a as follows.
Using the Jacobian matrix evaluated around the equilibrium point P a~ga , the linearized form of the system of Eqs (4) near this equilibrium point can be written as follows.
Here we have defined various matrix elements as follows.
g~m a =b a g a v a ; c~l a zk a zx a A a ð Þ =e a ; d~x a m a =b a g a e a ; A a~na m a l a zk a ð Þ1{b a g a ð Þ =g a The coefficients associated with the characteristic polynomial Y 4 zrY 3 zsY 2 ztY zu~0 (PI) of the Jacobian matrix defined in Eqs (10) can be written as follows.
From the Routh-Hurwitz criterion for a biquadratic polynomial [30] we find that the equilibrium point of the system P a~ga will be stable only when all the following inequality conditions are true. In other words there may be limit-cycle oscillations around the steady state only when any of these inequality conditions is not true.
Here the first inequality condition rw0 will be always true since g a !1=b a and subsequently we find A a w0. The second inequality condition rs{t ð Þw0 can be reversed only when either s,0 or rsvt for s.0. When the third inequality in Eqs (11) is not true then the biquadratic polynomial can have a complex root such that y Re +iy Im with positive real-part (here y Re §0 y Im §0) that results in the generation of limit-cycle oscillations of the concentration of TF protein A and the period of such oscillations [29][30][31][32] will be given by t p~2 p=y Im . This means that the period of GG oscillator can be modified by tuning the lifetime (c pa ) of the protein product of TF gene A (since t~c pa t) though the value of y Im is a function of other Group I parameters (w a , v a , e a ) which are in turn linearly depend on c pa . Since the Hill coefficient term n a presents only in the coefficient terms s, t and u, the third inequality in Eqs (11) can be reversed by increasing the value of n a for any set of Group I, II and III parameters. This means that the parameter space that is required to generate oscillations can be expanded by increasing the value of n a . Inequality conditions given by Eqs (11) for a stable motion of the dynamical system of Eqs (4) can be directly derived from the following Routh table (RT GG ) [30] corresponding to the biquadratic polynomial (PI).
When y Im~0 then the steady state solution will be either asymptotically stable or unstable depending on the values of realparts. From Eqs (10-11) we find that the system will be inconsistent near the fixed point both at very large as well as small values of Group I type parameters as v a ,w a ,e a ð Þ ?0 or ?. This means that there exists a critical range of these parameters to generate limit-cycle oscillations in the cellular levels of TF protein A. One should note that Group I parameters appear in the denominator of definitions of various coefficients of the characteristic polynomial (PI) which means that the period of oscillations will increase proportionately with respect to an increase in these parameters since we find y Im !1= v a ,w a ,e a ð Þ . Contrasting from Group I, there exist critical or threshold values of Group II type ordinary perturbation parameters x a ,s a ,k a ð Þ below which the oscillations occur. As in case of Group I type parameters, there exist a critical range of values of l a ,m a ð Þ Group III to generate limit-cycle oscillations. The critical value of Hill coefficient n a~C n a for a given set of parameters can be iteratively calculated by numerically solving the third inequality in Eqs (11) at various values of n a . When the dynamics of promoter state occupancy is much faster than the rate of change in the concentrations of other variables then one can set v a &0 and the system of Eqs (4) reduces to the following form.
Most of the earlier studies on GG oscillator consider Eqs (13) as the base model however with the conditions such that s a ,k a ð Þ~0. Using detailed numerical simulations we will show later that this assumption is reasonably invalid. The corresponding Jacobian matrix around the steady state P a~ga can be written as follows.
Here we have defined A 0 a~n a b a 1{b a g a ð Þl a zk a ð Þ and the characteristic polynomial associated with this equation is where the coefficients are defined as follows.
The Routh-Hurwitz [30] condition required by the system of Eqs (13)(14) to generate limit-cycle oscillations will be r 0 s 0 {t 0 ð Þ v0. Upon solving this inequality for the Hill coefficient n a , the expression for the critical value of n a that is required to generate limit cycle oscillations can be obtained as follows.
Here one should note that the term g a in the right hand side of this equation is still a function of n a and m a and the following limiting conditions exist.
Eqs (15)(16)(17) suggest that strong binding of the end-product (m a [ 0) at the promoter of TF gene A is required along with the conditions such as k a , s a ð Þ~0, and l a , w a , e a ð Þ1 to decrease the required critical Hill coefficient towards the minimum possible value as C n a & 9. When there is an additional dimerization step as in Eqs (4) and (5) then the resulting characteristic polynomial of the Jacobian matrix will be of fifth order as Y 5 zrY 4 zsY 3 ztY 2 zuY zm~0 (PIII) and the Routh criterion that is required to generate oscillations can be written as follows.
The physiological ranges of various parameters associated with Goodwin-Griffith oscillators are listed in Table 1.

One-to-one dual feedback oscillators
One can extend these scaling ideas for one-to-one negative feedback oscillator or toggle switches ( Figure 2B1) and repressilator models ( Figure 2C1). In case of one-to-one negative feedback oscillator n a number of end-product molecules of TF gene A bind with the cis-regulatory elements associated with the promoter of TF gene B and subsequently down-regulates whereas n b number of end-product molecules of TF gene B down-regulate the promoter of TF protein A upon binding with the corresponding cis-regulatory elements ( Figure 2B1). The set of differential rate equations associated with the two TFs one-to-one feedback system can be written in the dimensionless form as follows.
In these equations the subscripts will be such that when h~a,b then q~b,a where (a, b) denote the TF genes A and B respectively. One should note that the Hill coefficients associated with the binding of the end-product of TF A at the promoter of TF B and end-product of TF gene B at the promoter of TF gene A are n a and n b respectively and in general n a =n b . Here we have defined various other dimensionless variables and parameters as follows.
t~c pa t; P h~ph =p hs ; M h~mh =m hs ; Z h~zh =p hs ; In these definitions for h = a, b one needs to set q = b, a. The steady state solutions g b ,g a ð Þ to the coupled set of Eqs (18) corresponding to this two TFs system with respect to the scaled protein levels can be given as follows.
The steady state values of other dynamical variables can be calculated using the steady state values of the protein products g a and g b as follows.
Here y in Eqs (19) is the real root of y n a n 1{b a e y ð Þ À Á 0 where we have defined the function W as, Eqs (18) have three possible steady state solutions viz. g b~ga ð Þ , g b vg a ð Þ and g b wg a ð Þ . Under identical values of all the parameters such as (m a~mb , n a~nb and so on) we find the unstable steady state solution of the two TFs system as g b~ga where 0ƒ g b~ga ð Þ ƒ1. This means that the limit-cycle oscillations around this unstable steady state can occur only when the values of all the control parameters and initial conditions are identical with respect to both the TF genes A and B. Using the eighth order characteristic polynomial of the Jacobian matrix associated with the linearized form of Eqs (18) (Methods section) near the steady state values g b~ga ð Þ , one can numerically derive the conditions for the occurrence of oscillations from the Routh-Hurwitz criterion. When the values of the control parameters are such that (m a =m b or n a =n b and so on) or there is a transient perturbation in the values of these parameters or initial conditions, then the oscillating system will be unstable and driven to any one of the stable steady state solutions as either g b vg a or g b wg a through asymptotic spirals. For example when n b &n a or m b &m a then the stable steady-state solution will be (g a &0, g b &1). These results suggest that contrasting from GG oscillator model the identical two-TF feedback system cannot generate self-sustained oscillations in the presence of stochastic noise. One can also construct one-to-one feedback oscillator via coupling two independent GG oscillators. Here these independent TF oscillators A and B can be coupled via either A-OR-B ( Figure 2B2) or A-AND-B ( Figure 2B3) type logics. One can consider various types of regulatory combinations associated with these network architectures. The combinations in A-OR-B type coupling can be denoted as 'AsAc-BsBc' where 'As' and 'Bs' are the types of selfregulation of TF genes A and B respectively whereas 'Ac' and 'Bc' are the types of their cross-regulation on each other. Each type of regulation can be either 'P' or 'N' where 'P' denotes positive type and 'N' denotes the negative type regulation. Using these notations one can denote the configuration given in Figure 2B2 as NN-NN type, Figure 2B4 as NN-PP type and Figure 2B5 as NP-NP type. The configurations given in Figures 2B4 and 2B5 are the well-studied robust dual-feedback oscillators [25], [33][34][35]. Similarly one can consider various possibilities in A-AND-B type architectures. Noting the symmetry of regulation we find three possible types as P-P, N-N and N-P out of which only N-N will be a robust oscillator. The configuration given in Figure 2B3 is an N-N type. When the coupling is via A-OR-B type logic then both the promoters of TF genes A and B will be independently downregulated upon binding of protein end-products of both the TF genes A and B at the respective cis-regulatory elements associated with each promoter and the first and last equations in Eqs (18) will be modified as follows.
v hq dX hq dt~Z q n hq 1{X h ð Þ{m hq X hq ; X h~X m~a,b X hm ; In these equations for each value of subscript 'h' the subscript 'q' will take a, b and there are totally four equations associated with the overall promoter state dynamics. Various modified parameters in Eqs (21) Here Z hs is defined as in Eqs (20). When m hq~ m m and b h~1 then we can calculate the steady state protein levels from the set of polynomial Eqs (22) as g h~lh zk h ð Þ e y where y is the real root of y n a z1 ð Þ{ ln m m 1{e y ð Þ{e y n b z1 ð Þ À Á 0. When the GG oscillators A and B are coupled through A-AND-B type logic then the dimer (y d ) of both end-products (z a -z b ) will be the key regulating molecule that binds at the cis-acting elements associated with the promoters of both TF genes A and B however with different Hill coefficients (n h ) and subsequently down-regulate them. The respective modified differential rate equations corresponding to the dimerization and binding of dimer at the promoters of TF genes A and B can be written as follows.
The modified and new parameters and variables in Eqs (23) are defined as follows. In the definition of x h for h = a, b one needs to substitute q = b, a. Here l df (M 21 s 21 ) and l dr (s 21 ) are the forward and reverse rate constants associated with the diffusion limited dimerization reaction between the protein end-products of TF gene A and B. The corresponding steady state solutions to Eqs (24) can be written as follows.
Here g h is the solution to the set of following polynomial equations.
In this set of equations we need to set q = b for h = a and for h = b we need to set q = a. The parameters associated with dual feedback oscillators are summarized in Table 2.

Three gene repressilator type oscillators
Similar to Eqs (18) one can write the set of differential rate equations associated with the three TFs repressilator model as follows ( Figure 2C1).
Here the subscripts will be such that h~a,b,c; s~c,a,b; q~b,c,a ð Þ where (a, b, c) denotes respectively TF gene A, B and C in a cyclic (h, s, q) manner and the variables as well as various control parameters are defined as in case of Eqs (18) and generally n a =n b =n c . Similar to Eqs (19) one can derive the steady state solutions with respect to the scaled protein levels for three TFs system as follows.
In this equations y is the real root of y n a n b n c z1 ð Þ {n c n b ln Dzn c ln w{ ln y~0 where we have defined various other terms as, 1{b c e {y=nc y 1=nc : Using these steady state values of protein products g h , one can write the steady state solutions to other dynamical variables can be written similar to Eqs (20) as follows.
Under identical values of all the control parameters such as (m a~mb~mc , n a~nb~nc and so on), the system reaches the steady state as g a~gb~gc ð Þwhich is an unstable fixed point since even a small perturbation in the parameter values or initial conditions will drive the system towards a stable limit-cycle. As depicted in Figures 2C2 and 2C3 one can also construct the repressilator type model by cyclically coupling three independent GG oscillators A/B/C through -AND-or -OR-type logical gates as we have constructed in Figures 2B2-3. When the type of interaction is through -AND-type logic, then the z a -z b dimer down-regulates TF gene B, z b -z c dimer down-regulates TF gene C and the z c -z a dimer down-regulates TF gene A. The set of modified differential equations corresponding to the configuration that is given in Figure 2C3 can be written as follows.
k~b,c,a; q~c,a,b The steady state solutions can be obtained by numerical methods from the following set of algebraic equations.
In these equations for h = a, b, c one needs to set k = c, a, b and various new and modified parameters are defined as follows. When the type of interaction is through -OR-type logic as depicted in Figure 2C2 then the end products of both TF genes A and B can independently down-regulate TF gene B and the end products of TF genes B and C can independently down-regulate TF gene C and so on. The set of modified differential equations associated with such system will be similar to that of Eqs (21)(22)(23) where the indices will be extended for three TF genes A/B/C. One can also consider a fully interconnected network of TF genes A/B/C. In these configurations the self-regulated promoters of each TF gene A/B/C will be negatively regulated by the endproducts of the remaining two other TF genes. Here the mode of overall combinatorial interactions among these regulating endproducts and the corresponding promoters can be either through -AND-or -OR-type logics as represented by the dashed lines in Figures 2C2-3. In case of -OR-type logical gate, the negatively self-regulated promoter of TF gene A will also have cis-regulatory binding sites for the end-products of both TF genes B and C and so on. One can write the modified set of differential equations associated with such fully interconnected configuration (Figure 2C2) In the first one of Eqs (31), there will be three equations for each promoter and there are totally nine equations associated with the overall promoter state dynamics. In the second set of three equations as well as in the associated parameters for each value the subscript 'h' from the set (a, b, c), the subscripts q and k will take the remaining values. This means that when h = a then q = b and k = c and so on. Various modified parameters in Eqs (31)  In case of fully interconnected configuration through -ANDtype logic that is depicted in Figure 2C3 (with dashed lines), the complex z a -z b -z c will be the key regulating molecule that binds with the promoters of all the three TF genes A/B/C and down-regulate them. Similar to Eqs (23) one can write the modified set of differential equations corresponding to repressilator configuration that is fully interconnected through -AND-type logic as follows.
Here we have defined Y d~yd =p as . The steady state solutions to this equation can be obtained by solving the following set of polynomial equations.
Here various modified and new parameters are defined as follows.
e d~cpa p cs p bs l df ; x dh~ldf p qs p ks l hf ; w d~ldr l df p cs p bs ; x h~p n k {1 as d hz k hf . l df p bs p cs In these equations similar to Eqs (31) for each value the subscript 'h' from the set (a, b, c), the subscripts q and k will take the remaining values. This means that when h = a then q = b and k = c and so on for other values.

Perturbation-responses of various gene oscillators
Sample trajectories and phase portraits of GG oscillator for v a~0 ,2|10 {4 À Á are shown in Figures 3A1-3 and 4A1-3. Irrespective of the type of initial conditions and magnitude of the control parameters, the trajectories always start with an overshoot of protein production that is followed by asymptotic spirals towards a stable limit cycle. This seems to be an inherent property of negatively self-regulated loops [9]. Figures 3B1-4 and 4B1-4 suggest that there exists an optimum range of Group I parameters v a~4 {12 ð Þ|10 {4 and w a ,e a ð Þ[ 0:2,1:8 ð Þat which the critical Hill coefficient ( C n a ) that is required to generate selfsustained oscillations is a minimum which is in turn strongly dependent on the promoter state occupancy parameter m a . This optimum range is also dependent on the values of other Group II and III parameters. The optimum range of the conversion parameter seems to be l a [ 0:6{1:4 ð Þ . Results suggest that strong binding conditions m a2 2|10 {4 (v a =0) and m a v10 {12 (v a~0 ) are required to minimize the value of critical Hill coefficient with respect to changes in Group I type parameters. The minimum achievable values of critical Hill coefficients seems to be C n a~6 (v a =0) and C n a~9 (v a~0 ). When there is an additional dimerization step as described in Eqs (4-5) corresponding to the negative-feedback-only (NFO) model considered in reference [25], the minimum achievable critical Hill coefficient seems to be C n a~3 . One should note that in the Lac I oscillatory system the effective Hill coefficient is n a~4 since four dimers of lac I endproducts involved in the overall negative feedback. Numerical analysis of this NFO model system using the physiological range of parameters as given in Table 1 suggests that the period of oscillator can be well tuned by changing the promoter state affinity m a of the repressor without compromising the amplitude much as shown in Figure 3C1. Figure 4B4 shows the strong influence of s a on the critical C n a which means that the approximation (s a~0 ) as in case of most of the earlier studies on various genetic oscillators is not a valid one. At the critical Hill coefficient, the period as well as amplitude of oscillations are strongly dependent on the Group I parameters as shown in Figures 5A1-3. These results also demonstrate how the oscillator responds to temporal perturbations is Group I parameters. As we have shown in the theory section, the period of oscillations increases with increase in the Group I parameters whereas the amplitude seems to decrease as the value of Group I parameters increase. One should note that square of period of oscillation is inversely proportional to the total energy of an oscillator whereas the total energy is directly proportional to the square of amplitude. This means that the total energy of a GG oscillator can be fine-tuned by perturbing the Group I parameters. The Goodwin-Griffith oscillator seems to abruptly enter into the modified limit-cycle orbit upon introducing the perturbation and relax back much faster upon removal of perturbation in the parameter v a~0 rather than perturbations in other parameters w a ,e a ð Þ. In the latter cases, as shown in Figures 5A2-3 the relaxation of oscillator to the original orbit upon removal of perturbation seems to be through slow asymptotic spirals. Figure 5B1 suggest that the period of oscillations increases monotonically with respect to increase in the value of Group I parameters as we have predicted in the theory section. When v a =0 then there exists a range of w a [ 0:3,1 ð Þat which the period of limit cycle oscillations and the required C n a are almost independent of changes in w a . Figure 5B2 shows that when v a =0 then the period of oscillations linearly increases as s a increases whereas it linearly decreases with increase in s a when v a~0 .
The dual feedback motif ( Figures 2B1-3) is an adaptable one that can act as a toggle switch as well as an oscillator depending on the type of configuration. As we have shown in the theory section, the configuration depicted in Figure 2B1 requires identical values of all the control parameters as well as initial conditions to generate coupled as well as synchronized oscillations. Particularly this configuration can efficiently act as a toggle switch since the fixed point g a~gb is an unstable one and even small perturbations in the parameters or initial conditions is enough for the system to exit from the synchronized limit cycle oscillations around this unstable fixed point and subsequently move towards any one of the two stable steady states. The configurations given in Figures 2B2-3 can act as coupled oscillators. Sample trajectories and phase portraits of one-to-one coupled oscillators corresponding to the configuration given in Figure 2B1 are shown in Figures 6A1-6. The minimum achievable value of the critical Hill coefficient that is required to generate self-sustained oscillations around the unstable fixed point (g a~gb ) of dual feedback oscillator seems to be C n h = 5 that is closer ( C n h = 6) to the critical value corresponding to GG oscillator. Variations of critical Hill coefficient with respect to changes in combination of different groups of control parameters are shown in Figures 6B1-8 Figure 6B2 we find that C n h is also independent on the changes in the ordinary perturbation parameter x h . However it is strongly dependent on s h and the condition s h ƒ0:3 is required to achieve the minimum value of critical C n h . Results suggested that when there are identical perturbations in the given control parameter, then the one-to-one coupled oscillator ( Figure 2B1) behaves similar to that of GG oscillator. That is to say the period of oscillations increases and amplitude decreases with an increase in Group I parameters. Here the identical perturbations are such that for the parameter w h we have w h [ w h zd wh where w a~wb prior to perturbation and the magnitude of perturbation is such that d wa~dw b . When any of these two conditions fails, then the system will be driven towards the corresponding stable steady state.
Upon receiving a transient pulse of perturbation or imbalance in the control parameters the dual feedback oscillator exits from the limit-cycle orbit with a time-delay (t del ) and subsequently reaches one of the stable steady-states via asymptotic spirals. Here the target steady state is dependent on the type of disproportion in the parameter values among TF genes A and B. For example when the perturbation is from m a~mb towards m a vm b then the target steady-state will be g a wg b since the binding of TF endproduct B at the promoter of TF gene A is stronger than the binding of end-product A at the promoter of TF gene B. It seems that the value of this time-delay is dependent on the extent of disproportion (p k , where the subscript 'k' denotes the control parameter under consideration such as m h ) as well as duration of the perturbation (t w ) in control parameters or initial conditions associated with the TF genes A and B. Here the percentage of disproportion or imbalance with respect to the parameter k b that is associated with TF gene B is defined as p k~1 00Dk a {k b D=k b . Figure 6C shows the variation of the time-delay with respect to changes in the extent of disproportion (p m ) in the control parameter m h and duration of perturbation t w . It seems that t del approaches zero independently upon increase in both t w and p k . Further simulation results suggest that the time-delay t del is independent on the time (t pulse ) at which the perturbation in the control parameter is introduced into the system.
Contrasting from the configuration given in Figure 2B1, the limit-cycle orbits of the coupled oscillators depicted in Figures 2B2-3 are robust against transient imbalances in the control parameters. The minimum achievable value of the critical Hill coefficient seems to be C n h = 4 for the oscillator with A-OR-B type logic ( Figure 2B2) whereas C n h = 2 for the coupled oscillators with A-AND-B type logic ( Figure 2B3). Results suggest that the limit cycle orbit of coupled oscillators with A-AND-B and A-OR-B type logics are stable one. When there are temporal perturbations in Group I parameters associated with one of the Goodwin oscillators (TF gene A/B) then the other unperturbed oscillator responds to the changes in the behavior of the perturbed oscillator depending on the type of logical coupling between them. As shown in Figures 7A1-2, B1-2 and C1-2 in case of A-OR-B coupling an increase in the magnitude of Group I parameters associated with one of the oscillators A/B does not change the period of the entire system of oscillators (period-buffering) though there is a decrease in the amplitude of the oscillator that is perturbed in w h ,e h ð Þ. The decrease in the amplitude might be partially owing to the period-buffering effect. In case of A-AND-B type logical coupling, increase in the magnitude of Group I parameters w h ,e h ð Þ increases the period of oscillations and decreases the amplitude of the entire system of oscillators that includes both TF genes A/B.  (21). Above results corresponding to A-OR-B type logical coupling with respect to changes in the parameter v h are valid only when the temporal perturbations are the same for a given promoter of TF gene A/B. This means that for TF gene A (here we have subscript k~a) the extent of perturbation should be the same for both v aa and v ab while the set of parameters associated with the TF gene B (k~b) remains unperturbed. Here one should note that v hh controls the dynamics associated with the binding of end-product of TF gene 'H' at its own promoter whereas v hk controls the dynamics associated with the binding of end-product of TF gene 'K' at the promoter of TF gene 'H' as we have shown in Eqs (21). When there are perturbations in only one of these two split parameters (as v a [ v aa ,v ab ð Þ) then the coupled system of oscillators seems to be dynamically unstable and also produces modulated beats as shown in Figures 7A3-4. The period of such beats increases as the imbalance in the set of split parameters v hk increases as shown in Figures 7A5-6. These dynamical instabilities as well as beats abruptly disappear once the perturbations in v hk are removed. Whereas the system of coupled oscillators relaxes back to the initial unperturbed limit-cycle orbit through asymptotic spirals upon removal of perturbations in case of Group I control parameters w h ,e h ð Þ. Results from Figures 7 and 8 suggest that coupled oscillators with -AND-type logic are more robust against promoter state perturbations than the -OR-type coupling. Period of a network of oscillators can be easily fine-tuned by manipulating merely one of the oscillators when the mode of coupling is via -AND-type.
Sample trajectories and phase portraits of a repressilator configuration given by Figure 2C1 are shown in the Figures S1A1-5 in File S1. The minimum achievable value of the critical Hill coefficient for the repressilator seems to be C n h = 2 similar to that of a one-to-one feedback oscillator with A-AND-B type logical coupling. As we have shown in the theory section, the steady state fixed point is more stable when the parameters and initial conditions are identical for all the TF genes A/B/C. When there is a transient perturbation in the control parameters then the system leaves the steady state and enters into a stable limit-cycle orbit through asymptotic spirals with a time delay t del as shown in Figure S1A4 in File S1. Here the magnitude of this time delay seems to be directly proportional to the extent of imbalance or disproportions in the parameter values as shown in Figure S1A4 in File S1. Perturbation in the control parameter v h associated with any one of the TF genes A/B/C results in the decrease of amplitude of the perturbed as well as the one that regulates it. However perturbation in v h does not affect the period of oscillations of the entire system of TF genes as shown in Figures  S1B1-2 in File S1. This means that when v a is increased then the amplitudes of oscillations of TF genes A and C decrease whereas the amplitude of B is not affected. Perturbation in the control parameters w h ,e h ð Þassociated with any one of the TF genes A/B/ C decreases the amplitude of oscillations of the TF gene that is regulated by the end-product of the perturbed TF gene and increases the amplitude as well as width of oscillations of the TF gene that is regulating the perturbed gene. This means that when w a ,e a ð Þincreases then the amplitude of oscillations of TF genes A and B decreases whereas the width and amplitude of TF gene C increases. Further results show that an increase in w h ,e h ð Þ of any one of the oscillators increases the period of oscillations of the entire system of oscillators as shown in the supplementary information Figures S1B3-4 in File S1.
Sample trajectories and phase plane portraits associated with the configuration given in Figure 2C2 are shown in Figures  S2A1-3 and SB1-4 in File S1. Contrasting from three TF genes repressilator model ( Figure 2C1) the configurations given in Figures 2C2-3 do not require any asymmetry in the values of control parameters or initial conditions to trigger the stable oscillations. When the mode of coupling of TF genes A/B/C of GG oscillators is through -OR-type logic then in the presence of identical values of all the sets of control parameters the TF genes A/B/C oscillate in a synchronized manner with respect to period and amplitude. When there is a perturbation in set of the control parameters v hq (here v h will be split into v hh ,v hq À Á for each promoter) associated with any one of the TF genes A/B/C then there are at least three different phases of responses. In the first phase, as shown in Figures S2B1-2 in File S1 the system tries to resist the perturbation by keeping the synchronized limit-cycle orbit intact whereas in the second phase the system becomes unstable and chaotic whose magnitude depends on the extent of perturbation. Upon removal of perturbation, in the third phase the system enters into new asynchronous limit-cycle orbit with stable phase differences among the TF gene oscillators. When there is a perturbation in one of the split parameters v hh ,v hq À Á , then as shown in Figure S2B1 in File S1 the second phase will have several repeating elements of resistance and instability. Perturbation in the control parameters w h ,e h ð Þseems to have similar effects which are evident from Figures S2B3-4 in File S1. Contrasting from these results, the oscillator depicted in Figure 2C3 seems to be more robust against changes in the Group I control parameters and also they return back to the initial coherent type limit-cycle orbit upon removal of perturbations as shown in Figures S2C1-2 in File S1. Sample trajectories of fully interconnected configurations given in Figures 2C2-3 (with dashed lines) are shown in Figures S3A1-3 and SB1-3 in File S1. These results suggests that the fully interconnected three TF genetic oscillator will be more stable against perturbations in the critical control parameters when the mode of coupling is through -AND-type logic than the -OR-type logic.
Tuning capabilities of A-AND-B ( Figure 2B3) and A-OR-B ( Figure 2B2) type coupled dual feedback oscillators are demonstrated in Figures S4A-D and S4A-D in File S1. These results  (4). One needs to substitute Q = M (scaled concentration of mRNA) for red line, Q = X (promoter occupancy) for blue line and Q = Z (endproduct) for pink line. Simulation settings are m a~2 |10 {4 , v a~1 0 {3 , s a ,k a ,x a ð Þ0 and we set e a ,l a ,w a ð Þ1 which required a critical Hill coefficient of C n a = 6 to generate oscillations. Total simulation time is 100 (measured in terms of number of lifetimes of the protein product of TF gene A) and integration step is Dt~10 {5 . A2. Trajectories corresponding to the settings in A1. A3. Roots of the (biquadratic) characteristic polynomial (PI) associated with the Jacobian matrix for settings in A1. B1. Variation of critical Hill coefficient with the parameter set m a ,v a ð Þ. Minimum of this critical value seems to be achieved at m a *10 {4 , and v a *10 {3 . B2. Variation of critical Hill coefficient with the parameter set m a ,e a ð Þ. With the optimized settings in B1, the system seems to be robust when e a [ 0:2,2 ð Þ. B3. Variation of critical Hill coefficient with the parameter set m a ,w a ð Þ. With the optimized settings in B1, the system seems to be robust when w a [ 0:2,2 ð Þ. B4. Variation of critical Hill coefficient with the parameter set m a ,l a ð Þ. Default values of other parameters in B1-4 are as in A1. C1. Variation of period and amplitude of the negative-feedback-only model considered in reference [25] with respect changes in the promoter affinity parameter m a . Simulation settings are given in Table 1. Red solid line in the period and blue solid line is amplitude of the oscillator. doi:10.1371/journal.pone.0104328.g003 (35) show that A-AND-B type coupled oscillators can be tuned by changing the promoter state binding parameter m h efficiently without conceding on the amplitude of oscillations. Increase in m h monotonically decreases the amplitude (and increases the period) of A-OR-B type coupled oscillators. Whereas A-AND-B type oscillator shows two distinct regions of responses with respect to changes in m h namely a responsive region and nonresponsive region. For the settings in Figure S4B in File S1, the A-AND-B  25 . When m h .10 25 then changes in m h will not affect the period of oscillations much. On the other hand the amplitude (measured in terms of P h /P hs ) of A-AND-B oscillator is $1 in a wide range of m h as well as w h values. Further A-OR-B type coupled oscillators behaves similar to that of single GG oscillator with respect to changes in the perturbation of Group I parameter w h . In case of both types of coupled oscillators the amplitude seems to be inversely proportional to the critical C n h rather than the period of oscillations as it is apparent from Figures S4C and D in File S1 and there exists an optimum value of Group I parameter w h at which the amplitude of oscillations is a maximum. Figures S5A-D in File S1 shows that increase in Group I parameter v h increases the period of oscillations and decreases the critical C n h in both -AND-and -OR-type coupled oscillators. However the period of -AND-type oscillator seems to be more robust against changes in v h than -OR-type coupled dual feedback oscillator. Contrasting from these there exist a cutoff value of e h ,2 (for the simulation settings in Figure S5 in File S1) beyond which the amplitude of oscillations is practically zero in both types of coupled oscillators that is apparent from Figures S5C and D in File S1. The discontinuities in plots depicted in Figures S4 and S5 in File S1 mainly originate from the fact that each data point belongs to different values of critical Hill coefficients. For example upon iterating the perturbation parameter, we first find out the critical Hill value for the perturbed values and then obtain its amplitude and period of the perturbed system with newly calculated hill coefficient. This means that each data point in those plots comes from altogether different dynamical systems since the Hill coefficients are different. Therefore we may not expect continuous type plots either.

Discussion
Transcription factor gene oscillators play critical roles in driving cell-cycle to circadian rhythms. Here we have identified the critical control parameters associated with self-sustained oscillations of such oscillators and classified them into Groups I, II and III depending on their functionality. Group I parameters control the intracellular dynamics of synthesis and degradation of various molecules associated with regulated TF gene (Figure 1). The parameter w h of Group I describes the strength of coupling between the rate of degradation of mRNAs and the rate of degradation of corresponding TF proteins. We should note that transcription and translation of various TF genes of prokaryotes are simultaneously taking place well within the cytoplasm whereas in case of eukaryotes the transcription is taking place inside the nucleus and the synthesized pre-mRNA transcripts need to be spliced within the nucleus and then transported to cytoplasm after other post-transcriptional modifications through nuclear pores for translation. These differences in the cellular architecture demands higher lifetimes (1/c mh ) for eukaryotic mRNAs than the prokary-otic ones which results in the general observation that the values of w h associated with various genes in prokaryotes are lower than eukaryotes genes [36], [37]. It seems that in case of yeast, the genome-wide values of w h varies from 0.1 to 1 with a median of ,0.3 [36]. These observations suggest that the naturally occurring or synthetic oscillatory motifs in the transcription networks should operate well within this range of w h [ 0:1,1 ð Þ with different median values depending on the type of organism. For most of the bacterial genes we find that w h e 0:1 [8], [9]. The parameter v a describes the dynamics of binding-unbinding of the end-product molecule with the promoter of TF gene A. Large values of v a indicate slower changes in the residence state of promoters whereas small values indicate the faster dynamics of the promoter state towards the equilibrium. Earlier studies suggested [8], [9] a typical physiological value of v a as v a e 10 {4 for a prokaryotic selfregulatory systems (n a = 1) such as the one in E. coli. Whereas in case of eukaryotic systems such as yeast and human its value will be scaled up respectively to v a e 10 {3 and v a e 10 {2 owing to dilution in the local concentrations of various molecules upon increase in the nuclear volumes. One should note that yeast nucleus is ,10 1 times higher than E. coli cell whereas human cell nucleus is ,10 2 times higher than E. coli cell. The parameter e a describes the dynamics associated with of conversion of the protein product to end-product towards the equilibrium. This conversion reaction can be either a first-order conformational transition or pseudo first-order chemical modification of the protein product via an additional catalyst. Here e a also acts as an indirect delay parameter that in turn relates the protein decay rate with the conversion rate. One should note that the binding-parameter m a is the central one that connects the entire Group I type singular perturbation parameters. By definition m a is inversely proportional to the affinity of the end-product towards the promoter sequence. Along with m a the ordinary parameters s a ,l a ,k a ð Þ decide the steady state value of TF protein. Further one should note that m a (or m h in general) is the only parameter that can be externally modified in a working model and all the others are fixed system parameters which can be modified only at the time of designing the oscillatory module. In case of Lac based oscillators m a can be modified through changing the concentration of IPTG [25] which is a gratuitous inducer of lacI gene.
Understanding the design principles associated with the genetic oscillator is one of the central topics in synthetic and systems biology. An efficient synthetic oscillator should be robust against transient perturbations in the control parameters and fluctuations in the promoter state occupancies. The period of oscillators should be easily tunable without compromising the amplitude. Stricker et.al in reference [25] suggested that the robustness of the dualfeedback oscillators depicted in Figures 2B4-5 can be further increased by the explicit delay owing to the oligomerization steps associated with the end-products of TF genes A/B (here gene A can be lacI and gene B can be araC). Further studies by Tsai et.al in reference [35] suggested that unlike pure negative feedback systems, inclusion of a positive feedback loop within the negative feedback oscillators ( Figure 2B4-5) can increase the robustness and tunable nature of the original oscillator without compromising on the overall amplitude. Here the positively auto regulated TF gene acts as a booster for the overall protein level of the negatively self-regulated oscillatory TF protein. In this context one can also consider a positive coupling of a positively auto-regulated TF protein K (booster) with N-N or NN-NN type dual feedback oscillators as depicted in Figure 2B6. Upon analyzing various classes of oscillators depicted in Figure 2 one can conclude that the effects of perturbations in the control parameters and promoter state fluctuations will be minimum when the mode of coupling among network of oscillators is through -AND-logic.
From Figures S4 and S5 in File S1 we find that the A-AND-B coupled dual feedback oscillator can be externally tuned by modifying the promoter state binding m h without conceding on the amplitude of oscillations. The inherent chaotic response of the -OR-type coupled oscillators could be one of the reasons why the oscillations of most of the synthetic oscillators disappeared after certain period of time [18], [24] under in vivo conditions. Nevertheless the overall period of network of oscillators coupled through -OR-type logic are more resistant (period-buffering) to perturbations in the control parameters than the oscillators coupled through -AND-type logic. One should note that the temporal perturbations in the system parameters ultimately result in perturbations in the intra cellular concentrations of protein endproducts of coupled TF genes that in turn are looped back into the system. In this context the -OR-type coupled genes respond to perturbations in an additive way whereas those genes coupled through -AND-type logic respond in a multiplicative way. This could be the reason why the period of oscillations of -OR-type coupled oscillators are more resistant to perturbations than -ANDtype coupled oscillators. This means that the period of network of oscillators coupled through -AND-type logic can be more easily tuned by perturbing any one of the coupled oscillators. Depending on the circuit functionality one can also chose a combination of both the types of coupling. So far we have assumed a copy number of TF genes as one (d hz = 1). This may not be true in the bacterial based synthetic circuits constructed on plasmids since the plasmid copy number will be more than one in most of the times which can lead to further complexity. For example, a plasmid copy number of two for a NFO model that was constructed in reference [25] can mimic a combination of oscillators depicted in Figure 2B2 and 2B3 with TF gene A = B since it can mimic both -AND-as well as -OR-type coupled GG oscillators. This could be one of the possible reasons for observing stable and robust oscillations with such constructs [25] over several bacterial generations. In summary, in this paper we have considered various types of transcription factor oscillators by explicitly incorporating the promoter state dynamics and other chemical reaction balances in detail. Using our detailed model we have identified and classified various critical control parameters and numerically obtained their physiological and optimum ranges to generate self-sustained oscillations in the intracellular levels of mRNAs and transcription factor proteins. We further derived the basic design principles associated with robust and tunable gene oscillators. We have further demonstrated that by coupling two or more independent Goodwin-Griffith oscillators via -OR-or -AND-type logics one can construct genetic-oscillators which are fine-tunable and also robust against perturbations in the system parameters. When there is a perturbation in one of the -OR-type coupled oscillators, then the overall period of the system remains constant whereas in case of -AND-type of coupling the overall period of the system moves towards the perturbed oscillator. Though there is a periodbuffering, the oscillators coupled through -OR-type logic seems to be more sensitive to perturbations in the parameters associated with the promoter state dynamics than -AND-type.

Materials and Methods
We use Euler type numerical scheme [8], [9], [38] to integrate the set of differential rate equations corresponding to various types of oscillatory loops. For example in case of Eqs (4) which describe the Goodwin-Griffith model, the numerical integration scheme can be written as follows.
Initial and boundary conditions are X a ,M a ,P a ,Z a ð Þ0 and X a ,M a ,P a ,Z a ð Þ ƒ1 respectively. The scaled time step Dt should be chosen such that it captures the dynamics of all the variables including the dynamics of promoter state occupancies (first one in Eqs (34)). We divide the total scaled simulation time T into N equal intervals such that Dt = T/N. For simulation purpose we set Dt = 10 25 and the corresponding Dt = 0.03 s for a lifetime 1/ c pa ,60 mins. We use Newton's method [38] to find the fixed point solutions to steady state equations. Here to obtain the solution to a nonlinear algebraic equation f x ð Þ~0 one uses the iterative scheme x iz1~xi {f x i ð Þ=f 0 x i ð Þ. For example, the iterative numerical scheme to obtain the fixed point solution P a~ga to Eqs (9) particularly for n a w4 can be written as follows. In this numerical scheme we set the initial guess value of the fixed point solution as g a,i~0~0 and the tolerance limit for stopping iteration as Dg a,i {g a,iz1 Dƒ10 {5 . We can use the following computational workflow. For example, in case of oneto-one dual feedback oscillator one needs to construct the eighth dimensional Jacobian matrix associated with Eqs (18) as follows.
are as in A1-2. Q = X for red line, Q = M for green line and Q = P for blue line. A5. Phase portraits of TF genes A and B. Q = P for red line, Q = M for blue line and Q = X for green line. A6. Roots of the eighth dimensional characteristic polynomial derived from the Jacobian matrix associated with Eqs (18) for the parameter settings given in A1-2. B1-8. Variation of critical Hill coefficient that is required to generate oscillations with respect to changes in various control parameters. Default settings of other parameters are as in A1-2. C. Variation of t del with respect changes in percentage of disproportion p m and pulse width t w . doi:10.1371/journal.pone.0104328.g006 Here subscripts (a, b) denotes the TF genes A and B respectively. Various matrix elements are defined as follows.
In these definitions for h = a, b one needs to substitute k = b, a. To obtain the critical value of the Hill coefficient ( C n a ) one need to first solve the steady state equations. Upon substituting the steady state protein values into to the corresponding Jacobian matrix one can construct the characteristic polynomial and the Routh table. The corresponding inequality conditions for oscillations will be derived from this table. This procedure need to be carried out at various values of n a from n a = 1 in an iterative way. We set the default initial value of promoter state variable as X h~5 |10 {2 for any one of the TF genes A/B/C to trigger the limit-cycle oscillations in case of repressilator so that the effect of other Group I control parameters on the initial delay in oscillations can be studied. One should note that this initial delay is still a function of the disproportion in the initial conditions. We measure the concentrations of various molecules in terms of number of molecules inside the cell. Considering a typical E. coli bacterial cell (volume ,10 218 m 3 ) [39] we set d hz~1 , m hs ,10 2 molecules and p hs ,10 4 molecules [8], [9] where h = a, b and c respectively denotes TF genes A, B and C. Concentration of a single TF molecule inside a bacterial cell will be ,2 nM [40][41][42]. We assume a lifetime of TF protein as ,2 min and mRNA lifetime as ,0.2 min for E. coli that gives a decay rate of c pa e 10 {2 s {1 so that w a = 0.1 and we measure the time in terms of number of protein lifetimes. Since the dynamics of binding-unbinding of a single transcription factor protein (n a = 1) with its cis-regulatory site on DNA is a typical diffusion-controlled site-specific DNAprotein interaction, under in vivo conditions of a bacterial cell we find that k hf e 10 {7 molecules 21 (c p ) that will be scaled down to the interaction of n h number of TF molecules with sequentially located combinatorial binding sites as k hf [ k hf n h À Á molecules {na s {1 [26], [27]. The time associated with the unbinding reaction will be closer to that of the protein lifetime.
Here .
Since v a ,10 24 for a typical bacterial promoter [8], [9] we assume c ph p na as k hf À Á e 10 {4 for an arbitrary Hill coefficient n a . Using these values one can estimate the physiological ranges of various groups of control parameters as given in Table 1. With these settings for a bacterial lacI (n a = 4) system (Eqs. (4)(5)) the values of parameters will be v a v10 {4 and m a e 10 {5 . For simplification purpose we can assume identical values for similar group of parameters n h = n hq , m h~mhq and so on for other family of parameters such as x h and e h .

Supporting Information
File S1 This file contains Figure S1- Figure S5. Figure S1. Three gene repressilator model. A1-4. Phase portraits and trajectories of TF genes A, B and C of a repressilator. Simulation settings are s h ,k h ,x h ð Þ0, e h ,l h ,w h ,r h ð Þ1, m h~1 0 {4 and v h~1 0 {4 which required a critical Hill coefficient of C n a = 2. Total simulation time is 200 (number of lifetimes of the protein product of TF gene A) and integration step is Dt~10 {5 . To trigger the oscillations, we have introduced the asymmetry in the initial condition for the promoter state occupancy of TF gene A as X a~5 |10 {2 . Oscillations starts with a time delay t del whose value depends of the magnitude of this disproportion in the parameter values. A5. Roots of the twelfth degree characteristic polynomial associated with the Jacobian matrix of Eqs (26) for settings given in A1. B1-2. Effects of perturbation in v a that is raised to v a~1 0 {3 (v a~1 0 {2 in B2) in the time interval from 0 to 100. Increase in v a increases the period of oscillations of the entire system from t p *23 to 24.5 and reduces the amplitudes of TF genes A and C. The amplitudes of TF genes A/B/C are such that A,C,B. B3-4. Effects of perturbation in w a ,e a ð Þwhich are raised to w a~3 in the time interval from 0 to 100. Increase in w a increases the period of oscillation of the entire system from t p *23 to 30 and reduces the amplitudes of TF genes A and B and increases the amplitude of C and the amplitudes of TF genes are such that B,A,C (B3). Increase in e a increases the period of oscillation of the entire system as in B3 where the amplitudes of TF Figure 7. Dynamics of dual feedback Goodwin-Griffith oscillators coupled through OR gate. A1-2. Trajectories of protein-products of TF genes A and B which are two independent GG oscillators coupled through A-OR-B type logic as given in Figure 2B2. Simulation settings are s h ,k h ,x hq À Á 0, e h ,l h ,w h ,r h ð Þ1, m h~1 0 {5 and v h~1 0 {5 which required a critical Hill coefficient of C n a = 8. Total simulation time is 200 (number of lifetimes of the protein product of TF gene A) and integration step is Dt~10 {5 . For each promoter A/B the parameter v h will be split into v h [ v hh ,v hq À Á where v hh corresponds to self-regulation and v hq corresponds to cross regulation. Under identical values of all the parameters the system generates synchronized oscillations with a period of t p *4:5. Upon introduction of perturbation in v a from scaled time 0 to 100 (where v a~1 5|10 {5 ), the amplitude of TF gene A is reduced with a phase shift and the period of entire system that includes both TF genes A and B remains the same. A2 is a magnification of certain range of A1. A3-4. Effect of perturbation in only one of the split parameters v hh ,v hq À Á associated with TF gens A/B. Here v aa is perturbed to v aa~1 5|10 {5 . The system seems to be unstable and generates beats. A4 is a magnification of certain range of A3. A5-6. Here v aa is perturbed to v aa~3 0|10 {5 . Period of beats seems to increase as the disproportion among the split parameters increases. A6 is a magnification of certain range of A5. B1-2. Effect of perturbation in the parameter w a which is raised to w a~1 :5 from the default value in the interval from 0 to 100. Increase in w a reduces the amplitude of both the TF genes A and B with a phase shift and the period of oscillations of the entire system remains the same as in A1-2. B2 is a magnification of certain range of B1. C1-2. Effect of perturbations in the parameter e a which is raised to e a~1 :2 from the default value in the time interval from 0 to 100. Increase in e a reduces the amplitude of both the TF genes A and B without a phase shift and the period of oscillations of the entire system remains the same as in A1-2. C2 is a magnification of certain range of C1. doi:10.1371/journal.pone.0104328.g007 Figure 8. Dynamics of dual feedback Goodwin-Griffith oscillators coupled through AND gate. A1-4. Trajectories of protein-products of TF genes A, B which are two independent GG oscillators coupled through A-AND-B type logic as given in Figure 2B3. Simulation settings are s h ,k h ,x hq À Á 0, e h ,l h ,w h ,w d ,r h ð Þ1, m h~1 0 {5 and v h~1 0 {5 which required C n a = 2 (here we have set n h = 4 for clarity of results). Total simulation time is 200 (number of lifetimes of the protein product of TF gene A) and integration step is Dt~10 {5 . Under identical values of all the parameters the system generates synchronized oscillations with a period of t p *9. Upon introduction of perturbation in v a from scaled time 0 to 100 (where genes A/B/C are such that B,A,C (B4). Figure S2. Dynamics of three independent Goodwin-Griffith oscillators cyclically coupled. A1-3. Phase portraits of TF genes A/ B/C which are three independent GG oscillators cyclically coupled through -OR-type logic as given in Figure 2C2 (without dashed lines). Simulation settings are s h ,k h ,x h ð Þ0, e h ,l h ,w h ,r h ð Þ1, m h~1 0 {5 and v h~1 0 {4 which required a critical Hill coefficient of C n a = 5 (we have set this to 6 for clarity of results). Total simulation time is 500 (number of lifetimes of the protein product of TF gene A) and integration step is Dt~10 {5 . In this configuration the end-products of TF gene A and C will regulate TF genes A through A-OR-C type logic whereas the promoter of TF gene B will be regulated by the end-products of TF genes A and B through A-OR-B type logic and so on. Under identical values of all the parameters the system generates synchronized oscillations with a period of t p *7. B1-2. Effects of perturbation in v h . Upon introduction of perturbation in v aa from scaled time 100 to 400 (where v a~6 |10 {4 ) there three phases of responses. In the first phase, the system tries to resist the perturbation whereas the second phase consists of repeating elements of resistance and chaos. Upon removal of perturbation the system enters into new limit-cycle orbit in the third phase. In B2 both v aa and v ac are perturbed as in B1. B3-4. Effects of perturbation in w a ,e a ð Þwhich are raised to e a ,w a ð Þ~2 in the time interval from 100 to 400. System responds to the perturbation as in B1-2. C1-2. Trajectories of TF genes A/B/C which are three independent GG oscillators cyclically coupled through -AND-type logic as given in Figure 2C3 (without dashed lines) and their responses to perturbations in Group I control parameters. Simulation settings are s h ,k h ,x h ð Þ0, e h ,l h ,w h ,w d ,r h ð Þ1, m h~1 0 {5 and v h~1 0 {4 which required C n a = 2. Total simulation time is 500 and integration step is Dt~10 {5 . Identical values of all the parameters of the system generate synchronized oscillations. Introduction of perturbation in v a from scaled time 100 to 400 (where v a~3 |10 {4 , C1) affects only TF gene A whereas the orbit of other TF genes B/C seems to be stable and upon removal of the perturbation the system returns back to initial limit-cycle orbit. In C2 the parameter w a is perturbed to w a~2 as in C1. Figure S3. Dynamics of three Goodwin-Griffith oscillators which are fully interconnected. A1-3. Phase portraits of TF genes A/B/C which are three independent GG oscillators which are fully interconnected with -OR-type logic as given in Figure 2C2 (with dashed lines). Simulation settings are s h ,k h ,x h ð Þ0, e h ,l h ,w h ,r h ð Þ1, m h~1 0 {5 and v h~1 0 {4 which required a critical Hill coefficient of C n a = 6. Total simulation time is 500 (number of lifetimes of the protein product of TF gene A) and integration step is Dt~10 {5 . In this configuration all the endproducts of TF gene A/B/C will regulate all the three TFs through A-OR-B-OR-C type logic. Identical values of all the parameters of the system generate synchronized oscillations. Perturbation in v a from scaled time 100 to 300 (where v a~5 |10 {5 , A1) seems to make the system unstable. In A2 w a is perturbed to w a~3 as in A1 and in A3 e a is perturbed to e a~3 as in A1. B1-3. Phase portraits of TF genes A/B/C which are three independent GG oscillators fully interconnected with -ANDtype logic as given in Figure 2C3 (with dashed lines). Simulation settings are s h ,k h ,x dh ð Þ0, e h ,l h ,w h ,r h ð Þ1, m h~1 0 {5 and v h~1 0 {5 which required a critical Hill coefficient of C n a = 2. Total simulation time is 500 (number of lifetimes of the protein product of TF gene A) and integration step is Dt~10 {5 . In this configuration all the end-products of TF gene A/B/C will be regulated by their complex. Identical values of all the parameters of the system generate synchronized oscillations. Introduction of perturbation in v a from scaled time 100 to 300 (where v a~1 0 {4 in B1) seems to make the system unstable. In B2 w a is perturbed to w a~3 as in B1. In B3 e a is perturbed to e a~3 as in B1. Figure  S4. Tuning dynamics of A-OR-B and A-AND-B types of coupled oscillators with respect to perturbations in m and w. A, C. Tuning capability of GG oscillators coupled through A-OR-B type logic as given in Figure 2B2. Default Simulation settings are s h ,k h ,x hq À Á 0, e h ,l h ,w h ,r h ð Þ1 and v h~1 0 {5 . B, D. Tuning capability of GG oscillators coupled through A-AND-B type logic as given in Figure 2B3. Default Simulation settings are s h ,k h ,x hq À Á 0, e h ,l h ,w h ,w d ,r h ð Þ1 and v h~1 0 {5 . Plots A and B show the variation of period, critical C n h and amplitude with respect to changes in m h (iterated from 5610 27 to 10 23 with w h = 1) whereas plots C and D show the variation of these quantities with respect to changes in w h (iterated from 0.1 to 10 with m h~1 0 {5 ). Here period of oscillator is measured in the number of lifetimes of TF protein A (1/c a ) and amplitude is measured in terms of number of P h /P hs . Figure S5. Tuning dynamics of A-OR-B and A-AND-B types of coupled oscillators with respect to perturbations in v and e. A, C. Tuning capability of GG oscillators coupled through A-OR-B type logic as given in Figure 2B2. Default Simulation settings are s h ,k h ,x hq À Á 0, e h ,l h ,w h ,r h ð Þ1 and m h~1 0 {5 . B, D. Tuning capability of GG oscillators coupled through A-AND-B type logic as given in Figure 2B3. Default Simulation settings are s h ,k h ,x hq À Á 0, e h ,l h ,w h ,w d ,r h ð Þ1 and m h~1 0 {5 . Plots A and B show the variation of period, critical C n h and amplitude with respect to changes in v h (iterated from 5610 27 to 10 24 with e h = 1) whereas plots C and D show the variation of these quantities with respect to changes in e h (iterated from 0.7 to 8 with v h~1 0 {5 ). Here period of oscillator is measured in the number of lifetimes of TF protein A (1/c a ) and amplitude is measured in terms of number of P h /P hs . (ZIP)

Author Contributions
Conceived and designed the experiments: RM. Performed the experiments: RM. Analyzed the data: RM. Contributed reagents/materials/ analysis tools: RM. Wrote the paper: RM. v a~1 5|10 {5 in A1, v a~3 0|10 {5 in A3), the amplitude of TF gene A is reduced with a phase shift and the period of entire system (both TF genes A and B) increases to t p *10 in A1. A2 and A4 are magnifications of certain range of A1 and A3. B1-2. Effect of perturbation in the parameter w a which is raised to w a~2 from the default value in the time interval from 0 to 100. Increase in w a increases the period of oscillations of the entire system to t p *9:5 as in A1-2 and reduces the amplitude of TF genes A with a phase shift. B2 is a magnification of certain range of B1. C1-2. Effect of perturbations in the parameter e a which is raised to e a~1 0 from the default value in the time interval from 0 to 100. Increase in e a increases the period of oscillations of the entire system to t p *9:5 as in A1-2 and reduces the amplitude of both TF genes A and B without a phase shift. C2 is a magnification of certain range of C1. doi:10.1371/journal.pone.0104328.g008 Theory on the Dynamics of Oscillatory Loops PLOS ONE | www.plosone.org