Madalena Chaves
email: madalena.chaves@inria.fr

Diego A Oyarzún
email: d.oyarzun@ed.ac.uk

Jean-Luc Gouzé
email: jean-luc.gouze@inria.fr

Analysis of a genetic-metabolic oscillator with piecewise linear models

Keywords: Metabolic pathways, synthetic biology, gene regulatory networks, piecewise affine systems, periodic orbits, metabolic engineering

Interactions between gene regulatory networks and metabolism generate a diversity of dynamics, including multistability and oscillatory behavior. Here, we characterize a regulatory mechanism that drives the emergence of periodic oscillations in metabolic networks subject to genetic feedback regulation by pathway intermediates. We employ a qualitative formalism based on piecewise linear models to systematically analyze the behavior of gene-regulated metabolic pathways. For a pathway with two metabolites and three enzymes, we prove the existence of two co-existing oscillatory behaviors: damped oscillations towards a fixed point or sustained oscillations along a periodic orbit. We show that this mechanism closely resembles the "metabolator", a genetic-metabolic circuit engineered to produce autonomous oscillations in vivo.

Introduction

Metabolic activity is strongly linked with gene expression patterns. Alterations in metabolism are highly controlled by and coordinated with modifications in gene expression patterns [START_REF] Shlomi | A genome-scale computational study of the interplay between transcriptional regulation and metabolism[END_REF][START_REF] Bennett | Metabolic gene regulation in a dynamically changing environment[END_REF][START_REF] Liu | Dynamic metabolic control: towards precision engineering of metabolism[END_REF]. Recent experimental [START_REF] Bennett | Metabolic gene regulation in a dynamically changing environment[END_REF][START_REF] Fung | A synthetic gene-metabolic oscillator[END_REF][START_REF] Kotte | Bacterial adaptation through distributed sensing of metabolic fluxes[END_REF] and theoretical [START_REF] Baldazzi | The carbon assimilation network in Escherichia coli is densely connected and largely sign-determined by directions of metabolic fluxes[END_REF][START_REF] Oyarzún | Multistability and oscillations in genetic control of metabolism[END_REF][START_REF] Oyarzún | Noise propagation in synthetic gene circuits for metabolic control[END_REF] works have studied the synergy and emergent dynamics in interconnected genetic-metabolic networks, showing how systems can take advantage of the control and regulation options offered by such coupling.

Current techniques in synthetic biology allow re-wiring the genetic networks that control expression of metabolic enzymes, so as to produce new dynamic functions [START_REF] Liu | Dynamic metabolic control: towards precision engineering of metabolism[END_REF]. One of the first examples of such engineered systems is known as the metabolator and was developed by Fung et al. in 2005 [4]: it couples metabolic and genetic networks to exhibit autonomous oscillations. Based on the mechanisms that link glycolysis to the respiration cycle, the metabolator circuit makes use of the native pathways of the Escherichia coli bacterium that transform sugars into acetate via acetyl-coenzyme A and acetyl phosphate.

However, the mechanisms directly responsible for sustained oscillatory behavior are not easy to pinpoint and, so far, it is not clear whether the mathematical model originally proposed by Fung et al. [START_REF] Fung | A synthetic gene-metabolic oscillator[END_REF] is indeed capable of periodic oscillations. The study in [START_REF] Reznik | The dynamics of hybrid metabolic-genetic oscillators[END_REF], using bifurcation analysis and a generalized modeling technique, suggests some modifications to the model [START_REF] Fung | A synthetic gene-metabolic oscillator[END_REF] in order to obtain existence of periodic orbits: either adding some delay on enzyme transcription or decoupling one of the metabolites into two different pools. To further explore the mechanisms at play in this system, our analysis proposes an alternative mathematical model for the metabolator circuit, based on a reinterpretation of the enzymes' role and introducing metabolite-controlled gene expression.

Throughout this paper, our goal is to better understand the mechanism that lies at the core of periodic oscillations in genetic-metabolic models, such as the metabolator. Here we study an unbranched metabolic pathway where each enzyme is regulated by one of the metabolites, a schematic view studied in [START_REF] Oyarzún | Multistability and oscillations in genetic control of metabolism[END_REF] (Fig. 1(b)). Our theoretical formalism is based both on timescale separation leading to model reduction, and on approximation of continuous activity functions by piecewise linear functions. Genetic networks are often represented by piecewise linear systems, with the state space partitioned into a family of bounded regular hypercubes. This genetic-metabolic formalism uses a piecewise affine description (such as developed by L. Glass [START_REF] Glass | Stable oscillations in mathematical models of biological control systems[END_REF] or R. Thomas [START_REF] Thomas | Regulatory networks seen as asynchronous automata: A logical description[END_REF][START_REF] Snoussi | Logical identification of all steady states: the concept of feedback loop characteristic states[END_REF] and their collaborators) to provide a rich qualitative model that captures the general scheme of interactions between metabolic and genetic networks, with the advantage of offering an algorithmic view of the dynamics.

Using our formalism, we therefore suggest an alternative updated model for the metabolator and show it allows the co-existence of two types of oscillatory behavior: a periodic orbit and damped oscillations around a Filippov-type fixed point.

A gene-regulated metabolic pathway as a model for the metabolator

The system known as metabolator is a synthetically engineered biological network constructed by [START_REF] Fung | A synthetic gene-metabolic oscillator[END_REF], based on the interactions between a metabolic pathway and its genetic regulation. The pathway transforms sugars into acetate via acetyl-coenzyme A and acetyl phosphate, and synthetic regulation was introduced to have gene expression of both enzymes (P ta and Acs) involved in these transformations under the control of acetyl phosphate. The third enzyme, Ack, catalyzes the outgoing flow and is assumed constitutively expressed, hence constant.

As illustrated in Fig. 1(a),(b), this system can be viewed as a metabolic pathway where enzymatic gene expression is controlled by the metabolites. The analysis of the reduced enzymatic network Fig. 1(a1) (see Appendix A) shows that the main dynamics are governed by a positive loop with two variables (P ta and Acs), both with self-negative regulation, and a constant Ack. However, the pioneering work of René Thomas on the analysis of network topology shows that, in general, a positive loop leads to multistability while a negative loop with at least three variables is required to produce sustained oscillations (see, for instance, [START_REF] Thomas | Dynamical behaviour of biological regulatory networks -i. biological role of feedback loops and practical use of the concept of feedback loop-characteristic states[END_REF]).

In this paper, we therefore propose an alternative mathematical model for the metabolator network by re-interpreting the role of Acs and regulating the production of enzyme Ack by metabolite AcCoA.

With this new model, the 2-variable positive loop is replaced by a 3-variable negative loop through P ta Ack Acs, as in Fig. 1(b1). The enzyme reduced network now exhibits three negative loops, a configuration which most likely exhibits periodic oscillations, as discussed by R. Thomas and colleagues [START_REF] Thomas | Regulatory networks seen as asynchronous automata: A logical description[END_REF][START_REF] Thomas | Dynamical behaviour of biological regulatory networks -i. biological role of feedback loops and practical use of the concept of feedback loop-characteristic states[END_REF][START_REF] Kaufman | A new necessary condition on interaction graphs for multistationarity[END_REF]. A comparison between the mathematical terms describing the models of Fig. 1(a) and (b) is given in Table A.1. Accordingly, there are two significant updates from scheme (a) to scheme (b) of Fig. 1: 1. The influx into AcCoA depends only on S 0 (see Table A.1, the term g 1 (AcP )Acs is replaced by g 1 (s 0 )Acs).

Enzyme

Ack is now assumed to be time-varying and regulated by AcCoA (Ack was assumed constant in the model [START_REF] Fung | A synthetic gene-metabolic oscillator[END_REF]), with an associated consumption flow of the form g 3 (AcP )Ack as shown in Table A.1. Previous work [START_REF] Chaves | Dynamics of complex feedback architectures in metabolic pathways under gene regulation[END_REF] showed that negative regulation of Ack rules out oscillatory behaviour, while positive regulation allows for different types of oscillations. The existence of the positive regulation AcCoA → Ack is hypothetical in E. coli but it is apparently supported in other organisms (as indicated in [START_REF] Wendisch | Regulation of acetate metabolism in corynebacterium glutamicum: transcriptional control of the isocitrate lyase and malate synthase genes[END_REF] on the regulation of acetate metabolism in Corynebacterium glutamicum).

A qualitative model for the metabolator

In Fig. 1(b), the metabolator is represented by an irreversible metabolic pathway with one constant influx (s 0 ) and two metabolites (s 1 and s 2 ), with each transformation step catalyzed by one enzyme. This alternative model can be studied with a new methodology (developed by us in previous work [START_REF] Oyarzún | Multistability and oscillations in genetic control of metabolism[END_REF][START_REF] Chaves | Dynamics of complex feedback architectures in metabolic pathways under gene regulation[END_REF]) to study the interactions between metabolic pathways and the gene network that regulates enzymatic expression. Our formalism considers a self-contained regulatory system, composed by an irreversible metabolic pathway where the transcription of each enzyme is regulated by a metabolite and, conversely, every metabolite regulates expression of at least one enzyme.

To introduce our alternative model, note that the path AcP → LacI P ta can be simplified to AcP P ta, so we will consider the following five variables: x = (s 1 , s 2 , e 1 , e 2 , e 3 ) = (AcCoA, AcP, Acs, P ta, Ack) .

The pathway model is constructed according to two standard properties (see [START_REF] Oyarzún | Multistability and oscillations in genetic control of metabolism[END_REF] for all details): (i) the concentration of metabolites is determined by the difference between incoming and outgoing flows, g i , which are continuous increasing functions [START_REF] Heinrich | The regulation of cellular systems[END_REF]; (ii) the enzymes are produced by metabolite-controlled gene expression modeled by sigmoidal (Hill) functions [START_REF] Yagil | On the relation between effector concentration and the rate of induced enzyme synthesis[END_REF][START_REF] Jong | Modeling and simulation of genetic regulatory systems: A literature review[END_REF] and have linear degradation rates [START_REF] Ropers | Model reduction using piecewise-linear approximations preserves dynamic properties of the carbon starvation response in Escherichia coli[END_REF]:

ṡi = g i (s i-1 )e i -g i+1 (s i )e i+1 , i = 1, 2 (1)
ė1 = κ 0 1 + κ 1 1 h + (s 2 , θ 1 ) -γ 1 e 1 , ė2 = κ 0 2 + κ 1 2 h -(s 2 , θ 2 ) -γ 2 e 2 , (2)
ė3 = κ 0 3 + κ 1 3 h + (s 1 , θ 3 ) -γ 3 e 3 ,
where

g i (X) = k cat,i X/(K M,i + X), h + (X, θ) = X n /(X n + θ n ), and h -(X, θ) = 1 -h + (X, θ) for some integer n > 1.
Each enzyme is degraded at a rate γ i and synthesized at basal and regular rates κ 0 i and κ 1 i , respectively. The θ i represent the activity thresholds. For metabolites, k cat,i is the maximal turn-over rate and K M,i the concentration at half-maximal rate.

A numerical comparative analysis (Appendix A) lends strength to our chosen representation of the metabolator. To illustrate the differences between the two representations in Fig. 1 (a) and (b), in Appendix A we study the system

dx dt = F b (x) + α (F a (x) -F b (x)) (3)
where F a and F b denote the vector fields for the schemes in Fig. 1 (a) and (b), respectively. The case α = 0 corresponds exactly to our system (1)+(2), while the case α = 1 recovers the metabolator model of [START_REF] Fung | A synthetic gene-metabolic oscillator[END_REF], and the intermediate cases 0 < α < 1 represent linear combinations of the two systems. For α > 0.5 the contribution of the Fung et al. model is more important so that the combined system (3) has a higher "resemblance" to F a then to F b . Fig. A.7 shows simulations obtained as α successively increases, taking the values {0, 0.2, 0.4, 0.6, 0.8, 0.9}. A periodic orbit is obtained for the case α = 0, as expected [START_REF] Chaves | Dynamics of complex feedback architectures in metabolic pathways under gene regulation[END_REF]. Then, as α increases up to 0.8 (an 80% resemblance to F a ), a qualitatively similar periodic is still present, with a different shape but within the same region of the state space. For higher values of α, the orbit becomes very much distorted and actually disappears for α = 1. This study not only confirms the dynamical similarity between the two metabolator schemes of Fig. 1, but also suggests that some extra components may be required in the model of Fung et al. in order to guarantee sustained oscillations (in agreement with a previous analysis [START_REF] Reznik | The dynamics of hybrid metabolic-genetic oscillators[END_REF]).

A piecewise affine system of equations

Among the most notable contributions of René Thomas to the mathematical modeling of genetic regulatory networks is the idea that gene expression may be represented using only two levels [START_REF] Thomas | Boolean formalization of genetic control circuits[END_REF]: on or off, 0 or 1. This Boolean representation greatly simplified the analysis of the dynamics of genetic networks and led to a formalism simultaneously developed by L. Glass and co-authors [START_REF] Glass | The logical analysis of continuous, nonlinear biochemical control networks[END_REF][START_REF] Glass | Stable oscillations in mathematical models of biological control systems[END_REF], which has many applications [START_REF] Snoussi | Logical identification of all steady states: the concept of feedback loop characteristic states[END_REF][START_REF] Ropers | Model reduction using piecewise-linear approximations preserves dynamic properties of the carbon starvation response in Escherichia coli[END_REF][START_REF] Abou-Jaoudé | A theoretical exploration of birhythmicity in the p53-Mdm2 network[END_REF][START_REF] Edwards | A modelling framework for gene regulatory networks involving transcription and translation[END_REF]. Based on these ideas, the Hill functions are approximated by step functions σ + (r, θ) = 0 if r < θ and σ + (r, θ) = 1 if r > θ, and σ -(r, θ) = 1 -σ + (r, θ), to facilitate theoretical analysis of (1)- [START_REF] Bennett | Metabolic gene regulation in a dynamically changing environment[END_REF].

A second approximation is based on the fact that interconversion of metabolites is much faster than enzyme synthesis [START_REF] Heinrich | The regulation of cellular systems[END_REF][START_REF] Alon | An introduction to systems biology: design principles of biological circuits[END_REF]. This timescale separation justifies the simplification of metabolite dynamics by setting ṡi ≈ 0, which are the fast variables relative to e i . This quasi-steady state approximation is also at the basis of the Michaelis-Menten kinetic mechanism [START_REF] Segel | On the validity of the steady state assumption of enzyme kinetics[END_REF] and is regularly used in the analysis of systems with both fast and slow components [START_REF] Ropers | Model reduction using piecewise-linear approximations preserves dynamic properties of the carbon starvation response in Escherichia coli[END_REF][START_REF] Edwards | A modelling framework for gene regulatory networks involving transcription and translation[END_REF][START_REF] Kuntz | Model reduction of genetic-metabolic systems using timescale separation[END_REF].

It is useful to introduce the minimal and maximal concentrations for each enzyme:

E 0 i = E off i = κ 0 i γ i , E 1 i = E on i = (κ 0 i + κ 1 i ) γ i .
The enzyme equations imply that the intervals [E 0 i , E 1 i ] are forward-invariant for each enzyme e i , hence we can assume that e i (t) > 0 for all t.

At quasi-steady state, the metabolite equations in (1) yield g 1 (s 0 )e 1 = g i (s i-1 )e i for all i. Using these equalities and the monotonicity of g i , we can write:

σ + (s 2 , θ 1 ) = σ + (g 3 (s 2 ), g 3 (θ 1 )) := σ + e 1 e 3 , 1 β 1 σ -(s 2 , θ 2 ) = σ -(g 3 (s 2 ), g 3 (θ 2 )) := σ - e 1 e 3 , 1 β 2 σ + (s 1 , θ 3 ) = σ + (g 2 (s 1 ), g 2 (θ 3 )) := σ + e 1 e 2 , 1 β
3
where

β 1 = g 1 (s 0 )/g 3 (θ 1 ), β 2 = g 1 (s 0 )/g 3 (θ 2 )
, and

β 3 = g 1 (s 0 )/g 2 (θ 3 ).
These approximations lead to the following three-dimensional, piecewise affine system:

ė1 = κ 0 1 + κ 1 1 σ + (e 1 /e 3 , 1/β 1 ) -γ 1 e 1 ė2 = κ 0 2 + κ 1 2 σ -(e 1 /e 3 , 1/β 2 ) -γ 2 e 2 (4) ė3 = κ 0 3 + κ 1 3 σ + (e 1 /e 2 , 1/β 3 ) -γ 3 e 3 .
From the properties of step functions, system (4) can be viewed as a switching system with state space partitioned into six conic regions separated by the three planes:

e 1 = 1 β 1 e 3 , e 1 = 1 β 2 e 3 , e 1 = 1 β 3 e 2 . (5)
In each region, system (4) is linear and its solutions follow trajectories with φ a = (φ a 1 , φ a 2 , φ a 3 ) as a limit point called the focal point:

ėi = γ i (φ a i -e i ) (6)
where

φ i = E σi i , e 1 > 1 βi e i+1 E 1-σi i , e 1 < 1 βi e i+1 , (7)
σ i ∈ {0, 1}

State transition graph

Piecewise affine systems of the form introduced in [START_REF] Glass | Stable oscillations in mathematical models of biological control systems[END_REF][START_REF] Snoussi | Logical identification of all steady states: the concept of feedback loop characteristic states[END_REF] have the property that their state space is partitioned into a finite number of rectangular regions; in each of these, the vector field is linear and its solution can be explicitly computed. If a trajectory leaves one region to enter another, this implies a switch in the vector field and a corresponding change of solution. Each trajectory can then be constructed as a concatenation of these solutions.

For the piecewise affine system (4) the regions are instead defined in terms of the planes e 1 = e i /β, which leads to a partition of the state space in terms of conic regions. This introduces some differences in the treatment and analysis of steady states and their stability [START_REF] Oyarzún | Multistability and oscillations in genetic control of metabolism[END_REF]. In either case, the global dynamics of a piecewise linear system can be represented by a directed graph whose nodes are the set of regions and the edges correspond to switches between two regions: this is called the state transition graph of the system (see Fig. 3).

The state transition graph for (4) has six regions R a , defined by

R 0,0 = {(e 1 , e 2 , e 3 ) : e 1 < 1 β 3 e 2 , e 1 < 1 β 1 e 3 }, R 1,0 = {(e 1 , e 2 , e 3 ) : e 1 > 1 β 3 e 2 , e 1 < 1 β 1 e 3 }, R 0,1 = {(e 1 , e 2 , e 3 ) : e 1 < 1 β 3 e 2 , 1 β 1 e 3 < e 1 < 1 β 2 e 3 }, (8)
R 1,1 = {(e 1 , e 2 , e 3 ) : e 1 > 1 β 3 e 2 , 1 β 1 e 3 < e 1 < 1 β 2 e 3 }, R 0,2 = {(e 1 , e 2 , e 3 ) : e 1 < 1 β 3 e 2 , e 1 > 1 β 2 e 3 }, R 1,2 = {(e 1 , e 2 , e 3 ) : e 1 > 1 β 3 e 2 , e 1 > 1 β 2 e 3 },
where it is assumed that β 2 < β 1 . Consider the following set of parameters based on data from [START_REF] Fung | A synthetic gene-metabolic oscillator[END_REF] (see Fig. 2)

γ = [0.07, 0.07, 0.07], κ 0 = [0.06, 0.06, 0.06], κ 1 = [0.4, 0.65, 1], β = 1 0.6 , 1 0.7 , 1 0.9 . (9)
For this set of parameters, the transition graph is as given in Fig. 3(a), predicting the existence of two types of oscillatory behavior, along either a "short" (four regions) or a "long" (six regions) cycle. In fact, as illustrated in Fig. 3, there is co-existence of a periodic orbit following the long cycle and a

Filippov fixed point at the intersection of the four regions

R 0,1 ,R 0,2 ,R 1,1 ,R 1,2 .
In addition, this dynamical behavior is robust with respect to variations in the parameters. A recent method analyses robustness of globally stable equilibrium points of genetic networks, by treating the parameters as time-variant linear coefficients [START_REF] Chesi | Robustness analysis of genetic regulatory networks affected by model uncertainty[END_REF]. In our model, to ensure the co-existence of a periodic orbit and a locally stable equilibrium, we first identify the region of parameters that allow a transition graph of the form given in Fig. 3(a). Setting π = (γ, κ 0 , κ 1 , β), define

δ 1,r (π) = sign(φ r 1 - 1 β 1 φ r 3 ), δ 2,r (π) = sign(φ r 1 - 1 β 2 φ r 3 ), δ 3,r (π) = sign(φ r 1 - 1 β 3 φ r 2 ).
The algorithm for constructing the transition graph is given in [START_REF] Chaves | Dynamics of complex feedback architectures in metabolic pathways under gene regulation[END_REF], and depends on the expressions δ i,r (π). It follows that all sets of parameters π which satisfy δ i,r (π) = δ i,r (π) share the same transition graph. The dynamical behavior is thus maintained over a reasonably sized open neighborhood around [START_REF] Reznik | The dynamics of hybrid metabolic-genetic oscillators[END_REF].

In the next sections, a deeper analysis of system (4) establishes the co-existence of the two types of dynamical behavior, under some extra mild conditions on the parameters. In this sense, Theorems 1 and 2 below further characterize robustness of (4) with respect to parameters.

Remark. An interesting parallel can be drawn between genetic-metabolic and signaling networks through their transition graphs: the same dynamics can be generated either with a gene-protein signaling network or with a gene-regulated metabolic pathway. Indeed, the transition graph for our enzyme system (Fig. 3(a)) is similar to the one obtained for a piecewise affine model of the p53-Mdm2 network in [START_REF] Abou-Jaoudé | A theoretical exploration of birhythmicity in the p53-Mdm2 network[END_REF], where it is shown that the system admits two stable periodic orbits, that is the p53-Mdm2 system is birhythmic. The p53-Mdm2 network was initially developed in [START_REF] Abou-Jaoudé | From structure to dynamics: frequency tuning in the p53-Mdm2 network i. logical approach[END_REF], it is composed of 3 variables (p53, nuclear, and cytoplasmic Mdm2) and the interactions consist of a 3D negative feedback loop (p53→Mdm2c→Mdm2n p53) together with a double negative loop between p53 and nuclear Mdm2.

Thus, the transition graph (same in Fig. 3(a) and [START_REF] Abou-Jaoudé | A theoretical exploration of birhythmicity in the p53-Mdm2 network[END_REF]) implies that there are two types of oscillatory behavior, characterized by either large or short amplitudes. In our case, the short amplitude oscillations are damped, while for the p53-Mdm2 system there are two distinct periodic orbits. It may be conjectured that, under further conditions (for instance, a finer partition of the state space as in [START_REF] Abou-Jaoudé | A theoretical exploration of birhythmicity in the p53-Mdm2 network[END_REF]) our enzyme system might also exhibit birhythmicity.

Analysis of the qualitative model by change of variables

To perform a more detailed analysis of system ( 4) and show that it exhibits both a locally steady state and a periodic orbit, we will make a change of variables defined by the boundaries of the conic regions represented in Fig. 3(a) (as proposed in [START_REF] Chaves | Dynamics of complex feedback architectures in metabolic pathways under gene regulation[END_REF]):

w = e 1 - 1 β 1 e 3 , y = e 1 - 1 β 3 e 2 , z = e 1 - 1 β 2 e 3 (10)
and assume, for simplicity, that degradation rates are equal, γ i = γ for all i. We obtain the piecewise linear set of equations

ẇ = κ0 1 + κ 1 1 σ + (w, 0) - κ 1 3 β 1 σ + (y, 0) -γw ẏ = κ0 3 + κ 1 1 σ + (w, 0) - κ 1 2 β 3 σ -(z, 0) -γy (11)
ż = κ0 2 + κ 1 1 σ + (w, 0) - κ 1 3 β 2 σ + (y, 0) -γz.
where κ0 j = κ 0 1 -κ 0 j /β j . It can be observed that this system has a scheme of interactions composed only of negative cycles plus a positive self-loop on w (Fig. 1(b)). According a conjecture due to R. Thomas, this type of circuit can admit periodic oscillations [START_REF] Thomas | Regulatory networks seen as asynchronous automata: A logical description[END_REF][START_REF] Thomas | Dynamical behaviour of biological regulatory networks -i. biological role of feedback loops and practical use of the concept of feedback loop-characteristic states[END_REF][START_REF] Kaufman | A new necessary condition on interaction graphs for multistationarity[END_REF].

The regions represented in Fig. 3 are now defined through the new variables: R 0,0 = {(w, y, z) : w < 0, y < 0, z < 0},

R 1,0 = {(w, y, z) : w < 0, y > 0, z < 0}, R 0,1 = {(w, y, z) : w > 0, y < 0, z < 0}, (12)
R 1,1 = {(w, y, z) : w > 0, y > 0, z < 0}, R 0,2 = {(w, y, z) : w > 0, y < 0, z > 0}, R 1,2 = {(w, y, z) : w > 0, y > 0, z > 0},
and the new focal points can be found by substitution of φ R i from ( 7) into (10):

∆ R w = φ R 1 - 1 β 1 φ R 3 , ∆ R y = φ R 1 - 1 β 3 φ R 2 , ∆ R z = φ R 1 - 1 β 2 φ R 3 ,
where R is one of the six domains defined in [START_REF] Oyarzún | Noise propagation in synthetic gene circuits for metabolic control[END_REF]. The parameters of [START_REF] Thomas | Regulatory networks seen as asynchronous automata: A logical description[END_REF] should satisfy some general conditions regarding the relative positions of the focal points. Note that, in the new variables, the state space is centered at the origin, and trajectories oscillate around the origin. To guarantee the existence of oscillatory behavior (either damped or sustained), trajectories need to switch between the various regions and, hence, each focal point φ R must be located outside its own region R. In addition, focal points need to be at a sufficient distance from each other. Moreover, the form of system ( 6) induces pairwise equalities between the coordinates of the focal points. To summarize, we have for variable w:

∆ 1,0 w < ∆ 1,1 w = ∆ 1,2 w < -1 and 0 < ∆ 0,0 w < 1 < ∆ 0,1 w = ∆ 0,2 w ( 13
)
and for variable y:

∆ 0,0 y = ∆ 1,0 y < ∆ 0,1 y = ∆ 1,1 y < -1 and 1 < ∆ 0,2 y = ∆ 1,2 y . (14)
Define also

E on 1 := β 1 -β 2 β 1 E on 1 .
In addition, we have two technical sufficient conditions to establish the existence of a periodic orbit:

∆ 1,1 w ∆ 1,0 w < ∆ 1,1 y ∆ 1,0 y , (15)
and

∆ 0,2 w > ∆ 1,1 w > E on 1 , ∆ 1,1 y < ∆ 0,2 y = ∆ 1,2 y < ∆ 0,0 y . (16)
Note that all conditions above are satisfied by the set of parameters (9

): ∆ 1,0 w = -8.228, ∆ 1,1 w = ∆ 1,2 w = -2.514, ∆ 0,0 w = 0.343, ∆ 0,1 w = ∆ 0,2 w = 6.057, ∆ 0,2 y = ∆ 1,2 y = 5.800, ∆ 0,1 y = ∆ 1,1 y = -2.557, ∆ 0,0 y = ∆ 1,0 y = -8.
271, E on 1 = 0.939. Inequality (15) holds: 0.3056 < 0.3092.

One of the key points of the analysis is to understand the dynamics of (11) on R 1,1 , the only region where the transition graph splits between two pathways. Trajectories may escape R 1,1 through two boundary planes, corresponding to y = 0 (to follow the short cycle) and w = 0 (to follow the long cycle). So, we need to compare the times at which the solutions y(t) and w(t) hit 0. In each region, the solutions are easily computed from:

w(t) = ∆ 1,1 w + (w(0) -∆ 1,1 w )e -γt ,
(similarly for y(t)) where

∆ 1,1 w = E on 1 - 1 β 1 E on 3 and ∆ 1,1 y = E on 1 - 1 β 3 E on 2 .
The trajectory crosses the plane w = 0 at time T w given by w(T w ) = 0 and, similarly, it crosses the plane y = 0 at time T y given by y(T y ) = 0:

e γTw = ∆ 1,1 w -w(0) ∆ 1,1
w and e γTy = ∆ 1,1 y -y(0) ∆ 1,1 y . Therefore, if e γTw < e γTy , the trajectory follows the long cycle, otherwise it follows the short cycle.

Region R 1,1 can thus be partitioned into two parts and a separatrix:

R 1,1 long = (w, y, z) ∈ R 1,1 : w ∆ 1,1 w > y ∆ 1,1 y , R 1,1 short = (w, y, z) ∈ R 1,1 : w ∆ 1,1 w < y ∆ 1,1 y , (17)
R 1,1 sep = (w, y, z) ∈ R 1,1 : w ∆ 1,1 w = y ∆ 1,1 y .
A projection of these regions on the plane (e 1 , e 2 ) is shown in Fig. 4. In particular, under the above conditions on the parameters (satisfied by ( 9)), each of the regions R 1,1 long (blue + symbols) and R 1,1 short (blue dots) is invariant under the flow of the system, as shown in the next two sections.

Damped oscillations along the short cycle

To establish the existence of damped oscillations along the short cycle, the main observation is that trajectories starting in R 1,1 short will remain within the four regions that form the short cycle, R 0,1 ,R 0,2 ,R 1,1 , and R 1,2 . Within these four regions, it follows that system (11) can be reduced to a basic negative feedback loop involving variables y and z (see (18) below), whose behavior is well known [START_REF] Glass | Stable oscillations in mathematical models of biological control systems[END_REF]: it converges towards a fixed point, as illustrated in Fig 5. This result is summarized in Theorem 1. It implies that a steady state is thus always present in the qualitative metabolator corresponding to initial conditions near the axis Acs = Ack β2 = P ta β3 .

Theorem 1. Consider system [START_REF] Thomas | Regulatory networks seen as asynchronous automata: A logical description[END_REF] under assumptions ( 13)-( 14) and initial conditions of the form (w(0), y(0), z(0)) such that w(0)/∆ 1,1 w < y(0)/∆ 1,1 y . Then, the solution converges to the point ∆ F = ((1 -β2 β1 )E on 1 , 0, 0). The corresponding solution of the e i system converges to φ F = E on 1 (1, β 3 , β 2 ).

Observe that, inside the four regions R 0,1 ,R 0,2 ,R 1,1 , and R 1,2 , it always holds that e 3 < β 1 e 1 or w > 0 implying σ + (w, 0) ≡ 1 and hence variables y, z can be analyzed as an autonomous system, which is in fact a negative loop z → y z:

ẏ = κ0 2 + κ1 1 + κ 1 2 β 3 σ + (z, 0) -γy ż = κ0 3 + κ1 1 + κ 1 3 β 2 σ -(y, 0) -γz (18)
with κ0 j = κ 0 1 -κ 0 j /β j and κ1 j = κ 1 1 -κ 1 j /β j . This is a well known planar piecewise linear system: in [START_REF] Glass | Stable oscillations in mathematical models of biological control systems[END_REF],

the authors show that it has damped oscillations in the case of aligned focal points, i.e. placed at the vertices of a rectangle, as in [START_REF] Yagil | On the relation between effector concentration and the rate of induced enzyme synthesis[END_REF]. The first return map of (18) on the half segment {(y, z) : z = 0, y > 0} is of the form:

P (y) = ρy 1 + ry , r > 0, ρ = ∆ 0,2 z ∆ 0,1 y ∆ 1,1 z ∆ 1,2 y ∆ 0,2 y ∆ 0,1 z ∆ 1,1 y ∆ 1,2 z . ( 19
)
In our case, ρ = 1 due to the location of the focal points, which satisfy ∆ 0,2 y = ∆ 1,2 y , ∆ 0,1 y = ∆ 1,1 y and ∆ 0,2 z = ∆ 0,1 z , ∆ 1,1 z = ∆ 1,2 z . This implies P (y) < y for y > 0.

To apply this result, it is first necessary to guarantee that the trajectories cannot move away from the four regions R 0,1 ,R 0,2 ,R 1,1 , and R 1,2 whenever the initial condition is in R 1,1 short , so that ( 18) holds for all times; the proof is done in Appendix B.

A periodic orbit along the long cycle

To prove the existence of a periodic orbit, we will study the dynamics of variables w and y along the long cycle,

ẇ = κ0 1 + κ 1 1 σ + (w, 0) - κ 1 3 β 1 σ + (y, 0) -γw = g w (t, w, y) ẏ = κ0 3 + κ 1 1 σ + (w, 0) - κ 1 2 β 3 σ -(z, 0) -γy = g y (t, w, y). (20)
However, in contrast to [START_REF] Yagil | On the relation between effector concentration and the rate of induced enzyme synthesis[END_REF], the reduced system ( 20) is not autonomous, since the transitions R 0,1 → R 0,2 and R 1,2 → R 1,1 do not take place at fixed values of the w coordinate but instead depend on time. Nevertheless, variable z can be written in terms of w(t) in the form z(t) = w(t) + ( 1 β1 -1 β2 )e 3 (t), as follows from the definition of w and z (10), so those two types of transitions happen at instants t such that z(t) ≡ 0 or, equivalently, w(t) = β1-β2 β1 e 1 (t).

To study system [START_REF] Ropers | Model reduction using piecewise-linear approximations preserves dynamic properties of the carbon starvation response in Escherichia coli[END_REF] we will use an idea recently developed in [START_REF] Poignard | Periodic oscillations for non monotonic smooth negative feedback circuits[END_REF] and [START_REF] Chaves | 2D piecewise affine models approximate real continuous dynamics up to invariant sets[END_REF], which consists of bounding the original system by interior and/or exterior systems. These bounding systems are defined so that their trajectories circumscribe those of the original system and determine a region of the state space where trajectories of the original system evolve. For example, in [START_REF] Poignard | Periodic oscillations for non monotonic smooth negative feedback circuits[END_REF], a system with an uncertainty window is bounded by interior and exterior systems, both having periodic orbits which define a torus where trajectories of the original system evolve; in particular, this implies the existence of periodic behavior for the original system with uncertainties. The bounding systems are typically easier to study than the original system. Accordingly, we will construct a new simplified two-dimensional system that provides an inner bound for system [START_REF] Ropers | Model reduction using piecewise-linear approximations preserves dynamic properties of the carbon starvation response in Escherichia coli[END_REF] in the sense of Def. 1: Definition 1. Consider two parametrized bi-dimensional curves, φ = {(w(t), y(t)) : t ∈ [0, T ]} and φ = {( w(t), ỹ(t)) : t ∈ [0, T ]}, defined in the same orthant of the plane R 2 . The curve φ is said to be an inner bound on R to curve φ if:

∀t ∈ [0, T ], t ∈ [0, T ] : w(t) = w( t) ⇒ |y(t)| ≤ ỹ( t) .
The inner bound curve φ is thus closer to the origin at all points (see each orthant of Fig. 6). Construct a new system ( ẇ, ẏ) = h(w, y) as follows:

ẇ = h Q w (y) = γ ∆ Q w -w ẏ = h Q y (z) = γ ∆ Q y -y , Q ∈ {A, B, C, D} (21)
where the four regions Q are defined in terms of the R i,j domains defined in [START_REF] Snoussi | Logical identification of all steady states: the concept of feedback loop characteristic states[END_REF]:

Q A = R 0,0 , Q B = R 0,1 ∪ R 0,2 , Q C = {R 1,2 | w > E on 1 }, Q D = R 1,1 long ∪ R 1,0 ∪ {R 1,2 | w < E on 1 },
where R 1,2 |p denotes the restriction of R 1,2 to the property p, and the corresponding focal points are chosen as

∆ A = (∆ 0,0 w , ∆ 0,0 y ) = (φ 0,0 1 -φ 0,0 3 /β 1 , φ 0,0 1 -φ 0,0 2 /β 3 ) , ∆ B = (∆ 0,2 w , ∆ 0,2 y ) = (φ 0,2 1 -φ 0,2 3 /β 1 , φ 0,2 1 -φ 0,2 2 /β 3 ) , ∆ C = (∆ 1,2 w , ∆ 1,2 y ) = (φ 1,2 1 -φ 1,2 3 /β 1 , φ 1,2 1 -φ 1,2 2 /β 3 ) , ∆ D = (∆ 1,1 w , ∆ 1,1 y ) = (φ 1,1 1 -φ 1,1 3 /β 1 , φ 1,1 1 -φ 1,1 2 /β 3 ) . Define X := R 0,0 ∪ R 0,1 ∪ R 0,2 ∪ R 1,0 ∪ R 1,1 long ∪ R 1,2
and note that also

X = Q A ∪ Q B ∪ Q C ∪ Q D .
Proposition 5.1. Assume that conditions on the parameters ( 13)-( 14) hold. Consider the solutions of systems ( 20) and ( 21) in each region R ∈ X , respectively denoted by φR = ( w(t), ỹ(t)) (t ∈ [0, T ]) and φ R = (w(t), y(t)) (t ∈ [0, T ]), for the same initial condition (w 0 , y 0 ) ∈ R. If (15) holds, then curve φ R is an inner bound on R to curve φR , for each R.

In other words, the trajectory of system (21) always evolves closer to the origin and, starting from the same initial condition, the two trajectories cannot cross, with the curves corresponding to trajectories of ( 21) providing an inner bound to those of ( 18), as depicted in Fig. 6. The two-dimensional system ( 21) is easier to study than [START_REF] Ropers | Model reduction using piecewise-linear approximations preserves dynamic properties of the carbon starvation response in Escherichia coli[END_REF], and it has a periodic orbit:

Proposition 5.2. Under conditions on the parameters ( 13)-( 16), and initial conditions in X , system (21) has a unique, attractive, periodic orbit on X . The proofs of the two Propositions are given in Appendix C and Appendix D, respectively. The existence of periodic oscillations for the metabolator system now follows immediately as a corollary of Props. 5.1 and 5.2. The main idea is that, since ( 21) is an inner bound to [START_REF] Ropers | Model reduction using piecewise-linear approximations preserves dynamic properties of the carbon starvation response in Escherichia coli[END_REF], trajectories of (20) must thus remain outside the periodic orbit of system (21).

Theorem 2. Under conditions on the parameters ( 13)-( 16), and initial conditions in X , system [START_REF] Yagil | On the relation between effector concentration and the rate of induced enzyme synthesis[END_REF] has a periodic orbit on X .

The proof will follow from four observations:

(1) one can assume all trajectories of (20) remain within the rectangle

Q = [E off 1 -E on 3 /β 1 , E on 1 - E off 3 /β 1 ] × [E off 1 -E on 2 /β 3 , E on 1 -E off 2 /β 3 ]
, a natural forward-invariant region given by the vector field;

(2) by Prop. 5.1, system ( 21) is an inner bound for (20);

(3) by Prop. 5.2, system (21) admits a periodic orbit. Let C be the area enclosed by this orbit and y C be the intersection of the orbit with the semi-axis y > 0 (4) for each region R a , the focal point is in the "next" neighbor region, in a sequence defined by the long cycle in the state transition diagram in Fig. 3(a).

Taken together, (1)-( 3) imply that trajectories of the (non-autonomous) system (20) must evolve inside the region Q \ C, for all times. Point (4) implies that every trajectory eventually leaves its starting region and enters the "next" neighbor region, also following the sequence defined by the long cycle. Consider the interval I y = [0, y C ]: from points (1)-( 4), one can construct a first return map, F : I y → I y for system [START_REF] Ropers | Model reduction using piecewise-linear approximations preserves dynamic properties of the carbon starvation response in Escherichia coli[END_REF]. Then F is continuous and I y is invariant under F . By Brouwer's fixed point theorem, there exists ȳ such that F (ȳ) = ȳ, which in turn implies the existence of a periodic orbit for system [START_REF] Ropers | Model reduction using piecewise-linear approximations preserves dynamic properties of the carbon starvation response in Escherichia coli[END_REF].

These results conclude the theoretical analysis of the qualitative metabolator system to establish two co-existing types of asymptotic behaviors. Theorems 1 and 2 characterize the behavior of the reduced enzyme model [START_REF] Fung | A synthetic gene-metabolic oscillator[END_REF], which is itself an approximation of the full genetic-metabolic oscillator (1)- [START_REF] Bennett | Metabolic gene regulation in a dynamically changing environment[END_REF]. Typically, these approximations give a good qualitative description of the dynamial behavior of the original system, so we can expect co-existence of the damped and sustained oscillatory behaviors for system (1)-(2) (at least for high Hill coefficients). This is indeed verified numerically in Fig A .7 (α = 0) as well as in [START_REF] Chaves | Dynamics of complex feedback architectures in metabolic pathways under gene regulation[END_REF]. Another link between a continuous system and its qualitative piecewise linear description is provided, for instance, by [START_REF] Poignard | Periodic oscillations for non monotonic smooth negative feedback circuits[END_REF] and [START_REF] Chaves | 2D piecewise affine models approximate real continuous dynamics up to invariant sets[END_REF], which motivated Def. 1. As summarized above, under appropriate conditions, we can find bounds for a continuous system in terms of two piecewise linear systems; these bounds remain close to the original system and provide a faithful description of its dynamical behavior.

Discussion and conclusions

The goal of the paper was to study the interplay between genetic and metabolic networks. We have focused on the self-contained case of an irreversible metabolic pathway whose steps are catalyzed by a family of distinct enzymes. Transcription of the genes encoding for the enzymes is itself regulated by the metabolites and each metabolite is required to regulate at least one enzyme.

This class of models demonstrates the variety of dynamical behavior induced by genetic regulation in metabolic pathways. In particular, the same interconnecting architecture may lead to (multi) stability, damped or sustained oscillatory behavior and, in particular, co-existence of different regimes, such as convergence towards a steady state or a periodic orbit.

To deal with and analyze these models, we have developed a theoretical formalism based on two approximations: first, model reduction due to the different timescales of genetic and metabolic processes and, second, continuous activity functions are replaced by similar piecewise linear functions. Genetic networks are often represented by piecewise linear systems, with the state space partitioned into a family of bounded regular rectangles. It is interesting to note that the effect of metabolites on the transcriptional network leads to a similar system of piecewise linear equations but now the state space is partitioned into a family of cones, starting at the origin. This formalism captures the core processes at play between metabolic pathways and the genetic network associated to enzymatic regulation, such as in the lac operon [START_REF] Oyarzún | Multistability and oscillations in genetic control of metabolism[END_REF] and other examples. It is also useful in the study and design of new synthetic systems [START_REF] Oyarzún | Design of a bistable switch to control cellular uptake[END_REF]. As a special example, we have analyzed a qualitative model for the synthetic metabolator circuit, which predicts two types of oscillatory behavior, and is in agreement with the experimental results reported in [START_REF] Fung | A synthetic gene-metabolic oscillator[END_REF]. Our analysis shows that, in current models, the enzymatic reduced system is basically governed by a positive feedback circuit with two enzymes and two negative self-loops (Fig. 1(a1)), which may not lead to sustained oscillatory behavior. We propose an alternative model, by adding a single interaction (Fig. 1(b)). The new reduced system has a negative feedback loop involving three variables, much more likely to exhibit periodic oscillations [START_REF] Thomas | Regulatory networks seen as asynchronous automata: A logical description[END_REF][START_REF] Kaufman | A new necessary condition on interaction graphs for multistationarity[END_REF], as confirmed by our results.

The advantages of this formalism include a qualitative but still refined description of the system interactions, as well as a strong theoretical and algorithmic methodology that gives a global view of the system's trajectories. The piecewise linear formalism facilitates the analytic study of the system, in particular to prove the existence of a periodic orbit for the qualitative metabolator example, where we used comparison to an inner bounding system.

Table A.1: Simplified differential equations for the original metabolator system (middle column) and using the formalism for gene regulated metabolic pathways (right column). The functions ĝi , g i are strictly increasing. The main differing terms, one in each of the metabolite equations, are highlighted. The original metabolator system also contains degradation terms not included in the alternative formalism, but these introduce only minor effects.

Variable

Metabolator, Fung et al., Gene-regulated metabolic pathway, vector field F a vector field F b

s 1 = AcCoA S 0 -γ A s 1 + ĝ1 (s 2 )e 1 -g 2 (s 1 )e 2 g 1 (s 0 )e 1 -g 2 (s 1 )e 2 s 2 = AcP g 2 (s 1 )e 2 -ĝ1 (s 2 )e 1 -k e3 (s 2 -H ext ) g 2 (s 1 )e 2 -g 3 (s 2 )e 3 e 1 = Acs V 10 + V 1 h + (s 2 , θ 1 ) -γ 1 e 1 V 10 + V 1 h + (s 2 , θ 1 ) -γ 1 e 1 e 2 = P ta V 20 + V 2 h -(s 2 , θ 2 ) -γ 2 e 2 V 20 + V 2 h -(s 2 , θ 2 ) -γ 2 e 2 e 3 = Ack - V 30 + V 3 h + (s 1 , θ 3 ) -γ 3 e 3
network for the model in [START_REF] Fung | A synthetic gene-metabolic oscillator[END_REF], we assume that the metabolites are at quasi-steady state and that the strictly increasing functions ĝ1 , g 2 are linear, of the form ĝ1 (s 2 ) = k 1 s 2 and g 2 (s 1 ) = k 2 s 1 , to obtain:

S 0 -γ A s 1 + k 1 s 2 • e 1 -k 2 s 1 • e 2 = 0 k 2 s 1 • e 2 -k 1 s 2 • e 1 -k e3 • s 2 = 0.
Solving these equations with respect to s 1 and s 2 leads to

s 1 = S 0 γ A k 1 e 1 + k e3 k 1 e 1 + k e3 + k e3 k2 γ A e 2 and s 2 = S 0 k2 γ A e 2 k 1 e 1 + k e3 + k e3 k2 γ A e 2 .
Now, letting k e3 ≈ k 3 e 3 denote e 3 activity and substituting these dependencies into the functions h + (s 2 ) or h -(s 2 ), yields a 3-dimensional system with an interaction graph as shown in Fig. 1(a1). The differences between the two enzyme networks are discussed in Section 2.

The second step in model comparison aims at showing that the dynamics generated by systems F a and F b are qualitatively similar, by numerically studying system (3):

dx dt = F b (x) + α(F a (x) -F b (x))
with 0 < α < 1. As α increases from 0 to 1, system (3) is transformed from our gene-regulated metabolic pathway into the metabolator of Fung et al. For α > 0.5, system (3) has a closer resemblance to F a .

From Table A.1, the difference between the two models is given by:

F a (x) -F b (x) =       S 0 -γ A s 1 + ĝ1 (s 2 )e 1 -g 1 (s 0 )e 1 -ĝ 1 (s 2 )e 1 + g 3 (s 2 )e 3 -k e3 s 2 0 0 -(V 30 + V 3 h + (s 1 , θ 3 ) -γ 3 e 3 )       (A.1)
The simulations in Fig. A.7(a) show that the system always exhibits a periodic orbit as α increases from 0 to 0.8. Although the orbit remains located within the same region of the state space, its shape is increasingly distorted and finally the periodic orbit disappears when the combined system resembles F a by more than 80%. s 2 . This means that the difference function would be:

F a (x) -F b (x) =       S 0 -γ A s 1 + ĝ1 (s 2 )e 1 -g 1 (s 0 )e 1 -ĝ 1 (s 2 )e 1 0 0 0       (A.2)
The solutions corresponding to this scenario are depicted in Fig. A.7(b). Once again, the combined system exhibits a periodic orbit for α values in the interval [0, 0.8]. However, in this case both the location and shape of the orbit remain qualitatively similar, and the orbit disappears only for α > 0.9.

In both scenarios, a system whose vector field is composed of the original metabolator by up to 80% and as low as 20% of the alternative metabolator still exhibits a periodic orbit. Altogether, these results support our formalism as a reasonable qualitative model to capture the scheme of interactions between the metabolic and genetic networks within the metabolator. The fact that oscillations disappear for α > 0.9 indicates that the original metabolator model may require extra terms (such as a dynamically evolving enzyme e 3 ), in order to guarantee sustained oscillations.

For the numerical simulations in

Appendix B. Proof of Theorem 1

To prove the Theorem, we need to show invariance of R 1,1 short for the flow of system [START_REF] Thomas | Regulatory networks seen as asynchronous automata: A logical description[END_REF]. In other words, we need to show that, starting from an initial condition (w(0), y(0), z(0)) in R 1,1 short , the trajectory returns to R 1,1 short after a full turn.

To do this, assume that (w(0), y(0), z(0)) ∈ R 1,1 short is also at the boundary between R 1,2 and R 1,1 , that is: w(0) > 0, y(0) > 0, z(0) = 0, w(0) ∆ where the last inequality follows from β 2 < β 1 and A > 0. This means w(T ) > w(0) > 0. Analysis of [START_REF] Yagil | On the relation between effector concentration and the rate of induced enzyme synthesis[END_REF] shows that, at t = T the map [START_REF] Jong | Modeling and simulation of genetic regulatory systems: A literature review[END_REF] applies and yields: y(T ) = P (y(0)) < y(0).

Finally, since both ∆ 1,1 w < 0 and ∆ 1,1 y < 0: implying that the trajectory still satisfies condition (B.1). In conclusion, the sub-region R 1,1 short is invariant under the flow of system [START_REF] Thomas | Regulatory networks seen as asynchronous automata: A logical description[END_REF]. Therefore, trajectories do no leave those four regions and e 1 (t) → E on 1 for all t > 0. The variables y, z evolve according to [START_REF] Yagil | On the relation between effector concentration and the rate of induced enzyme synthesis[END_REF]. By the results in [START_REF] Glass | Stable oscillations in mathematical models of biological control systems[END_REF], we conclude that trajectories converge towards the point at the intersection of the four regions: (0, 0). This implies: e 1 (t) -1 β2 e 3 (t) → 0 and so e 3 (t) → β 2 E on 1 and w(t) → (1 -β2 β1 )E on 1 .

In terms of the variables e i the (locally) stable steady state φ F is thus given by φ F 1 = E on 1 and:

∆ F = (0, 0) ⇔ y = 0 and z = 0 implying φ F 2 = β 3 E on 1 and φ F 3 = β 2 E on 1 .

Figure 1 :

 1 Figure 1: Two models for the metabolator network, with solid bold arrows representing the metabolic network and solid light arrows representing genetic regulation. There are two metabolites, acetyl-coenzyme A (AcCoA) and the acetate pool represented by acetyl phosphate (AcP ), and three enzymes, acetyl-CoA synthetase (Acs), phosphate acetyltransferase (P ta), and acetate kinase (Ack). (a) The metabolator network analyzed by Fung et al., and (b) a qualitative view as an irreversible pathway with two metabolites and three enzymes, where all production sources for each metabolite are merged into one arrow and all the enzymes are regulated by a metabolite in the pathway; (a1), (b1) show a corresponding reduced network involving only the enzymes, for the schemes in (a), (b), respectively.

  denotes the sign of the corresponding step function, and the indexes ( 1 , 2 , 3 ) = (3, 3, 2) follow from (4).Fig.2clearly shows that model (4) has the capacity to reproduce the oscillatory behavior observed in the synthetic metabolator, since it simultaneously captures the shape of the oscillations, the relative minimal and maximal amplitudes, as well as the period (around 46 minutes).

Figure 2 :

 2 Figure 2: Model (4) fits synthetic metabolator data (blue dots), taken from Figure 1(a) in [4], with GFP arbitrary units normalized to 5. Parameters are γ i = 0.07, κ 0 i = 0.06 (for all i), κ 1 = [0.4, 0.65, 1], and β = [1/0.6, 1/0.7, 1/0.9]. The period of the model solution is approximately 45.88 minutes.

Figure 3 :

 3 Figure 3: (a) Transition graph of system (4) with parameters (9); (b) Dynamics of the model (4) among the regions R a . Solutions are shown (red curves) for two sets of initial conditions (cyan dots). Three planes are pictured: e 3 = β 2 e 1 , e 3 = β 1 e 1 , and e 2 = β 3 e 1 (vertical). One trajectory follows the long transition cycle in (a) and converges to a periodic orbit. The other follows the short cycle, spiraling towards an equilibrium at the intersection of the four regions (blue dot). Parameters are as in (9).

Figure 4 :

 4 Figure 4: The basins of attraction for the reduced enzyme network (4) with parameters (9), projected on the plane (e 1 , e 2 ). Initial conditions are taken on the plane e 3 = β 2 e 1 with e 2 < β 3 e 1 . The solid black line represents the plane e 2 = β 3 e 1 . The dashed black line represents the separatrix region R 1,1 sep , between the two basins of attraction. The red curve is the projection of the periodic orbit on the plane (e 1 , e 2 ). The blue + symbols represent initial conditions leading to the periodic orbit along the long cycle. The blue dots represent initial conditions leading to convergence along the short cycle towards the fixed point φ F = E on 1 (1, β 3 , β 2 ) (red star), at the intersection of the planes e 3 = β 2 e 1 and e 2 = β 3 e 1 .

Figure 5 :

 5 Figure 5: Convergence of a trajectory towards the fixed point φ F = E on 1 (1, β 3 , β 2 ), along the short cycle, projected on the plane (z, y).

Figure 6 : 1 E

 61 Figure 6: Convergence of a trajectory towards the periodic orbit, projected on the plane (w, y). The focal points for the six regions are labeled ∆ R . The blue curve represents a trajectory of the (original) system (18) and the red curve represents a trajectory of the inner bound system (21). The dashed vertical line represents the value w = β 1 -β 2 β 1 E on 1 and the solid black line on the orthant w > 0, y > 0 represents the separatrix R (1,1) sep .

Figure A. 7 :

 7 Figure A.7: The periodic solutions of linear combination system (3), as parameter α increases from 0 (the metabolic pathway formalism) to 1 (the simplified metabolator). Case (a) illustrates this transformation, using (A.1), while case (b) illustrates the role of enzyme e 3 = Ack through (A.2). In (a), the orbit for α = 0.8 is truncated for better comparison, as it extends to s 1 ≈ 25.

  Fig. A.7, the functions ĝi , g i are of the form:g i (s) = k cat,i s s + k M,iwith the following parameters:kcat,1 = k cat,1 = 60, k cat,2 = 180, k cat,3 = 480, kM,1 = 0.002, k M,i = 1, i = 1, 2, 3, s 0 = 1.5, S 0 = 20s 0 , γ A = 20γ 1 , k e3 = 20, H ext = 0.The parameters for the enzyme equations are:γ i = 0.0116, i = 1, 2, 3, [θ 1 , θ 2 , θ 3 ] = [0.047, 0.055, 0.22], n = 10, [V 10 , V 20 , V 30 ] = 0.0116 × [1, 4, 1.5], [V 1 , V 2 , V 3 ] = 0.0116 ×[14, 16, 28.5].

1 )Furthermore, note that e 1

 11 Then the corresponding trajectory ψ(t; w(0), y(0), z(0)) will perform at least one turn along the short cycle. Let T be the first return to the boundary between R 1,2 and R 1,1 : z(T ) = 0 and w(T ) > 0, y(T ) > 0. Thenz(0) = z(T ) = 0 ⇒ e 1 (T )is an increasing function while w > 0, since ė1 = γ(E on 1 -e 1 ), so e 1 (T ) = e 1 (0) + A for some A > 0, which implies e 3 (T ) = e 3 (0) + β 2 A. Thereforee 1 (T ) -1 β 1 e 3 (T ) = e 1 (0) + A -1 β 1 (e 3 (0) + β 2 A)= e 1 (0) -1 β 1 e 3 (0) + A(1 -

Acknowledgments M.C. and J.-L.G. were partly supported by the French agency for research through projects ICycle ANR-16-CE33-0016-01 and Maximic ANR-17-CE40-0024-01. D.O. was partly supported by the Human Frontier Science Program through a Young Investigator Grant (RGY0076-2015).

Appendix A. Comparison of two Metabolator models

To better understand how closely the two models are related, this section compares our suggested metabolator model with the one studied in [START_REF] Fung | A synthetic gene-metabolic oscillator[END_REF]. The equations for the two models are summarized in Table A. [START_REF] Shlomi | A genome-scale computational study of the interplay between transcriptional regulation and metabolism[END_REF].

The first step is to compare the topology of the two reduced enzyme networks (shown in Fig. 1(a1) and (b1)). The network for our model is directly obtained from (4). To deduce the corresponding Appendix C. Proof of Proposition 5.1

The proof follows by checking that, in each region R a , the slope of the trajectory of the inner system is always steeper than that of the original system. This is clear from Fig. 6, and means that trajectories of system [START_REF] Ropers | Model reduction using piecewise-linear approximations preserves dynamic properties of the carbon starvation response in Escherichia coli[END_REF] always enclose those of system [START_REF] Thomas | Boolean formalization of genetic control circuits[END_REF].

For either system [START_REF] Ropers | Model reduction using piecewise-linear approximations preserves dynamic properties of the carbon starvation response in Escherichia coli[END_REF] or [START_REF] Thomas | Boolean formalization of genetic control circuits[END_REF], observe that the w and y have affine and decoupled equations in each region so can easily be solved explicitly. To show that ( 21) is an inner bound for system [START_REF] Ropers | Model reduction using piecewise-linear approximations preserves dynamic properties of the carbon starvation response in Escherichia coli[END_REF] in the sense of Def. 1, we will check each region separately:

(a) In regions R 0,0 , R 0,2 , R 1,1 long , and

for all t ∈ min{T, T }. Definition 1 is satisfied since w(t) = w(t) and |y(t

for all t ∈ min{T, T }. Definition 1 is satisfied since w(t) = w(t) and |y(t

for all t ∈ min{T, T }. To check Definition 1, find t < t such that w(t) = w( t):

And now compute y(t) -ỹ( t) and replace 1 -e -γ t by (C.1):

From ∆ 1,0 y < 0 and condition (15) the first term is negative. From t > t, y 0 > 0, ∆ 1,0 y < 0 and w 0 /∆ 1,0 w > 0 the second term is also negative. In conclusion, in all regions w(t) = w( t) implies |y(t)| < ỹ( t) , which establishes the necessary property in Def. 1.

Appendix D. Proof of Proposition 5.2

To establish the existence of a periodic orbit for system [START_REF] Thomas | Boolean formalization of genetic control circuits[END_REF], we will proceed by first constructing a Poincaré map (P ) which represents the intersection of the orbits with a lower-dimensional section transversal to the flow. For any initial condition on this section (say y 0 ), the Poincaré map gives the first return to the section after a full orbit (P (y 0 )). The system admits a periodic orbit if the Poincaré map has a fixed point.

Construction of the Poincaré map. It is not difficult to check that the half-line

at the boundary between Q C and Q D (see dashed vertical line in Fig. 6), is transversal to the flow of system ( 21), so we will define P : S → S.

Starting at S, trajectories of (21) will successively cross the regions Q D ,Q A ,Q B , and

with y 0 , y 4 > 0, w 1 < 0, y 2 < 0, and w 3 > 0. The Poincaré map P can be constructed as the composition of the impact maps p i corresponding to the crossing of each region Q, from a starting point in one boundary to the next boundary:

w -w 3

.

Using ( 13)-( 14) and ( 16), |∆ 1,1 w | > E on 1 , for maps p i , i = 1, 2, 3, 4, it is easy to check the signs of first and second derivatives:

By composition, the derivatives of map P satisfy: P > 0 and P > 0, so that the Poincaré map is a, negative, strictly increasing and strictly convex function.

Invariance of the interval [y * 0 , +∞) under map P . First, we need to check that trajectories remain within the boxes that constitute the long cycle for appropriate initial conditions. From the proof of Prop. 5.1, it follows that we need (w 0 , y 0 , z 0 ) ∈ R 1,1 long and w 0 > E on 1 :

We will consider the limiting point

y and show that it satisfies P (y * 0 ) > y * 0 . Then, because P is increasing and y 0 is positive, y 0 > y * 0 implies y * 0 < P (y * 0 ) < P (y 0 ) and, together with ∆ 1,1 y < 0, this guarantees that trajectories do not leave the long cycle. To compute P (y * 0 ) proceed by substitution:

where

The signs and magnitudes of the focal points [START_REF] Thomas | Dynamical behaviour of biological regulatory networks -i. biological role of feedback loops and practical use of the concept of feedback loop-characteristic states[END_REF] and assumptions ( 16) imply

We want to check that P (y * 0 ) -y * 0 = y * 4 -y * 0 > 0. Substitution yields:

where

The denominator of y * 4 -y * 0 is positive and the term α 1 is positive, from the signs of q 3 , r 3 and the focal points. The term α 2 is also positive since E on 1 / ∆ 1,1 w < 1 and ∆ 1,2 y / ∆ 1,1 y > 1 using assumptions [START_REF] Wendisch | Regulation of acetate metabolism in corynebacterium glutamicum: transcriptional control of the isocitrate lyase and malate synthase genes[END_REF].

The desired result follows.

Fixed point of map P . The system (21) admits a periodic orbit if the Poincaré map P (y) has a fixed point: P (y p ) = y p > 0. To show this, recall that P is a continuous, increasing, and strictly convex map, defined on the interval [y * 0 , +∞). We know that P (y * 0 ) > y * 0 . In the positive limit, we have:

w which is a finite quantity and implies P (y ∞ ) < y ∞ for some sufficiently large y ∞ . By convexity of P , there is exactly one fixed point in the interval (y * 0 , y ∞ ), and hence there is a unique periodic orbit for the system [START_REF] Thomas | Boolean formalization of genetic control circuits[END_REF].

In addition, since P is a convex function, P (x) ≤ (P (v) -P (x))/(v -x). Letting v = y p and x = y *