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

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 switch9 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.8 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 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 (G1 and G2) interconvert two substrates (M1 and M2). Gene G1 codes for the enzyme which converts M1 into M2, and gene G2 codes for the enzyme which converts M2 into M1. We explicitly assume that the rate of gene expression is substantially faster than protein translation, so that the abundances of the mRNA transcripts encoding G1 and G2 may be assumed to be in quasi-steady-state. Metabolite M2 represses the expression of gene G1 and activates the expression of gene G2. Furthermore, we assume that the rate of inflow of M1 is constant, and the rate of outflow of M2 is in proportion to its concentration. Because of the symmetry in the dynamics of M1 and M2 (both are produced and degraded by the same metabolic reactions), we replace M1 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

dg1dt=11+m22g1,
(1a)
dg2dt=m221+m22g2,
(1b)
dcdt=ρ(1m2),
(1c)
dm2dt=ρm2+η1g1(cm2)η2g2m2,
(1d)

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.

FIG. 1.

Several different metabolator designs. (a) The core metabolator, studied in Secs. II and III. Metabolite M2 represses the transcription of gene G1 and activates the transcription of gene G2. (b)Core metabolator with intermediary delay gene G3, and (c) core metabolator with un-lumped metabolites, both studied in Sec. IV. (d) Reverse metabolator, studied in Sec. V. In contrast to the core metabolator, M2 activates the transcription of G1 and represses the transcription of G2.

FIG. 1.

Several different metabolator designs. (a) The core metabolator, studied in Secs. II and III. Metabolite M2 represses the transcription of gene G1 and activates the transcription of gene G2. (b)Core metabolator with intermediary delay gene G3, and (c) core metabolator with un-lumped metabolites, both studied in Sec. IV. (d) Reverse metabolator, studied in Sec. V. In contrast to the core metabolator, M2 activates the transcription of G1 and represses the transcription of G2.

Close modal

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

[1001201012000ρ2ρ+η2η2η122ρη1η22].
(2)

The characteristic equation for J is

(1+λ)(λ3+λ2(2ρ+η1+η2+22)+λ(4ρ+3η2+η1+η1ρ2)+η1ρ2)=0.
(3)

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 ci 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

η1ρ2+η1+η2+2ρ=η1+3η2+4ρ+η1ρ2.
(4)

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

(1+ρ)η12+((1+ρ)(2+2ρ+η2)+2ρ+3η2)η1+(2+2ρ+η2)(4ρ+3η2)=0.
(5)

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

dXdt=I(X,Y)O(X,Y),dYdt=F(X,Y)G(X,Y).
(6)

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)

1X0dXdt=1X0(I(X0,Y0)I(X0,Y0)I(X,Y)O(X0,Y0)O(X0,Y0)O(X,Y)).
(7)

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

dxdt=1X0(I(X0,Y0)I(X0,Y0)I(xX0,yY0)O(X0,Y0)O(X0,Y0)O(xX0,yY0)).
(8)

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 hx|(1,1) and hy|(1,1), we find

hx|(1,1)=1X0(I(X0,Y0)θxIO(X0,Y0)θxO),hy|(1,1)=1X0(I(X0,Y0)θyIO(X0,Y0)θyO),
(9)

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 analysis7) 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 G1 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))

dG1dt=α1K12+M22β1G1.
(10)

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

dg1dt=1G1,0(α1K12+(m2M2,0)2β1g1G1,0).
(11)

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

μ(m2)=α1K12+m22M2,02α1K12+M2,02=K12+M2,02K12+m22M2,02.

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

dμdm2|m2=1=(K12+M2,02)(2m2M2,02)(K12+m22M2,02)2|m2=1=(2M2,02)(K12+M2,02)(K12+M2,02)2=2M2,02K12+M2,02.

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

dM1dt=I+G2R2(M2)G1R1(M1),dM2dt=G1R1(M1)G2R2(M2)R3(M2),dG1dt=P1(M2)D1(G1),dG2dt=P2(M2)D2(G2).
(12)

The terms R1 and R2 correspond to the reactions catalyzed by G1 and G2, respectively, while I and R3 correspond to the inflow and outflow of mass from the metabolator. The Pi and Di terms correspond to the production and degradation of proteins. Importantly, the Pi terms capture the feedback regulation of metabolite M2 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 Di are generic, monotonically increasing functions of their arguments (i.e., R1M1>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,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 Eq. (12), we find that

G1,0R1(M1,0)G2,0R2(M2,0)=I,G1,0R1(M1,0)G2,0R2(M2,0)=R3(M2,0),P1(M2,0)=D1(G1,0),P2(M2,0)=D2(G2,0).
(13)

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 M1 to M2) to be precisely balanced by the rate of reaction G2R2 (the conversion of M2 to M1) and I, the rate of inflow of M1 to the system. Similarly, G2R2 must be precisely equal to the sum of G1R1 and R3 (the outflow of M2 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)

dm1dt=1M1,0(αv+vG2,0R2(M2,0)g2G2,0R2(m2M2,0)(1+α)vG1,0R1(M1,0)g1G1,0R1(m1M1,0)),dm2dt=1M2,0((1+α)vG1,0R1(M1,0)g1G1,0R1(m1M1,0)vG2,0R2(M2,0)g2G2,0R2(m2M2,0)αvR3(M2,0)R3(m2M2,0)),dg1dt=1G1,0(L1P1(M2,0)P1(m2M2,0)L1D1(G1,0)D1(g1G1,0)),dg2dt=1G2,0(L2P2(M2,0)P2(m2M2,0)L2D2(G2,0)D2(g2G2,0)),
(14)

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

The Jacobian matrix J of the core metabolator is then

J=[θ1(1+α)vM1,0θ2vM1,0(1+α)vM1,0vM1,0θ1(1+α)vM2,0θ2vαθ3vM2,0(1+α)vM2,0vM2,00L1θf1G1,0L1θd1G1,000L2θf2G2,00L2θd2G2,0].
(15)

The parameters θi in Eq. (15) correspond to the elasticities of metabolic reactions described earlier

θ1=1G1,0R1(M1,0)((g1G1,0R1(m1M1,0))m1)|m1=g1=1,θ2=1G2,0R2(M2,0)((g2G2,0R2(m2M2,0))m2)|m2=g2=1,θ3=1R3(M2,0)(R3(m2M2,0)m2)|m2=1.
(16)

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

θf1=1P1(M2,0)(P1(m2M2,0)m2)|m2=1,θf2=1P2(M2,0)(P2(m2M2,0)m2)|m2=1,θd1=1D1(G1,0)(D1(G1,0)g1)|g1=1,θd2=1D2(G2,0)(D2(G2,0)g2)|g2=1.
(17)

For the generalized model of the core metabolator, θf1<0,θf2>0 since M2 inhibits G1 and activates G2. 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.2 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 G3 and representing LacI from Ref. 8, might act as an intermediary delay in the negative feedback from M2 to G1 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 M2 into its appropriate constitutive molecules (now labeled M2 and M3) in Figure 1(c). In the revised topology, the enzyme transcribed from gene G2 reconverts M3 to M1, but the regulation of G1 and G2 remains dependent on M2. 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 103 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 103). 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 G3 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 M2 and the two genes are reversed. Thus, M2 now activates G1, and hence its own production, and represses G2, and hence its own degradation (see Figure 1(d)).

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

dg1dt=m221+m22g1,
(18a)
dg2dt=11+m22g2,
(18b)
dcdt=ρ(1m2),
(18c)
dm2dt=ρm2+η1g1(cm2)η2g2m2.
(18d)

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

[1001201012000ρ2ρ+η2η2η122ρη1η22].

The characteristic equation of the RM is now

(λ+1)(λ3+λ2(2+2ρ+η1+η22)+λ(η1η2+η1ρ2)+η1ρ2)=0.
(19)

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 ci is the coefficient of the λith term of the cubic polynomial in Eq. (19)). Assuming that λ1,2=±iω, we find two conditions

ω2=η1ρ2+2ρ+η1+η2,ω2=η1(1+ρ)η22.
(20)

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

η12(1+ρ)+η1(ρ(η2+2ρ+2)+2)η2(η2+2ρ+2)=0.
(21)

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,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 2(c)). Second, for large values of ρ, unstable periodic solutions showed characteristically sudden changes in amplitude, reminiscent of relaxation oscillations commonly found in neuroscience.16 

FIG. 2.

(a) Continuation of a limit cycle from a Hopf bifurcation in the RM (18). For small values of ρ, two limit cycles exist (a large amplitude unstable cycle, and a small amplitude stable cycle), while for large values of ρ, a single, unstable limit cycle exists. (b) For ρ = 5.5 (black), 5.6 (red), 5.8 (green), and 6.0 (blue), the unstable limit cycle exhibits increasingly sudden changes in abundance of M2. (c) For ρ=5.15, the two limit cycles (stable, black; unstable, red) exhibit smooth changes in amplitude. For all cases, η1=1.27,η2=7.92.

FIG. 2.

(a) Continuation of a limit cycle from a Hopf bifurcation in the RM (18). For small values of ρ, two limit cycles exist (a large amplitude unstable cycle, and a small amplitude stable cycle), while for large values of ρ, a single, unstable limit cycle exists. (b) For ρ = 5.5 (black), 5.6 (red), 5.8 (green), and 6.0 (blue), the unstable limit cycle exhibits increasingly sudden changes in abundance of M2. (c) For ρ=5.15, the two limit cycles (stable, black; unstable, red) exhibit smooth changes in amplitude. For all cases, η1=1.27,η2=7.92.

Close modal

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.

FIG. 3.

(a) Two-parameter bifurcation diagram of the RM (18). As ρ and η1 are varied, a curve of Hopf bifurcations (red) and a fold of limit cycles (black) collide at a Bautin point. Insets (B–D) correspond to the dynamics of the RM when η1=2.435, in between the two Bautin points. Phase portraits correspond to the dynamics of metabolite M2 (vertical axis) and total metabolite C=M1+M2 (horizontal axis), with green dots indicating initial conditions and red dots indicating final positions. All simulations were run for a time interval of 1000. Black curves indicate trajectories approaching a fixed point, blue curves correspond to trajectories approaching a stable limit cycle, and red curves correspond to trajectories spiralling away from an unstable limit cycle. (b) For small ρ, the RM exhibits dampled, spiralling oscillations to a stable fixed point. (c) For ρ=2.55, two distinct limit cycles exist. (d) For large values of ρ, the stable limit cycle in (c) undergoes a Hopf bifurcation and creates a stable fixed point, surrounded by a large amplitude unstable limit cycle.

FIG. 3.

(a) Two-parameter bifurcation diagram of the RM (18). As ρ and η1 are varied, a curve of Hopf bifurcations (red) and a fold of limit cycles (black) collide at a Bautin point. Insets (B–D) correspond to the dynamics of the RM when η1=2.435, in between the two Bautin points. Phase portraits correspond to the dynamics of metabolite M2 (vertical axis) and total metabolite C=M1+M2 (horizontal axis), with green dots indicating initial conditions and red dots indicating final positions. All simulations were run for a time interval of 1000. Black curves indicate trajectories approaching a fixed point, blue curves correspond to trajectories approaching a stable limit cycle, and red curves correspond to trajectories spiralling away from an unstable limit cycle. (b) For small ρ, the RM exhibits dampled, spiralling oscillations to a stable fixed point. (c) For ρ=2.55, two distinct limit cycles exist. (d) For large values of ρ, the stable limit cycle in (c) undergoes a Hopf bifurcation and creates a stable fixed point, surrounded by a large amplitude unstable limit cycle.

Close modal

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.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 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.3 

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 M2 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 M1

d(δm1)dt=θ1(1+α)vM1,0δm1,0+θ2vM1,0δm2+(1+α)vM1,0δg1+vM1,0δg2,
(22)

where δm1,δm2,δg1, and δg2 correspond to perturbations in m1,m2,g1, and g2, 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 M1. 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.

FIG. 4.

Dependence of stability of the RM on metabolic and genetic time scales (red = lower stability). Instability in the reverse metabolator increases as the genetic time scale becomes faster. For simplicity, we assumed that M1=M2,β1=β2, and v = 1, so that there existed only two time scales, one metabolic and one genetic. Additional numerical investigations indicate results are identical when this assumption is relaxed (results not shown).

FIG. 4.

Dependence of stability of the RM on metabolic and genetic time scales (red = lower stability). Instability in the reverse metabolator increases as the genetic time scale becomes faster. For simplicity, we assumed that M1=M2,β1=β2, and v = 1, so that there existed only two time scales, one metabolic and one genetic. Additional numerical investigations indicate results are identical when this assumption is relaxed (results not shown).

Close modal

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.18 

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 G3 (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,8 effectively “un-lumping” M2. 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 M2 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 M2 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 M2 now activates, rather than represses, the gene G1 and hence M2 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,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.

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

dG1dt=α1K12+M22β1G1F1,
(A1a)
dG2dt=α2M22K22+M22β2G2F2,
(A1b)
dM1dt=Vinv1G1M1+v2G2M2F3,
(A1c)
dM2dt=VoutM2+v1G1M1v2G2M2F4.
(A1d)

Lets begin nondimensionalizing F1. Let m2=M2K1. We can rewrite Eq. (A1a) as

1β1dG1dt=α1β1K12(1+m2)G1,
(A2)
1K1β1dG1dt=α1β1K13(1+m22)G1K1,
(A3)
dg1dτ=γ11+m22g1,
(A4)

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

Then, the nondimensionalized equation for G2 immediately as follows:

c1dg2dτ=γ2c22m221+c22m22g2,
(A5)

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 M1 so that m1=M1K1. Doing so, we get

β1K1dm1dτ=Vinv1g1m1K12+v2g2m2K12,
(A6)
dm1dτ=ρ1η1g1m1+η2g2m2,
(A7)

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

Finally, the nondimensionalized equation for M2 follows simply and is

dm2dτ=ρ2m2+η1g1m1η2g2m2,
(A8)

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̇=VinVoutM2. We can non-dimensionalize this system as well to find that

dcdτ=ρ1ρ2m2,
(A9)

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)

JλI=[θ1(1+α)vM1,0λθ2vM1,0(1+α)vM1,0vM1,0θ1(1+α)vM2,0θ2vαθ3vM2,0λ(1+α)vM2,0vM2,00β1θr1β1λ00β2θr20β2λ],
(B1)

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

JλI=[θ1(1+α)vM1,0λθ2vM1,0(1+α)vM1,0vM1,0λM1,0M2,0αθ3vM2,0λ000β1θr1β1λ00β2θr20β2λ].
(B2)

Now expanding the determinant of Eq. (B2) along the first column, we find

|JλI|=(θ1(1+α)vM1,0λ)(αθ3vM2,0λ)(β1λ)(β2λ)+M1,0M2,0λ|θ2vM1,0(1+α)vM1,0vM1,0β1θr1β1λ0β2θr20β2λ|.
(B3)

Expanding again, we find

|JλI|=(θ1(1+α)vM1,0λ)(αθ3vM2,0λ)(β1λ)(β2λ)+M1,0M2,0λ(vM1,0(β2θr2)(β1+λ)+(β2λ)(θ2vM1,0(β1λ)+β1θr1(1+α)vM1,0)).
(B4)

Expanding this expression a final time and grouping terms, we find the coefficients of the characteristic equation for the Jacobian J of the core metabolator are (where ai corresponds to the coefficient of the λith term)

a4=1,a3=β1+β2+(1+α)vθ1M1,0+vθ2+αvθ3M2,0,a2=β1β2+α(1+α)v2θ1θ3M1,0M2,0+(β1+β2)((1+α)vθ1M1,0+vθ2+αvθ3M2,0)+β2vM2,0(θr2)+β1vM2,0((1+α)θr1),a1=β1β2v((1+α)θ1M1,0+αθ3+θr2+θ2(1+α)θr1M2,0)+(β1+β2)α(1+α)v2θ1θ3M1,0M2,0,a0=α(1+α)β1β2v2θ1θ3M1,0M2,0.
(B5)

The Routh Hurwitz criterion states that all of the eigenvalues of a 4 × 4 matrix have negative real parts when the following conditions on its characteristic polynomial are satisfied

an>0,n=0,1,2,3,4,a2a3>a1a4,a1a2a3>a12a4+a0a32.
(B6)

First, note that the condition an>0 is satisfied. Next, we calculate the quantity a2a3

a2a3=(β1β2)(β1+β2+(1+α)vθ1M1,0+vθ2+αvθ3M2,0)+α(1+α)v2θ1θ3M1,0M2,0(β1+β2+(1+α)vθ1M1,0+vθ2+αvθ3M2,0)+(β1+β2)((1+α)vθ1M1,0+vθ2+αvθ3M2,0)(β1+β2+(1+α)vθ1M1,0+vθ2+αvθ3M2,0)+(β2vθr2M2,0(1+α)vβ1θr1M2,0)(β1+β2+(1+α)vθ1M1,0+vθ2+αvθ3M2,0).
(B7)

Now, a1a4 = a1 because a4=1. Explicitly calculating a2a3a4

a2a3a1a4=β1β2(β1+β2)+α(1+α)v2θ1θ3M1,0M2,0((1+α)vθ1M1,0+vθ2+αvθ3M2,0)+(β1+β2)((1+α)vθ1M1,0+vθ2+αvθ3M2,0)(β1+β2+(1+α)vθ1M1,0+vθ2+αvθ3M2,0)+(β2vθr2M2,0(1+α)vβ1θr1M2,0)((1+α)vθ1M1,0+vθ2+αvθ3M2,0)+β22vθr2M2,0(1+α)vβ12θr1M1,0.
(B8)

All the terms are positive, so the condition a2a3>a1a4 is satisfied.

Now, we rewrite the final condition as a1(a2a3a1a4)a0a32>0. Calculating a0a32, we find

a0a32=(β1+β2)2α(1+α)β1β2v2θ1θ3M1,0M2,0+2(β1+β2)((1+α)vθ1M1,0+vθ2+αvθ3M2,0)(α(1+α)β1β2v2θ1θ3M1,0M2,0)+((1+α)vθ1M1,0+vθ2+αvθ3M2,0)2(α(1+α)β1β2v2θ1θ3M1,0M2,0).
(B9)

Also, rewriting a1 and highlighting the two key components

a1=β1β2v((1+α)θ1M1,0+αθ3+θr2+θ2(1+α)θr1M2,0)Term1+(β1+β2)α(1+α)v2θ1θ3M1,0M2,0Term2.
(B10)

Now, note that the first line of Eq. (B8) multiplied by term 2 of Eq. (B10) gives the first term of Eq. (B9). Similarly, we can isolate a term like 2(β1+β2)(α(1+α)v2θ1θ3M1,0M2,0) from the third line of Eq. (B8) and multiply it by term 1 of Eq. (B10) to obtain the second term of Eq. (B9), so that all the remaining terms of the product of the two lines are positive. Finally, we can multiply the second line of Eq. (B8) by term 1 of Eq. (B10) to give the third term of Eq. (B9) and have all leftover terms be positive. The remaining terms are all positive, proving that the Routh-Hurwitz criteria are satisfied and all of the eigenvalues of the Jacobian have negative real part.

For the metabolator in Figure 1(c), we used the GM

dM1dt=I+G2R2(M2)G1R1(M1),dM2dt=G1R1(M1)G2R2(M2)R3(M2),dG1dt=P1(G3)D1(G1),dG2dt=P2(M2)D2(G2),dG3dt=P3(M2)D3(G3).
(C1)

The Jacobian of Eq. (C1) is

J=[θ1(1+α)vM1,0θ2vM1,0(1+α)vM1,0vM1,00θ1(1+α)vM2,0θ2vαθ3vM2,0(1+α)vM2,0vM2,0000β10β1θr10β2θr20β200β3θr300β3].
(C2)

We assume again that θr1<0 and all other elasticities are greater than zero.

For the metabolator in Figure 1(d), we use the GM

dM1dt=I+G2R2(M3)G1R1(M1),dM2dt=G1R1(M1)R4(M2),dM3dt=R4(M2)G2R2(M3)R3(M3),dG1dt=P1(M2)D1(G1),dG2dt=P2(M2)D2(G2).
(C3)

The Jacobian of Eq. (C3) is

J=[θ1(1+α)vM1,00θ2vM1,0(1+α)vM1,0vM1,0θ1(1+α)vM2,0θ4(1+α)vM2,00(1+α)vM2,000θ4(1+α)vM3,0θ2vαθ3vM3,00vM3,00β1θr10β100β2θr200β2].
(C4)

Above, θ4 corresponds to the elasticity of reaction R4(M2). We again assume again that θr1<0 and all other elasticities are greater than zero.

1.
U.
Alon
,
An Introduction to Systems Biology: Design Principles of Biological Circuits, Chapman & Hall/CRC Mathematical and Computational Biology
(
Chapman and Hall/CRC
,
2006
).
2.
K. J.
Astrom
and
R. M.
Murray
,
Feedback Systems: An Introduction for Scientists and Engineers
(
Princeton University Press
,
2008
).
3.
A.
Bar-Even
,
E.
Noor
,
Y.
Savir
,
W.
Liebermeister
,
D.
Davidi
,
D. S
Tawfik
, and
R.
Milo
, “
The moderately efficient enzyme: Evolutionary and physicochemical trends shaping enzyme parameters
,”
Biochemistry
50
(
21
),
4402
4410
(
2011
).
4.
K. A.
Dill
,
K.
Ghosh
, and
J. D.
Schmit
, “
Physical limits of cells and proteomes
,”
Proc. Natl. Acad. Sci. U.S.A.
108
(
44
),
17876
17882
(
2011
).
5.
E. J.
Doedel
,
A. R.
Champneys
,
T.
Fairgrieve
,
Y.
Kuznetsov
,
B.
Oldeman
,
R.
Paffenroth
,
B.
Sandstede
,
X.
Wang
, and
C.
Zhang
, “
Auto-07p: Continuation and bifurcation software for ordinary differential equations
,” (
2007
).
6.
M. B.
Elowitz
and
S.
Leibler
, “
A synthetic oscillatory network of transcriptional regulators
,”
Nature
403
(
6767
),
335
338
(
2000
).
7.
D. A.
Fell
, “
Metabolic control analysis: A survey of its theoretical and experimental development
,”
Biochem. J.
286
(
Pt 2
),
313
330
(
1992
).
8.
E.
Fung
,
W. W.
Wong
,
J. K.
Suen
,
T.
Bulter
,
S.-g.
Lee
, and
J. C.
Liao
, “
A synthetic gene-metabolic oscillator
,”
Nature
435
(
7038
),
118
122
(
2005
).
9.
T. S.
Gardner
,
C. R.
Cantor
, and
J. J.
Collins
, “
Construction of a genetic toggle switch in Escherichia coli
,”
Nature
403
(
6767
),
339
342
(
2000
).
10.
E.
Gehrmann
,
C.
Gläßer
,
Y.
Jin
,
B.
Sendhoff
,
B.
Drossel
, and
K.
Hamacher
, “
Robustness of glycolysis in yeast to internal and external noise
,”
Phys. Rev. E
84
(
2
),
021913
(
2011
).
11.
D.
Girbig
,
J.
Selbig
, and
S.
Grimbs
, “
A MATLAB toolbox for structural kinetic modeling
,” in
Bioinformatics
(Oxford, England,
2012
).
12.
S.
Grimbs
,
J.
Selbig
,
S.
Bulik
,
H.-G.
Holzhutter
, and
R.
Steuer
, “
The stability and robustness of metabolic states: Identifying stabilizing sites in metabolic networks
,”
Mol. Syst. Biol.
3
, Article No.
146
(
2007
).
13.
T.
Gross
and
U.
Feudel
, “
Generalized models as a universal approach to the analysis of nonlinear dynamical systems
,”
Phys. Rev. E
73
,
016205
016214
(
2006
).
14.
T.
Gross
,
L.
Rudolf
,
S. A.
Levin
, and
U.
Dieckmann
, “
Generalized models reveal stabilizing factors in food webs
,”
Science
325
(
5941
),
747
750
(
2009
).
15.
M.
Heinemann
and
U.
Sauer
, “
Systems biology of microbial metabolism
.”
Curr. Opin. Microbiol.
13
(
3
),
337
343
(
2010
).
16.
E. M.
Izhikevich
,
Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting
(
MIT
,
2006
).
17.
N.
Jamshidi
and
B.
ØPalsson
, “
Formulating genome-scale kinetic models in the post-genome era
,”
Mol. Syst. Biol.
4
, Article No.
171
(
2008
).
18.
E.
Klipp
, “
Timing matters
,”
FEBS Lett.
583
(
24
),
4013
4018
(
2009
).
19.
J. D.
Orth
,
I.
Thiele
, and
B.
ØPalsson
, “
What is flux balance analysis?
Nat. Biotechnol.
28
(
3
),
245
248
(
2010
).
20.
J. A.
Papin
,
T.
Hunter
,
B. O.
Palsson
, and
S.
Subramaniam
, “
Reconstruction of cellular signalling networks and analysis of their properties
,”
Nat. Rev. Mol. Cell Biol.
6
(
2
),
99
111
(
2005
).
21.
F.
Picard
,
C.
Dressaire
,
L.
Girbal
, and
M.
Cocaign-Bousquet
, “
Examination of post-transcriptional regulations in prokaryotes by integrative biology
,”
C. R. Biol.
332
(
11
),
958
973
(
2009
).
22.
O.
Purcell
,
N. J.
Savery
,
C. S.
Grierson
, and
M.
di Bernardo
, “
A comparative analysis of synthetic genetic oscillators
,”
J. R. Soc. Interface
7
(
52
),
1503
1524
(
2010
).
23.
E.
Reznik
and
D.
Segrè
, “
On the stability of metabolic cycles
,”
J. Theor. Biol.
266
(
4
),
536
549
(
2010
).
24.
B.
Schwanhäusser
,
D.
Busse
,
Na.
Li
,
G.
Dittmar
,
J.
Schuchhardt
,
J.
Wolf
,
W.
Chen
, and
M.
Selbach
, “
Global quantification of mammalian gene expression control
,”
Nature
473
(
7347
),
337
342
(
2011
).
25.
R.
Steuer
,
T.
Gross
,
J.
Selbig
, and
B.
Blasius
, “
Structural kinetic modeling of metabolic networks
,”
Proc. Natl. Acad. Sci. U.S.A.
103
,
11868
11873
(
2006
).
26.
R.
Steuer
,
A.
Nunes Nesi
,
A. R.
Fernie
,
T.
Gross
,
B.
Blasius
, and
J.
Selbig
, “
From structure to dynamics of metabolic pathways: Applications to the plant mitochondrial TCA cycle
,”
Bioinformatics
23
,
1378
1385
(
2007
).