The repressilator is an artificial synthetic gene network created and named by Elowitz and Leibler  . The simplicity of the network’s structure and design makes the repressilator an ideal candidate for studies on the dynamic features of synthetic gene networks   . The network exhibits negative feedback and nonlinear interactions, both of which are important features of dynamical systems exhibiting oscillations    . The biochemical details of the design principles for this oscillator are complicated, yet for the purposes of this work we present a brief background on the biology of the network: A gene undergoes a process called transcription, which consists in making a “copy” of the gene called messenger RNA (mRNA). Once mRNA is produced, it diffuses out of the nucleus into the cytoplasm and attaches to a ribosome, where it is “read” through a process called translation and thus a protein is created. Through this gene-mRNA-protein process, the cell has the opportunity to carry out all the various tasks encoded in its DNA.
The regulation of the process described above is an important enterprise that the cell needs to control and fine-tune constantly  . Fortunately, there are various control mechanisms that the cell uses to regulate mRNA and protein concentrations. One of these is exemplified by the repressilator and consists in autoregulating mRNA production at the transcriptional level. This autoregulation mechanism happens through a feedback loop between three mRNAs and three proteins. Figure 1 describes the network structure of the repressilator, where it shows a system of three mRNAs coupled to their associated protein products in a cyclic way. The first protein represses production of the second protein, the second represses production of the third, and the third represses production of the first (for more details see   ). These three cyclic repression mechanisms happen at the transcriptional level and are considered to be nonlinear processes modeled via a Hill function. The mathematical model of the repressilator  is given by the following six coupled first order ordinary differential equations (ODEs)
where and are defined as , , and . Here α represents dimensionless transcription rate in the absence of repressor, the background production rate of protein in the presence of saturation, β is the dimensionless ratio between protein decay and mRNA decay, and n is the Hill coefficient representing the degree of cooperation of repression (for more biological details see  ).
The biological and experimental applications of this model are well documented    . The biochemical construction of the repressilator was carried out in the bacteria Escherichia coli (E. coli) and the design was carefully crafted so that oscillations would be exhibited by a green fluorescent reporter gene  . The outcome of this genetically programmed network was a periodic green flashing of E. coli cells which demonstrated that oscillations can be synthetically constructed. Interestingly, in the original experimental measurements by Elowitz and Leibler  oscillations were found to be somewhat noisy and erratic. However, recent new measurements by Potvin-Trottier et al.  show that “some of the erratic behaviour originally reported was due to the limited imaging platforms available at the time.” In their study, Potvin-Trottier et al. show that the repressilator exhibits “highly regular and robust oscillations,” which suggests that the system is stable with respect to perturbations and thus resembles a natural biochemical network. These experimental findings are confirmed in this
Figure 1. Repressilator network diagram. The first protein, , represses production of the second protein, , the second represses production of the third, , and the third represses production of the first. These repressions occur at the transcriptional level, where mRNA is repressed by protein for and defined as , , and . The feedback loop is characterized by nonlinear transcriptional interactions and the six-dimensional ODE model associated to the network is given in Equations (1) and (2). Here the arrow ( ) represents production and the perpendicular symbol ( ) represents repression.
theoretical work by studying the direction and stability of the periodic solutions of the model.
In this work we study the periodic solutions of the repressilator by means of a bifurcation analysis. In Section 2 we compute the steady state solutions and show that these undergo a transition from stable to unstable giving rise to a limit-cycle born in a Hopf bifurcation. In Section 3 we present a nonlinear analysis of the equations, which involves a center manifold reduction on the six-dimensional system. The latter yields closed form expressions for the steady state, frequency, and amplitude of the oscillation, all of which are analyzed through a parameter study in Section 5 and confirmed in a numerical continuation analysis in Section 6. In Section 7 we present our conclusions and discuss the biological significance of our results.
2. Linear Stability Analysis of the Equilibrium Solution
We start our analysis by showing that Equations (1)-(2) have at least one biologically significant equilibrium solution. The equilibria are found by setting
. This gives
where represents the steady state solution. Substituting (4) into (3) we obtain
which can be transformed into a nonlinear algebraic system of equations. The solutions to the latter can be approximated using a numerical root finding technique, such as Newton’s method for systems of nonlinear equations. However, in this work we are interested in biologically significant solutions where . Substituting the latter into (5) gives the following polynomial
and since we are only interested on the positive real roots, by Descartes’ rule of signs the polynomial (6) has either one or three real positive roots for all . This shows that the system has at least one positive biologically significant equilibrium solution when . Simulations and plots are provided in Section 5.
To find the stability of the steady state, , we define and to be deviations from equilibrium and . Substituting these into Equations (1) and (2) results in the following nonlinear system
where and , , and . Expanding for small values of , Equation (7) becomes
where the Taylor coefficients A, , and are given by
Next we analyze the linearized system coming from Equations (8) and (9)
which were obtained by truncating the nonlinear terms in Equation (9). The Jacobian at the origin for system (13)-(14) is
and the associated characteristic equation is found by setting , which gives the following equation
Since the steady state is stable then the real parts of the eigenvalues are all negative, and as we vary β there is a critical value, , where the first pair of complex conjugate eigenvalues cross the imaginary axis. Thus substituting into (16) gives
which is the condition for a change in stability and a Hopf bifurcation. Notice that for (i.e. ) the system (13) and (14) will exhibit solutions of the form
where and are the amplitudes of the and oscillations, and where is a phase angle. As we add a small detuning off of the critical value, , the nonlinear system (1) and (2) is expected to exhibit periodic solutions, which come with a change in stability for the steady state. This change in stability is given by Equation (17), where the stable region is
given by and the unstable by . In Section 5
we compute and plot the associated stability diagrams for some parameter values.
3. Center Manifold Analysis
We use a center manifold reduction to determine the amplitude and direction of the limit cycle bifurcation. This will be accomplished by studying the system at the critical parameter values and . Solving Equation (17) for we obtain
where A is given by Equation (10). Substituting condition (17) into (16) for , setting , and solving for gives
which consists of two branches associated with the in Equation (20) for . More details presented in Section 5.
We start the analysis by expressing system (8) and (9) as follows
where so that , J is given by Equation (15), and
For and we know that J has a complex conjugate pair of eigenvalues on the imaginary axis, , and thus we let be the associated eigenvectors corresponding to and , respectively. These eigenvectors will satisfy the following three conditions
where is the standard scalar product in . This yields
and where we have chosen the scaling factor so that . This completes the linear part of the analysis.
To study the nonlinear part of the center manifold approximation we start by rewriting as follows
where and the multilinear functions and are given by the following expressions
By the Center Manifold Theorem  we know there is a locally defined smooth 2-dimensional invariant manifold (surface) that is tangent to the 2-dimensional eigenspace generated by q and p. Thus we express the solution of (22) as
where is the coordinate of q on the (flat) two-dimensional eigenspace or center subspace spanned by q and p. Here w is the “rest of the solution” which does not lie in the center subspace, but rather on the center manifold. Notice that the center subspace is tangent to the center manifold at the origin, and since y is the coordinate of q such that , then it will be the projection of x onto the center subspace. This gives the following expression for the time derivative of
where and since then it may written as so that Equation (34) can be expressed as follows
Equations (35) and (36) represent the flow on the flat 2-dimensional eigenspace, which will yield closed form expressions for the limit cycle oscillations. To calculate and we start by solving Equation (33) in terms of w and substitute to obtain
Taking the derivative of the latter and substituting Equation (34) we obtain
The normal form for the local parametrization of the truncated center manifold (neglecting cubic and higher order terms) is given by (see   )
which we compute by taking the derivate and equating to Equation (38). The latter yields the following expressions for the parametrization coefficients
and where we omit the expressions of Equations (47)-(49) for brevity. Equation (43) represents the center manifold approximation, , which we substitute into (34) to obtain the flow on the 2-dimensional center subspace
Since then Equation (50) can be expressed as follows
Using the expressions for Equations (56) and (57) we may express the results in terms of polar coordinates and use a near-identity transformation to change the flow to the following equations
where the stability coefficient is given as follows
We refer the reader to   for the standard derivation of Equation (60). Substituting the expressions for J, w0, q, and p from Equations (15), (21), (27), and (28), respectively, yields the stability coefficient for the system at the critical parameter values. For brevity, we omit here the full expression for Q and present the corresponding numerical plots and results in Section 5.
4. Unfolding the Center
In this section we use the center manifold computation to approximate the amplitude of the limit cycle born at the Hopf bifurcation. Our efforts in the previous chapter gave the expressions of the Hopf point at the critical parameter values, which provided the “reduced” system at and . In this chapter we take the second step to compute the normal form of a system with bifurcating parameters, which consists in using a perturbation method to add “unfolding” terms to the reduced normal form found previously. We begin the perturbation approach by adding a small detuning off of the critical value
so that the nonlinear system (1)-(2) is expected to exhibit periodic solutions. This occurs due to the transversality condition of the eigenvalues which will yield a small increase in the real and imaginary parts of the . We thus assume the resulting expression is of the form , where R and are given by and . Thus Equations (35) and (36) will take the approximate form
which can be transformed into polar coordinates (and by means of a near-identity transformation) into the following flow on the center subspace
Setting Equation (64) to zero and solving for r gives the limit cycle amplitude as
where Q is given by Equation (60) and R is found by analyzing the linearized system (13) and (14) which has solutions of the form
where . Substituting the latter two equations and (61) into (13) and (14), linearizing for small , and solving for R and we obtain expressions of the form
where and are long and complicated expressions in terms of β, α, n, and A omitted here for brevity. Substituting (69) and (60) into (66) gives the closed form expression for the limit cycle amplitude, r, which can be computed numerically for different parameter values. In the next chapter we present our results.
5. Parameter Study
In this section we use our closed form expressions to obtain more information on the dynamics of the system. We start by computing the steady state concentration, , for different parameter values. For our model, a biologically significant range for the Hill coefficient, n, is given by (see  ), which also agrees with the equilibria condition for positive real roots of the polynomial in Equation (6). Thus the steady state solution, , is determined by solving Equation (6) for given values of α, α0, and n, where we note that β does not play a role in the values of . Figure 2 shows displayed as a function of α for , and . Both of these figures were plotted by numerically solving Equation (6) and confirmed using MATLAB’s built-in function ode 45.m. Figure 2 for shows one numerical simulation for and , denoted on the plot with an asterisk (*), where we can see that on both cases.
Next we use Equation (17) to obtain the α-β stability diagram for the system. To find an analytical expression involving only α we solve Equation (6) for when and to obtain the real solution
which we use to substitute into the Hopf bifurcation condition (17) where A is given by Equation (10). Alternatively, we may use Equation (20) to compute
Figure 2. Equilibrium solution as a function of α for and (left) and (right). Solution is computed from Equation (6) for given α, α0, and n. One MATLAB numerical simulation is presented on the (left) plot, where and .
which will give us two branches: and . Figure 3(a) shows the α-β parameter
space divided into stable and unstable regions for and . The dividing curve is formed by a lower β1-branch and an upper β2-branch with the two meeting at , which is where the Hopf bifurcation occurs and the system exhibits periodic solutions. In addition, for Figure 3(b) we use Equation (21) and set to obtain the critical frequencies, , as a function of α on both β-branches.
To study the direction of the Hopf bifurcation we find the stability coefficient, Q, for the system. Using Equation (60) we obtain Figure 3(c) where we have plotted Q as a function of α for and . Here we see that for all on both β-branches, which shows that the system exhibits a surface of periodic solutions on the center manifold with stable limit cycles  . The latter together with the results of Figure 3(a) implies that that the stable limit cycle only exists inside the unstable region bounded by the β-curves, which shows that the Hopf bifurcation is supercritical.
From Equations (66) and (69) we see that the amplitude, r, of the oscillation is
given by the product . Figure 3(d) shows r as a function of α for
, , , and . Notice that the amplitude of the oscillation grows larger as α increases, which shows how the limit cycle’s size changes as we follow the bifurcation curve β for a some fixed . For we also obtain growing oscillations as α increases to 5.4426 and after which the amplitudes start decreasing (figure not presented here). Figure 3(d) also presents one MATLAB ode45.m numerical simulation represented with an asterisk (*) in (a)-(d). Thus the numerical results obtained with our closed form expressions for , , and are the following: from Equation (6) we find that
Figure 3. Numerical results as a function of α for and . (a) α-β parameter space divided into stable and unstable regions. The critical curve is formed by an upper β2-branch and a lower β1-branch; (b) Values of associated with a dividing curve formed by an upper β1-branch and a lower β2-branch; (c) Stability coefficient, Q, showing that for all on both β-branches (upper , lower ); (d) Limit cycle amplitude for where we have fixed the detuning as off of the Hopf bifurcation. Here we present one MATLAB simulation labeled (*) on (a)-(d) for and .
and from (10) we obtain which gives and using Equations (21) and (20), respectively. Substituting these results into the center manifold analysis we obtain from Equation (60) and from Equation (66). We confirm these results via MATLAB’s ode45.m by numerically simulating the system for the appropriate parameter values to obtain the time course in Figure 3(d) where the numerical amplitude of the oscillation is , in good agreement with our center manifold results. The biological significance of knowing the amplitude and frequency of the oscillation is mainly exemplified by the experimental biochemical measurements in  where the authors confirm an amplitude that varies roughly between concentrations of 2 and 3 with a frequency of about t = 10 (generations).
6. Continuation Analysis
In this section we complete the Hopf bifurcation analysis by computing the bifurcation diagram using the numerical continuation software package AUTO. We start by setting as the starting point and define α as the primary continuation parameter. For this particular set of parameters, the system exhibits stable steady state behavior with as predicted in Figure 2(a). Continuing in the positive α-direction we encounter a Hopf bifurcation at as predicted in Figure 3(a). We then perform a second α-continuation starting from the Hopf point, which constructs the bifurcation diagram presented in Figure 4 (solid) for .
The second bifurcation curve, presented in Figure 4 (dotted), corresponds to the parameter value and it shows how a small increase in the value of the Hill coefficient, n, dramatically changes the dynamic behavior of the solution within the same α-β parameter region. The curve shows that for the system will always exhibit periodic solutions, regardless of the value of β. To confirm our numerical continuation results, we used MATLAB’s built-in function ode45.m to compute the three time course simulations presented in Figure 4: (left, circled asterisk), (right top, asterisk), and (right bottom, asterisk). These plots exhibit the correct change in dynamic behavior along when and . Other interesting dynamic behavior can be studied using our numerical continuation results, for example, increasing will shift both bifurcation curves to the right (not presented here). This shift accounts for the effects that “leakiness’’  might have on the system dynamics, because it increases the minimum α for which oscillations are possible while keeping β unchanged. In the following section we present our conclusions and discussions on the biological interpretation of these computational results.
Figure 4. Bifurcation diagram for the repressilator model when . Hopf bifurcation curves computed with AUTO for (solid) and (dotted). Presented here are three MATLAB simulations showing stable behavior for (left, circled asterisk), periodic solutions for (right top, asterisk), and stable behavior (right bottom, asterisk).
This work provides a Hopf bifurcation study of the repressilator model. The linear analysis of the model yields analytical expressions for the steady state solution and critical parameter values. These were used to categorize stability diagrams separating α-β parameter regions where the bifurcation occurs. Setting our parameters close to their critical values allows the system to be close to the bifurcation point, which provides the set up to carry out a center manifold reduction on the system. The final outcome of the center manifold reduction allowed us to find closed form approximate expressions for the amplitude of the limit cycle born at the Hopf, frequency of the oscillation, and the stability coefficient. These expressions, along with our linear analysis results, are then used in a parameter study to construct a more comprehensive picture of the system’s dynamic behavior. Finally, we confirmed our theoretical results in two ways: 1) by numerically simulating appropriate time courses via MATLAB’s built-in function ode45.m and 2) by numerical continuation of the full nonlinear system using AUTO.
From a biological perspective, the two main parameters (α and β) have “opposite” effects on the system. The parameter α roughly represents the maximum production or transcription rate of mRNA in the absence of repression
[1, 2]. Based on the mathematical structure of the Hill term, , we can
see that an increase in protein concentration, , makes the term smaller and
thus decreases its influence on the production rate of mRNA, . On the
other hand, the parameter β roughly represents the ratio between mRNA and protein degradation. Thus, if β = (degradation of protein)/(degradation of mRNA) then increasing β means that either degradation of protein becomes faster and/or degradation of mRNA becomes slower. These rough descriptions of α and β explain how their associated opposite effects give rise to negative feedback, which is an essential feature of any biological oscillator. These conclusions can be verified with our computational results presented in Figure 3 and Figure 4. By considering fixed degradation rates for protein and mRNA (i.e. β = constant), if repression is small (i.e. α small) then there is not enough opposing force to counterbalance the cell’s natural tendency to stay at equilibrium. The latter is represented in Figure 5(a) where the system “behaves” as a network with positive feedback (due to small repression) with one stable steady state. This is also confirmed in Figure 4 where the α-β parameter region to the left of the bifurcation curves shows stable state behavior. A numerical simulation for and confirms the system’s tendency to reach equilibrium regardless of the cell’s initial conditions.
Increasing α for a fixed β creates an opposing reaction, which will offset the cell’s tendency for equilibrium. At the critical value, the system will switch from positive to negative feedback as represented in Figure 5(b). The numerical
Figure 5. Dynamic behavior of the repressilator for different parameter values. (a) Small repression (α small) yields a system that behaves as a network with positive feedback and exhibits a stable solution for fixed β. Numerical simulation for (bottom left); (b) System exhibits oscillations when repression is strong enough ( ) producing a network with negative feedback and nonlinear interactions. Numerical simulation for and (bottom right).
simulation for high α confirms the cyclic clockwise repression effect in the network’s protein-mRNA concentration dynamics. In Figure 4 this occurs as we cross the bifurcation curves and in Figure 3 as we cross the β-branches, which confirms that for any β there exists a threshold on α for the system to exhibit oscillations. From Figure 4 we can deduce that this threshold becomes smaller as n becomes larger, as was exemplified when we compared the and bifurcation curves. Furthermore, this shows how a small increase in the Hill coefficient, n, dramatically changes the biological significance of the bifurcation diagram. The latter is due to the nonlinear effects of the Hill term, giving α a stronger influence over the system’s switch from positive to negative feedback.
 Garcia-Ojalvo, J., Elowitz, M.B. and Strogatz, S.H. (2004) Modeling a Synthetic Multicellular Clock: Repressilators Coupled by Quorum Sensing. PNAS, 101, 10955-10960.
 Pett, J.P., Korencic, A., Wesener, F., Kramer, A. and Herzel, H. (2016) Feedback Loops of the Mammalian Circadian Clock Constitute Repressilator. PLOS Computational Biology, 12, e1005266.
 Buse, O., Kuznetsov, A. and Perez, R.A. (2009) Existence of Limit Cycles in the Repressilator Equations. International Journal of Bifurcation and Chaos, 19, 4097-4106.