A Appendix: Formulation and Non-dimensionalization of the Core Metabolator

The Dynamics of Hybrid Metabolic-Genetic Oscillators

Abstract

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

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

I Introduction

Synthetic biology is in need of theoretical tools to design circuits with a high probability of exhibiting specific, user-defined dynamical behaviors. Conventionally, efforts in this regard have focused on explicit kinetic formulations of the dynamics of a circuit. Indeed, deterministic models of the toggle switch [9] and the repressilator [6], two early examples of successful synthetic circuits, directly predicted the dynamical behavior of these constructs and likely influenced the final design of the circuits themselves. Perhaps the most prominent challenge in assembling these models was a fundamental lack of knowledge regarding kinetic details, including the mechanistic form of the rate laws, as well as the precise values of rate constants themselves.

The primary focus of this work is to develop a modeling framework for understanding the dynamics of a biological circuit for which we have little detailed knowledge of the true kinetics. In particular, we will try to derive generic results dependent exclusively on the “topology” of the circuit, rather than its detailed kinetics. To do so, we will focus our efforts on a particular system known as the metabolator, a synthetic circuit constructed in Escherichia coli K12 which exhibits limit cycle oscillations. At a coarse level of description, the metabolator consists of two metabolites which are interconverted between each other by two enzymes. One of the metabolites in the circuit acts as a regulator, repressing the expression of the enzyme catalyzing the metabolite’s production and activating the expression of the enzyme catalyzing its degradation. In Ref. [8], the authors observed that in conditions of sufficiently high metabolic inflow to the metabolator, a gene encoding one of the regulated enzymes exhibited sustained oscillations. The metabolator is particularly interesting because it remains (to our knowledge) the only example of a synthetic circuit containing a combination of both metabolic and genetic (transcriptional/translational) components. In all other cases, the building blocks of circuits have been entirely transcriptional and translational. The “hybrid” nature of the metabolator is particularly interesting because metabolism is well-known to take place at a faster rate than transcription and translation. The role such separation of time scales actually plays in the dynamics of the metabolator is unresolved, although it has been observed to be quite in important in other metabolic contexts [17].

One of our primary approaches to studying the metabolator is based on a non-dimensionalization technique for studying dynamical systems commonly referred to as generalized modeling (GM) [25]; [13]; [12]; [26]; [10]. In GM, a change of variables is applied to a dynamical system so that many of the otherwise unconstrained parameters in the system (e.g., rate constants) are replaced by “elasticity” parameters with well-defined ranges. Then, the dynamics of the system can be studied around an arbitrary and potentially unknown steady-state. GM’s conclusions are invariant to the particular choice of rate law, and depend only on the sensitivity of the rate law at the system’s equilibrum. This enables the study of a topological model of the metabolator, in which no assumptions are generally made on the explicit form of the rate law, nor the true values of the parameters in the system. In practice, many GM studies have focused on using large-scale numerical simulations to tease out statistical properties relating elasticities to stability properties of the system. In this work, we supplement such numerical work with complementary analytical studies, and highlight the additional insight that can be gained through explicit calculations. We complement our use of GM with explicit dynamical modeling of the metabolator, and highlight the complementary conclusions that can be drawn from the two techniques (e.g. local vs. global bifurcation phenomena when applying GM or explicit dynamical modeling, respectively). Importantly, we demonstrate how the results from GM can prompt further study with explicit modeling, and vice versa.

The main results of the work are as follows. First, a simple, “core” model of the metabolator (two genes, two metabolites) originally illustrated in Ref. [8] is investigated using a conventional dynamical model. We find that the system lacks the capability to pass through a Hopf bifurcation. Second, using GM, we analytically prove that any dynamical model with a topology analogous to the core metabolator is incapable of Hopf bifurcations. Concluding that additional topological elements are necessary for a Hopf bifurcation, we use numerical simulations with GM to progressively restore elements of the metabolator model from Ref. [8] which we omitted in our core model. Doing so, we reveal the crucial topological components of the metabolator which endow it with the ability to oscillate. Third, we investigate the possibility of oscillations in an explicit model of a simply “re-wired” version of the metabolator in which the regulatory connections are reversed (the reverse metabolator, “RM”). We show that the RM is capable of sustained oscillations by both a local Hopf bifurcation as well as global fold of limit cycles. Fourth, prompted by findings in the explicit model and using the GM tools we developed earlier, we study the role that metabolic and genetic time scales play in generating instability and potential oscillations in the RM in vivo.

Figure 1: Several different metabolator designs. (a) The core metabolator, studied in Sections II-III. Metabolite represses the transcription of gene and activates the transcription of gene . (b) Core metabolator with intermediary delay gene , and (c) Core metabolator with un-lumped metabolites, both studied in Section IV. (d) Reverse metabolator, studied in Section V. In contrast to the core metabolator, activates the transcription of and represses the transcription of .

Ii The Core Metabolator

To begin to understand the fundamental components which generate sustained oscillations in the metabolator, we developed a four-dimensional dynamical model of the system. Our particular choice of design was motivated by the schematic provided in Ref. [8] and re-illustrated in 1 A. Two enzymes ( and ) interconvert two substrates ( and ). Gene codes for the enzyme which converts into , and gene codes for the enzyme which converts into . We explicitly assume that the rate of gene expression is substantially faster than protein translation, so that the abundances of the mRNA transcripts encoding and may be assumed to be in quasi-steady-state. Metabolite represses the expression of gene and activates the expression of gene . Furthermore, we assume that the rate of inflow of is constant, and the rate of outflow of is in proportion to its concentration. Because of the symmetry in the dynamics of and (both are produced and degraded by the same metabolic reactions), we replace with the total amount of substrate in the system . 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

(1a)
(1b)
(1c)
(1d)

where the parameters are assumed to be strictly positive. See Appendix A for a derivation of (1). We emphasize that the core metabolator serves as a jumping off-point: observing oscillations here would suggest the presence of oscillations in more detailed models capturing the true topological elements of the in vivo metabolator. On the other hand, not observing oscillations in this simplified model implies that one needs to examine a model with more components in order to observe Hopf bifurcations, see Section 4.

The fixed point of (1) is and . The Jacobian , evaluated at the fixed point, becomes

(2)

The characteristic equation for is

(3)

In order to identify potential candidates for a Hopf bifurcation, we make use of a constraint relating the coefficients of (3) which must be satisfied in order for (3) to have a pair of purely imaginary roots. First, we note that (3) has a root at . Then, we can substitute into the cubic polynomial in (3) we find that these conditions are , where is the coefficient of the th term of the cubic polynomial in parentheses in (3). After substituting the appropriate coefficients into these conditions and some algebra, we find

(4)

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

(5)

Now the outcome becomes plainly clear: because (5) only contains positive coefficients, Descartes’ Rule of Signs immediately proves that it cannot have positive real roots. Thus, there exists no combination of which will yield a pair of purely imaginary eigenvalues, and the system does not exhibit a Hopf bifurcation.

Iii Generalized Modeling of the Core Metabolator

Can we say anything more general about the oscillatory capabilities of the core metabolator illustrated in Figure 1A? In the previous section, we showed that one particular realization of this schematic model did not exhibit a Hopf bifurcation, but it remains unclear whether this was a quality of our choice of model. One may naturally wonder if choosing different kinetic rate laws for the metabolic reactions, or different regulatory kinetics for the feedback regulation, would yield qualitatively different behavior like oscillations. In this section, we propose a method for addressing these concerns by applying a technique known as generalized modeling (GM).

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.

iii.1 Generalized Modeling

First, we present the normalization procedure upon which GM is based, known as generalized modeling, see for example Ref. [13]. Consider a system of ordinary differential equations,

(6)

Here, and are both functions of time. The stability of (6) around its equilibrium is determined by the eigenvalues of the Jacobian . Here, we approach the problem of calculating after first executing a particular change of variables central to GM. Assume that there exists an equilibrium and normalize (6) by its steady state concentration (assuming ):

(7)

Letting so that the variables are nondimensionalized and the equilibrium is at , we can write

(8)

To find the components of the Jacobian for , we simply differentiate each term by each independent variable and evaluate at the equilibrium . If we let and write the components of the Jacobian as and , we find

(9)

where and . By repeating the calculations outlined above for , we can calculate the remaining entries of the Jacobian.

We argue that these ’s (herein referred to as elasticities and similar to the elasticities from metabolic control analysis [7]) have well-defined and limited ranges which makes performing calculations with them much easier than with similar calculations with the original parameters. To demonstrate the utility of elasticities, we will reformulate the dynamics of 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 (23a)),

(10)

We let and (where and are the steady-state concentrations of and , respectively) . Then, we can write

(11)

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

To calculate the elasticity, we take the derivative of with respect to , evaluated at ,

Upon inspection, this elasticity can only take values in . When , the elasticity approaches negative two and the production of becomes very sensitive to small changes in the abundance . In contrast, when , the elasticity approaches zero and the production of is insensitive to changes in .

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.

iii.2 Generalized Model of the Core Metabolator

We can use the aforementioned normalization procedure to develop a generalized model of the core metabolator. As we will show, the generalized model enables us to draw conclusions about a topological model of the core metabolator, while relaxing many of the assumptions we made in our explicit model of the core metabolator in Section 2. We begin with the following “generalized” form of the dynamics:

(12)

The terms and correspond to the reactions catalyzed by and , respectively, while and correspond to the inflow and outflow of mass from the metabolator. The and terms correspond to the production and degradation of proteins. Importantly, the terms capture the feedback regulation of metabolite on the two enzymes.

We make a few assumptions regarding the dynamics of the core metabolator in our generalized model. First, we assume and are generic, monotonically increasing functions of their arguments (i.e. ). Many frequently used kinetic rate laws, including Hill, Michaelis-Menten, and mass-action kinetics satisfy this assumption, among many other possible forms of kinetics. We also retain the earlier assumption that the rates of metabolic reactions are linear with respect to the amount of gene product. The structure of (12) corresponds precisely to the topological structure of the core metabolator in Figure 1A.

Before we proceed to calculating the Jacobian of (12), we briefly describe a simple method for calculating the prefactors (those parameters which were not elasticities, i.e. ) which appear in the Jacobian and were described in the prior section. In prior work on GM [25], it has been shown that for metabolic networks, these prefactors correspond precisely to the rates of metabolic reactions at equilibrium. As we show below, these rates can be calculated by enforcing that at equilibrium the fluxes of reactions flowing into and out of each metabolite are balanced [19]. Furthermore, for the non-metabolic components of the metabolator, we can apply a similar balancing procedure by noting that at equilibrium, the production and degradation rates for each protein must be equal. Applying these balancing conditions to (12), we find that

(13)

Here, and indicate the steady state concentrations of metabolites 1 and 2, and gene products 1 and 2, respectively. The first two constraints in (13) force the rate of reaction (corresponding to the conversion of to ) to be precisely balanced by the rate of reaction (the conversion of to ) and , the rate of inflow of to the system. Similarly, must be precisely equal to the sum of and (the outflow of from the system). The final two conditions impose that, for each gene, the rates of production and degradation must be balanced.

Because of the interdependency of all the prefactors in (13), it becomes useful to explicitly introduce the parameter . We also introduce a free parameter by setting . From the equilibrium condition (13), it follows that . Finally, we let . We can now rescale (12) in the same manner as (8):

(14)

where .

The Jacobian matrix of the core metabolator is then

(15)

The parameters in (15) correspond to the elasticities of metabolic reactions described earlier:

(16)

Here, follows directly from the earlier assumption that all are monotonically increasing functions of their arguments. Similarly, the parameters in (15), corresponding to elasticities of the regulatory reactions, can be calculated by:

(17)

For the generalized model of the core metabolator , since inhibits and activates . Furthermore, since an increase in the abundance of a gene product leads to an increase in its rate of degradation.

The utility of reformulating the core metabolator (1) with the GM in (14) is now plainly clear. Many of the prior assumptions made on the rate laws in our explicit dynamical model (1) have been relaxed. Our sole assumption is that the rate laws are monotonic, so that it is always the case that the signs of the elasticities are known. Now, by analytically studying the eigenvalues of , we can understand the stability properties of the core metabolator across the entire spectrum of possible parameter values and rate laws.

In Appendix B, we show that (15) 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 (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 1A, 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, that 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 the following section.

Iv Hopf Bifurcations in More Complete Generalized Models of the Metabolator

The analysis in Sections 2 and 3 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 Figure 1B-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 and representing LacI from Ref. [8], might act as an intermediary delay in the negative feedback from to and consequently induce oscillations (see Figure 1B). 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 into its appropriate constitutive molecules (now labeled and ) in Figure 1C). In the revised topology, the enzyme transcribed from gene reconverts to , but the regulation of and remains dependent on . 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 and for all other elasticities. Since the parameters and remained unconstrained, we decided to sample these parameters across seven orders of magnitude (from to 1/sec, assuming ). For each instance of our simulation, the eigenvalues of were calculated. The process was repeated 20,000 times for each choice of and . In total, over 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 ). 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 had no bearing on the appearance of a Hopf bifurcation. In fact, the model remained entirely stable, failing to exhibit a single instance of an eigenvalue with positive real part. On the other hand, the refined three-metabolite model did exhibit instability and Hopf bifurcations. These results suggest that the precise nature of the regulation is crucial for the metabolator to exhibit oscillations. In order for a Hopf bifurcation to be possible, the regulatory metabolite could not be repressing the expression of its own degrading enzyme. This delayed positive feedback appears to be a crucial element endowing the metabolator with the potential to undergo a Hopf bifurcation.

V Reverse Metabolator

The difficulty in constructing synthetic circuits with poorly studied components often motivates synthetic biologists to construct several alternative designs. Frequently, these alternative designs are simply composed of rewired versions of the same regulatory components. In light of this, we investigated the dynamical behavior of one such re-wired version of the core metabolator, which we dub the reverse metabolator (RM). In the RM, the regulatory connections found in the core metabolator between and the two genes are reversed. Thus, now activates , and hence its own production, and represses , and hence its own degradation (see Figure 1D). We selected this design precisely because synthetic biologists frequently build a number of alternative constructs in the course of circuit assembly, in order to increase the odds of identifying a succesful design.

In nondimensional form, the RM is identical to (1) but with the regulation of and reversed:

(18a)
(18b)
(18c)
(18d)

v.1 Hopf Bifurcations in the Reverse Metabolator

For the RM, we repeated the explicit dynamical analysis outlined in Section 2. The fixed points of (18) are identical to those of the core metabolator ( and ), but the Jacobian becomes

The characteristic equation of the RM is now

(19)

Next, we calculate the necessary condition for (19) to have a pair of purely imaginary roots (i.e. that , where is the coefficient of the th term of the cubic polynomial in (19)). Assuming that , we find two conditions:

(20)

Equating these two expressions for , we arrive at a quadratic polynomial in

(21)

Using Descartes’ rule of signs, (21) must contain positive, real roots since and the term is negative. Therefore, there exists an open set of points 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 , and can be identified.

Figure 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 . (C) For , the two limit cycles (stable, black; unstable, red) exhibit smooth changes in amplitude. For all cases, .

We used AUTO to study the behavior of the RM’s limit cycles in different regions of parameter space. As shown in Figure 2, the limit cycles exhibit two notable features depending on the magnitude of inflow into the RM, characterized by the parameter . First, for small , two distinct periodic orbits exist: a small amplitude, stable limit cycle and a large amplitude unstable limit cycle. These two cycles annihilate at a fold of limit cycles. Both limit cycles display relatively smooth changes in amplitude (Figure 2C). Second, for large values of , unstable periodic solutions showed characteristically sudden changes in amplitude, reminiscent of relaxation oscillations commonly found in neuroscience [16].

v.2 Criticality of the Hopf Bifurcation

Figure 3: (A) Two-parameter bifurcation diagram of the RM (18). As and are varied, a curve of Hopf bifurcations (red) and a fold of limit cycles (black) collide at a Bautin point. Insets (B), (C), and (D) correspond to the dynamics of the RM when , in between the two Bautin points. Phase portraits correspond to the dynamics of metabolite (vertical axis) and total metabolite (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 , 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.

We investigate the stability of the limit cycles generated by the Hopf bifurcations in the RM. Limit cycles arising from supercritical Hopf bifurcations are stable to small perturbations, while those generated by subcritical Hopf bifurcations are not. To explore the criticality of the observed Hopf bifurcations, a two-parameter bifurcation analysis of the RM was conducted in AUTO, shown in Figure 3A. As the two parameters and were varied, the criticality of the Hopf bifurcation was tracked using the magnitude of the leading Floquet multiplier. We observed that as the value of was decreased, the Hopf bifurcation switched from being subcritical to being supercritical. As the value of was decreased further, a second transition in criticality was observed, and the Hopf bifurcations returned to being subcritical.

The point at which a Hopf bifurcation switches criticality is called a Bautin bifurcation. This higher-order bifurcation is associated with the collision of a fold of limit cycles with a manifold of Hopf bifurcations. Using AUTO, we were able to follow the fold of limit cycles as it branched off from one Bautin bifurcation, tracked closely with the manifold of Hopf bifurcations, and reintersected with the Hopf curve at the second Bautin point (Figure 3A). The fold of limit cycles is the boundary between two regions of parameter space where distinct qualitative dynamical behaviors of occur. On one side of the fold (in this case, for Region B in Figure 3), no sustained oscillations exist. On the other side of the fold (Region C in Figure 3), a pair of limit cycles appears. For the RM, this pair contained a small amplitude stable limit cycle and a large amplitude unstable limit cycle.

The existence of a Bautin bifurcation indicates there are two distinct routes to oscillation in the reverse metabolator. First, by decreasing for a fixed (moving from point D to point C in Figure 3A), the reverse metabolator passes through a supercritical Hopf bifurcation. In this case, oscillations appear at infinitesimal amplitude, and their amplitude grows as is decreased further. This type of bifurcation is identical to the one observed in the in vivo metabolator [8]. A second, and potentially more interesting, route to oscillation appears if the reverse metabolator passes through a fold of limit cycles. As the value of is increased for fixed (moving from point B to C in Figure 3A), a finite-amplitude stable limit cycle appears, surrounded by a larger amplitude unstable limit cycle. Points perturbed from the unstable steady state will immediately be attracted to oscillations that are much larger in magnitude than those observed in the immediate vicinity of a Hopf bifurcation. Perturbations outside of the basin of attraction of the stable limit cycle will lead to spiraling instability, and eventually to the extinction of one of the metabolites. We examine the biological implications of alternative modes of oscillation in the RM in the Discussion, Section 7.

Vi Generalized Model of the Reverse Metabolator

The appearance of relaxation oscillations inspired us to investigate the role time scales play in the dynamics of the RM. The existence of such divergent time scales is well-known in the regulation of metabolism. The allosteric feedback of small metabolites and the post-translational modification of enzymes target the activity of an enzyme pool at time-scales on the order of microseconds [4]; [1]; [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]; [24]; [15], and affects the actual size of the enzyme pool directly. Remarkably, the rate of an enzymatically catalyzed reaction itself, which the above processes ultimately regulate, is known to take place at characteristic rates spanning seven orders of magnitude, from microseconds to minutes [3].

In order to study the inherent time scales of the RM, we constructed a generalized model. This model was identical in structure to (14), except that the signs of the elements corresponding to the transcriptional regulation of were reversed (). Our intent was to use the entries in the Jacobian to identify time-scales which characterized the dynamics of metabolites and genes in the RM. To understand the connection between and time-scales, recall that each row in describes the linearized dynamics of a small perturbation away from the equilibrium of one state variable. For example, the first row of corresponds to the dynamics of a perturbation in :

(22)

where , and correspond to perturbations in , and , respectively. Upon inspection of (22), all terms on the right hand side share a common divisor which has units of time. We can treat this shared prefactor as a natural measure of the time scale of . Repeating this process for all of the rows in (15), we identified the characteristic time scale for each state space variable. Entries in rows one and two of , corresponding to the dynamics of the two metabolites, shared common divisors and , respectively. Similarly, entries in the final two rows of , corresponding to the genetic elements of the metabolator, shared common divisors and (again with units of time).

Figure 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 , and , 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).

We studied how the stability of depends on the ratio of the metabolic and genetic time scales. To do so, we generated a large number of random instances of and calculated the real component of the leading eigenvalue of (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 30x30 grid of and values. For each point on the grid, we sampled 1000 instances of all other ’s uniformly from the interval for each point in the grid. The results of this simulation are shown in Figure 4.

Two major trends regarding the stabiilty of the RM emerged from our simulations. First, we found that as the dynamics of the genetic elements in the reverse metabolator grew slower, the stability of the circuit as a whole increased. As the value of grew smaller and eventually became faster than the metabolic time scale, the stability of the circuit was observed to fall significantly. Second, our simulations also revealed that small regulatory elasticities contributed to the stability of the circuit, while large elasticities significantly reduced the probability of a stable circuit. This observation qualitatively agrees with our first observation relating high values of to instability. Large magnitude regulatory elasticities increase the same terms in the Jacobian as . Together, our observations suggest that exceptionally slow regulation in the metabolator stabilizes the circuit, a conclusion that parallels commonly held notions on the speed of regulation and the role it plays in stability [18].

Vii Conclusions and Discussion

Contemporary biology is charged with understanding how fundamentally different cellular processes such as signalling, metabolism, and transcriptional regulation, interact and function as a coherent whole. Our work here explored the different types of conclusions which may be drawn from detailed kinetic models and more generic generalized models in the study of a hybrid metabolic-genetic oscillator. We first showed that the core metabolator (1) has a unique fixed point which is stable for all parameter values (recall the end of Section 2). Hence, the core metabolator cannot undergo Hopf bifurcation to oscillations. Then, we developed a GM model of the core metabolator (12), and we showed that this system also cannot undergo Hopf bifurcations by applying the Routh-Hurwitz criterion to the Jacobian (15) (see Section 3.2 and Appendix B).

These analytical results for the core metabolator (1) and its generalized model (12) suggested that the network topology and the number of components in the network play an important role in determining whether or not oscillations are possible. As a result, in Section 4, we explored two more detailed versions of the core metabolator. The first of these has an intermediary delay gene (Figure 1B), but did not exhibit Hopf bifurcations in any of the large number of numerical simulations with randomly chosen system parameters we carried out. The second of these consists of an extra metabolite, which is also present in the synthetically constructed circuit [8], effectively “un-lumping” . 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 induced by the network topologies for the core metabolator, generalized model of the core metabolator, and the first “enhanced” network (Figures 1A and 1B), is responsible for the absence of oscillations in those systems.

Theory and simulation play important roles in synthetic biology. As we have shown, they can highlight generic designs which have certain desirable behaviors, e.g. oscillations, and they can help exclude those that do not. In this respect, the analysis of the enhanced metabolator model with a third metabolite (Figure 1C) suggested that we explore some other four-component models with related designs that do exhibit oscillations. In particular, we studied system (18), the reverse metabolator, in which the wiring from to the two genes is reversed compared to the core metabolator, as shown in Figure 1D. This fundamentally altered the system dynamics so that now activates, rather than represses, the gene and hence is now self-activating. The analysis of Section 5 of the RM reveals that it exhibits both local and global routes to oscillations, via Hopf bifurcations and saddle-node bifurcations of limit cycles, respectively, as shown in Figure 3.

Most of the experimentally observed oscillations in synthetic circuits have arisen through local Hopf bifurcations [6]; [8]; [22]. This is very likely to be due to the comparable ease of predicting local bifurcations. Indeed, the GM methods which are used throughout our work are only capable of identifying global bifurcations in very special and exceedingly rare cases (i.e. when they are known to occur near the intersection of two local bifurcations, like a double Hopf bifurcation). In contrast to the Hopf (where oscillations appear at infinitesimal amplitude), oscillations arising from a fold of limit cycles (as seen in the RM) immediately appear at a finite, non-zero amplitude. This feature may make such a “global” route to oscillation desirable in synthetic biology, which is still tackling the problem of designing robust cellular oscillators. The similarity between the RM and the original metabolator (ultimately, simply just a switch of promoters in the architecture of the circuit) suggests that constructing the RM and investigating its dynamics in vivo may be relatively easy, and may lead to the experimental observation of a global bifurcation.

Finally, we propose that the construction of large generalized models of biological circuits may enable the resolution of several open problems regarding the dynamics of large, complex systems. We investigated one such question here (recall Section 6): the role that time scales play in determining the stability of the system. While we found that stability generally increased as the time scale of regulation became longer (Figure 4), it remains unclear how general this conclusion is among biological networks. A closely related question is whether strongly stabilizing components exist within biological networks. Such components have been argued for in ecological networks [14], as well as in prior work on metabolism using GM [25]. However, the extent to which components of distinct types of biological networks (e.g. signalling, transcriptional regulation, post-translational modifications) act to stabilize or destabilize metabolism is unknown. The investigation of such questions, aided by generalized modeling, seems potentially tractable numerically; computational toolboxes for the automated assembly and inspection of generalized models already exist [11], and their extension to GM with regulation appears straightforward.

Acknowledgements.
ER and DS 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. TJK was supported in part by NSF grant DMS-1109587. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. ER would like to acknowledge the kind support of the Santa Fe Institute while studying the metabolator at the 2011 Summer School.

Appendix A Appendix: Formulation and Non-dimensionalization of the Core Metabolator

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

(23a)
(23b)
(23c)
(23d)

Let’s begin nondimensionalizing . Let . We can rewrite (23a) as

(24)
(25)
(26)

where and .

Then, the nondimensionalized equation for immediately follows as

(27)

where and .

For (23c), we rewrite it by substituting those nondimensionalized variables and parameters which we have already defined and choose to scale so that . Doing so, we get

(28)
(29)

where and .

Finally, the nondimensionalized equation for follows simply and is

(30)

where .

We also note that a simpler way to define the system is in terms of “total” metabolite in the system . Then, . We can non-dimensionalize this system as well to find that

(31)

where .

For the purpose of simplifying the analysis in the article, we let and . Then, the system formed by (26),(27),(31), and (30) is the nondimensionalized core metabolator we study in this article.

Appendix B Appendix: Proof that Generalized Model of the Core Metabolator Does Not Admit Hopf Bifurcations

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

(32)

where we have let and . This means that 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 (32) by and add it to row 2 of (32) to find