Birhythmicity and Hard Excitation from Coupled Synthetic Feedback Loops

Synthetic biology opens up the possibility of creating circuits that would not survive in the natural world and studying their behaviors in living cells, expanding our notion of biology. Based on this, we analyze on a synthetic biological system the effect of coupling between two instability-generating mechanisms. The systems considered are two topologically equivalent synthetic networks. In addition to simple periodic oscillations and stable steady state, the system can exhibit a variety of new modes of dynamic behavior: coexistence between two stable periodic regimes (birhythmicity) and coexistence of a stable periodic regime with a stable steady state (hard excitation). Birhythmicity and hard excitation have been proved to exist in biochemical networks. Through bifurcation analysis on these two synthetic cellular networks, we analyze the function of network structure for the collapse and revival of birhythmicity and hard excitation with the variation of parameters. The results have illustrated that the bifurcation space can be divided into four subspaces for which the dynamical behaviors of the system are generically distinct. Our analysis corroborates the results obtained by numerical simulation of the dynamics.


Introduction
The successful construction of the first synthetic oscillator in 2000 signaled the entry into the new era of artificial cellular rhythms.Known as the repressilator, this oscillator, expressed in E. coli, consists of a set of three repressors coupled cyclically [1].Having been shown to produce sustained oscillations in neural networks, such regulatory structure is reminiscent of recurrent cyclic inhibition [2].Following this first success, a series of synthetic oscillatory networks expressed in bacteria or mammalian cells had been developed, mostly based on genetic regulation [3][4][5][6][7].These synthetic networks exhibit oscillations with tunable frequencies covering a wide range, different from tens of minutes up to 24 h.All these synthetic oscillators are based on one form or another of negative feedback, the realization of which is often highly sophisticated.Since then, negative feedback has proved the necessity of occurrence of oscillations.What interests researchers the most in synthetic oscillators is that they allow finely tuning the parameters that control the oscillatory dynamics of the regulatory network.This contributes to a detained understanding of the conditions in which periodic bahavior arises in biological systems.Such endeavor could be facilitated by the construction of synthetic oscillators uncoupled from natural cellular rhythms, while the latter are often intertwined.
Birhythmicity, namely, the coexistence between two stable regimes of limit cycle oscillations, had been reported in the various types of complex Ca 2+ oscillations which can arise in a model based on the mechanism of Ca 2+ -induced Ca 2+ release (CICR) that takes into account the Ca 2+ -stimulated degradation of inositol 1,4,5-trisphosphate (InsP3) by a 3kinase [21].Besides simple periodic behavior, this model for cytosolic Ca 2+ oscillations in nonexcitable cells shows complex oscillatory phenomena like bursting or chaos.They showed that the model also admitted a coexistence between two stable regimes of sustained oscillations (birhythmicity).And what is more, birhythmicity also had been discovered in Drosophila; circadian oscillations in the levels of two proteins, PER and TIM, result from the negative feedback exerted by a PER-TIM complex on the expression of the PER and TIM genes which code for these two proteins.On the basis of these experimental observations, they recently proposed a theoretical model for circadian oscillations of the PER and TIM proteins in Drosophila.Here they showed that for constant environmental conditions this model is capable of displaying birhythmicity.Birhythmicity, even multirhythmicity, had been reported in synthetic multicellular systems induced by cell-to-cell communication [25].And induced by communication, the two coupled synthetic repressilators could also show inhomogeneous limit cycles [26,27].Inhomogeneous limit cycle, which is also named as hard excitation, is defined as coexistence of a stable periodic regime with a stable steady state.
Moreover, it is ubiquitous for coupled feedback loops occurring in many contexts (such as metabolism, signalling, and development) to control important aspects of cell physiology, for example, circadian rhythms [16][17][18], DNA synthesis, mitosis, and the development of somites in vertebrate embryos.For the sake of interactions between numerous intracellular or extracellular biomolecules, cells exhibit complex behaviors.Therefore, in order to understand complex cellular behaviors, it is important to investigate the topology of cellular circuits and corresponding dynamical characteristics.Recently, structures and functions of network motifs have been investigated numerously.For example, negative feedback induces oscillation, positive feedback induces tunability, and so on.The present researches on natural network or synthetic networks all proved that feedback loops play an important role in maintaining cellular homeostasis, producing bistability, multistability, sustained oscillation, and birhythmicity [24], and making critical decisions such as cell fate and cell development decisions.
However, such feedback loops are usually found as a coupled structure rather than a single isolated form in various cellular circuits [28][29][30].Although there have been some studies on the coupled feedback loops for particular cases, few unified investigations have been reported.Then, questions appear.What are the advantages of such coupled feedback loops, since they must have evolved to achieve specific regulatory functions in natural cellular circuits.To find the answer, we explore the dynamic characteristic of the coupled feedback structures, which are composed of one central negative feedback loop (three-node feedback loop just like the repressilator) and one additional feedback loop.In this paper, we consider the coupled feedback loops that add only one node and two regulations to the central three-node negative feedback loops.The all possible twelve coupled network structures are divided into three types, and the regular dynamics had been discussed [31].According to detailed investigation, we found two topologically equivalent synthetic networks shown in Figure 1 could illustrate enormous dynamical behaviors.
The central negative feedback loop consisting of three repressors (, , and ) has the potential function to generate sustained oscillations.However, many biological oscillators also have an additional positive or negative feedback loop, raising the question of what advantages the extra feedback loop imparts.For simplicity, the extra feedback loops are composed of the original components of the central negative feedback loop and one additional component, , with two regulations.Through computational analysis, we investigate the dynamics behavior of all the twelve coupled feedback structures.We analyze the birhythmicity and hard excitation by means of bifurcation diagrams and locate the different domains of complex oscillatory behavior in parameter space.Through the investigation of the three types of coupled feedback loops [31], only type (2) NN (a negative feedback loop + a negative feedback loop) could show birhythmicity and hard excitation in proper parameter intervals.
Based on our modeling method, we gave the mathematical equation about the network shown in Figure 1(c) as above.In the context of biological networks, we take  as the main bifurcation parameter.Here we choose [0, 200] as the considering interval of parameter  during the single parameter bifurcation.Somewhat surprisingly we have found that the two parameters bifurcation model incorporating the formation of the two negative feedback loops not only can account for regular oscillations in proper parameter intervals, but also may give rise in other parameter intervals to more complex oscillatory phenomena including birhythmicity and hard excitation, the collapse and revival of which depend on the parameters differently.

Mathematical Modeling
Since it is common for cellular networks composed of complicated interconnections among components, those subnetworks of particular functioning are often identified as network motifs.Intriguingly, among such network motifs, feedback loops are very often found as a coupled structure in cellular circuits, which are thought of playing important dynamical roles.Among these synthetic integrated genetic circuits, we mainly focus on one negative feedback loop with three repressors plus one additional feedback loop by adding another regulator.The coupled feedback loops with different topology structures are shown in Figure 1.In this paper we investigate the function and dynamical behaviors of these two topologically equivalent genetic circuits of coupled negative feedback loops.
In order to investigate the potential dynamical behaviors of cellular circuits, which are in general quite complicated due to the nonlinear interaction among the components, we model the interconnected negative feedback loops in the present transcriptional regulatory network with Hill's kinetics without considering the kinetic equation of genes.When component  activates , we describe the dynamics as (/  )  /(1 + (/  )  ); otherwise, when component  inhibits , the kinetics equation is written as 1/(1+(/  )  ).Among the synthetic networks, the transcription processes of gene  are regulated by  or .The whole network is composed of two negative feedback networks, where one consists of components , , and  and the other one consists of components , , , and .In Figure 1, the black and the white circles show opposite regulating functions.
When the regulations among components , , and  are "activation → + repression ⊣, " namely,  activates  and  represses , based on the above assumption, the dimensionless ordinary differential equation (ODE) model for the coupled negative feedback networks in Figure 1(c) is described in (1).Then we use "CA1" (one central negative feedback loop plus one additional feedback loop of components , , , and , where regulations among components , , and  are "activation → + repression ⊣") to denote these coupled negative feedback loops: (1) Otherwise, when the regulations among these three components are "repression ⊣ + activation → , " namely  represses  and  activates , the ODE model for Figure 1(d) is described in (2).And these coupled negative feedback loops are denoted by "CA2" in this paper:

𝑑𝑍
In the above two mathematical models,  ,,, ,  ,,, , , and  ,,, denote basal, activation/inhibition, and degradation rates of the regulatory factors, respectively.The values used in the simulation of parameters are   =   =   =   = 216,   =   =   =   = 0,   =   =   =   = 1,   =   =   = 10, and   =   = 100. indicates the sensitivity of one component with respect to its regulators.  , (,  = , , , ) denote the threshold of  inducing a significant response of .These values are taken based on those in [32,33].To see how the additional negative feedback loop affects the dynamical behaviors of the central negative feedback loop, we solve the differential equations numerically and plot the bifurcation curves with the help of software Matlab and XPPAUT.

Results
By using XPPAUT, we identify the ranges of the parameter, , over which the system exhibits limit cycle oscillations, birhythmicity, hard excitation, and so forth, by locating the bifurcations point at which corresponding dynamics appear and disappear.What is more, how these various dynamical behaviors are influenced by two different parameters is investigated in detail with two parameters bifurcation analysis.(1).The mathematic equations for networks in Figures 1(c) and 1(d) are described as model (1) and model (2).And values of parameters are just as those given in Section 1.In examining the dynamics of these two models, we choose  as a bifurcation parameter and then identify the range of the parameter for corresponding dynamical behaviors.To clearly investigate the evolving of Notations are SLC for single stable limit cycle, SSS for stable steady state, BR for birhythmicity, and HE for hard excitation.Black solid lines stand for limit point bifurcation (LP); dot-dashed lines denote Hopf bifurcation (HB).All these notations represent the same meaning in the following bifurcation diagrams about two parameters throughout the whole paper.The diagrams have been established by numerical integration of (1) with the parameter values listed in Section 2.

Bifurcation Analysis for Model
the system's dynamics, the coupled negative feedback loops in Figure 1(c) are analyzed at first.From Figure 2, it is easy to observe that simple stable steady state (SSS), single limit circle (SLC), birhythmicity (BR), and hard excitation (HE) exist by changing the values of parameter .For the sake of being observed clearly, the main part, [12,30], of the considered parameter interval, [0, 200], is shown in Figure 2 with bifurcation processing exhibited.Through numerical analysis with the help of XPPAUT, we discover that the considered parameter interval, [0, 200], is divided into five subintervals: (I-V).During the subintervals I = [0, 14.91] and V = [21.32, 200], there are one unstable steady state (the dashed lines shown in Figure 2) and one stable limit cycle (the solid blue circles) about the system.Points  = 14.91 and  = 21.32 are limit points, at which bifurcations appear.In detail, in subintervals II = [14.91, 15.85] and IV = [20.85,21.32], the system exhibits various dynamical behaviors: (1) an unstable steady state (the dashed lines); (2) a small-amplitude stable periodic solution (the red solid circles); (3) a middle-amplitude unstable periodic solution (the black circles); and (4) a large-amplitude stable periodic solution (the black solid circles).That is to say, birhythmicity exists in subintervals (II) and (IV).Depending on the initial values taken, one of the two stable periodic solutions is to appear.Points  = 15.85 and  = 20.85 are Hopf bifurcation points.Lastly, in the subinterval III = [15.85, 20.85], the system shows hard excitation, namely, inhomogeneous limit cycles: (1) one stable steady state (the black solid line); (2) one unstable periodic solution (the black circles); and (3) one stable large-amplitude periodic solution.
According to the investigation above, there are four bifurcation points: two limit points 14.91 and 21.32; two Hopf bifurcation points 15.85 and 20.85.As  becomes at the beginning of low values, the system shows single limit cycle; then it shows birhythmicity as  climbs to and crosses the limit point 14.91, without jumping over the Hopf bifurcation What is more, in Figure 3, the time courses of state variable  are illustrated for three examples of these different dynamical behaviors: single stable limit cycle (SLC) at  = 14; two different amplitudes oscillations at  = 15.5, achieved by using different initial values; and hard excitation, namely, one stable steady state and one stable periodic orbit at  = 18.Phase diagrams of the two stable periodic orbits in Figure 3  Models of synthetic genetic applets usually either consist of single synthetic units [1,34] or exhibit different oscillatory behaviors [35,36].There are many study results reported recently describing different types of oscillation in networks of synthetic genetic oscillators.In particular, their models exhibit multistability, oscillation death, inhomogeneous limit cycles, and inhomogeneous steady state.All of these have been proved to exist for biologically realistic parameter ranges.In this paper we present a detailed bifurcation analysis that allows us to determine the origin of the different solutions and the scenarios of transitions between them, therefore providing deeper qualitative and quantitative conclusions about the structure and dynamical behavior of the system.
The further bifurcation analysis about two parameters is performed using the XPPAUT and Matlab packages for two topological equivalent coupled genetic oscillators and shows that already two negative feedback loops provide a large variety of possible regimes.From Figure 4, it is easy to observe that the phenomena of birhythmicity and hard excitation exist in U-shape domains, dominated by the limit point bifurcation lines.Then the Ushape domains are divided into several small domains by two Hopf bifurcation lines.Nevertheless, the values of parameters   and  on the U-shape lines are just the threshold for the system illustrating birhythmicity and hard excitation.If the values of   ,   ,   , and   are lower than the threshold values (the values on the bottom lines of the U-shape), the system just shows single stable limit cycle only or stable steady state only.Moreover, no matter the values of  are smaller than the left threshold lines or larger than the right threshold lines of the U-shape, the system just shows single stable oscillations.All the threshold values for   ,   ,   , and   are Journal of Applied Mathematics lower limit.However, there are minimum and maximum of  for the system showing BR and HE.Moreover, the limit point bifurcation lines and Hopf bifurcation lines will become superposition as the values of   ,   ,   , and   increase.
Figure 5 shows the bifurcation diagrams of parameters  versus   ,   ,   , and   separately in four subpanels.However, the situations of bifurcation domains for BR and HE have the reverted U-shape compared to those above except that the -  bifurcation domains still have the similar Ushape to -  .The U-shape and reversed U-shape lines also represent the limit point bifurcation lines.Then the U-shape and the reversed U-shape domains are also divided into several small domains by Hopf bifurcation lines.Nevertheless, the values of parameters   and  on the U-shape and the reversed U-shape lines are just the threshold for the system illustrating birhythmicity and hard excitation.If the values of   ,   , and   are lower than the threshold values (the values on the top lines of the reversed U-shape), the system just shows single stable limit cycle only or stable steady state only.Moreover, no matter the values of  are smaller than the left threshold lines or larger than the right threshold lines of the reversed U-shape, the system just shows single stable oscillations.All the threshold values for   ,   ,   , and   are lower limit.However, there are minimum and maximum of  for the system showing BR and HE.The exceptional situations in Figure 5(a) are similar to those in Figure 4.
Figure 6 shows the bifurcation diagrams of parameters  versus   ,   ,   , and   separately in four subpanels.From Figure 6, it is easy to observe that the domains of birhythmicity and hard excitation in Figures 6(c) and 6(d) are all closed.In the closed domains dominated by the limit point bifurcation lines, they are also divided into several small domains by Hopf bifurcation lines.Nevertheless, the values of   ,   ,   ,   , and  on the U-shape lines are just the threshold for the system illustrating birhythmicity and hard excitation.There are upper and lower limits for   ,   , and   and only upper limits for   since the system could show BR and HE.
Figure 7 shows the bifurcation diagrams of parameters  versus   ,   ,   ,   ,   , and  separately in six subpanels.From Figure 7, it is easy to observe that the domains of BR and HE are almost closed, where the shapes of the closed domains are different from those in Figure 6.The superposition of limit point bifurcation lines and Hopf bifurcation lines also appears.(2).With the same analysis as above about the network shown in Figure 1(d), we discover that the considered parameter interval [0, 200] of  is also divided into five subintervals (the single parameter bifurcation of  is similar to that in Figure 2 and not shown here): (I-V).During the subintervals I = [0, 16.06] and V = [23.41, 200], there are one unstable steady state and one stable limit cycle about the system.Points  = 16.06 and  = 23.41 are limit points, where bifurcations are present.In detail, in subintervals II = [16.06, 16.62] and IV = [22.88, 23.41], the system exhibits existence of many kinds of dynamic behaviors: (1) an unstable steady state; (2) a small-amplitude stable periodic solution; (3) a middle-amplitude unstable periodic solution; and (4) a large-amplitude stable periodic solution.Namely, birhythmicity exists in subintervals (II) and (IV).Depending on the initial values, one of the two stable periodic solutions appears.Points  = 16.62 and  = 22.88 are Hopf bifurcation points.Lastly, in the subinterval III = [16.62, 22.88], the system shows hard excitation: (1) one stable steady state; (2) one unstable periodic solution; and (3) one stable large-amplitude periodic solution.The different time courses of different rhythmicity are easily obtained and very similar to those in Figure 3. Therefore, they are not shown here.

Bifurcation Analysis for Model
In model (1), the additional regulations among , , and  are "activation → + repression ⊣" shown in Figure 1(c) while in model (2) the additional regulations are "repression ⊣ + activation → " shown in Figure 1(d).As the central negative feedback loop shows oscillation, the additional loops should show oscillation or stable steady state for "CA1" and "CA2" to exhibit birhythmicity and hard excitation.From Figure 8, it is easy to observe that the phenomena of BR and HE exist in U-shape domains in Figures 8(c) and 8(d), dominated by the limit point bifurcation lines.These are the same phenomena as those in model (1).However, in Figure 8(a), the threshold shaped by limit point bifurcation for parameters  and   is different from that in Figure 4(a) for model (1).To illustrate BR and HE, there is a lower limit of parameter   for model (1) and an upper limit of   for model (2).Then the reversed U-shape domains in Figure 8 Nevertheless, the values of parameters   and  on the Ushape lines are just the threshold for the system illustrating birhythmicity and hard excitation.If the values of   are higher than the threshold values (the values on the upper lines of the reversed U-shape), the system just shows single stable limit cycle only or stable steady state only.Moreover, no matter the values of  are smaller than the left threshold lines or larger than the right threshold lines of the reversed U-shape, the system just shows single stable oscillations.All the threshold values for   are upper limit.However, there are minimum and maximum of  for the system showing BR and HE.Except for  versus   , the bifurcation situations for  versus   ,   , and   are similar to those in model (1) shown in Figure 4.Moreover, the limit point bifurcation lines and Hopf bifurcation lines will also become superposition as the values of   ,   , and   increase and   decreases.
From Figure 9, it is easy to observe that the bifurcations for  versus   ,   , and   are similar to those in model (1) shown in Figure 5 with BR and HE existing in reversed Ushape domains.However, the existing domains of BR and HE for  versus   are U-shape in model (1) while being reversed U-shape in model (2).In Figure 10, the bifurcation of parameter  versus   is also different from that in model ( 1) shown in Figure 6, with the other three panels having the similar phenomena.
From Figure 11, we obtain the same results as that compared in models (1) and (2).Only the bifurcation parameters about the additional component  have contrary phenomena between the two models.

Conclusion and Comment
Feedback loops are omnipresent in natural cellular circuits.There have been a lot of functions reported for single feedback loops in recent years; that is, negative feedback loops could reduce response signal amplitude and response time, maintain homeostasis, and are sufficient for oscillation, especially the recent research synthesized oscillator, namely, repressilator, in E. coli, whereas positive feedback loops amplify signals, cause bistability or hysteresis, and elongate response time.However, in natural cellular circuits, feedback loops are usually coupled together rather than existing in isolation.
In this paper, we have investigated the dynamics behaviors of coupled feedback loops.The coupled feedback loops considered here are composed of one central negative feedback loop pluse another additional feedback loops.There could be twelve topology structures for the coupled feedback loops that are added by only one node and two regulations.Based on the components' number of the additional feedback loops, we divide the twelve systems into three types.With the help of XPPAUT and Matlab, we obtained the bifurcation diagrams for every system.Through numerical analysis, we find that the complex oscillatory phenomena including birhythmicity and hard excitation exist in both of these two models.By comparative analysis of two parameters bifurcation, we proved that there are almost the same bifurcation phenomena in models (1) and (2) except that the bifurcation parameters about the additional component  have contrary bifurcation domain shapes between the two models.The present results on the occurrence of complex oscillatory phenomena in our synthetic model are of particular interest for understanding the conditions in which birhythmicity may arise in biological systems, such as Ca + oscillation and circadian rhythm.However, the relative smallness of these domains raises doubts about the possible physiological significance of birhythmicity in regard to circadian rhythm generation and inhomogeneous limit cycles in complex biological systems.Beyond the particular context of synthetic regulatory coupled feedback loops we discuss the results in the light of other mechanisms underlying birhythmicity and inhomogeneous limit cycles in regulated biological systems.
Birhythmicity requires stringent conditions both on the kinetics and on the parameter values.Thus, it is probably less frequent than its well-known stationary counterpart, bistability, in which two stable steady states coexist for a given set of experimental conditions, as demonstrated for several biochemical systems such as the peroxidase reaction.Birhythmicity provides a new mode of physiological regulation as it allows for a switch between two periodic regimes upon suitable perturbation.It would be of interest to search for this phenomenon not only in chemical or metabolic oscillatory systems but also in the many rhythmic processes occurring in the brain, which arise precisely from multiple regulatory interactions between neurons.
Although the relative smallness of these domains raises doubts about the possible physiological significance of birhythmicity in regard to circadian rhythm generation, beyond the particular context of circadian rhythms, we discuss the results in the light of other mechanisms underlying birhythmicity and inhomogeneous limit cycles in regulated biological systems.Furthermore, the development of artificial cellular oscillators opens the way to pharmacological applications such as pulsatile drug delivery.the Natural Science Foundation of the Education Department of Henan, China (no.2011B110003), and the Excellent Young Scientific Talents Cultivation Foundation of Henan University (no.0000A40351).

Figure 1 :
Figure 1: Schematic diagrams of the coupled two negative feedback genetic circuits.Panel (a) is the transcriptional regulation networks for two different coupled negative feedback loops.panel (b) is the simplified and separated schematic diagram of panel (a).In panel (b), the black circle and the white circle are of opposite regulating functions, namely, "activation → + repression ⊣" shown in panel (c) or "repression ⊣ + activation → " shown in panel (d).Arrows denote activation and "⊣" symbols denote repression.The central negative feedback loop consists of three components (, , and ); meanwhile, the additional feedback loop consists of components , , , and .

Figure 4 :
Figure 4: Bifurcation diagrams of two parameters showing the domains of different dynamical behaviors for model (1).(a)  and   , (b)  and   , (c)  and   , and (d)  and   .Notations are SLC for single stable limit cycle, SSS for stable steady state, BR for birhythmicity, and HE for hard excitation.Black solid lines stand for limit point bifurcation (LP); dot-dashed lines denote Hopf bifurcation (HB).All these notations represent the same meaning in the following bifurcation diagrams about two parameters throughout the whole paper.The diagrams have been established by numerical integration of (1) with the parameter values listed in Section 2.

Figure 8 :
Figure 8: Bifurcation diagrams of two parameters showing the domains of different dynamical behaviors for model (2).(a)  and   , (b)  and   , (c)  and   , and (d)  and   for model (2).
(a) or U-shapes in Figures 8(b), 8(c), and 8(d) domains are divided into several small domains by Hopf bifurcation lines.