The dynamics of hybrid metabolic-genetic oscillators

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⁹ and the repressilator,⁶ 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.¹⁷ 

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 steady-state. 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.

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 Fig. 1(a). Two enzymes (G₁ and G₂) interconvert two substrates (M₁ and M₂). Gene G₁ codes for the enzyme which converts M₁ into M₂, and gene G₂ codes for the enzyme which converts M₂ into M₁. 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₁ and G₂ may be assumed to be in quasi-steady-state. Metabolite M₂ represses the expression of gene G₁ and activates the expression of gene G₂. Furthermore, we assume that the rate of inflow of M₁ is constant, and the rate of outflow of M₂ is in proportion to its concentration. Because of the symmetry in the dynamics of M₁ and M₂ (both are produced and degraded by the same metabolic reactions), we replace M₁ with the total amount of substrate in the system $C=M1+M2$⁠. 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 Eqs. (1a). 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 Sec. IV.

The fixed point of Eqs. (1a) is $g1*=g2*=12,m2*=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 Eq. (3) which must be satisfied in order for Eq. (3) to have a pair of purely imaginary roots. First, we note that Eq. (3) has a root at $λ=−1$⁠. Then, we can substitute $λ1,2=±iω$ into the cubic polynomial in Eq. (3) we find that these conditions are $c0=ω2c2 and c1=ω2c3$⁠, where c_i is the coefficient of the $λith$ term of the cubic polynomial in parentheses in Eq. (3). After substituting the appropriate coefficients into these conditions and some algebra, we find

It is not clear if Eq. (4) can be satisfied with real, positive values of the system parameters. Rewriting Eq. (4) as a quadratic polynomial in $η1$⁠, we obtain

Now the outcome becomes plainly clear: because Eq. (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.

Can we say anything more general about the oscillatory capabilities of the core metabolator illustrated in Figure 1(a)? In Sec. II, 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 GM.

GM is a nondimensionalization procedure for reformulating the kinetics of a dynamical system, enabling the study of local dynamics near an equilibrium in an efficient way. To do so, parameters in the system are rewritten in terms of normalized parameters known as elasticities. The elasticities have a direct connection to the original kinetic parameters, but are much easier to work with. Most importantly, elasticities typically have well-defined and limited ranges (e.g., [0,1]), and sampling them across this range effectively samples all possible values of the original kinetic parameters (which may otherwise span several potentially unknown orders of magnitude). Furthermore, these elasticities are not necessarily tied to any single choice of kinetic rate law; instead, they simply represent the sensitivity of any rate law to small changes in the equilibrium value of a variable, normalized by the steady-state magnitude of that rate law. In the context of metabolism, elasticities are similar to logarithmic derivatives of the rate law with respect to a change in substrate concentration.^7,23

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.

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 Eq. (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 $(X0,Y0)$ and normalize Eq. (6) by its steady state concentration (assuming $I(X0,Y0),O(X0,Y0)≠0$⁠)

Letting $x=XX0,y=YY0$ so that the variables are nondimensionalized and the equilibrium is at (1, 1), we can write

To find the components of the Jacobian for $dxdt$⁠, we simply differentiate each term by each independent variable and evaluate at the equilibrium (x, y) = (1, 1). If we let $h=dxdt$ and write the components of the Jacobian as $∂h∂x|(1,1)$ and $∂h∂y|(1,1)$⁠, we find

where $θxI=1I(X0,Y0)(∂I(xX0,yY0)∂x)$ and $θxO=1O(X0,Y0)(∂O(xX0,yY0)∂x)$⁠. By repeating the calculations outlined above for $dydt$⁠, 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⁷) 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₁ 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 Eq. (A1a))

We let $g1=G1G1,0$ and $m2=M2M2,0$ (where $G1,0$ and $M2,0$ are the steady-state concentrations of G₁ and M₂, respectively). Then, we can write

We will calculate the elasticity of the first term in Eq. (11) corresponding to the production of g₁. To do so, we normalize this production term by its magnitude at steady state

To calculate the elasticity, we take the derivative of $μ(m2)$ with respect to m₂, evaluated at $m2=1$

Upon inspection, this elasticity can only take values in [–2, 0]. When $M2,0≫K1$⁠, the elasticity approaches negative two and the production of G₁ becomes very sensitive to small changes in the abundance $M2,0$⁠. In contrast, when $M2,0≪K1$⁠, the elasticity approaches zero and the production of G₁ is insensitive to changes in $M2,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.

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 Sec. II. We begin with the following “generalized” form of the dynamics

The terms R₁ and R₂ correspond to the reactions catalyzed by G₁ and G₂, respectively, while I and R₃ 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₂ on the two enzymes.

We make a few assumptions regarding the dynamics of the core metabolator in our generalized model. First, we assume $Ri,Pi$⁠, and D_i are generic, monotonically increasing functions of their arguments (i.e., $∂R1∂M1>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 Eq. (12) corresponds precisely to the topological structure of the core metabolator in Figure 1(a).

Before we proceed to calculating the Jacobian of Eq. (12), we briefly describe a simple method for calculating the prefactors (those parameters which were not elasticities, i.e., $O(X0,Y0)$⁠) which appear in the Jacobian and were described in Sec. III A. In prior work on GM,²⁵ 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.¹⁹ 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 Eq. (12), we find that

Here, $M1,0,M2,0,G1,0$⁠, and $G2,0$ indicate the steady state concentrations of metabolites 1 and 2, and gene products 1 and 2, respectively. The first two constraints in Eq. (13) force the rate of reaction $G1R1$ (corresponding to the conversion of M₁ to M₂) to be precisely balanced by the rate of reaction $G2R2$ (the conversion of M₂ to M₁) and I, the rate of inflow of M₁ to the system. Similarly, $G2R2$ must be precisely equal to the sum of $G1R1$ and R₃ (the outflow of M₂ 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 Eq. (13), it becomes useful to explicitly introduce the parameter $v=G2,0R2(M2,0)$⁠. We also introduce a free parameter $α>0$ by setting $αv=R3(M2,0)$⁠. From the equilibrium condition (13), it follows that $I=αv and G1,0R1(M1,0)=(1+α)v$⁠. Finally, we let $L1=P1(−M2,0)=D1(G1,0),L2=P2(M2,0)=D2(G2,0)$⁠. We can now rescale Eq. (12) in the same manner as Eq. (8)

where $m1=M1M1,0,m2=M2M2,0,g1=G1G1,0,g2=G2G2,0$⁠.

The Jacobian matrix J of the core metabolator is then

The parameters $θi$ in Eq. (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 $Pi,Di,Ri$ are monotonically increasing functions of their arguments. Similarly, the parameters $θf1,θf2,θd1,θd2$ in Eq. (15), corresponding to elasticities of the regulatory reactions, can be calculated by

For the generalized model of the core metabolator, $θf1<0,θf2>0$ since M₂ inhibits G₁ and activates G₂. 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 Eq. (1a) with the GM in Eq. (14) is now plainly clear. Many of the prior assumptions made on the rate laws in our explicit dynamical model Eq. (1a) 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 Eq. (15) cannot have eigenvalues with positive real part for any admissible value of parameters. The proof in the Appendix follows by an application of the Routh-Hurwitz criterion, a well-known tool from classical control theory.² Since Eq. (15) does not have any eigenvalues with positive real part, we have proven that all equilibria of the GM of the core metabolator are exponentially stable to small perturbations. This directly implies that the GM of the core metabolator cannot undergo a Hopf bifurcation.

The main result of this section bears some elaboration: the GM of the core metabolator encompasses all possible explicit realizations of the topological model shown in Figure 1(a), under the general hypotheses about the rate laws stated above. With elementary techniques, it was possible to show that none of these potential realizations could ever pass through a Hopf bifurcation. We note that there may still be global mechanisms, such as a saddle node of limit cycles, which introduce oscillations. This naturally suggests that some additional component of the topology of the in vivo metabolator, which has been removed in our model of the core metabolator, must be crucial to the existence of a Hopf bifurcation. This idea is explored further in Sec. IV.

The analysis in Secs. II and III is suggestive: which topological simplifications in the core metabolator and in the GM of the core metabolator eliminated the possibility of a Hopf bifurcation? To address this question, we used GM to inspect which components would lead to oscillations when reintroduced into the core metabolator. To do so, successively more complete versions of the in vivo metabolator were created, as shown in Figures 1(b) and 1(c). We developed two hypotheses for the crucial topological features of the metabolator which would be necessary for sustained oscillations. First, we reasoned that a third gene, labeled G₃ and representing LacI from Ref. 8, might act as an intermediary delay in the negative feedback from M₂ to G₁ and consequently induce oscillations (see Figure 1(b)). Such delays in feedback have been implicated in the literature as a major source for the onset of oscillations. Second, we separated the lumped metabolite M₂ into its appropriate constitutive molecules (now labeled M₂ and M₃) in Figure 1(c). In the revised topology, the enzyme transcribed from gene G₂ reconverts M₃ to M₁, but the regulation of G₁ and G₂ remains dependent on M₂. For each of our new topological models of the metabolator, the corresponding GM was assembled (see Appendix  C) and numerical simulations completed to sample the space of parameters in the model. In each instance of the simulation, a random sample from the entire possible space of θ parameters was drawn. In order to efficiently sample the numerical space, we sampled uniformly from the ranges $θr2=(0,2],θr1=[−2,0)$ and $θi=(0,1]$ for all other elasticities. Since the parameters $M1,0v,M2,0v,$ and $1β$ remained unconstrained, we decided to sample these parameters across seven orders of magnitude (from $10−3$ to $103$ 1/s, assuming $M1,0=M2,0$⁠). For each instance of our simulation, the eigenvalues of J were calculated. The process was repeated 20 000 times for each choice of $M1,0v,M2,0v$⁠, and $1β$⁠. In total, over $106$ simulations were conducted. Putative Hopf bifurcations were detected by examining whether a pair of complex conjugate eigenvalues existed which were sufficiently close to the imaginary axis (in our case, the magnitude of their real component was required to be less than $10−3$⁠). 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₃ 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.

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,” or RM. In the RM, the regulatory connections found in the core metabolator between M₂ and the two genes are reversed. Thus, M₂ now activates G₁, and hence its own production, and represses G₂, and hence its own degradation (see Figure 1(d)).

In nondimensional form, the RM is identical to Eq. (1a) but with the regulation of g₁ and g₂ reversed

For the RM, we repeated the explicit dynamical analysis outlined in Sec. II. The fixed points of Eq. (18a) are identical to those of the core metabolator (⁠$g1*=g2*=12,m2*=1$⁠, and $c*=η1 + η2 + 2ρη1$⁠), but the Jacobian becomes

The characteristic equation of the RM is now

Next, we calculate the necessary condition for Eq. (19) to have a pair of purely imaginary roots (i.e., that $c0=ω2c2and c1=ω2c3$⁠, where c_i is the coefficient of the $λi$th term of the cubic polynomial in Eq. (19)). Assuming that $λ1,2=±iω$⁠, we find two conditions

Equating these two expressions for $ω2$⁠, we arrive at a quadratic polynomial in $η1$

Using Descartes' rule of signs, Eq. (21) must contain positive, real roots since $η2>0$ and the $η10$ 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,⁵ 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 2(c)). Second, for large values of ρ, unstable periodic solutions showed characteristically sudden changes in amplitude, reminiscent of relaxation oscillations commonly found in neuroscience.¹⁶ 

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 3(a). 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 3(a)). 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 3(a)), 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.⁸ 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 3(a)), 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 Sec. VII.

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 post-translational 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.³ 

In order to study the inherent time scales of the RM, we constructed a generalized model. This model was identical in structure to Eq. (14), except that the signs of the elements corresponding to the transcriptional regulation of M₂ 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₁

where $δm1,δm2,δg1$⁠, and $δg2$ correspond to perturbations in $m1,m2,g1$⁠, and g₂, respectively. Upon inspection of Eq. (22), all terms on the right hand side share a common divisor $M1,0v$ which has units of time. We can treat this shared prefactor as a natural measure of the time scale of M₁. Repeating this process for all of the rows in Eq. (15), we identified the characteristic time scale for each state space variable. Entries in rows one and two of J, corresponding to the dynamics of the two metabolites, shared common divisors $M1,0v$ and $M2,0v$⁠, respectively. Similarly, entries in the final two rows of J, corresponding to the genetic elements of the metabolator, shared common divisors $1β1$ and $1β2$ (again with units of time).

We studied how the stability of J depends on the ratio of the metabolic and genetic time scales. To do so, we generated a large number of random instances of J and calculated the real component of the leading eigenvalue $λmax$ of J (a positive value indicated instability, while a negative value indicated stability). The sampling was completed by selecting a particular choice of the ratio of metabolic and genetic time scales, and then creating a 30 × 30 grid of $θr1$ and $θr2$ values. For each point on the grid, we sampled 1000 instances of all other θ's uniformly from the interval [0, 1] for each point in the grid. The results of this simulation are shown in Figure 4.

Two major trends regarding the stability 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.¹⁸ 

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 Eq. (1a) has a unique fixed point which is stable for all parameter values (recall the end of Sec. II). Hence, the core metabolator cannot undergo Hopf bifurcation to oscillations. Then, we developed a GM model of the core metabolator Eq. (12), and we showed that this system also cannot undergo Hopf bifurcations by applying the Routh-Hurwitz criterion to the Jacobian Eq. (15) (see Sec. III B and Appendix  B).

These analytical results for the core metabolator Eq. (11) and its generalized model Eq. (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 Sec. IV, we explored two more detailed versions of the core metabolator. The first of these has an intermediary delay gene G₃ (Figure 1(b)), 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,⁸ effectively “un-lumping” M₂. 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₂ induced by the network topologies for the core metabolator, generalized model of the core metabolator, and the first “enhanced” network (Figures 1(a) and 1(b)), 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 1(c)) suggested that we explore some other four-component models with related designs that do exhibit oscillations. In particular, we studied system Eq. (18a), the reverse metabolator, in which the wiring from M₂ to the two genes is reversed compared to the core metabolator, as shown in Figure 1(d). This fundamentally altered the system dynamics so that M₂ now activates, rather than represses, the gene G₁ and hence M₂ is now self-activating. The analysis of Sec. V 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 Sec. VI): 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,¹⁴ as well as in prior work on metabolism using GM.²⁵ 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,¹¹ and their extension to GM with regulation appears straightforward.

E.R. and D.S. were supported by grants from the Office of Science (BER), U.S. Department of Energy (DE-SC0004962), and National Science Foundation NSF DMS-0602204 EMSW21-RTG Biodynamics at Boston University. T.J.K. was supported in part by NSF Grant No. DMS-1109587. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. E.R. would like to acknowledge the kind support of the Santa Fe Institute while studying the metabolator at the 2011 Summer School.

We begin by writing out the equations for the metabolator. In dimensionalized form, they are

Lets begin nondimensionalizing F₁. Let $m2=M2K1$⁠. We can rewrite Eq. (A1a) as

where $γ1=α1β1K13,g1=G1K1$⁠, and $τ=β1t$⁠.

Then, the nondimensionalized equation for G₂ immediately as follows:

where $g2=G2K1,c1=β1β2,c2=K1K2$⁠, and $γ2=α2β2K1$⁠.

For Eq. (A1c), we rewrite it by substituting those nondimensionalized variables and parameters which we have already defined and choose to scale M₁ so that $m1=M1K1$⁠. Doing so, we get

where $ρ1=Vinβ1K1,η1=v1K1β1$⁠, and $η2=v2K1β1$⁠.

Finally, the nondimensionalized equation for M₂ 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=M1+M2$⁠. Then, $Ċ=M1̇+M2̇=Vin−VoutM2$⁠. We can non-dimensionalize this system as well to find that

where $c=CK1$⁠.

For the purpose of simplifying the analysis in the article, we let $c1=c2=γ1=γ2=1$ and $ρ1=ρ2=ρ$⁠. Then, the system formed by Eqs. (A4), (A5), (A9), and (A8) is the nondimensionalized core metabolator we study in this article.

We work with the Jacobian of the generalized model of the core metabolator (15)

where we have let $βi=Liθd,iGi,0$ and $θr,i=θf,iθd,i$⁠. This means that $θr1<0,θr2>0$ and removes two otherwise free parameters from the analysis.

Adding multiples of one row to another does not change the value of the determinant. Thus, to simplify matters, we multiply row one of Eq. (B1) by $M1,0M2,0$ and add it to row 2 of Eq. (B1) to find