The Dynamics of Hybrid Metabolic-Genetic Oscillators

The synthetic construction of intracellular circuits is frequently hindered by a poor knowledge of appropriate kinetics and precise rate parameters. Here, we use generalized modeling (GM) to study the dynamical behavior of topological models of a family of hybrid metabolic-genetic circuits known as"metabolators."Under mild assumptions on the kinetics, we use GM to analytically prove that all explicit kinetic models which are topologically analogous to one such circuit, the"core metabolator,"cannot undergo Hopf bifurcations. Then, we examine more detailed models of the metabolator. Inspired by the experimental observation of a Hopf bifurcation in a synthetically constructed circuit related to the core metabolator, we apply GM to identify the critical components of the synthetically constructed metabolator which must be reintroduced in order to recover the Hopf bifurcation. Next, we study the dynamics of a re-wired version of the core metabolator, dubbed the"reverse"metabolator, and show that it exhibits a substantially richer set of dynamical behaviors, including both local and global oscillations. Prompted by the observation of relaxation oscillations in the reverse metabolator, we study the role that a separation of genetic and metabolic time scales may play in its dynamics, and find that widely separated time scales promote stability in the circuit. Our results illustrate a generic pipeline for vetting the potential success of a potential circuit design, simply by studying the dynamics of the corresponding generalized model.

The engineering of biological circuits, synthetically constructed in the laboratory and designed to exhibit novel behaviors, has matured rapidly as a discipline over the past decade. A major challenge for synthetic biology is understanding how different cellular subsystems, composed of distinct components and operating at widely different speeds, interface and work together to generate coherent cellular behaviors. In particular, metabolic pathways can rapidly harvest energy and nutrients for reproduction, while genetic regulatory networks slowly control these pathways in response to environmental changes. Here, we study a generalized version of the metabolator, a unique synthetic circuit composed of both metabolic and regulatory components, previously constructed experimentally and shown to exhibit oscillations. In addition to exploring multiple alternative network topologies, our generalized approach overcomes the inherent uncertainty associated with precise kinetic details of biological circuits. Through this analysis we identify circuit components that are essential for producing sustained oscillations. Furthermore, we find that certain dynamical regimes of a rewired metabolator display unexpected oscillatory properties, such as the capacity to suddenly transition to high amplitude oscillations. This novel behavior could be tested experimentally, and lead to novel synthetic biology modules and applications.

I. INTRODUCTION
Synthetic biology is in need of theoretical tools to design circuits with a high probability of exhibiting specific, user-defined dynamical behaviors. Conventionally, efforts in this regard have focused on explicit kinetic formulations of the dynamics of a circuit. Indeed, deterministic models of the toggle switch 9 and the repressilator 6 , two early examples of successful synthetic circuits, directly predicted the dynamical behavior of these constructs and likely influenced the final design of the circuits themselves. Perhaps the most prominent challenge in assembling these models was a fundamental lack of knowledge regarding kinetic details, including the mechanistic form of the rate laws, as well as the precise values of rate constants themselves.
The primary focus of this work is to develop a modeling framework for understanding the dynamics of a biological circuit for which we have little detailed knowledge of the true kinetics. In particular, we will try to derive generic results dependent exclusively on the "topology" of the circuit, rather than its detailed kinetics. To do so, we will focus our efforts on a particular system known as the metabolator, a synthetic circuit constructed in Escherichia coli K12 which exhibits limit cycle oscillations. At a coarse level of description, the metabolator consists of two metabolites which are interconverted between each other by two enzymes. One of the metabolites in the circuit acts as a regulator, repressing the expression of the enzyme catalyzing the metabolite's production and activating the expression of the enzyme catalyzing its degradation. In Ref. 8 , the authors observed that in conditions of sufficiently high metabolic inflow to the metabolator, a gene encoding one of the regulated enzymes exhibited sustained oscillations. The metabolator is particularly interesting because it remains (to our knowledge) the only example of a synthetic circuit containing a combination of both metabolic and genetic (transcriptional/translational) components. In all other cases, the building blocks of circuits have been entirely transcriptional and translational.
The "hybrid" nature of the metabolator is particularly interesting because metabolism is well-known to take place at a faster rate than transcription and translation. The role such separation of time scales actually plays in the dynamics of the metabolator is unresolved, although it has been observed to be quite in important in other metabolic contexts 17 .
One of our primary approaches to studying the metabolator is based on a non-dimensionalization technique for studying dynamical systems commonly referred to as generalized modeling (GM) 10,12,13,25,26 . In GM, a change of variables is applied to a dynamical system so that many of the otherwise unconstrained parameters in the system (e.g., rate constants) are replaced by "elasticity" parameters with well-defined ranges. Then, the dynamics of the system can be studied around an arbitrary and potentially unknown steady-state. GM's conclusions are invariant to the particular choice of rate law, and depend only on the sensitivity of the rate law at the system's equilibrum. This enables the study of a topological model of the metabolator, in which no assumptions are generally made on the explicit form of the rate law, nor the true values of the parameters in the system. In practice, many GM studies have focused on using large-scale numerical simulations to tease out statistical properties relating elasticities to stability properties of the system. In this work, we supplement such numerical work with complementary analytical studies, and highlight the additional insight that can be gained through explicit calculations. We complement our use of GM with explicit dynamical modeling of the metabolator, and highlight the complementary conclusions that can be drawn from the two techniques (e.g. local vs. global bifurcation phenomena when applying GM or explicit dynamical modeling, respectively). Importantly, we demonstrate how the results from GM can prompt further study with explicit modeling, and vice versa.
The main results of the work are as follows. First, a simple, "core" model of the metabolator (two genes, two metabolites) originally illustrated in Ref. 8 is investigated using a conventional dynamical model. We find that the system lacks the capability to pass through a Hopf bifurcation. Second, using GM, we analytically prove that any dynamical model with a topology analogous to the core metabolator is incapable of Hopf bifurcations. Concluding that additional topological elements are necessary for a Hopf bifurcation, we use numerical simulations with GM to progressively restore elements of the metabolator model from Ref. 8 which we omitted in our core model. Doing so, we reveal the crucial topological components of the metabolator which endow it with the ability to oscillate. Third, we investigate the possibility of oscillations in an explicit model of a simply "re-wired" version of the metabolator in which the regulatory connections are reversed (the reverse metabolator, "RM"). We show that the RM is capable of sustained oscillations by both a local Hopf bifurcation as well as global fold of limit cycles. Fourth, prompted by findings in the explicit model and using the GM tools we developed earlier, we study the role that metabolic and genetic time scales play in generating instability and potential oscillations in the RM in vivo.

II. THE CORE METABOLATOR
To begin to understand the fundamental components which generate sustained oscillations in the metabolator, we developed a four-dimensional dynamical model of the system.
Our particular choice of design was motivated by the schematic provided in Ref. 8 and re- illustrated in I A. Two enzymes (G 1 and G 2 ) interconvert two substrates (M 1 and M 2 ). Gene G 1 codes for the enzyme which converts M 1 into M 2 , and gene G 2 codes for the enzyme which converts M 2 into M 1 . We explicitly assume that the rate of gene expression is substantially faster than protein translation, so that the abundances of the mRNA transcripts encoding G 1 and G 2 may be assumed to be in quasi-steady-state. Metabolite M 2 represses the expression of gene G 1 and activates the expression of gene G 2 . Furthermore, we assume that the rate of inflow of M 1 is constant, and the rate of outflow of M 2 is in proportion to its concentration.
Because of the symmetry in the dynamics of M 1 and M 2 (both are produced and degraded by the same metabolic reactions), we replace M 1 with the total amount of substrate in the system C = M 1 + M 2 . Furthermore, we assume bilinear mass-action rate laws for the dynamics of the metabolic reactions. We dub our simplified model the core metabolator. The nondimensional equations describing the dynamics of the core metabolator are where the parameters ρ, η 1 , η 2 are assumed to be strictly positive. See Appendix A for a derivation of (1). We emphasize that the core metabolator serves as a jumping off-point: observing oscillations here would suggest the presence of oscillations in more detailed models capturing the true topological elements of the in vivo metabolator. On the other hand, not observing oscillations in this simplified model implies that one needs to examine a model with more components in order to observe Hopf bifurcations, see Section 4.
The fixed point of (1) is g * 1 = g * 2 = 1 2 , m * 2 = 1 and c * = η 1 +η 2 +2ρ η 1 . The Jacobian J, evaluated at the fixed point, becomes The characteristic equation for J is In order to identify potential candidates for a Hopf bifurcation, we make use of a constraint relating the coefficients of (3) which must be satisfied in order for (3) to have a pair of purely imaginary roots. First, we note that (3) has a root at λ = −1. Then, we can substitute λ 1,2 = ±iω into the cubic polynomial in (3) we find that these conditions are c 0 = ω 2 c 2 and c 1 = ω 2 c 3 , where c i is the coefficient of the λ i th term of the cubic polynomial in parentheses in (3). After substituting the appropriate coefficients into these conditions and some algebra, we find It is not clear if (4) can be satisfied with real, positive values of the system parameters.
Rewriting (4) as a quadratic polynomial in η 1 , we obtain Now the outcome becomes plainly clear: because (5) only contains positive coefficients, Descartes' Rule of Signs immediately proves that it cannot have positive real roots. Thus, there exists no combination of (η 1 , η 2 , ρ) which will yield a pair of purely imaginary eigenvalues, and the system does not exhibit a Hopf bifurcation.

III. GENERALIZED MODELING OF THE CORE METABOLATOR
Can we say anything more general about the oscillatory capabilities of the core metabolator illustrated in Figure 1A? In the previous section, we showed that one particular realization of this schematic model did not exhibit a Hopf bifurcation, but it remains unclear whether this was a quality of our choice of model. One may naturally wonder if choosing different kinetic rate laws for the metabolic reactions, or different regulatory kinetics for the feedback regulation, would yield qualitatively different behavior like oscillations. In this section, we propose a method for addressing these concerns by applying a technique known as generalized modeling (GM). This section is organized as follows. First, we introduce a motivating example of the GM nondimensionalization procedure. Next, we illustrate the role this non-dimensionalization plays in the formulation of generalized models. Then, we apply this procedure to generate a GM of the core metabolator. Finally, we analyze the stability properties of the GM, and analytically prove that a topological model of the core metabolator with fairly general rate laws is incapable of a Hopf bifurcation.

A. Generalized Modeling
First, we present the normalization procedure upon which GM is based, known as generalized modeling, see for example Ref. 13 . Consider a system of ordinary differential equations, Here, X and Y are both functions of time. The stability of (6) around its equilibrium is determined by the eigenvalues of the Jacobian J. Here, we approach the problem of calculating J after first executing a particular change of variables central to GM. Assume that there exists an equilibrium (X 0 , Y 0 ) and normalize (6) by its steady state concentration Letting x = X X 0 , y = Y Y 0 so that the variables are nondimensionalized and the equilibrium is at (1, 1), we can write To find the components of the Jacobian for dx dt , we simply differentiate each term by each independent variable and evaluate at the equilibrium (x, y) = (1, 1). If we let h = dx dt and write the components of the Jacobian as ∂h ∂x | (1,1) and ∂h ∂y | (1,1) , we find where . By repeating the calculations outlined above for dy dt , we can calculate the remaining entries of the Jacobian. We argue that these θ's (herein referred to as elasticities and similar to the elasticities from metabolic control analysis 7 ) have well-defined and limited ranges which makes performing calculations with them much easier than with similar calculations with the original parameters. To demonstrate the utility of elasticities, we will reformulate the dynamics of G 1 in the core metabolator. To do so, we will begin with the first equation of the dimensionalized dynamics of the metabolator, reproduced from Appendix A (see (A1a)), We let g 1 = G 1 G 1,0 and m 2 = M 2 M 2,0 (where G 1,0 and M 2,0 are the steady-state concentrations of G 1 and M 2 , respectively) . Then, we can write We will calculate the elasticity of the first term in (11) corresponding to the production of g 1 . To do so, we normalize this production term by its magnitude at steady state: To calculate the elasticity, we take the derivative of µ(m 2 ) with respect to m 2 , evaluated Upon inspection, this elasticity can only take values in [−2, 0]. When M 2,0 K 1 , the elasticity approaches negative two and the production of G 1 becomes very sensitive to small changes in the abundance M 2,0 . In contrast, when M 2,0 K 1 , the elasticity approaches zero and the production of G 1 is insensitive to changes in M 2,0 .
The well-defined and limited range of elasticities becomes extremely powerful in studying the dynamics of a system near its equilibrium. As we show in the ensuing section, these elasticities feature prominently in the Jacobian of a dynamical system reformulated with GM. We will demonstrate how elasticities can be used to prove the generic stability of a topological model of a dynamical system, with only mild assumptions on the explicit rate laws governing the dynamics.

B. Generalized Model of the Core Metabolator
We can use the aforementioned normalization procedure to develop a generalized model of the core metabolator. As we will show, the generalized model enables us to draw conclusions about a topological model of the core metabolator, while relaxing many of the assumptions we made in our explicit model of the core metabolator in Section 2. We begin with the following "generalized" form of the dynamics: The terms R 1 and R 2 correspond to the reactions catalyzed by G 1 and G 2 , respectively, while I and R 3 correspond to the inflow and outflow of mass from the metabolator. The P i and D i terms correspond to the production and degradation of proteins. Importantly, the P i terms capture the feedback regulation of metabolite M 2 on the two enzymes.
We make a few assumptions regarding the dynamics of the core metabolator in our generalized model. First, we assume R i , P i and D i are generic, monotonically increasing functions of their arguments (i.e. ∂R 1 ∂M 1 > 0). Many frequently used kinetic rate laws, including Hill, Michaelis-Menten, and mass-action kinetics satisfy this assumption, among many other possible forms of kinetics. We also retain the earlier assumption that the rates of metabolic reactions are linear with respect to the amount of gene product. The structure of (12) corresponds precisely to the topological structure of the core metabolator in Figure 1A.
Before we proceed to calculating the Jacobian of (12), we briefly describe a simple method for calculating the prefactors (those parameters which were not elasticities, i.e. O(X 0 , Y 0 )) which appear in the Jacobian and were described in the prior section. In prior work on GM 25 , it has been shown that for metabolic networks, these prefactors correspond precisely to the rates of metabolic reactions at equilibrium. As we show below, these rates can be calculated by enforcing that at equilibrium the fluxes of reactions flowing into and out of each metabolite are balanced 19 . Furthermore, for the non-metabolic components of the metabolator, we can apply a similar balancing procedure by noting that at equilibrium, the production and degradation rates for each protein must be equal. Applying these balancing conditions to (12), we find that Here, M 1,0 , M 2,0 , G 1,0 and G 2,0 indicate the steady state concentrations of metabolites 1 and 2, and gene products 1 and 2, respectively. The first two constraints in (13) force the rate of reaction G 1 R 1 (corresponding to the conversion of M 1 to M 2 ) to be precisely balanced by the rate of reaction G 2 R 2 (the conversion of M 2 to M 1 ) and I, the rate of inflow of M 1 to the system. Similarly, G 2 R 2 must be precisely equal to the sum of G 1 R 1 and R 3 (the outflow of M 2 from the system). The final two conditions impose that, for each gene, the rates of production and degradation must be balanced.
Because of the interdependency of all the prefactors in (13), it becomes useful to explicitly introduce the parameter v = G 2,0 R 2 (M 2,0 ). We also introduce a free parameter α > 0 by setting αv = R 3 (M 2,0 ). From the equilibrium condition (13), it follows that I = αv and G 1,0 R 1 (M 1,0 ) = (1 + α)v. Finally, we let L 1 = P 1 (−M 2,0 ) = D 1 (G 1,0 ), L 2 = P 2 (M 2,0 ) = D 2 (G 2,0 ). We can now rescale (12) in the same manner as (8): where The Jacobian matrix J of the core metabolator is then The parameters θ i in (15) correspond to the elasticities of metabolic reactions described earlier: Here, θ i > 0, i = 1, 2, 3 follows directly from the earlier assumption that all P i , D i , R i are monotonically increasing functions of their arguments. Similarly, the parameters (15), corresponding to elasticities of the regulatory reactions, can be calculated by: For the generalized model of the core metabolator , θ f 1 < 0, θ f 2 > 0 since M 2 inhibits G 1 and activates G 2 . Furthermore, θ d1 , θ d2 > 0 since an increase in the abundance of a gene product leads to an increase in its rate of degradation.
The utility of reformulating the core metabolator (1) with the GM in (14) is now plainly clear. Many of the prior assumptions made on the rate laws in our explicit dynamical model (1) have been relaxed. Our sole assumption is that the rate laws are monotonic, so that it is always the case that the signs of the elasticities are known. Now, by analytically studying the eigenvalues of J, we can understand the stability properties of the core metabolator across the entire spectrum of possible parameter values and rate laws.
In Appendix B, we show that (15)  The size of the dynamical system (greater than four state-space variables) made confirming these results using the Routh-Hurwitz criterion analytically intractable. The results of our numerical studies were unexpected. We found that reintroducing G 3 had no bearing on the appearance of a Hopf bifurcation. In fact, the model remained entirely stable, failing to exhibit a single instance of an eigenvalue with positive real part. On the other hand, the refined three-metabolite model did exhibit instability and Hopf bifurcations. These results suggest that the precise nature of the regulation is crucial for the metabolator to exhibit oscillations. In order for a Hopf bifurcation to be possible, the regulatory metabolite could not be repressing the expression of its own degrading enzyme. This delayed positive feedback appears to be a crucial element endowing the metabolator with the potential to undergo a Hopf bifurcation.

V. REVERSE METABOLATOR
The difficulty in constructing synthetic circuits with poorly studied components often motivates synthetic biologists to construct several alternative designs. Frequently, these alternative designs are simply composed of rewired versions of the same regulatory components. In light of this, we investigated the dynamical behavior of one such re-wired version of the core metabolator, which we dub the reverse metabolator (RM). In the RM, the regulatory connections found in the core metabolator between M 2 and the two genes are reversed.
Thus, M 2 now activates G 1 , and hence its own production, and represses G 2 , and hence its own degradation (see Figure 1D). We selected this design precisely because synthetic biologists frequently build a number of alternative constructs in the course of circuit assembly, in order to increase the odds of identifying a succesful design.
In nondimensional form, the RM is identical to (1) but with the regulation of g 1 and g 2 reversed:
Using Descartes' rule of signs, (21) must contain positive, real roots since η 2 > 0 and the η 0 1 term is negative. Therefore, there exists an open set of points (ρ * , η * 1 , η * 2 ) at which the reverse metabolator undergoes a Hopf bifurcation. Using the numerical continuation software package AUTO 5 , a number of Hopf bifurcations spanning several orders of magnitude in η 1 , η 2 , and ρ can be identified.
We used AUTO to study the behavior of the RM's limit cycles in different regions of parameter space. As shown in Figure 2, the limit cycles exhibit two notable features depending on the magnitude of inflow into the RM, characterized by the parameter ρ. First, for small ρ, two distinct periodic orbits exist: a small amplitude, stable limit cycle and a large amplitude unstable limit cycle. These two cycles annihilate at a fold of limit cycles.
Both limit cycles display relatively smooth changes in amplitude ( Figure 2C). Second, for

B. Criticality of the Hopf Bifurcation
We investigate the stability of the limit cycles generated by the Hopf bifurcations in the RM. Limit cycles arising from supercritical Hopf bifurcations are stable to small perturbations, while those generated by subcritical Hopf bifurcations are not. To explore the criticality of the observed Hopf bifurcations, a two-parameter bifurcation analysis of the RM was conducted in AUTO, shown in Figure 3A. As the two parameters ρ and η 1 were varied, the criticality of the Hopf bifurcation was tracked using the magnitude of the leading Floquet multiplier. We observed that as the value of ρ was decreased, the Hopf bifurcation switched from being subcritical to being supercritical. As the value of ρ was decreased further, a second transition in criticality was observed, and the Hopf bifurcations returned to being subcritical.
The point at which a Hopf bifurcation switches criticality is called a Bautin bifurcation.
This higher-order bifurcation is associated with the collision of a fold of limit cycles with a manifold of Hopf bifurcations. Using AUTO, we were able to follow the fold of limit cycles as it branched off from one Bautin bifurcation, tracked closely with the manifold of Hopf bifurcations, and reintersected with the Hopf curve at the second Bautin point ( Figure   3A). The fold of limit cycles is the boundary between two regions of parameter space where distinct qualitative dynamical behaviors of occur. On one side of the fold (in this case, for Region B in Figure 3), no sustained oscillations exist. On the other side of the fold (Region C in Figure 3), a pair of limit cycles appears. For the RM, this pair contained a small amplitude stable limit cycle and a large amplitude unstable limit cycle.
The existence of a Bautin bifurcation indicates there are two distinct routes to oscillation in the reverse metabolator. First, by decreasing ρ for a fixed η 1 , η 2 (moving from point D to point C in Figure 3A), the reverse metabolator passes through a supercritical Hopf bifurcation. In this case, oscillations appear at infinitesimal amplitude, and their amplitude grows as ρ is decreased further. This type of bifurcation is identical to the one observed in the in vivo metabolator 8 . A second, and potentially more interesting, route to oscillation appears if the reverse metabolator passes through a fold of limit cycles. As the value of ρ is increased for fixed η 1 , η 2 (moving from point B to C in Figure 3A), a finite-amplitude stable limit cycle appears, surrounded by a larger amplitude unstable limit cycle. Points perturbed from the unstable steady state will immediately be attracted to oscillations that are much larger in magnitude than those observed in the immediate vicinity of a Hopf bifurcation. Perturbations outside of the basin of attraction of the stable limit cycle will lead to spiraling instability, and eventually to the extinction of one of the metabolites. We examine the biological implications of alternative modes of oscillation in the RM in the Discussion, Section 7.

VI. GENERALIZED MODEL OF THE REVERSE METABOLATOR
The appearance of relaxation oscillations inspired us to investigate the role time scales play in the dynamics of the RM. The existence of such divergent time scales is well-known in the regulation of metabolism. The allosteric feedback of small metabolites and the posttranslational modification of enzymes target the activity of an enzyme pool at time-scales on the order of microseconds 1,4,20,21 . In contrast, the transcriptional regulation of enzyme levels (which occurs through the binding of protein transcription factors to DNA sequences) takes place at a substantially slower rate, on the order of minutes 1,15,24 , and affects the actual size of the enzyme pool directly. Remarkably, the rate of an enzymatically catalyzed reaction itself, which the above processes ultimately regulate, is known to take place at characteristic rates spanning seven orders of magnitude, from microseconds to minutes 3 .
In order to study the inherent time scales of the RM, we constructed a generalized model. This model was identical in structure to (14), except that the signs of the elements corresponding to the transcriptional regulation of M 2 were reversed (θ r1 > 0, θ r2 < 0). Our intent was to use the entries in the Jacobian J to identify time-scales which characterized the dynamics of metabolites and genes in the RM. To understand the connection between J and time-scales, recall that each row in J describes the linearized dynamics of a small perturbation away from the equilibrium of one state variable. For example, the first row of J corresponds to the dynamics of a perturbation in M 1 : where δm 1 , δm 2 , δg 1 , and δg 2 correspond to perturbations in m 1 , m 2 , g 1 , and g 2 , respectively.
Upon inspection of (22), all terms on the right hand side share a common divisor M 1,0 v which has units of time. We can treat this shared prefactor as a natural measure of the time scale Two major trends regarding the stabiilty of the RM emerged from our simulations. First, we found that as the dynamics of the genetic elements in the reverse metabolator grew slower, the stability of the circuit as a whole increased. As the value of 1 β grew smaller and eventually became faster than the metabolic time scale, the stability of the circuit was observed to fall significantly. Second, our simulations also revealed that small regulatory elasticities contributed to the stability of the circuit, while large elasticities significantly reduced the probability of a stable circuit. This observation qualitatively agrees with our first observation relating high values of β to instability. Large magnitude regulatory elasticities increase the same terms in the Jacobian as β. Together, our observations suggest that exceptionally slow regulation in the metabolator stabilizes the circuit, a conclusion that parallels commonly held notions on the speed of regulation and the role it plays in stability 18 .

VII. CONCLUSIONS AND DISCUSSION
Contemporary biology is charged with understanding how fundamentally different cellular processes such as signalling, metabolism, and transcriptional regulation, interact and function as a coherent whole. Our work here explored the different types of conclusions which may be drawn from detailed kinetic models and more generic generalized models in the study of a hybrid metabolic-genetic oscillator. We first showed that the core metabolator (1) has a unique fixed point which is stable for all parameter values (recall the end of Section 2). Hence, the core metabolator cannot undergo Hopf bifurcation to oscillations. Then, we developed a GM model of the core metabolator (12), and we showed that this system also cannot undergo Hopf bifurcations by applying the Routh-Hurwitz criterion to the Jacobian (15) (see Section 3.2 and Appendix B).
These analytical results for the core metabolator (1) and its generalized model (12) suggested that the network topology and the number of components in the network play an important role in determining whether or not oscillations are possible. As a result, in Section 4, we explored two more detailed versions of the core metabolator. The first of these has an intermediary delay gene G 3 ( Figure 1B), but did not exhibit Hopf bifurcations in any of the large number of numerical simulations with randomly chosen system parameters we carried out. The second of these consists of an extra metabolite, which is also present in the synthetically constructed circuit 8 , effectively "un-lumping" M 2 . This enhanced network exhibited Hopf bifurcations and instability for many of the parameter values we simulated.
The delayed positive feedback created by the presence of this third metabolite is responsible for the induced oscillations. Moreover, this result also suggests that the self-repression of M 2 induced by the network topologies for the core metabolator, generalized model of the core metabolator, and the first "enhanced" network ( Figures 1A and 1B), is responsible for the absence of oscillations in those systems.
Theory and simulation play important roles in synthetic biology. As we have shown, they can highlight generic designs which have certain desirable behaviors, e.g. oscillations, and they can help exclude those that do not. In this respect, the analysis of the enhanced metabolator model with a third metabolite ( Figure 1C) suggested that we explore some other four-component models with related designs that do exhibit oscillations. In particular, we studied system (18), the reverse metabolator, in which the wiring from M 2 to the two genes is reversed compared to the core metabolator, as shown in Figure 1D. This fundamentally altered the system dynamics so that M 2 now activates, rather than represses, the gene G 1 and hence M 2 is now self-activating. The analysis of Section 5 of the RM reveals that it exhibits both local and global routes to oscillations, via Hopf bifurcations and saddle-node bifurcations of limit cycles, respectively, as shown in Figure 3.
Most of the experimentally observed oscillations in synthetic circuits have arisen through local Hopf bifurcations 6,8,22 . This is very likely to be due to the comparable ease of predicting local bifurcations. Indeed, the GM methods which are used throughout our work are only capable of identifying global bifurcations in very special and exceedingly rare cases (i.e. when they are known to occur near the intersection of two local bifurcations, like a double Hopf bifurcation). In contrast to the Hopf (where oscillations appear at infinitesimal amplitude), oscillations arising from a fold of limit cycles (as seen in the RM) immediately appear at a finite, non-zero amplitude. This feature may make such a "global" route to oscillation desirable in synthetic biology, which is still tackling the problem of designing robust cellular oscillators. The similarity between the RM and the original metabolator (ultimately, simply just a switch of promoters in the architecture of the circuit) suggests that constructing the RM and investigating its dynamics in vivo may be relatively easy, and may lead to the experimental observation of a global bifurcation.
Finally, we propose that the construction of large generalized models of biological circuits may enable the resolution of several open problems regarding the dynamics of large, complex systems. We investigated one such question here (recall Section 6): the role that time scales play in determining the stability of the system. While we found that stability generally increased as the time scale of regulation became longer (Figure 4), it remains unclear how general this conclusion is among biological networks. A closely related question is whether strongly stabilizing components exist within biological networks. Such components have been argued for in ecological networks 14 , as well as in prior work on metabolism using GM 25 . However, the extent to which components of distinct types of biological networks (e.g. signalling, transcriptional regulation, post-translational modifications) act to stabilize or destabilize metabolism is unknown. The investigation of such questions, aided by generalized modeling, seems potentially tractable numerically; computational toolboxes for the automated assembly and inspection of generalized models already exist 11 , and their extension to GM with regulation appears straightforward.
where ρ 1 = V in β 1 K 1 , η 1 = v 1 K 1 β 1 and η 2 = v 2 K 1 β 1 . Finally, the nondimensionalized equation for M 2 follows simply and is where ρ 2 = Vout β 1 . We also note that a simpler way to define the system is in terms of "total" metabolite in the system C = M 1 + M 2 . Then,Ċ =Ṁ 1 +Ṁ 2 = V in − V out M 2 . We can non-dimensionalize this system as well to find that where c = C K 1 . For the purpose of simplifying the analysis in the article, we let c 1 = c 2 = γ 1 = γ 2 = 1 and ρ 1 = ρ 2 = ρ. Then, the system formed by (A4),(A5),(A9), and (A8) is the nondimensionalized core metabolator we study in this article.
to the coefficient of the λ i th term): The Jacobian of (C1) is We assume again that θ r1 < 0 and all other elasticities are greater than zero.
For the metabolator in Figure 1D, we use the GM: