Abstract
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.
1. 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–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.
As a result of interactions between numerous intracellular or extracellular biomolecules, complex cellular behaviors could be easily found, such as bistability [8–12], oscillation [13–18], birhythmicity [19–24], hard excitation, and even chaotic phenomena. Interestingly, originating from different mechanisms based on distinct modes of cell regulation, several rhythms of different periods in certain cells may coexist.
Birhythmicity, namely, the coexistence between two stable regimes of limit cycle oscillations, had been reported in the various types of complex oscillations which can arise in a model based on the mechanism of -induced release (CICR) that takes into account the -stimulated degradation of inositol 1,4,5-trisphosphate (InsP3) by a 3-kinase [21]. Besides simple periodic behavior, this model for cytosolic 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–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–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.
(a)
(b)
(c)
(d)
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.
2. 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 ; otherwise, when component inhibits , the kinetics equation is written as . 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:
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 , , , , and . 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.
3. 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.
3.1. Bifurcation Analysis for Model (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 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 and 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 and 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 point 15.85. Between two Hopf bifurcation points 15.85 and 20.85, it shows hard excitation. With the values of exceeding the Hopf bifurcation point 20.85, the system provides almost symmetrical dynamical phenomena, birhythmicity at first, then returning back to single limit cycles.
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 ; two different amplitudes oscillations at , achieved by using different initial values; and hard excitation, namely, one stable steady state and one stable periodic orbit at . Phase diagrams of the two stable periodic orbits in Figure 3(b) are shown in Figure 3(d). Figure 3(c) shows stable steady state with , , , and and oscillation with , , , and as .
(a)
(b)
(c)
(d)
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. Figures 4–7 and Figures 8–11 illustrate the two parameter bifurcation domains for model (1) and model (2), respectively.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(e)
(f)
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 U-shape 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 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 U-shape 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.
3.2. Bifurcation Analysis for Model (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 and 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 and 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.
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(a) or U-shapes in Figures 8(b), 8(c), and 8(d) domains are divided into several small domains by 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 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 U-shape 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.
4. 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 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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors acknowledge the support from the Natural Science Foundation (nos. 11105040 and 61104138) of China, 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).