Bifurcation Analysis and Probabilistic Energy Landscapes of Two-Component Genetic Network

Genetic oscillators play an important role in cell physiology. In this paper, we investigate the global kinetics of an activator-repressor genetic oscillator through the bifurcation analysis and probabilistic energy landscape. The results reveal that the bifurcation can induce the rich dynamics with the variation of the parameters, including monostable state, oscillations, and bistable state. These dynamics are further investigated by the probabilistic energy landscapes in which the force from the negative gradient of the energy landscapes attract the system down to the oscillation path, and the flux is the driving force of the oscillations. Besides, the probabilistic energy landscapes show that the dynamics exhibit distinct transition features when the system is close to the critical points of the oscillation region. Based on the theory of energy landscape, we also explore how the system parameters and noise affect the stability of the activator-repressor oscillations through the physical methods of barrier heights and energy dissipation. The simulation results show that the system parameters and noise intensity have a negative relation with the stability of the oscillations, and the control parameter can play a positive role in maintaining the stability of oscillations. Our methods are quite general and can be utilized to explore the potential global dynamics of more complicated gene regulatory systems.


I. INTRODUCTION
It is no coincidence that the gene regulatory network (GRN) is of such general importance: they determine the fate of cells by gene expression mechanism. Recently, there has been a growing focus on gene regulatory networks due to its fundamental role in living organisms [1]- [4]. In particular, oscillation dynamics play an important role in many aspects such as cell division, cell movement and cell differentiation [5], [6]. In the GRNs, different gene motifs represent different biological function, such as promotion (positive) and repression (negative). One gene unit controls the expression of others, and different control relationships form gene regulatory networks.
The associate editor coordinating the review of this manuscript and approving it for publication was Jianquan Lu .
Thus, the feedback loops in GRN are the basic architectures to induce oscillations [7]. The cell-cycle network of Xenopus laevis embryos, which contains a negative feedback loop between CcB and APC, and a positive feedback loop between CcB and Cdc25, leading to a sustained oscillation [8]. The p53-Mdm2 oscillatory network with multiple positive and negative feedback loops has been widely studied [9]- [12]. Besides, the circadian clock with negative feedback loops forms an oscillation and the period of which is close to 24 h [13]. Due to these oscillations are often based on the transcriptional or translational feedback loops, we regard them as genetic oscillators. More information about the genetic oscillators can refer to the review article [6].
It is not sufficient to rely on biological experiments to study the gene regulation networks. Therefore, mathematical models must be introduced to the gene regulatory networks to study their dynamic behaviors. In general, the feedback loops can be formulated as a mathematical description such as the nonlinear ODE equations with Hill function or DDE equations [14], [15]. The theoretical analysis of ODEs can help us find the steady-state solutions of the network and understand the dynamics of genetic oscillators more intuitively [16]. However, in a cell, the number of molecules is sparse. So, the fluctuations of intrinsic or external noise can be important for dynamics [17]- [19]. Besides, some studies have pointed out that the process of promoter and DNA binding at the transcription level is random [6]. Thus, instead of the deterministic ODEs, master equations or stochastic differential equations should be introduced to describe the corresponding regulate process [20], [21].
Bifurcation analysis can be introduced to explore the local properties of the deterministic system. However, how to explore the global dynamics of the system is still a challenge. Motivated by this, we will study the global properties of the genetic oscillator through the methods of the probabilistic energy landscape and flux [22]- [25]. Furthermore, once we get the energy landscape of the genetic oscillator, we can measure the robustness of the oscillation from the physical point of view.
The genetic oscillator is a nonequilibrium network, which can exchange the energies and information with its surroundings. Thus, the global dynamics of the genetic oscillator are determined by both the probabilistic energy landscape and flux [26]- [28]. The energy landscape is closely related to the steady-state probability which can be obtained by calculating the corresponding stochastic differential equation. In the probabilistic energy landscapes of a nonlinear system, different energy weights represent different states. For example, the stable state corresponds to higher steady state probability and lower energy weight.
It is significant to study the dynamics of the twocomponent genetic oscillator due to its fundamental role in the GRNs [29]. Here, according to the genetic oscillator model proposed by Guantes et al. in 2006 [5], we will investigate the global kinetics and the robustness of the activator-repressor oscillator through the methods of bifurcation analysis and probabilistic energy landscapes.
The rest of this paper is organized as follows. In Section II, we give the circuits of the two-component genetic oscillator and its mathematical representation. And then, the global bifurcation analyses of the two-dimensional system are performed. In Section III, we firstly give the theoretical derivation of the relationship between energy landscape and steady state probability. Then, the numerical simulation results of the probabilistic energy landscape are given to explore the global dynamics of the corresponding stochastic system. Finally, we explore the effects of system parameters and noise intensity on the stability of the oscillations through the physical methods of barrier height and energy dissipation. We conclude this paper in Section IV.

II. MODEL AND BIFURCATION ANALYSIS A. TWO-COMPONENT GENETIC OSCILLATOR MODULE
In this section, we will introduce the construction principle of the two-component genetic oscillator model, which is proposed by Guantes and his colleagues [5]. Fig. 1 shows the key features of the minimal network architecture and gene transcription circuit of activator x 1 and repressor x 2 . The two-component genetic oscillator model can be described by a set of ODEs: where x i (i = 1, 2) are the concentrations of the activator A and the repressor R, respectively. a represents the ratio of decline rates between activator and repressor. k denotes the protein accumulation through the interaction of activator and the promoter, r represents the reaction intensity of repressor, m represent the effective basal rate of activator and φ is the ratio of effective basal rate between activator and repressor, i.e., n = φm. For convenience, we named parameters k, m, r and φ as the system parameter, and parameter a as the control parameter. Please refer to [5] for detailed information about this model.

B. BIFURCATION ANALYSIS
Hopf bifurcation is a relatively common bifurcation phenomenon which can induce oscillations in nonlinear dynamic systems. In this section, we give the theoretical analysis of Hopf bifurcation. Suppose (A 0 , R 0 ) is an equilibrium point of the system (1), let x(t) = x 1 (t) − A 0 , y(t) = x 2 (t) − R 0 , Then, linearizing system (1) at its equilibrium point (A 0 , R 0 ) yields: (2) VOLUME 8, 2020 where The matrix J have the form as follows: then, we get the characteristic equation of the system (2) at the equilibrium point (0,0): that is, . Then, we get the theoretical results as follows: Theorem When a two-dimensional system satisfies the following conditions, the Hopf bifurcation occurs.
When a = a 0 , system (1) have a equilibrium point (A 0 , R 0 ). The First Lyapunov coefficient a 0 has the following form [30]: To explore the effects of the parameters on the dynamics of the genetic oscillator, we performed the bifurcation analysis with the help of MATCONT and XPPAUT software. Numerical simulation results are presented in Fig. 2-Fig. 7. Fig. 2 illustrates bifurcation diagrams under different parameters. From Fig. 2(a), one can see that the blue line divides the(a, k) plane into two parts: the monostable region and the oscillation region. The blue line is divided into three parts by the generalized Hopf(GH ) bifurcation point and the Bogdanov-Takens(BT ) bifurcation point. The bifurcation point GH is the critical point between the curve for supercritical Hopf bifurcation points and the curve for the subcritical Hopf bifurcation points. The bifurcation point BT is the critical point between the curve for the Hopf bifurcation (supercritical Hopf bifurcation or subcritical Hopf bifurcation) points and the curve for the fold bifurcation points. The first part is the curve above the bifurcation point BT and represents the fold bifurcation points, the second part is the , and (g)k = 6.0. The red line represents stable equilibrium, and the black dashed line is the unstable equilibrium. Full green circles represent the maxima and minima of the stable limit cycles, and blue open circles are the maxima and minima of the unstable limit cycles. The black solid circles are the bifurcation points. The labels suph and suph i represent the supercritical Hopf bifurcation points, subh represents the subcritical Hopf bifurcation point, f i represents the fold bifurcation point, hc represents the homoclinic bifurcation point, lpc represents the fold bifurcation of the limit cycles, and snic represents the saddle-node homoclinic bifurcation point.
curve between the bifurcation point BT and the bifurcation point GH , which is the curve for the subcritical Hopf bifurcation points, and the third part is the rest of the blue line, which is the curve for supercritical Hopf bifurcation points. The curve for the fold bifurcation points is labeled F, and the curves for the subcritical Hopf bifurcation points and supercritical Hopf bifurcation points are labeled H . The lowest point of curve is located at (a, k) = (12.90, 2.46), the GH point is located at (a, k) = (10.98, 3.31), and the BT point is located at (a, k) = (10.27, 5.85). Thus we can conclude that there is the monostable region with a single stable equilibrium point when k < 2.46; When k increase, the system enters the oscillating region with two supercritical Hopf bifurcation points when 2.46 < k < 3.31, indicating the stable limit cycles will be bifurcated from the corresponding bifurcation points. When 3.31 < k < 5.85, the system will go to another region with a subcritical Hopf bifurcation point and a supercritical Hopf bifurcation point, indicating the emerges of both  the unstable limit cycles and the stable limit cycles. When k > 5.85, there is the oscillation region with a supercritical Hopf bifurcation point. The corresponding codimension-one bifurcation diagrams under different system parameter k are described below:   (i) When k = 2.0( Fig. 2(b)), there is a stable fixed point in the network when the system parameter k belongs to the region (0, 2.46). And with the increase of the bifurcation VOLUME 8, 2020 parameter a, the concentrations of A changed from high to low.
(ii) When k = 3.0( Fig. 2(c)), which belongs to the region (2.46, 3.31), the two supercritical Hopf bifurcation points suph 1 and suph 2 appear at a = 11.249610 and a = 14.766415, respectively. Thus the curve of a single stable steady state is divided into three parts: stable equilibria points with the higher concentration of A, unstable equilibria points around by the stable limit cycles, and stable equilibria points with the lower concentration of A.
(iii) When k increases to k = 4.0( Fig. 2(d)), which belongs to the region (3.31, 5.85), a subcritical Hopf bifurcation point sub and a supercritical Hopf bifurcation point sup appear at a = 10.629802 and a = 15.949615, respectively. As the increase of the bifurcation parameter a, the unstable limit cycles are bifurcated from the subcritical Hopf bifurcation point sub and become stable at the fold bifurcation of the limit cycles lpc. The stable limit cycles are vanishing at the supercritical Hopf bifurcation point sup.
(iiii) For large k values k = 5.0( Fig. 2(e)(f)), which also belongs to the region (3.31, 5.85), a subcritical Hopf bifurcation point sub and a supercritical Hopf bifurcation point sup appear at a = 10.383334 and a = 16.767119. Besides, the two fold bifurcation f 1 and f 2 appear such that the unstable branch is divide into three parts. Moreover, instead of the lpc, the unstable limit cycles ended with the homoclinic bifurcation point hc due to the appears of fold bifurcation point.
(iiiii) When further increase k to k = 6.0 ( Fig. 2(g)), a supercritical Hopf bifurcation point sup appears at a = 17.442166, and a saddle-node homoclinic bifurcation point snic arises due to the disappearance of the subcritical Hopf bifurcation. With the increases of a, the stable limit cycle is developing until it meets the sup point.
Next, we take the subcritical Hopf bifurcation in case (iiii) as an example to verify the Hopf bifurcation theory. Fixed the parameter values a = 10.383334, k = 5.0, m = 1.58, r = 1.0, φ = 0.05, respectively. And the corresponding equilibrium point is (A 0 , R 0 ) = (3.498835, 3.853638). Substituted the parameter values in equation (6), we have: Hence, Hopf bifurcation conditions are satisfied. In addition, according to (7) we obtain the First Lyapunov coefficient a 0 = 2.979859 > 0. Therefore, according to the Hopf bifurcation theorem, the supercritical Hopf bifurcation occurs when a = 10.383334. Fig. 3 performs the bifurcation analysis of the system parameter m and the control parameter a. From Fig. 3, one can see that the system parameter k can induce similar bifurcation diagrams with system parameter m. Thus we omit the description of the bifurcation diagram in Figure 3.
In Fig. 4 and Fig. 5, we plot the bifurcation diagrams of system parameter r and φ with respect to the control parameter a, respectively. Compared Fig. 4 with Fig. 5, one can see that the system parameter r and φ can induce similar bifurcation diagrams. Here, we take Fig. 4 an example to illustrate the bifurcation phenomena. From Fig. 4(a), one can see that the curve on the left side is divided into three parts by two BT points. The three parts (from top to bottom) are the curve for the supercritical Hopf bifurcation points, the curve for the fold bifurcation points, and the curve for the subcritical Hopf bifurcation points, respectively. The curve on the right side is for the supercritical Hopf bifurcation points. The two TB points are located at (a, r) = (6.76, 2.37) and (a, r) = (2.38, 18.55), respectively. There will be a subcritical Hopf bifurcation point and a supercritical Hopf bifurcation point in the region 0 < r < 2.37, a supercritical Hopf bifurcation point in the region 2.37 < r < 18.55, and two supercritical Hopf bifurcation points in the region 18.55 < r < 27.86. The corresponding codimension-one bifurcation diagrams under different system parameter r are described below: (i) When r = 1.5 ( Fig. 4(b)), which belongs to the region (0, 2.37), a subcritical Hopf bifurcation point sub and a supercritical Hopf bifurcation point sup appear at a = 13.533403 and a = 17.07, respectively. The unstable limit cycle is generate from the subcritical Hopf bifurcation point sub and ended with the homoclinic bifurcation point hc due to the appears of fold bifurcation points f 1 and f 2 . With the increases of a, the stable limit cycle is developing until it meets the sup point at a = 17.07.
(ii) With increases r to r = 10.0 ( Fig. 4(c)), which belongs to the region (2.37, 18.55), the supercritical Hopf bifurcation point is disappear with the rise of a saddle-node homoclinic bifurcation point snic and a fold bifurcation f 1 . A stable limit cycle is bifurcated from the subcritical Hopf bifurcation point sub and extends left to the point snic.
(iii) Further increases r to k = 25 ( Fig. 4(d)), which belongs to the region (18.55, 27.86), the two supercritical Hopf bifurcation points suph 1 and suph 2 appear at the lower a values a = 2.075413 and a = 2.443606. The saddle-node homoclinic bifurcation point snic disappears due to the maxima of the stable limit cycle becomes lower than the fold bifurcation f 1 .
Besides, in Fig. 6 and Fig. 7, we also plot the bistable bifurcation diagrams of system parameter m and k, respectively. In Fig. 6(a) and Fig. 7(a), we plot the codimension-two bifurcation diagram of system parameters m and k with respect to the control parameter a when φ is relatively large (φ = 0.5). From Fig. 6(a) and Fig. 7(a), one can see that the two curves for the fold bifurcation points collide and disappear at the cusp bifurcation point (CP) with the variation of system parameters m and k, respectively. The bifurcation point CP is the critical point between the two curves for the fold bifurcation points. The blue line divides the plane into two parts: the monostable region and the bistable region. In Fig. 6(a), the CP point is located at (a, m) = (1.04, 1.58). In Fig. 2(f), the CP point is located at (a, k) = (1.05, 4.48). Thus we conclude that the bistability appears in the system when m ≥ 1.58 and k ≥ 4.48 when φ is relatively large. Here, we describe the corresponding bistable codimension-one bifurcation in Fig. 6.
(i) When m = 1.0 ( Fig. 6(b)), which belongs to the region of monostability, there is a stable fixed point in the system. And with the increase of the bifurcation parameter a, the concentrations of activator changed from high to low.
(ii) When increased m to m = 3.0 (Fig. 6(c)), which belongs to the region of bistability, the equilibria exhibit the bistability of a Z-shaped bifurcation diagram. The Z-shaped curve is divided into three parts by the two fold bifurcation points f 1 and f 2 : the stable branch with higher concentrates of x 1 , the stable branch with lower concentrates of x 1 , and the unstable branch between the two stable branches. The abscissa values of the two fold bifurcation points f 1 and f 2 are a = 0.8115 and a = 1.009, respectively.
(iii) Further increases m to m = 15 ( Fig. 6(d)), the abscissa values of the two fold bifurcation points f 1 and f 2 become a = 0.3858 and a = 1.0. This indicates that the bistability region is increasing with the increase of system parameter m.
To have a clear insight into the codimension-one bifurcation diagrams, we consider the corresponding phase analysis and time evolution analysis are as shown in Fig. 8 and Fig. 9, respectively. Here, we take the bifurcation diagrams Fig. 2(e) and Fig. 6(c) as an example. From Fig. 8(a) one can see that there is a stable equilibrium when the control parameter a = 10.0, which corresponds to the stable state in the bifurcation diagrams Fig. 2(e). When a = 10.34( Fig. 8(b)), the system has three equilibria(two unstable and one stable). In addition to a stable equilibrium, a stable limit cycle begins to appear when a = 10.37 (Fig. 8(c)). When a = 10.38 (Fig. 8(d)), an unstable limit cycle separates the stable equilibrium from the stable limit cycle. When a = 10.39 (Fig. 8(e)), there are three unstable equilibria around by a stable limit cycle. When a = 10.40 (Fig. 8(f)), there is a stable equilibrium point around by a stable limit cycle. Fig. 8(g) shows that there is one unstable equilibrium between the two stable equilibria, which corresponds to the bistable bifurcation diagram in Fig. 6(c). Fig. 9 illustrate the time evolution of the system under different bifurcation parameters a. From Fig. 9, one can see that with the increase of parameter a, the system gradually goes from the monostability to oscillations. Further increase parameter a, The system returns to monostability. Also, we find that a small parameter value a corresponds to a high stable concentration of activator and a big parameter value a corresponds to a low stable concentration of activator. Fig. 10 shows the period and amplitude of the oscillations under different control parameter a.

III. PROBABILISTIC ENERGY LANDSCAPE A. THEORETICAL ANALYSIS
In this section, we give the theoretical analysis of the relationship between probabilistic energy landscape and steady-state probability distribution [31]. In general, we cannot write the potential function of a nonequilibrium system directly. Here we will directly map out the probabilistic energy landscape U by solving the corresponding FokkerĺCPlanck equation to get the steady-state probability P.
As mentioned in the introduction, the fluctuations of intrinsic or external noises are inevitable in cellular environments. To analyze the genetic oscillators more accurately, the deterministic system (1) should to be changed to the corresponding Langevin equations: that is,ẋ VOLUME 8, 2020  where ζ i , (i = 1, 2) is the noise term. Here, we assume that the random term is Gaussian white noise which can be described as: Here, δ(t) represents the Dirac delta function and D ij is the random matrix which can measure the intensity of the random force.
According to [31], there is a unique transformation of the system equation (9): Here S(x) and A(x) are represent the semi-positive symmetric matrix and anti-symmetric matrix, respectively. The scalar function U (x) is the potential function and we named it the probabilistic energy landscape due to it has a direct relationship with the steady-state probability distribution. ξ is the stochastic force. It's easy to see that the semi-positive symmetric matrix is satisfies:ẋ T S(x)ẋ ≥ 0, that is, matrix S(x) is 'dissipative'; the anti-symmetric matrix is satisfies: is 'non-dissipative'. Therefore, one can conclude that matrix S(x) represents the dissipate force of the dynamic system, matrix A(x) represents the transverse force of the dynamic system.
In addition, we find that there are two independent terms F and ζ in system (9) while the system (10) has four independent terms S, A, U and ξ . Hence, we should add some additional constraints to ensure that system (10) has a unique form. One constraint condition is the stochastic-dissipation relation [31]: Next let's prove the uniqueness of equation (10). Substitute equation (9) into equation (10) and eliminate theẋ, we have: Since the random term and the state variables are independent of each other, we can decompose equation (11) into the following two forms: one is that for the deterministic force: and the other is that for the stochastic force: From equation (13) and equation (14), one can see that there is a 'rotation' from the force F to the probabilistic energy landscape U and a 'rotation' from the stochastic force ζ to the stochastic force ξ at every point of x. And then, we substitute equations (14), we have: This indicates a generalized Einstein relation to transverse matrix A.
Next we need to construct an auxiliary function Here''-1'' represents the inverse of S(x)+A(x). Thus we have Based on the properties of the probabilistic energy landscape U : (13) becomes: A(x), we obtain:
If the auxiliary function G is uniquely determined, according to the relationships (S(x) , the probabilistic energy landscape U , matrix S(x) and matrix A(x) can be only determined: Next, we will analyze the relationship between the probabilistic energy landscape U and the steady-state probability distribution P ss in the non-equilibrium genetic oscillator. The Fokker-Planck equation can be written as: For the sake of simplicity, we regard the diffusion coefficient matrix as a constant (D 11 = D 22 = D, D 12 = D 21 = 0). Rewrite the equation (20) as the following form: which indicates a conservation law of probability. The flux vector J can be defined as: If the steady-state probability distribution exists, we have ∂P ss ∂t = 0 and ∇ · J (x, t) = 0. There are two different cases. One case is in the equilibrium system, which corresponds to the condition of detailed equilibrium, that is At the steady-state, the flux becomes Using equations (23) and (24), we have which indicates the system can be written as the negative gradient of a potential function under the condition of the detailed balance. The other case is that in a non-equilibrium system, the detailed equilibrium is broken, that is Using equations (24) and (26), we obtain: The above equation implies that the driving force F in a non-equilibrium system can be decompose into two parts: the probabilistic energy landscape −D · ∂ ∂x U , and the flux J ss /P ss . Hence, the relationship between probabilistic energy landscape U and steady-state probability distribution P ss is: U = −lnP ss . This implies that the probabilistic energy landscape has a negative correlation with steady-state probability distribution, that is, the lower landscape corresponds to higher probability.
In summary, for an equilibrium system, the dynamics can be completely determined by the probabilistic energy landscape U ; for a non-equilibrium open system (genetic oscillator, in this case), the global kinetics can be determined by the probabilistic energy landscape U and the flux J . In other words, once we obtain the steady-state probability distribution P ss , the global stationary dynamics of the system can be characterized by mapping out the corresponding probabilistic energy landscape U and the flux J .

B. RESULTS
As discussed above, the steady-state probability distribution can be obtained by calculating the Fokker-Planck equation. However, it is a challenging work to calculate the FP equations directly. Here, we impose the basic finite difference method to calculate the Fokker-Planck equation with the help of the Virtual Cell [32]. The drawing region is a square with 0 ≤ x 1 ≤ 10 and 0 ≤ x 2 ≤ 10, the initial condition is P = 0.01 and the time step is 0.05s. The Neumann boundary condition is used in the simulation process. After we obtain the steady-state probability distribution P ss , we can obtain the probabilistic energy landscape U by the relationship U = −lnP ss . The numerical results of the probabilistic energy landscapes are shown in Fig. 11.
(I ) When a = 10, which is close to the critical pointa = 10.383334), the system is in the monostable region. From Fig. 11(a), one can see that the landscape shows there is a global minimum value which corresponds to the stable equilibrium point with a higher concentration under the small noise intensity (D = 1.0 × 10 −8 ). When the noise intensity increases (D = 0.05, Fig. 11(p)), the landscape shows a closed ring valley, which represents the stable limit cycle of VOLUME 8, 2020 the deterministic oscillation path. Compared Fig. 11(a) and Fig. 11(p), we can conclude that increasing the noise intensity can induce the system from the monostability to oscillation when the control parameter a is close to the critical point. (II ) When a = 10.34, which is closer to the bifurcation value, the system is still in the monostable region. From Fig. 11(a), one can see that the landscape implies a clear closed ring valley, which represents the stable limit cycle of the deterministic oscillation path. Compared Fig. 11(a), Fig. 11(b) and Fig. 11(p), we can conclude that the closer the control parameter a is to its bifurcation points, the less noise intensity is required to induce the system to oscillate.
(III ) When a = 10.37, the system is in the bistable region. From Fig. 11(c), one can see that in addition to a deep funnel toward the local minimum value (marked by P), there is a closed ring valley in the landscape.
(IV ) When a = 11, 15, 16, respectively, the system is in the oscillation region. From Fig. 11(d)(e)(f), one can see that both the landscapes exhibit a distinct closed ring valley shape, which corresponds to the stable limit cycle of the deterministic oscillation path. When a = 17, 18, 20, respectively, which is close to the bifurcation value(a = 16.767119), the system is in the monostable region. However, the corresponding landscapes (Fig. 11(g)(h)(i))still exhibit a distinct closed ring valley shape. This phenomenon further indicates that the noise can induce the system from the monostability to oscillation under certain conditions.
(V ) When a = 21, 23, 24, respectively, the system is in the monostable region. From Fig. 11(j)(k)(l), one can see that the closed ring valley is gradually disappear and both the landscapes have a global minimum value which represents the stable equilibrium point with a lower concentration.
(VI ) Fig. 11(m) illustrates the energy landscape when φ is relatively large. When a = 0.9, the system is in the bistable region. From Fig. 11(m), one can see that the landscape has two local minimum values which represent the stable equilibrium point with a higher concentration and the stable equilibrium point with a lower concentration, respectively. When increases D (D = 0.01, Fig. 11(n)), the landscape shows that the stable fixed point with a lower concentration is gradually disappear and completely disappear when d = 0.05 (Fig. 11(o)).
In summary, we can draw the following conclusions from Fig. 11: Noise can extend the oscillation region of the system when a is close to the critical points; The closer the control parameter is to the critical points, the less noise intensity is required to induce the system to oscillate; Noise can induce the dynamics to exhibit distinct transition features, including the transition from the monostability with a single stable equilibrium point to the oscillation with a stable limit cycle and the transition from the bistability with two stable fixed points to the monostability with a single stable equilibrium point. Finally, compared Fig. 11 and Fig. 2(e) we find that the global evolution trend of the energy landscapes ( Fig. 11 (a-l)) is completely consistent with the bifurcation diagram in Fig. 2(e).
Oscillation dynamics of the genetic oscillator are the most important and interesting properties in many aspects of cell physiology. To further investigated the probabilistic energy landscape when the system in the stable oscillation region, we calculate the force from the negative gradient landscape (−D ∂U ∂x ) and the force from flux (J ss /P ss ) when a = 10.40, which belongs to the stable oscillation region (see Fig. 12). From Fig. 12, one can see that the force from the negative gradient landscape is primarily vertical to the deterministic oscillation trajectory, and the force from flux is primarily parallel along the oscillation trajectory. This phenomenon indicates that the force from the negative gradient of the energy landscape attracts the function states walk toward the oscillation path, and the flux is the driving force of the oscillations.

C. BARRIER HEIGHT, ENERGY DISSIPATION, AND SENSITIVITY ANALYSIS
To investigate the stability of the oscillations, we introduce the barrier heights and the energy dissipation to measure the stability of the limit cycle [22].

1) BARRIER HEIGHT
The definition of the barrier height is: barrier height=U max − U min , where U max is the local maximum point of the land-scape along inside the limit cycle circle, U min is the minimum point of the landscape along with the limit cycle attractor. The higher the barrier height, the more stable the limit cycle.

2) ENERGY DISSIPATION
As mentioned in the introduction, the genetic oscillator is a non-equilibrium gene regulatory network, which can exchange the energies and information with its surroundings. This property can lead to the energy dissipation, which can be a global measurement of the non-equilibrium gene regulatory network. Furthermore, the energy dissipation has a close relation with the entropy production rate if the steady state exists. The entropy S has the following form [33]: Then, we obtain the increase in entropy S at temperature T by differentiating equation (29): in the above equation, − (k B T ∇lnP − F) · Jdx represents the entropy production rate (e p ), and F · Jdx represents the energy dissipation (h d ). At steady state, we haveṠ = 0, and the entropy production rate is equal to the energy dissipation. Hence, we can use the expression F · Jdx to compute the energy dissipation in steady state. The higher the energy dissipation, the less stable the limit cycle. Fig. 13 illustrates the influence of the different parameters on the barrier height under the same noise intensity. From Fig. 13, one can see that the parameters m, a and k have the top effects on the barrier height among all parameters. Therefore, here we will explore the effects of parameters m, a and k on the stability of the oscillations as shown in Fig. 14. From Fig. 14, we find that the barrier height increases with the increase of parameter a, while decreases with the increase of parameter m and k. We thus conclude that the parameter a (m, k )has the positive (negative) effects on the stability of the oscillations, that is, when we increase the control parameter a in the oscillation region, the limit cycle becomes more (less) stable. Besides, we also study the effects of the noise intensity on the stability of the oscillations by using the physical methods of barrier height and energy dissipation. From Fig. 15, we can see that the barrier height decrease, the energy dissipation increase as the noise intensity increase. These results indicate that the noise intensity has a negative contribute effects on the stability of the oscillations. To further investigate the effects of the important parameters m, a and k on the system, we give the sensitivity analysis about them as shown in Fig. 16. Here, we introduce the index k/k to measure the percentage increase in parameters. From Fig. 16, we find that the parameter a has a significant  growth in barrier height as the index k/k increase, while the parameters m and k stayed almost unchanged. This shows that the parameter a is more sensitive to influence the stability of the system.

IV. CONCLUSION
Genetic oscillators based on the interaction of several gene units have been found to be involved in the regulation of the cell cycles, the 24-hour circadian clocks, or enzyme synthesis. Therefore, uncovering the fundamental properties of such oscillators becomes more important for understanding many aspects of cell physiology.
Here, we investigate the global dynamics of a twocomponent genetic oscillator by using the methods of bifurcation analyzes and probabilistic energy landscapes. Firstly, we perform the local bifurcation analysis at equilibria, including codimension-two bifurcation and codimension-one bifurcation. The results show that the system undergoes the monostable state with higher concentrations of A, oscillations from the Hopf bifurcation and the monostable state with lower concentrations of A with the increase of the control parameter a. In particular, the system can display the bistable state of the coexistence of a stable equilibrium point and a stable limit cycle at a certain region of parameter a. Besides, the system can exhibit the bistable state of the coexistence of two stable equilibrium points when the system parameter φ is relatively large. In addition, compared the bifurcation diagrams we find an interesting phenomenon that the system parameter k can induce similar bifurcation dynamics with system parameter m, and the system parameter r can induce similar bifurcation dynamics with system parameter φ. By analyzing the system equations, we find that system parameter k (r) and m (φ) have the same effect on protein concentration. Then, we further demonstrate these dynamics through the probabilistic energy landscapes in which the force from the negative gradient of the energy landscape attract the system down to the oscillation path, and the flux is the driving force of the oscillations. Also, based on the landscapes, we have the following results: 1)Noise can extend the oscillation region of the system when the control parameter a is close to the critical points. 2) The closer the control parameter is to the critical points, the less noise intensity is required to induce the system to oscillate. 3) Noise can induce the dynamics to exhibit distinct transition features, including the transition from the monostability with a single stable fixed point to the oscillation with a stable limit cycle, and the transition from the bistability with two stable equilibrium points to the monostability with a single stable equilibrium point. Finally, we study the global stability of the oscillations utilizing the physical methods of barrier height and energy dissipation. Our results reveal that the system parameters m, k and the noise parameter D play a negative role in maintaining the stability of the oscillations, and there is a positive relationship between the control parameter a and the robustness of the oscillations. In addition, the sensitivity analysis of the parameters indicates that the parameter a is more sensitive to influence the stability of the system. Our approach is quite general to provide a realistic pathway for studying the global stability, robustness and kinetics of a more complex gene network.