Hopf Bifurcation Analysis of the Repressilator Model

The repressilator is a genetic network that exhibits oscillations. The network is formed of three genes, each of which represses each other cyclically, creating a negative feedback loop with nonlinear interactions. In this work we present a computational bifurcation analysis of the mathematical model of the repressilator. We show that the steady state undergoes a transition from stable to unstable giving rise to a stable limit-cycle in a Hopf bifurcation. The nonlinear analysis involves a center manifold reduction on the six-dimensional system, which yields closed form expressions for the frequency and amplitude of the oscillation born at the Hopf. A parameter study then shows how the dynamics of the system are influenced for different parameter values and their associated biological significance.


Introduction
The repressilator is an artificial synthetic gene network created and named by Elowitz and Leibler [1].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 [2] [3].The network exhibits negative feedback and nonlinear interactions, both of which are important features of dynamical systems exhibiting oscillations [4] [5] [6].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 nuc-A.Verdugo DOI: 10.4236/ajcm.2018.82011138 American Journal of Computational Mathematics leus 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 [3].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 [1] [2]).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 [1] is given by the following six coupled first order ordinary differential equations (ODEs) ( ) where 1, 2, 3 i = and i j are defined as 1 3 j = , 2 1 j = , and 3 2 j = .Here α represents dimensionless transcription rate in the absence of repressor, 0 α 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 [1]).
The biological and experimental applications of this model are well documented [1] [7] [8].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 [1].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 [1] oscillations were found to be somewhat noisy and erratic.However, recent new measurements by Potvin-Trottier et al. [8] 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 American Journal of Computational Mathematics 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.

Linear Stability Analysis of the Equilibrium Solution
We start our analysis by showing that Equations ( 1  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 2 n ≥ .This shows that the system has at least one positive biologically significant equilibrium solution when 2 n ≥ .Simulations and plots are provided in Section 5.
To find the stability of the steady state, * , we define i ξ and i η to be deviations from equilibrium these into Equations ( 1) and ( 2) results in the following nonlinear system ( ) ( ) ( ) where 1, 2,3 i = and 1 3 j = , 2 1 j = , and 3 2 j = .Expanding for small values of j η , Equation ( 7) becomes where the Taylor coefficients A, 2 K , and 3 K 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 = ± ) the system (13) and ( 14) will exhibit solutions of the form ( ) ( ) ( ) ( ) where i A and i B are the amplitudes of the ( ) η oscillations, and where i φ is a phase angle.As we add a small detuning off of the critical value, cr β β = + ∆ , the nonlinear system (1) and ( 2) is expected to exhibit peri- odic 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 ( ) we compute and plot the associated stability diagrams for some parameter values.

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 where A is given by Equation (10).Substituting condition (17) into (16) for More details presented in Section 5.
To study the nonlinear part of the center manifold approximation we start by rewriting ( ) By the Center Manifold Theorem [9] 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 6 x ∈  of (22) as where y ∈  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 , y p x = , then it will be the projection of x onto the center subspace.This gives the following expression for the time derivative of , , y p x p yq yp w where ( ) The normal form for the local parametrization of the truncated center manifold (neglecting cubic and higher order terms) is given by (see [9] [10]) ( ) ( ) ( ) where , , , H B q q G q p B q q q = − − (47) , , , H B q q G q p B q q q = − − and where we omit the expressions of Equations ( 47)-(49) for brevity.Equation (43) represents the center manifold approximation, ( ) , w y y , which we substitute into (34) to obtain the flow on the 2-dimensional center subspace ( ) , .

H y y T y T y T y y T yy T y
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 Q p C q q q p B q J B q q p B q i I J B q q We refer the reader to [9] [10] for the standard derivation of Equation ( 60).
Substituting the expressions for J, w 0 , 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.

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 cr β β = and 0 i λ ω = ± .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 , 1 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 0 i λ ω = ± .We thus assume the resulting expression is of the form R i λ = ± Ω , where R and Ω are given by ( ) ( ) 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 R i λ = ± Ω .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 1 R and 1 ω 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.

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, ( ) , , , , , p p p p p p , for different parameter values.For our model, a biologically significant range for the Hill coefficient, n, is given by 2 10 n ≤ ≤ (see [1]), which also agrees with the equilibria condition for positive real roots of the polynomial in Equation (6).Thus the steady state solution, * p , 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 * p . Figure 2 shows * p displayed as a function of α for 0 0, . 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 2 n = shows one numerical simulation for 5 α = and 0 1 α = , denoted on the plot with an asterisk (*), where we can see that * 2 p ≈ 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 * p when 0 0 α = and 2 n = 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 American Journal of Computational Mathematics cr β which will give us two branches: ω , 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 2 n = and 0 0 α = .Here we see that 0 Q < for all 4.2426 α > on both β-branches, which shows that the system exhibits a surface of periodic solutions on the center manifold with stable limit cycles [10].
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 0.1 ∆ = .Notice that the amplitude of the oscilla- tion grows larger as α increases, which shows how the limit cycle's size changes as we follow the bifurcation curve β for a some fixed  are the following: from Equation ( 6) we find that

Conclusions
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 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, .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 100 β = and 10 α = 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  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 2.0 n = and 2.1 n = 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.

Figure 1 .
Figure 1.Repressilator network diagram.The first protein, 1 p , represses production of the second protein, 2 p , the second represses production of the third, 3 p , and the third represses production of the first.These repressions occur at the transcriptional level, where mRNA i m is repressed by protein condition for a change in stability and a Hopf bifurcation.
two branches associated with the ± in Equation (20) for cr β .

=
we know that J has a complex conjugate pair of eigenvalues on the imaginary axis, by taking the derivate and equating to Equation (38).The latter yields the following expressions for the parametrization coefficients

Figure 2 .
Figure 2. Equilibrium solution as a function of α for 0 0,1, 2 α = and 2 n = (left) and 5 n = (right).Solution * p is computed from Equation (6) for given α, α 0 , and n.One MATLAB numerical simulation is presented on the 2 n = (left) plot, where 5 α = and a) shows the α-β parameter space divided into stable and unstable regions for 0 0 α = and 2 n = .The di- viding curve is formed by a lower β 1 -branch and an upper β 2 -branch with the two meeting at 4.2426 α = , 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 cr β β = to obtain the critical frequencies, 0

Figure 3 (
d) shows r as a function of α for 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 0 0 α = , 2 n = , and 30 α = increase in protein concentration, j p , makes the term smaller and thus decreases its influence on the production rate of mRNA, (b).The numerical A. Verdugo DOI: 10.4236/ajcm.2018.82011151 American Journal of Computational Mathematics

Figure 5 .
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 100 α β = = (bottom left); (b) System exhibits oscillations when repression is strong enough ( cr α α > ) producing a network with negative feedback and nonlinear interactions.Numerical simulation for 10 α = x as follows