ABSTRACT

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.

Keywords:

Hopf Bifurcation, Repressilator, Center Manifold Analysis

1. 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 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 [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)

$\frac{\text{d}{m}_i}{\text{d}t}=-m_i+\frac{\alpha }{1+p_{j_i}^n}+{\alpha }_0$ (1)

$\frac{\text{d}{p}_i}{\text{d}t}=-\beta \left(p_i-m_i\right)$ (2)

where $i=1,2,3$ and $j_i$ are defined as $j_1=3$ , $j_2=1$ , and $j_3=2$ . Here α represents dimensionless transcription rate in the absence of repressor, ${\alpha }_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

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

$\frac{\text{d}{p}_i}{\text{d}t}=\frac{\text{d}{m}_i}{\text{d}t}=0$ . This gives

$0=-m_i^{*}+\frac{\alpha }{1+{\left(p_j^{*}\right)}^n}+{\alpha }_0$ (3)

$0=p_i^{*}-m_i^{*},$ (4)

where $\left(m_1^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{m}_2^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{m}_3^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{p}_1^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{p}_2^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{p}_3^{\mathrm{<mo>\; *\; </mo>}}\right)$ represents the steady state solution. Substituting (4) into (3) we obtain

$p_i^{*}=\frac{\alpha }{1+{\left(p_j^{*}\right)}^n}+{\alpha }_0,$ (5)

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 $m_i^{*}=p_i^{*}=p^{*}$ . Substituting the latter into (5) gives the following polynomial

${\left(p^{*}\right)}^{n+1}-{\alpha }_0{\left(p^{*}\right)}^n+p^{*}-\left(\alpha +{\alpha }_0\right)=\mathrm{<mi>\; </mi>0},$ (6)

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 $n\ge 2$ . This shows that the system has at least one positive biologically significant equilibrium solution when $n\ge 2$ . Simulations and plots are provided in Section 5.

To find the stability of the steady state, $p_i^{\mathrm{<mo>\; *\; </mo>}}=m_i^{\mathrm{<mo>\; *\; </mo>}}=p^{\mathrm{<mo>\; *\; </mo>}}$ , we define ${\xi }_i$ and ${\eta }_i$ to be deviations from equilibrium ${\xi }_i=m_i-p^{*}$ and ${\eta }_i=p_i-p^{*}$ . Substituting these into Equations (1) and (2) results in the following nonlinear system

$\frac{\text{d}{\xi }_i}{\text{d}t}=-\left({\xi }_i+p^{*}\right)+\frac{\alpha }{1+{\left({\eta }_{j_i}+p^{*}\right)}^n}+{\alpha }_0$ (7)

$\frac{\text{d}{\eta }_i}{\text{d}t}=-\beta \left({\eta }_i-{\xi }_i\right)$ (8)

where $i=1,2,3$ and $j_1=3$ , $j_2=1$ , and $j_3=2$ . Expanding for small values of ${\eta }_j$ , Equation (7) becomes

$\frac{\text{d}{\xi }_i}{\text{d}t}\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}-{\xi }_i+A{\eta }_{j_i}+K_2{\eta }_{j_i}^2+K_3{\eta }_{j_i}^3+\cdots$ (9)

where the Taylor coefficients A, $K_2$ , and $K_3$ are given by

$A=-\frac{\alpha n{\left(p^{*}\right)}^{n-1}}{{\left(1+{\left(p^{*}\right)}^n\right)}^2}$ (10)

$K_2=\frac{\alpha n{\left(p^{*}\right)}^{n-2}\left(n{\left(p^{*}\right)}^n+{\left(p^{*}\right)}^n-n+1\right)}{2{\left(1+{\left(p^{*}\right)}^n\right)}^3}$ (11)

$K_3=-\frac{\alpha n{\left(p^{*}\right)}^{n-3}\left(\left(n+1\right)\left(n+2\right)p^{2n}+4\left(1-n^2\right)p^n+\left(n-2\right)\left(n-1\right)\right)}{6{\left(1+{\left(p^{*}\right)}^n\right)}^4}.$ (12)

Next we analyze the linearized system coming from Equations (8) and (9)

$\frac{\text{d}{\xi }_i}{\text{d}t}=-{\xi }_i+A{\eta }_{j_i}$ (13)

$\frac{\text{d}{\eta }_i}{\text{d}t}=-\beta {\eta }_i+\beta {\xi }_i$ (14)

which were obtained by truncating the nonlinear terms in Equation (9). The Jacobian at the origin for system (13)-(14) is

$J\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}\left[\begin{array}{cccccc}-1& 0& 0& 0& 0& A\\ 0& -1& 0& A& 0& 0\\ 0& 0& -1& 0& A& 0\\ \beta & 0& 0& -\beta & 0& 0\\ 0& \beta & 0& 0& -\beta & 0\\ 0& 0& \beta & 0& 0& -\beta \end{array}\right]$ (15)

and the associated characteristic equation is found by setting $det\left(\lambda I-J\right)=0$ , which gives the following equation

${\left(\lambda +1\right)}^3{\left(\lambda +\beta \right)}^3-A^3{\beta }^3=\mathrm{<mi>\; </mi>0.}$ (16)

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, $\beta ={\beta }_{cr}$ , where the first pair of complex conjugate eigenvalues cross the imaginary axis. Thus substituting $\lambda =\pm i{\omega }_0$ into (16) gives

$\frac{{\left({\beta }_{cr}+1\right)}^2}{{\beta }_{cr}}\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}\frac{3A^2}{4+2A}$ (17)

which is the condition for a change in stability and a Hopf bifurcation. Notice that for $\beta ={\beta }_{cr}$ (i.e. $\lambda =\pm i{\omega }_0$ ) the system (13) and (14) will exhibit solutions of the form

${\xi }_i\left(t\right)=A_{i}cos\left({\omega }_{0}t+{\varphi }_i\right)$ (18)

${\eta }_i\left(t\right)=B_i\cos \left({\omega }_{0}t\right)$ (19)

where $A_i$ and $B_i$ are the amplitudes of the ${\xi }_i\left(t\right)$ and ${\eta }_i\left(t\right)$ oscillations, and where ${\varphi }_i$ is a phase angle. As we add a small detuning off of the critical value, $\beta ={\beta }_{cr}+\Delta$ , 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 $\frac{{\left(\beta +1\right)}^2}{\beta }>\frac{3A^2}{4+2A}$ and the unstable by $\frac{{\left(\beta +1\right)}^2}{\beta }<\frac{3A^2}{4+2A}$ . 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 $\beta ={\beta }_{cr}$ and $\lambda =\pm i{\omega }_0$ . Solving Equation (17) for ${\beta }_{cr}$ we obtain

${\beta }_{cr}\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}\frac{3A^2-4A-8}{4A+8}\pm \frac{A\sqrt{9A^2-24A-48}}{4A+8}$ (20)

where A is given by Equation (10). Substituting condition (17) into (16) for $\beta ={\beta }_{cr}$ , setting $\lambda =\pm i{\omega }_0$ , and solving for ${\omega }_0$ gives

${\omega }_0=\frac{\sqrt{3}{\beta }_{cr}A}{2\left({\beta }_{cr}+1\right)}$ (21)

which consists of two branches associated with the $\pm$ in Equation (20) for ${\beta }_{cr}$ . More details presented in Section 5.

We start the analysis by expressing system (8) and (9) as follows

$\frac{\text{d}x}{\text{d}t}\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}Jx+F\left(x\right)$ (22)

where $x\in {ℝ}^6$ so that $x={\left({\xi }_1\mathrm{,}{\xi }_2\mathrm{,}{\xi }_3\mathrm{,}{\eta }_1\mathrm{,}{\eta }_2\mathrm{,}{\eta }_3\right)}^{\text{T}}$ , J is given by Equation (15), and

$F\left(x\right)=\left[\begin{array}{c}{K}_2{\eta }_3^2+K_3{\eta }_3^3\\ K_2{\eta }_1^2+K_3{\eta }_1^3\\ K_2{\eta }_2^2+K_3{\eta }_2^3\\ 0\\ 0\\ 0\end{array}\right]+\mathcal{O}\left({\Vert x\Vert }^4\right)\mathrm{.}$ (23)

For $\beta ={\beta }_{cr}$ and $\omega ={\omega }_0$ we know that J has a complex conjugate pair of eigenvalues on the imaginary axis, $\lambda =\pm i{\omega }_0$ , and thus we let $q\mathrm{,}p\in {ℂ}^6$ be the associated eigenvectors corresponding to $i{\omega }_0$ and $-i{\omega }_0$ , respectively. These eigenvectors will satisfy the following three conditions

$Jq=i{\omega }_{0}q$ (24)

$J^{*}p=-i{\omega }_{0}p$ (25)

$\langle p,q\rangle =1$ (26)

where $\langle p,q\rangle ={\sum }_{k=1}^6{\overline{p}}_k{q}_k$ is the standard scalar product in ${ℂ}^6$ . This yields

$q={\left(a_1\mathrm{,}{a}_1^2\mathrm{,1,}{a}_1{a}_2\mathrm{,}{a}_1^2{a}_2\mathrm{,}{a}_2\right)}^{\text{T}}$ (27)

$p=k_1{\left({\overline{a}}_1,1,{\overline{a}}_1^2,{\overline{a}}_1{\overline{a}}_3,{\overline{a}}_3,{\overline{a}}_1^2{\overline{a}}_3\right)}^{\text{T}}$ (28)

where

$a_1=\frac{A\beta }{\left(1+i\omega \right)\left(\beta +i\omega \right)},\mathrm{<mi>\; </mi>}{a}_2=\frac{\beta }{\beta +i\omega },\mathrm{<mi>\; </mi>}{a}_3=\frac{1+i\omega }{\beta },\mathrm{<mi>\; </mi>}{k}_1=\frac{1}{3a_1^2\left(1+a_2{a}_3\right)}$ (29)

and where we have chosen the scaling factor $k_1$ so that $\langle p,q\rangle =1$ . This completes the linear part of the analysis.

To study the nonlinear part of the center manifold approximation we start by rewriting $F\left(x\right)$ as follows

$F\left(x\right)=\frac{1}{2}B\left(x\mathrm{,}x\right)+\frac{1}{6}C\left(x\mathrm{,}x\mathrm{,}x\right)+\mathcal{O}\left({\Vert x\Vert }^4\right)$ (30)

where $x={\left(x_1\mathrm{,}{x}_2\mathrm{,}{x}_3\mathrm{,}{x}_4\mathrm{,}{x}_5\mathrm{,}{x}_6\right)}^{\text{T}}$ and the multilinear functions $B\left(x\mathrm{,}y\right)$ and $C\left(x\mathrm{,}y\mathrm{,}z\right)$ are given by the following expressions

$B\left(x\mathrm{,}y\right)=2K_2{\left(x_6{y}_6\mathrm{,}{x}_4{y}_4\mathrm{,}{x}_5{y}_5\mathrm{,0,0,0}\right)}^{\text{T}}$ (31)

$C\left(x,y,z\right)=6K_3{\left(x_6{y}_6{z}_6,x_4{y}_4{z}_4,x_5{y}_5{z}_5,0,0,0\right)}^{\text{T}}.$ (32)

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 $x\in {ℝ}^6$ of (22) as

$x\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}yq+\overline{y}p+w$ (33)

where $y\in {ℂ}^1$ 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=\langle p\mathrm{,}x\rangle$ , then it will be the projection of x onto the center subspace. This gives the following expression for the time derivative of $y=\langle p\mathrm{,}x\rangle =\langle p\mathrm{,}yq+\overline{y}p+w\rangle$

$\frac{\text{d}y}{\text{d}t}=i{\omega }_{0}y+\langle p\mathrm{,}F\rangle$ (34)

where $F=F\left(yq+\overline{y}\overline{q}+w\right)$ and since $y\in {ℂ}^1$ then it may written as $y=y_1+i{y}_2$ so that Equation (34) can be expressed as follows

$\frac{\text{d}{y}_1}{\text{d}t}=-{\omega }_0{y}_2+\mathrm{Re}\langle p,F\rangle$ (35)

$\frac{\text{d}{y}_2}{\text{d}t}={\omega }_0{y}_1+\mathrm{Im}\langle p,F\rangle .$ (36)

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 $y_1$ and $y_2$ we start by solving Equation (33) in terms of w and substitute $y=\langle p\mathrm{,}x\rangle$ to obtain

$w=x-\langle p,x\rangle q-\langle \overline{p},x\rangle \overline{q}.$ (37)

Taking the derivative of the latter and substituting Equation (34) we obtain

$\frac{\text{d}w}{\text{d}t}=Jw+F-\langle p,F\rangle q-\langle \overline{p},F\rangle \overline{q}$ (38)

where

$\langle p,F\rangle =\frac{1}{2}{G}_{20}{y}^2+G_{11}y\overline{y}+\frac{1}{2}{G}_{02}{\overline{y}}^2+y\langle p,B\left(q,w\right)\rangle +\overline{y}\langle p,B\left(\overline{q},w\right)\rangle +\text{h}\text{.o}\text{.t}$ (39)

and

$G_{20}=\langle p,B\left(q,q\right)\rangle =2a_1{a}_2^2{\overline{k}}_1\left(a_1^2+a_1+1\right)\left(a_1^3-a_1^2+1\right)$ (40)

$G_{11}=\langle p,B\left(q,\overline{q}\right)\rangle =2a_1{a}_2{\overline{a}}_2{\overline{k}}_1\left(a_1^3{\overline{a}}_1^2+{\overline{a}}_1+1\right)$ (41)

$G_{02}=\langle p,B\left(\overline{q},\overline{q}\right)\rangle \mathrm{<mi>\; </mi>}=2{\overline{a}}_2^2{\overline{k}}_1\left(a_1^2{\overline{a}}_1^4+{\overline{a}}_1^2+a_1\right).$ (42)

The normal form for the local parametrization of the truncated center manifold (neglecting cubic and higher order terms) is given by (see [9] [10] )

$w\left(y,\overline{y}\right)=\mathrm{<mi>\; </mi>}{w}_{20}{y}^2+w_{11}y\overline{y}+w_{02}{\overline{y}}^2$ (43)

which we compute by taking the derivate and equating to Equation (38). The latter yields the following expressions for the parametrization coefficients

$w_{20}=\frac{1}{2}{\left(2i{\omega }_{0}I-J\right)}^{-1}{H}_{20}$ (44)

$w_{11}=-J^{-1}{H}_{11}$ (45)

$w_{02}=-\frac{1}{2}{\left(2i{\omega }_{0}I+J\right)}^{-1}{H}_{02}$ (46)

where

$H_{20}=B\left(q,q\right)-G_{20}q-\langle \overline{p},B\left(q,q\right)\rangle \overline{q}$ (47)

$H_{11}=B\left(q,\overline{q}\right)-G_{11}q$ (48)

$H_{02}=B\left(\overline{q},\overline{q}\right)-G_{02}q-\langle \overline{p},B\left(\overline{q},\overline{q}\right)\rangle \overline{q}$ (49)

and where we omit the expressions of Equations (47)-(49) for brevity. Equation (43) represents the center manifold approximation, $w\left(y\mathrm{,}\overline{y}\right)$ , which we substitute into (34) to obtain the flow on the 2-dimensional center subspace

$\frac{\text{d}y}{\text{d}t}=i{\omega }_{0}y+T_0{\overline{y}}^2+T_1{y}^3+T_2{y}^2\overline{y}+T_{3}y{\overline{y}}^2+T_4{\overline{y}}^3$ (50)

where

$T_0=\frac{1}{2}{G}_{02}$ (51)

$T_1=\frac{1}{2}\langle p,B\left(q,w_{20}\right)\rangle$ (52)

$T_2=\langle p,B\left(q,w_{11}\right)\rangle +\frac{1}{2}\langle p,B\left(\overline{q},w_{20}\right)\rangle$ (53)

$T_3=\langle p,B\left(\overline{q},w_{11}\right)\rangle +\frac{1}{2}\langle p,B\left(q,w_{02}\right)\rangle$ (54)

$T_4=\frac{1}{2}\langle p,B\left(\overline{q},w_{02}\right)\rangle .$ (55)

Since $y=y_1+i{y}_2$ then Equation (50) can be expressed as follows

$\frac{\text{d}{y}_1}{\text{d}t}=-{\omega }_0{y}_2+\mathrm{Re}\left\{H\left(y_1,y_2\right)\right\}$ (56)

$\frac{\text{d}{y}_2}{\text{d}t}=\mathrm{<mi>\; </mi>}{\omega }_0{y}_1+\mathrm{Im}\left\{H\left(y_1,y_2\right)\right\}$ (57)

where

$H\left(y_1,y_2\right)=T_0{\overline{y}}^2+T_1{y}^3+T_2{y}^2\overline{y}+T_{3}y{\overline{y}}^2+T_4{\overline{y}}^3.$ (58)

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

$\frac{\text{d}r}{\text{d}t}\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}Q{r}^3+\mathcal{O}\left(r^5\right)\mathrm{<mi>\; </mi>},\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}\frac{\text{d}\theta }{\text{d}t}\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}{\omega }_0+\mathcal{O}\left(r^2\right)$ (59)

where the stability coefficient is given as follows

$\begin{array}{l}Q=\frac{1}{2{\omega }_0}\mathrm{Re}[\langle p,C\left(q,q,\overline{q}\right)\rangle -2\langle p,B\left(q,J^{-1}B\left(q,\overline{q}\right)\right)\rangle \\ +\langle p,B\left(\overline{q},{\left(2i{\omega }_{0}I-J\right)}^{-1}B\left(q,q\right)\right)\rangle ].\end{array}$ (60)

We refer the reader to [9] [10] for the standard derivation of Equation (60). Substituting the expressions for J, w₀, 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 $\beta ={\beta }_{cr}$ and $\lambda =\pm i{\omega }_0$ . 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

$\beta ={\beta }_{cr}+\Delta ,\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}\left|\Delta \right|\ll 1$ (61)

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 $\lambda =\pm i{\omega }_0$ . We thus assume the resulting expression is of the form $\lambda =R\pm i\Omega$ , where R and $\Omega$ are given by $R=R_1\Delta +\mathcal{O}\left({\Delta }^2\right)$ and $\Omega ={\omega }_0+{\omega }_1\Delta +\mathcal{O}\left({\Delta }^2\right)$ . Thus Equations (35) and (36) will take the approximate form

$\frac{\text{d}{y}_1}{\text{d}t}=R{y}_1-\Omega y_2+\mathrm{Re}\langle p,F\rangle$ (62)

$\frac{\text{d}{y}_2}{\text{d}t}=R{y}_2+\Omega y_1+\mathrm{Im}\langle p,F\rangle .$ (63)

which can be transformed into polar coordinates (and by means of a near-identity transformation) into the following flow on the center subspace

$\frac{\text{d}r}{\text{d}t}=Rr+Q{r}^3+\mathcal{O}\left(r^5\right)$ (64)

$\frac{\text{d}\theta }{\text{d}t}=\Omega +\mathcal{O}\left(r^2\right)\mathrm{.}$ (65)

Setting Equation (64) to zero and solving for r gives the limit cycle amplitude as

$r^2=-\frac{R}{Q}$ (66)

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

${\xi }_i=A_i{\text{e}}^{\lambda t}$ (67)

${\eta }_i=B_i{\text{e}}^{\lambda t}$ (68)

where $\lambda =R\pm i\Omega$ . Substituting the latter two equations and (61) into (13) and (14), linearizing for small $\Delta$ , and solving for R and $\Omega$ we obtain expressions of the form

$R=R_1\Delta$ (69)

$\Omega ={\omega }_0+{\omega }_1\Delta$ (70)

where $R_1$ and ${\omega }_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.

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, $\left(p^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{p}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{p}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{p}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{p}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{p}^{\mathrm{<mo>\; *\; </mo>}}\right)$ , for different parameter values. For our model, a biologically significant range for the Hill coefficient, n, is given by $2\le n\le 10$ (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^{\mathrm{<mo>\; *\; </mo>}}$ , is determined by solving Equation (6) for given values of α, α₀, and n, where we note that β does not play a role in the values of $p^{\mathrm{<mo>\; *\; </mo>}}$ . Figure 2 shows $p^{\mathrm{<mo>\; *\; </mo>}}$ displayed as a function of α for ${\alpha }_0=0,1,2$ , and $n=2,5$ . 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 $n=2$ shows one numerical simulation for $\alpha =5$ and ${\alpha }_0=1$ , denoted on the plot with an asterisk (*), where we can see that $p^{\mathrm{<mo>\; *\; </mo>}}\approx 2$ 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^{\mathrm{<mo>\; *\; </mo>}}$ when ${\alpha }_0=0$ and $n=2$ to obtain the real solution

$p^{*}\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}{\left(\frac{\sqrt{27{\alpha }^2+4}}{6\sqrt{3}}+\frac{\alpha }{2}\right)}^{\frac{1}{3}}-\frac{1}{3}{\left(\frac{\sqrt{27{\alpha }^2+4}}{6\sqrt{3}}+\frac{\alpha }{2}\right)}^{-\frac{1}{3}}$ (71)

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

${\beta }_{cr}$ which will give us two branches: ${\beta }_1=\frac{3A^2-4A-8}{4A+8}+\frac{A\sqrt{9A^2-24A-48}}{4A+8}$ and ${\beta }_2=\frac{3A^2-4A-8}{4A+8}-\frac{A\sqrt{9A^2-24A-48}}{4A+8}$ . Figure 3(a) shows the α-β parameter

space divided into stable and unstable regions for ${\alpha }_0=0$ and $n=2$ . The dividing curve is formed by a lower β₁-branch and an upper β₂-branch with the two meeting at $\alpha =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 $\beta ={\beta }_{cr}$ to obtain the critical frequencies, ${\omega }_0$ , 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 $n=2$ and ${\alpha }_0=0$ . Here we see that $Q<0$ for all $\alpha >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 $\sqrt{-\frac{R_1}{Q}}\cdot \sqrt{\Delta }$ . Figure 3(d) shows r as a function of α for

${\alpha }_0=0$ , $n=2$ , $\beta ={\beta }_1$ , and $\Delta =0.1$ . 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 $\left|\Delta \right|\ll {\beta }_{cr}$ . For $\beta ={\beta }_2$ 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 ${\alpha }_0=0$ , $n=2$ , and $\alpha =30$ are the following: from Equation (6) we find that

$p^{*}=3$ and from (10) we obtain $A=-1.8$ which gives ${\omega }_0=-1.492$ and ${\beta }_{cr}=22.255$ using Equations (21) and (20), respectively. Substituting these results into the center manifold analysis we obtain $Q=-0.151$ from Equation (60) and $r=2.996$ 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 $r=\left(7.384-1.395\right)/2=2.995$ , 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 [8] 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 $\left(\beta \mathrm{,}\alpha \mathrm{,}n\mathrm{,}{\alpha }_0\right)=\left(10,1,2,0\right)$ 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 $p^{\mathrm{<mo>\; *\; </mo>}}\approx 0.682$ as predicted in Figure 2(a). Continuing in the positive α-direction we encounter a Hopf bifurcation at $\alpha \approx 12.929$ 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 $n=2$ .

The second bifurcation curve, presented in Figure 4 (dotted), corresponds to the parameter value $n=2.1$ 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 $n=2.1$ curve shows that for $\alpha >100$ 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(\beta \mathrm{,}\alpha \mathrm{,}n\right)=\left(100,10,2.1\right)$ (left, circled asterisk), $\left(\beta \mathrm{,}\alpha \mathrm{,}n\right)=\left(100,100,2.1\right)$ (right top, asterisk), and $\left(\beta \mathrm{,}\alpha \mathrm{,}n\right)=\left(100,100,2.0\right)$ (right bottom, asterisk). These plots exhibit the correct change in dynamic behavior along $\beta =100$ when $n=2.0$ and $n=2.1$ . Other interesting dynamic behavior can be studied using our numerical continuation results, for example, increasing ${\alpha }_0>0$ will shift both bifurcation curves to the right (not presented here). This shift accounts for the effects that “leakiness’’ [1] 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.

7. 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 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, $\frac{\alpha }{1+p_j^n}$ , we can

see that an increase in protein concentration, $p_j$ , makes the term smaller and

thus decreases its influence on the production rate of mRNA, $\frac{\text{d}{m}_i}{\text{d}t}$ . 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 $\beta =100$ and $\alpha =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(b). The numerical

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 $n=2.0$ and $n=2.1$ 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.

Cite this paper

Verdugo, A. (2018) Hopf Bifurcation Analysis of the Repressilator Model. American Journal of Computational Mathematics, 8, 137-152. https://doi.org/10.4236/ajcm.2018.82011

References