Alesia Nazim Fatès
email: nazim.fates@inria.fr

A pedagogical example: a family of stochastic cellular automata that plays

Alesia is a two-player zero-sum game which is quite similar to the rock-paper-scissors game: the two players simultaneously move and do not know what the opponent plays at given round. The simultaneity of the moves implies that there is no deterministic good strategy in this game, otherwise one would anticipate the moves of the opponent and easily win the game. We explore how to build a family of one-dimensional stochastic cellular automata to play this game. The rules are built in an iterative way by progressively increasing the complexity of the transitions. We show the possibility to construct a family of rules with interesting results, including a good performance when confronted to the Nash-equilibrium strategy.

Introduction

The purpose of this note is to present a sketch on how stochastic cellular automata can be used to play a simple strategy game in which randomness has a central role. This game, named Alesia after the battle that opposed Gallic tribes and the army of J. Caesar 52 BC, has simple rules [?]: (1) The two players initially have the same number of soldiers, say 50. (2) Each round, the two players fight a battle by simultaneously engaging a given number of soldiers. The soldiers are then lost, whatever the outcome. Players must engage at least one soldier at each round (if they have not lost all their soldiers). (3) The winner of battle is simply the player who has engaged more troops; the front moves by one step in the direction of his opponent. (4) The winner of the game is the player who succeeds to reach his opponent's camp, that is, to make the front advance more than W steps, where W is fixed in advanced (here we take W = 2). (5) The game ends in a draw if none of the players reaches the opponent's camp.

To illustrate these rules, an example of game is given on Fig. 1 After playing a few games, one discovers that any good strategy should find a balance between playing too much or too little soldiers at each battle. Indeed, a large number of soldiers ensures victory but depletes the reserves for the next rounds. On the contrary, a small number of soldiers saves the reserves but increases the risks to loose the battle.

The other important point is that a good strategy is necessarily stochastic. Indeed, if for instance a player always plays 4 at the first round, the opponent can foresee this decision and decide to play 5 and thus win the first round. If the player always plays 20, the opponent's interest is to play only 1, and to deliberately loose the first round to gain an advantage in the number of soldiers. Of course, this argument can be recursively applied to the next rounds... Our aim is to analyse whether simple cellular automata could calculate a good strategy for this game. Our motivation is to explore how "noisy" components with an elementary behaviour can cooperate to perform an interesting computation, as in a biological organisms. We are aware that cellular automata constitute a Turing-universal model of computation, and can thus compute any function that is computable by a classical machine, but our goal here is to make a decision emerge from simple non-deterministic mechanisms.

The links between cellular automata and game theory have been mainly explored with iterated two-player games on grids or graphs: these models distribute n players on a grid (or on a graph) and make the players interact by pairs with their neighbours. The players generally update their strategy according to the payoff received and the goal is to observe how the strategies dynamically evolve on the graph [?]. Interestingly enough, it can be observed that even small variations in the updating method (e.g., from synchronous to partially synchronous) may radically change the ultimate evolution of the system [?].

However, the use of cellular automata for the design of strategies is almost absent from the literature. Fraenkel has presented a pioneering work of how cellular automata may compute a particular strategy [?]. Cook, Larsson and Neary also made interesting connections between (deterministic) one-dimensional cellular automata and a generalisation of the game of Nim called the Blocking Wythoff Nim [?].

The purpose of this paper is to present some other simple examples by taking the Alesia game as a starting point. It can thus be considered as a pedagogical example in order to solve a problem in a distributed way. Our goal is not to reach optimality, but we will nevertheless compare our solutions to an optimal strategy in order to have a quantitative estimation of the quality of our cellular automata strategy.

Definitions

Formalisation of the game

We now introduce the formal definitions of the game. We assume that the game starts with N soldiers and that the "arena" is such that a player needs to make W (cumulative) steps in the direction of his opponent to win.

For the sake of simplicity we model a game by an infinite sequence of moves. Let a t , b t be the number of soldiers of player A and B, respectively, at time t and α t , β t the number of soldiers that they respectively engage at time t. We have:

∀t ∈ N, a t+1 = a t -α t and ∀t ∈ N, b t+1 = b t -β t .
The position of the front evolves according to w 0 = 0 and:

w t+1 = w t +      1 if α t > β t -1 if α t < β t 0
otherwise. The game stops when one of the players hits the camp of his adversary or when there are no more soldiers. As the outcome of the game depends on the sequence of moves s = (a t , b t ) t∈N , we call this sequence the game for the sake of simplicity.

Since each player is forced to play at least one soldier at each time step, the number of rounds is finite; we denote it by T (s).

The gain G(s) that results from a game s is 1, -1 or 0, if player A wins, or if player B wins, or if there is a draw, respectively. Formally, we have:

T (s) = min t∈N {w t = W + 1 or w t = -W -1 or (a t , b t ) = (0, 0)}, and : G(s) =      1 if w T = W + 1, -1 if w T = -W -1, 0 otherwise.
The rules of the game impose a 0 = b 0 = N and: ∀t < T , a t > 0 =⇒ 1 ≤ α t ≤ a t , and b t > 0 =⇒ 1 ≤ β t ≤ b t and, for t ≥ T , we set α t = β t = 0.

Stochastic cellular automata

For the sake of simplicity, we will use one-dimensional cellular automata with nearest-neighbours interaction. To describe our model, we simply take an infinite line of cells : Z. Note however that for the simulations, it is more convenient to use periodic boundary conditions (Z/nZ with a large value for n ).

We denote by Q the set of states the cells can hold. A configuration represents the global state of the system; it is an element of Q Z .

A stochastic cellular automaton is defined with a local transition function

ϕ : Q 3 × Q → [0, 1]
and we write ϕ (x, y, z), q = p to express that a cell with a neighbourhood state (x, y, z) has a probability p to update to the state q. We require: ∀(x, y, z) ∈ Q 3 , q∈Q ϕ (x, y, z), q = 1. The global transition function thus maps a configuration x t to the configuration x t+1 such that:

∀i ∈ Z, Pr[x t+1 i = q] = ϕ (x t i-1 , x t i , x t i+1 ), q ,
and where all the probabilities are drawn independently for each cell.

In order to build a strategy that uses a cellular automaton, we will translate the state of the game into a configuration, make this configuration evolve, and then interpret the resulting configuration as an action of the player. Given two players A and B, we will measure their relative strength with the expected gain that player A has against player B. This statistical estimation is obtained by repeating N games opposing A and B. If N A and N B are the number of games that A or B won, respectively, then the expected gain of A against B is given by (N A -N B )/N .

We now present how we build a cellular automaton player by progressively improving the strategy players. This improvements are made by defining families of rules and making each new family compete against the previously found players. This method can be seen as a "bootstrapping" technique because each level of complexity emerges from a previously defined level of complexity.

Cellular automata players

All our players are Markovian: they base their decision on the current state of the game only. For the sake of readability, we will assume that the player we are describing is player A, while its opponent is player B. We recall our notations: a, b denote the number of soldiers of player A and B, respectively, and w denotes the position of the front at time t, the number of soldiers played by a given player Π is denoted by Π(a, w, b). We also assume that the game is played in "standard" conditions, that is, with W = 2 and N = 50 [?].

Table 1: Expected gains of uni-M players opposed to other uni-M player. The deviation from an antisymmetrical form are due to the statistical variations (10 5 samples).

uni-10 uni-12 uni-14 uni-16 uni-18 uni-20 uni-22 uni-24 uni-10 -0.

Uniform distribution players

The first player we can consider, denoted by uniform, plays a uniform number between 1 and a. It will serve as the basis of our bootstrapping process. Formally: uni-M(a, w, b) = U{a}, where U{0} = 0 and U{k} draws a random number uniformly in {1, . . . , k} for k > 0. This strategy is easy to defeat with the following opponent: beatUnif always plays 1 unless the front is in position -W (danger of loosing the game), in which case it plays b (as many soldiers as the opponent has). Formally: beatUnif(a, w, b) = b if w = -W and beatUnif(a, w, b) = 1 otherwise.

Our simulations show that beatUnif has an expected gain of 0.93 against unif: it wins the game with a probability greater than 95%.

Uniform distribution players with saturation

The weakness of uniform comes from the fact that it is too "generous": it is easy to "exhaust" simply by using a defensive strategy. A straightforward improvement of this player is to limit the maximum number of soldiers engaged at each time step. We define uni-M as the player which draws a number uniformly between 1 and a and then plays this value if it is lower than a threshold value M , or plays M otherwise. Formally: uni-M(a, w, b) = max{U{a}, M }.

The question then comes to know what is the best setting for M . Again, as there is no "absolute" good player; we thus simply oppose players with various settings of M and observe how they perform one against the other.

The results are presented in Tab. 1. The data represents the expected gain of uni-M players with different settings of M . These experiments indicate that a good setting of M is in the interval 16-18. Indeed uni-16 and uni-18 have a positive expected gain defeat when opposed to the other players and when opposed to each other, the difference in expected gain is not significant.

We are now in position to continue our bootstrapping process: our next objective is to build a CA player that performs better than uni-18. Before going on, let us observe that "coding" the uniform or uni-M players with a onedimensional cellular automaton is not that easy. Indeed, if the input is coded in the form of a configuration that has n cells in a given state, it is not clear how cells could interact locally in order to produce every possible output between 0 and n -1 (or n) with an equal probability.

Another simple player: the binomial player

If we have a set of cells which can hold a state with a given probability, the most intuitive "computation" is to draw a number according to a binomial distribution. We define the binomial players as the strategy where each soldier of the player has a probability ρ "to be played". For a given ρ, we set R = 100 * ρ and denote by bin-R the binomial with parameter ρ = R/100. To encde our player, we simply use the binary alphabet Q = {0, 1} and map a game state (a, w, b), to the initial condition x ∈ Q Z formed by a consecutive 1's on a background of 0's. Then each cell independently applies the rule where a 0 remains a 0 and a 1 remains a 1 with probability ρ and becomes a 0 with probability 1 -ρ. We have: binomial(a, f, b) = card{i ∈ Z, y i = 1}, where y denotes the configuration obtained by a one-step transformation of x.

Table 2 shows the expected gain of various binomial players against various uniform-M players. We observe that for each value of M , there is a different value of ρ which maximizes the expected gain. The player which has the highest minimal expected gain is bin-24. It is interesting to note that if we take ρ = 0.28, the "best" opponent of this binomial players is uni-12 and not uni-18, as expected from what was seen by making uniform-M players together. This illustrates the fact that it is not possible to compare the expected gain of rules with a total order. Given three players A,B,C, we can observe that B performs better than C versus A, and nevertheless B is beaten by C. When opposing the binomial players one against another, the best player is bin-22 (Tab. 3). Our objective is now to find a player that beats this new challenger.

4 Taking into account the situation: the q2h2 player

The previous model was not really a cellular automaton since there was no interaction between cells. The next improvement we can do is to take into account the position of the front and number of soldiers of the opponent. We request that the initial state of the game to be translated in an initial condition of the cellular automaton with a simple method. Typically, the number of soldiers and the position of the front should be coded with a "unary" code, that is, each number of soldiers should correspond to the same number of cells in a given state, plus or minus some constants. Let us now present a player which respects these constraints, we name it q2h2.

The model employs 8 states: Q = {E, A, V, V*, R, F, F*, B}; its elements respectively represent the following states: empty, A-soldier, voluntary, voluntary-star, reluctant, front, front-star, B-soldier. The 'star' represents an information that travels from right to left in order to transmit an influence from the B-soldier cells to the different cells which represent the strength of player A. The following rules describe how this influence is transmitted.

We associate to a game position (a, w, b) an initial configuration x ini built as follows : the a first cells are in state A. Next, we put δ cells in state F, where δ = W + w + 1 represents the state of the front. The next following b cells are in state B. We thus have:

x ini (a, w, b) = ..EE AAAAA a times F F F δ times BBBBBBB b times EEE...
The evolution of the cellular automaton can be described with the following scenario: soldiers of player A (state A) turn to the voluntary state (V) or to the reluctant state (R) in one step. After that, some reluctant soldiers of A may turn to voluntary according to the "danger" they feel. This danger is evaluated as a combination of b and w: a) the more soldiers the opponent has, the greater the danger; b) the smaller the front is, the greater the danger. In practice, each cell in state B has a given probability (p T ) to initiate a signal that will travel to the left until it eventually reaches a cell in state R it then turns this cell to a V. This signal can also be absorbed with probability p A .

As the information travels from the right to the left, we can define the local function ϕ (x, y, z), q with the use of the probabilistic function ξ(q, q ), which takes as an input the state of the cell itself q and the state of right neighbour q and outputs a state in Q with a given probability. This function, which depends on three probabilities p V , p T and p A , is defined as follows.

• An empty cell remains empty: ξ(E, •) = E.

• Each A immediately decides if it turns to a voluntary state or to a reluctant state:

ξ(A, •) = V with probability p V , R with probability 1 -p V .
• The behaviour of front cells and front-star cells depend on what they see on their right: a) a B cell: this corresponds to the case where the "stars" are initiated by the soldiers of the opponent. b) another F or F* cell : this case corresponds to the transmission of the star to the left. We set a probability p A to be absorbed, i.e., the star is not transmitted. This reads:

ξ(F, B) = ξ(F*, B) = F* with probability p T , F
with probability 1 -p T , and:

ξ(F, F*) = ξ(F*, F*) = F* with probability 1 -p A , F with probability p A .
In all other cases, the front cells remains stable, the front-star cells become front cells: ξ(F, q) = F(F*, q) = F for q / ∈ {B, F, F*}. (Note that in a normal behaviour, the only useful case is q = E.)

• A voluntary cell or voluntary-star cell simply transmits the star from right to left. This reads: ξ(V, q ) = ξ(V*, q ) = V* if q ∈ {V*, F*} and ξ(V, q ) = ξ(V*, q ) = V otherwise.

• A refractory cell remains refractory unless it sees a star on its right. This is translated by: ξ

(R, V*) = ξ(R, F*) = V and ξ(R, q ) = R if q / ∈ {V*, F*}.
• Cells in state B simply disappear at the rate of one cell per time step. This reads: ξ(B, E) = E and ξ(B, q

) = B if q = B.
An illustration of this behaviour can be seen on Fig. 1-right. (As it can be easily guessed, the front-star are drawn in red and the voluntary-star cells are in green.)

To analyse this rule, first, let us simply set p A = 0, in other words, we do not take into account the state of the front. We ask how we can tune p V and p T in order to beat the binomial players. Recall that so far our best binomial player is bin-22 (ρ = 0.22). Table 4-left shows the result of this player against bin-22 for different values of p V and p T . It can be seen that this player is easy to defeat: for example, for p V = 0.40 and p T = 0, one obtains an expected gain greater than 0. A zero value indicates that in fact, in this case it is sufficient to take into account only the strength of the opponent to obtain good results. Note that the value of (p V , p T ) which maximises the expected gain against binomial players varies greatly with ρ. For example, for bin-35, the best expected gain is obtained for p V = 0.20 and p T = 0.10, and this gain (∼0.15) is much lower than for bin-22. This may seem paradoxical but remember that bin-22 is the player which has the best result against the other binomial players, and that we know nothing for the other rules. All we can say is that we are sure that even with p A = 0 the family of q2h2 players dominates the family of binomial players: for every binomial player, there exists a q2h2 player which can defeat it. This is a direct consequence of the fact that the q2h2 players include all the binomial players simply by setting p T = 0.

Table 5 show how various binomial players can be defeated by fixing p V = 0.22 and setting an appropriate value of p T . Similar results can be obtained for other values of p V , which confirms the great advantage of the q2h2 players against the binomial players.

Despite these encouraging results, we could not find any setting of (p V , p T ) which would dominate all the other binomial players. Once again, the difficulty stems from the impossibility to establish a total order between rules. For example the rule with the setting (p V , p T ) = (0.5, 0.8) defeats bin-22 (with an expected gain of 0.60) and is defeated by bin-35 (with an expected gain of 0.67)... but, as seen above, bin-22 defeats bin-35!

In a second step, to demonstrate how setting a positive value of p A can be useful, simply examine a precise situation and leave a more systematic study for future work. In the paragraph above we showed how setting p A = 0 and varying p V and p T allows us to defeat every binomial player. However, against bin-30, only a small region of (p V , p T ) has a positive expected gain, and the maximum positive gain one can obtain is around 0.1.

Table 6-left shows the expected gain of q2h2 against bin-30 when we set p T = 0.25 and vary p A and p V . It can be seen that allowing p A > 0 gives much better results: in particular for p V = 0.15 and p A = 0.8, the expected again is above 0.4, which is quite impressive.

We have seen that it is not possible to totally order the players according to their respective scores. There exists however a strategy which is never defeated (on average). We denote this player by Nash, as it corresponds to what is called the Nash equilibrium: Nash gives a guarantee to have a positive expected gain against any player P, but against a third player P', it may be that P performs better than Nash.

We have tested the performance of q2h2 against Nash, with a setting of p T = 0.25 and different values of p V and p A . The results are displayed on Tab. 6-right: we see that in some cases, the average expected loss is close to zero and of the same amplitude as the noise on the measures. In practice, this implies that one can hardly distinguish between the Nash player and the q2h2 player with the proper settings.

These first results are rather encouraging and there are many directions in which they can be deepened. For example, it is interesting to examine the scaling properties of our cellular automata: obtaining the optimal Nash-equilibrium player demands more time as the size of the game increases, our models can easily be applied to larger sizes, maybe with an adjustment of the three probabilities. Of course, the behaviour of this cellular automaton can also be obtained with classical mathematical functions but here we wanted to examine how simple interacting elements would play this game. Our goal was to show that a non-trivial behaviour can be obtained by progressively increasing the complexity of the rules. Another research direction would be to make the system evolve autonomously and see if it can discover new levels of complexity without an external aid.

acknowledgements

We express our sincere gratitude to Bruno Scherrer for introducing us to the game of Alesia and for providing us with an optimal Nash player program. We thank Irne Marcovici for her valuable comments on the manuscript.

Table 2 :

 2 Expected gains of binomial players against uniform players (10 5

	samples).
		uni-8 uni-10 uni-12 uni-14 uni-16 uni-18 uni-20
	bin-22	0.554	0.409	0.370	0.375	0.368	0.377	0.391
	bin-24	0.589	0.428	0.400	0.411	0.430	0.441	0.463
	bin-26	0.608	0.433	0.399	0.417	0.451	0.475	0.501
	bin-28	0.622	0.428	0.380	0.405	0.448	0.485	0.513
	bin-30	0.628	0.420	0.345	0.368	0.420	0.468	0.504
	bin-32	0.627	0.398	0.314	0.319	0.374	0.431	0.480

Table 3 :

 3 Expected gains of binomial players against other binomial players (10 5 samples).

		bin-15 bin-20 bin-25 bin-30
	bin-15		-0.182	-0.327	-0.417
	bin-20	0.180		-0.002	-0.021
	bin-25	0.328	0.003		0.098
	bin-30	0.419	0.017	-0.099
		bin-20 bin-22 bin-24 bin-26
	bin-20		-0.004	-0.005	0.000
	bin-22	0.006		0.015	0.039
	bin-24	0.005	-0.018		0.035
	bin-26	0.005	-0.030	-0.029

Table 4 :

 4 Expected gains of q2h2 players with p A = 0 against bin-22 (left) and bin-35 (right). The two parameters p T (columns) and p V (lines) are varied (2.10 4 samples).

			p T
	p V	0.10	0.20	0.30	0.40	0.50
	0.0 -0.986 -0.189 0.470	0.637 -0.240
	0.5 -0.565	0.133 0.530	0.346 -0.656
	0.10 -0.070	0.264 0.507 -0.076 -0.864
	0.15	0.021	0.342 0.343 -0.500 -0.953
	0.20	0.062	0.328 0.023 -0.774 -0.989
		0.10	0.20	0.30	0.40	0.50	0.60
	0.0 -0.997 -0.870 -0.212 -0.026 -0.246 -0.844
	0.10 -0.415	0.007 -0.111 -0.311 -0.776 -0.984
	0.20	0.107 -0.171 -0.382 -0.753 -0.971 -0.999
	0.30 -0.205 -0.436 -0.735 -0.960 -0.998 -1.000
	0.40 -0.501 -0.743 -0.948 -0.996 -1.000 -1.000
	0.50 -0.778 -0.949 -0.997 -1.000 -1.000 -1.000

Table 5 :

 5 Expected gains of binomial players (lines) against q2h2 players (columns) with p V = 0.22, p A = 0 and a varying value of p T (10 4 samples).

	p T	0.0	0.5	0.10	0.15	0.20	0.25	0.30
	bin-15 -0.330 -0.541 -0.650 -0.709 -0.607 -0.276 0.277
	bin-20	0.007 -0.114 -0.269 -0.389 -0.376 -0.102 0.323
	bin-25 -0.008	0.086	0.058 -0.020 -0.048	0.083 0.387
	bin-30 -0.082	0.015	0.157	0.206	0.203	0.286 0.478
	bin-35 -0.192 -0.143	0.009	0.188	0.326	0.408 0.554

Table 6 :

 6 Expected gains of q2h2 players with p T = 0.25 against bin-30 (left) and against Nash (right). The two parameters p A (lines) and p V (columns) are varied (10 5 samples).

	0.10	0.15	0.20	0.25	0.0	0.2	0.4	0.6	0.8	0.10
	0.20 0.111 -0.065 -0.256 -0.370	0.20 -0.050 -0.058 -0.075 -0.095 -0.115 -0.148
	0.40 0.254	0.208 -0.010 -0.213	0.40 -0.039 -0.033 -0.043 -0.056 -0.063 -0.088
	0.60 0.280	0.348	0.178 -0.027	0.60 -0.035 -0.029 -0.027 -0.029 -0.039 -0.048
	0.80 0.276	0.395	0.272	0.081	0.80 -0.093 -0.051 -0.038 -0.029 -0.027 -0.038
	0.100 0.267	0.403	0.280	0.108	0.100 -0.183 -0.139 -0.081 -0.046 -0.043 -0.033


Madalena Chaves
email: madalena.chaves@inria.fr

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

Dynamics of complex feedback architectures in metabolic pathways

Keywords: Metabolic pathways, synthetic biology, stability analysis, piecewise linear models, gene regulatory networks

Cellular metabolism contains intricate arrays of feedback control loops. These are a key mechanism by which cells adapt their metabolism to survive environmental disturbances. Gene regulation, in particular, typically displays complex architectures that combine positive and negative feedback loops between metabolites and enzymatic genes. Yet because of strong nonlinearities and high-dimensionality, it is challenging to determine the closed-loop dynamics of a given feedback architecture. Here we present a novel technique for the analysis of metabolic pathways under gene regulation. Our theory blends ideas from timescale separation and piecewise affine dynamical systems, applied to a wide class of unbranched metabolic pathways under steep nonlinear feedback. We propose a systematic method to construct a state transition graph for a given regulatory architecture, from where candidate closed-loop dynamics can be singled out for further analysis. The method recasts a high-dimensional nonlinear system into a piecewise affine system defined on a polytopic partition of the state space. In its most general setup, our theory allows to characterize the dynamics of pathways of arbitrary length, with any number of regulators, and under any combination of positive and negative feedback loops. We illustrate our results on an exemplar system that displays a stable limit cycle and bistable dynamics. We also discuss the implications of our results for the design of gene circuits in synthetic biology and metabolic engineering.

Introduction

Cellular metabolism adapts to environmental fluctuations through networks of regulatory feedback loops [START_REF] Chubukov | Coordination of microbial metabolism[END_REF]. In natural systems, such feedbacks control the expression of catalytic enzymes and match the activity of specific metabolic pathways to extracellular nutrients and the cellular energy budget. Last decades have witnessed striking examples on the construction of synthetic feedback control loops for metabolism, with the goal of increasing production of valuable chemicals [START_REF] Farmer | Improving lycopene production in Escherichia coli by engineering metabolic control[END_REF][START_REF] Keasling | Synthetic biology and the development of tools for metabolic engineering[END_REF][START_REF] Xu | Improving fatty acids production by engineering dynamic pathway regulation and metabolic control[END_REF] and the generation of entirely new metabolic dynamics such as oscillations [START_REF] Fung | A synthetic genemetabolic oscillator[END_REF] or cell-to-cell communication [START_REF] Gupta | Dynamic regulation of metabolic flux in engineered bacteria using a pathwayindependent quorum-sensing circuit[END_REF][START_REF] Soma | Self-induced metabolic state switching by a tunable cell density sensor for microbial isopropanol production[END_REF].

Regulatory motifs in metabolic systems include architectures with various combinations of feedback and feedforward loops, whose signs and interactions determine the overall pathway dynamics. A fundamental design problem is the characterization of control architectures that achieve a prescribed system response [START_REF] Stevens | Designing RNAbased genetic control systems for efficient production from engineered metabolic pathways[END_REF][START_REF] Oyarzún | Design of a bistable switch to control cellular uptake[END_REF]. However, there are no general stability analysis methods that can deal with the type of nonlinearities and connectivity of metabolic systems. A number of recent reviews in metabolic engineering have highlighted the need for model-based approaches for design [START_REF] Irene | Dynamic metabolic engineering: New strategies for developing responsive cell factories[END_REF][START_REF] He | Synthetic biology and regulatory networks: where metabolic systems biology meets control engineering[END_REF][START_REF] Liu | Dynamic metabolic control: towards precision engineering of metabolism[END_REF], while the use of control theoretical principles has gained substantial traction in the field [START_REF] Oyarzún | Synthetic gene circuits for metabolic control: design trade-offs and constraints[END_REF][START_REF] Del | Control theory meets synthetic biology[END_REF].

In this paper we provide two contributions on the analysis of control architectures in metabolic pathways under gene regulation. First, we show that the interaction graph of a large family of feedback systems can be classified into few distinct classes of interaction graphs. Our findings suggest a catalogue of architectures in the form of interactions graphs with nested positive or negative loops, where the number of distinct classes depends only on pathway length (Section 3). This catalogue links each feedback architecture to a systematized feedback control network and its corresponding dynamical behav-ior. Second, we provide an algorithm to build a state transition graph for a given feedback architecture, from which candidate asymptotic dynamics can be identified and further analyzed (Section 4). The analysis relies on the combination of timescale separation and qualitative, piecewise-affine, models for gene regulation. Here we extend a technique originally described in [START_REF] Oyarzún | Multistability and oscillations in genetic control of metabolism[END_REF] for pathways with one regulator, to a much larger class of control systems with any number of regulators and any combination of positive and negative feedback loops. We demonstrate our theory with an example motivated by [START_REF] Fung | A synthetic genemetabolic oscillator[END_REF] that illustrates the dynamical diversity within each class of feedback architectures (Section 5).

Metabolic pathways with multiple regulators

Metabolic pathways interconvert chemicals (metabolites) through sequences of biochemical reactions. Such reactions are catalyzed by specific proteins (enzymes) synthesized by the cell. By controlling the expression of genes that code for the enzymes, cells activate and deactivate the activity of pathways depending on the environmental conditions. Because enzyme synthesis requires energy, this control mechanism avoids the inefficient production of enzymes that are not needed in a specific environment. This is the case, for example, when key metabolites for growth are already present in the growth media and do not need to be produced by cells themselves.

Here we consider a general class of unbranched metabolic pathways with n metabolites and n + 1 enzymes, as shown in Fig. 1. The conversion of each metabolite s i into s i+1 is catalyzed by an enzyme e i+1 . The substrate s 0 represents an extracellular nutrient and we assume it to remain constant in the timescale of the pathway. This assumption corresponds to cases in which, for example, s 0 represents a pool of extracellular nutrient in continuous culture. Under a constant substrate, the pathway reaches a nonzero steady state flux and remains outside thermodynamic equilibrium, which is the typical operating regime of metabolism.

We further consider that enzyme synthesis is controlled by some of the metabolites in the pathway. We consider a large class of feedback architectures satisfying the following assumptions: A1. Each enzyme is regulated by a single metabolite; A2. There are 1 ≤ k ≤ n regulatory metabolites.

Assumptions A1-A2 define a large class of feedback architectures that may include complex combinations of positive and negative feedback loops. This is a significant extension to previous work [START_REF] Oyarzún | Multistability and oscillations in genetic control of metabolism[END_REF] that considered the case of only one regulatory metabolite (k = 1). In particular, Assumption A2 allows for any number of regula- tory metabolites to control enzyme expression. In most natural instances, e. g. the lactose operon or amino acid metabolism [START_REF] Chubukov | Coordination of microbial metabolism[END_REF], pathways have one regulatory metabolite (k = 1), with the case of k > 1 being particularly interesting for the design of regulatory architectures for metabolic engineering [START_REF] Liu | Dynamic metabolic control: towards precision engineering of metabolism[END_REF].

If s i and e i denote the concentration of metabolites and enzymes respectively, the system in Fig. 1 can be described by the nonlinear model:

ṡi = g i (s i-1 )e i -g i+1 (s i )e i+1 , i = 1, . . . , n, (1) ėi
= κ 0 i + κ 1 i σ i (s i , θ i ) -γe i , i = 1, . . . , n + 1, (2)
where s i is the metabolite that regulates e i , g i : R + ≥0 → [0, M i ] are continuous and strictly increasing functions representing the kinetics of each enzyme [START_REF] Heinrich | The regulation of cellular systems[END_REF], and σ i ∈ {σ + , σ-} are increasing or decreasing sigmoidal functions describing the nonlinear regulation from metabolites to enzyme synthesis [START_REF] Mannan | Fundamental Design Principles for Transcription-Factor-Based Metabolite Biosensors[END_REF][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]. The positive parameters (κ 0 i , κ 1 i ) represent the minimal and maximal rates of enzyme synthesis, and γ > 0 describes the dilution of protein concentrations caused by cellular growth.

Our objective in this paper is to characterize the asymptotic dynamics of nonlinear feedback systems described by equations ( 1)- [START_REF] Batt | Temporal logic analysis of gene networks under parameter uncertainty[END_REF]. To this end, we will use two techniques for simplifying the models [START_REF] Oyarzún | Multistability and oscillations in genetic control of metabolism[END_REF]: first, we approximate enzyme synthesis by piecewise linear functions; and second, we exploit the inherent separation of timescales and the shape of the nonlinearities of the enzyme kinetics.

Piecewise linear models for gene regulation

Gene expression can often be approximated by two levels, corresponding to "high" and "low" levels of transcription. In such cases, the sigmoidal functions that describe enzyme synthesis can be simplified to step functions, following the approaches suggested by Thomas [START_REF] Thomas | Boolean formalization of genetic control circuits[END_REF] and Glass and Kauffman [START_REF] Glass | The logical analysis of continuous, nonlinear biochemical control networks[END_REF]:

σ + (r, θ) = 1 -σ -(r, θ) = 1, r > θ, 0, r < θ,
where σ + (resp., σ -) corresponds to an activation (resp., inhibition), and the parameter θ is a threshold of regulation. In this model, enzyme concentrations eventually converge to a bounded interval, E off i ≤ e i (t) ≤ E on i for all t ≥ t 0 > 0, with

E off i = κ 0 i γ , E on i = (κ 0 i + κ 1 i ) γ . (3)
The step function approximation leads to a qualitative description of the dynamics in terms of piecewise linear systems which, in turn, allows for more detailed analytical results (see examples for Escherichia coli metabolic networks [START_REF] Oyarzún | Multistability and oscillations in genetic control of metabolism[END_REF][START_REF] Ropers | Model reduction using piecewise-linear approximations preserves dynamic properties of the carbon starvation response in Escherichia coli[END_REF] and genetic transcription and translation networks [START_REF] Edwards | A modelling framework for gene regulatory networks involving transcription and translation[END_REF][START_REF] Hudson | Dynamics of transcriptiontranslation networks[END_REF]). Piecewise linear systems induce a partition of the state space into well defined domains which can be naturally mapped onto the vertices of a state transition graph [START_REF] Batt | Temporal logic analysis of gene networks under parameter uncertainty[END_REF][START_REF] Casey | Piecewise-linear models of genetic regulatory networks: equilibria and their stability[END_REF]; this graph provides a global view of the dynamics of the system and will be constructed in Section 4.

Timescale separation

Metabolic timescales are much faster than those of enzyme synthesis [START_REF] Heinrich | The regulation of cellular systems[END_REF][START_REF] Ropers | Model reduction using piecewise-linear approximations preserves dynamic properties of the carbon starvation response in Escherichia coli[END_REF], so the system (1)-( 2) can be simplified by assuming the metabolites to be in quasi-steady state. The variables s i evolve fast towards a slow manifold of the form s = G(e), where s and e are the vectors of metabolite and enzyme concentrations, respectively. By substitution of the metabolites s i in the ODEs for the enzymes, the slow dynamics in (2) can be simplified to depend only on the enzyme concentrations e i .

Such timescale separation ideas have been used for analysis of various enzymatic mechanisms [START_REF] Baldazzi | Importance of metabolic coupling for the dynamics of gene expression following a diauxic shift in Escherichia coli[END_REF][START_REF] Kuntz | Model reduction of genetic-metabolic systems using timescale separation[END_REF], genetic transcription and translation networks [START_REF] Edwards | A modelling framework for gene regulatory networks involving transcription and translation[END_REF], and justify the well known Michaelis-Menten kinetic mechanism [START_REF] Segel | On the validity of the steady state assumption of enzyme kinetics[END_REF]. Briefly, the goal is to rewrite (1)-(2) in the form: ε ṡ = G(s, e, ε), ė = F (s, e, ε), where ε 1. Under suitable conditions (see [START_REF] Khalil | Nonlinear systems[END_REF] for more details), G(s, e, 0) = 0 defines the slow manifold s = G(e), and Tikhonov's Theorem guarantees that solutions of the slower dynamics ė = F ( G(e), e, 0) remain close to those of the original system.

In the particular case where all enzymes follow Michalis-Menten kinetics, that is, g i (s i ) = k cat i s i-1 /(K M,i + s i-1 ) (see example in Section 5), the timescale separation can be applied to (1)-(2) provided that

ε i = max γK M,i E on i k cat i , γK M,i+1 E on i k cat i 1, (4)
for i = 2, . . . , n. These conditions can be obtained by applying an appropriate change of variables: ẽi = e i /E on i , si = s i /K M,i , and t = γt, similar to those used in [START_REF] Kuntz | Model reduction of genetic-metabolic systems using timescale separation[END_REF][START_REF] Segel | On the validity of the steady state assumption of enzyme kinetics[END_REF]. The conditions typically hold for metabolic systems, since dilution by growth is much slower than enzyme turnover rates, and hence γ k cat i . The quasi-steady state assumption is: A3. We assume that ṡi ≈ 0 for each t ≥ 0, leading to g 1 (s 0 )e 1 = g i (s i-1 )e i , for i = 2, . . . , n + 1. [START_REF] Chaves | Analysis of a genetic-metabolic oscillator with piecewise linear models[END_REF] Using assumption A3 we can substitute

s i = g -1 i+1 g 1 (s 0 ) e 1 e i+1 (6)
into equation [START_REF] Batt | Temporal logic analysis of gene networks under parameter uncertainty[END_REF] for enzyme e i . To guarantee the existence of solutions of (6), a sufficient condition is:

E off i ≥ g 1 (s 0 ) M i+1 E on i , for i = 2, . . . , n + 1. (7)
Now, let I reg = {1, 1 +1, . . . , k +1} be the set of indices corresponding to the enzymes that catalyze the outgoing flux for each of the regulator metabolites together with index 1, the enzyme which catalyzes the constant input substrate; set I nrg = {1, 2, . . . , n} \ I reg . From [START_REF] Chubukov | Coordination of microbial metabolism[END_REF], it becomes clear that the subset of k + 1 equations corresponding to {e i : i ∈ I reg }, does not depend on the variables {e j , i ∈ I nrg }, so it can be studied separately. Therefore, without loss of generality we will thus consider a modified version of Assumption A2: A2'. Each metabolite in the pathway regulates at least one enzyme, or equivalently, n ≡ k.

Taken together points A1 and A2' imply that there is only one metabolite (denoted as s p )that regulates two enzymes and all other metabolites regulate a single enzyme. It is then useful to establish an index correspondence:

I : {1, . . . , n + 1} → {1, . . . , n} (8)
with I(i) = i and I(p 1 ) = I(p 2 ) = p. Its inverse is defined as: I -1 (j) = q j and I -1 (p) = {p 1 , p 2 }.

Enzyme dynamics and regulatory motifs

The unbranched metabolic pathway dynamics can thus be reduced to the enzyme equations. In order to capture this genetic regulation, we will introduce, with a slight abuse of notation, the following definition of a regulatory motif :

σ = (σ 1 , . . . , σ n ), σ i = 1, σ i (r, θ) = σ + (r, θ) -1, σ i (r, θ) = σ -(r, θ). (9)
In addition, the quasi steady state assumption ( 6) and the fact that g i are non-decreasing imply

s i > θ i ⇔ g 1 (s 0 )e 1 > g i+1 (θ i )e i+1 ⇔ e 1 > 1 β i e i+1 . ( 10
)
where we have defined the parameter

β i = g 1 (s 0 )/ g i+1 (θ i ).
Under the relations in [START_REF] Farmer | Improving lycopene production in Escherichia coli by engineering metabolic control[END_REF], the step functions satisfy

σ i (s i , θ i ) = σ i (s i -θ i , 0) = σ i (e 1 - 1 β i e i+1 , 0), (11)
and therefore the enzyme dynamics in ( 2) can be rewritten as

ėi = κ 0 i + κ 1 i σ i e 1 - 1 β i e i+1 , 0 -γe i . (12)
Equation ( 12) is a much simpler version of the original model in ( 1)-( 2) and, as we show in the next sections, it is amenable to detailed analysis.

General structure of the interaction graph

Given a metabolic pathway such as the one in Fig. 1, one would like identify and catalogue the possible feedback networks induced by its specific regulatory motif σ. In this section we provide such classification in terms of the possible interactions graphs associated to a regulatory motif. The pattern of feedback loops and their signs determine the dynamical behavior of the system: for instance, it has been shown that a necessary condition for multistability is the existence of a positive loop, while a negative loop is required for stable oscillations [START_REF] Gouzé | Positive and negative circuits in dynamical systems[END_REF][START_REF] Soulé | Graphical requirements for multistationarity[END_REF][START_REF] Richard | Necessary conditions for multistationarity in discrete dynamical systems[END_REF][START_REF] Remy | Graphics requirement for multistability and attractive cycles in a Boolean dynamical framework[END_REF]].

An interaction graph describes the dependencies between model variables according to the sign pattern of the Jacobian matrix [START_REF] Jong | Modeling and simulation of genetic regulatory systems: A literature review[END_REF]. Although the Jacobian of ( 12) has a non-constant sign pattern, next we describe a suitable change of variables under which the interaction graph can be readily built. We first note that the right hand side of ( 12) depends on sign of e 1 -1 βi e i+1 , and thus we introduce the following definition:

Definition 1 [Cone membership function] Let β > 0, r ∈ R n+1
>0 , and r i the i th component of r. Define the cone membership function by

ν i (β, r) : R >0 × R n+1 >0 -→ R (β, r) -→ r 1 - 1 β r i+1 .
In other words, sign of ν i determines whether or not the (n + 1)-dimensional vector r belongs to the planar cone defined by a half-line with slope β.

The new variables can now be defined as:

x j = ν j-1 (β qj-1 , e), j ∈ {2, . . . , n + 1} \ {p + 1}, z 1 = ν p (β p1 , e), z 2 = ν p (β p2 , e),
From the model in [START_REF] Glass | The logical analysis of continuous, nonlinear biochemical control networks[END_REF], it follows that:

ẋj = κ0 1 + κ 1 1 σ 1 (x 1+1 , 0) - κ 1 j β qj-1 σ j (x j +1 , 0) -γx j , j ∈ {2, . . . , n + 1} \ {p + 1} (13) żj
= κ0 pj + κ 1 1 σ 1 (x 1+1 , 0) - κ 1 p+1 β pj σ p+1 (x p+1+1 , 0) -γz j , j = 1, 2.
where κ0

j = κ 0 1 -κ 0 j /β qj-1 and κ0 pj = κ 0 1 -κ 0 p+1 /β pj .
The Jacobian of ( 13) is sign-constant and therefore we can build the interaction graph for the new variables by adding a directed arrow from x j to x i if ∂( ẋi )/∂x j = 0, and (for simplicity) omitting all arrows related to the natural degradation terms -γx j or -γz j . This interaction graph has some immediate properties: R1 all nodes have exactly two incoming arrows; R2 node x 1+1 has n + 1 outgoing arrows, signed σ 1 ; R3 node x p+1+1 has exactly two outgoing arrows, both signed -σ p+1 , towards z 1 and z 2 ; R4 all other nodes x j +1 have exactly one outgoing arrow, signed -σ j , towards x j . In particular, node z j has an outgoing arrow, signed -σ pj , towards x pj .

A self arrow counts as both an incoming and an outgoing arrow for the node. If a node has an outgoing arrow to x p1 , then it has necessarily an outgoing arrow to x p2 , hence that node is either x 1+1 or x p+1+1 .

The above properties suggest that the large family of feedback systems described by equations ( 1)-( 2) can be classified into relatively few distinct interaction graphs, based on the identification of four special nodes:

i the node w = x 1+1 with n + 1 outgoing arrows; ii the node y = x p+1+1 with two outgoing arrows; iii the two nodes z 1 and z 2 (which can coincide with either w or y, when p j = p+1 + 1, for some j ∈ {1, 2}); iv finally, list the remaining nodes, x j +1 for j / ∈ {p + 1, 1}, which can be re-labeled x a , x b , etc.

To obtain all the distinct interaction graphs, consider all combinations of z 1 and z 2 according to (iii), and add an edge between two of these nodes according to rules (R1)-(R4). The possible interaction graphs are listed in Fig. 2, for pathways with three or four enzymes. Next we consider a specific example to illustrate the construction of the interaction graphs. Example 2 Consider the pathway in Fig. 1 (top), which illustrates assumption A2'. There are only two regulators, s 1 and s 2 , so it follows from ( 6) that e 4 will not influence the dynamics of e i , i = 1, 2, 3. Since s 2 regulates two enzymes, it follows that p = 2 and, from the definition in (8), the index correspondence is

-σ p +1 y z 2 w -σ p +1 σ 1 σ 1 -σ p +1 -σ p 2 y z 2 w -σ p +1 σ 1 -σ p 2 y z 2 w y z 1 w z 2 y x a w z 2 y x a w z 2 -σ p 2 x b x a w y -σ a y x a w z 2 y x a w z 2 B A σ 1 σ 1 σ 1 σ 1 σ 1 σ 1 -σ p +1 -σ p +1 -σ p 2 -σ p +1 -σ p +1 -σ p +1 -σ p +1 -σ p +1 -σ p +1 -σ p +1 -σ p +1 -σ p +1 -σ p +1 -σ p +1 -σ p +1 σ 1 σ 1 σ 1 σ 1 σ 1 σ 1 σ 1 σ 1 σ 1 σ 1 σ 1 σ 1 σ 1 σ 1 σ 1 σ 1 σ 1 σ 1 σ 1 σ 1 σ 1 σ 1 σ 1 σ 1 -σ a -σ a -σ b -σ a -σ p 2 -σ p 2 -σ a -σ p 2 -σ p 2 -σ p 1
( 1 , 2 , 3 ) = (1, 2, 2), (q 1 , p) = (1, {2, 3}).
The change of variables is as follows:

x 2 = e 1 - 1 β 1 e 2 , z 1 = e 1 - 1 β 3 e 3 , z 2 = e 1 - 1 β 2 e 3 ,
with w = x 2 and y = z 1 , yielding an interaction graph as in Fig. 2A(middle).

To classify the interaction graphs, a first observation is that all (strongly connected) graphs with n nodes exhibit n nested loops of growing lengths, starting with the self loop on w, each loop adding only one new variable. For instance, for the first case in Fig. 2A:

w ; w y; w → z 2 → y → w,
and for the first case in Fig. 2B w ; w z

1 ; w → y → z 1 → w; w → z 2 → y → z 1 → w.
In fact, this result can be generalized as follows:

Proposition 3 [Nested loops] Consider the n + 1 dimensional system [START_REF] Gupta | Dynamic regulation of metabolic flux in engineered bacteria using a pathwayindependent quorum-sensing circuit[END_REF]. If the interaction graph is formed of a single strongly connected component, and if either w = z j or y = z j , for all j ∈ {1, 2}, then it contains n + 1 nested loops, {L } n+1 =1 , defined by:

L 1 = {w}, L 2 = {w, v 1 : w → v 1 , v 1 → w}, L +1 = {w, v 1 , . . . , v : w → v , v j+1 → v j , v 1 → w},
where j = 1, . . . , -1 in L +1 . If the interaction graph contains more than one strongly connected component, than one of these components has j+1 nodes and contains nested loops L 1 ,. . .,L j+1 .

Proof. This structure follows directly from the properties (R1)-(R4) and (i)-(iv). First, note that w has an arrow toward every node, including itself, so w forms a trivial first loop, L 1 . Second, choose the node v 1 = w which has an edge towards w: by point (R1), each node has exactly two incoming arrows (with one of them w), hence there is a unique choice for v 1 . This gives a 2-node circuit L 2 . Now, to use an induction argument, suppose we have chosen j ≥ 1 distinct nodes, v 1 , . . . , v j , defining loops L 1 up to L j+1 . Note that, by construction, each of the nodes v i , i < j has exactly two incoming arrows, w and v i+1 . If j = n, all the loops have been found and the interaction graph is formed of a single strongly connected component. If j < n, choose next the node v j+1 = w which has an edge towards v j : again, by point (R1) there can be only one such node. If v j+1 = v j , then we have found loop L j+2 and can continue. Otherwise, v j+1 = v j , which means that v j = y since it has two outgoing arrows. But it also has two incoming arrows, w and v j . This means that the set of nodes w, v 1 , . . . , v j (with j < n) cannot have any further incoming arrows and thus forms a strongly connected component, which contains the nested loops L 1 ,. . .,L j+1 . Finally, the remaining nodes (v j+1 , . . . , v n ) can be organized into at least one more strongly connected component. This finishes the proof of the Proposition.

Remark 4 In the case both w = z 1 and y = z 2 it follows that w y, as well as self-arrows for both w and y, thus preventing the generation of n + 1 nested loops. This is illustrated in Fig. 2B bottom left.

As a corollary, it follows that an n + 1 dimensional system [START_REF] Gupta | Dynamic regulation of metabolic flux in engineered bacteria using a pathwayindependent quorum-sensing circuit[END_REF] whose interaction graph has more than one strongly connected component can in fact be studied as a (lower) (j + 1)-dimensional system with nested loops L 1 ,. . .,L j+1 . Indeed, the dynamics of this sub-system w, v 1 , . . . , v j can be studied as an autonomous system, whose behaviour is an input to the system formed by the remaining variables v j+1 , . . . , v n . Examples are shown in Fig. 2B (bottom middle and right) which can both be decomposed into a 3-dimensional system (w, y, z 2 ) which controls x a . Interestingly, these 3-dim systems share the same graphs of Fig. 2A (left and middle). In other words, the systems whose graphs are represented in Fig. 2B(bottom middle and right) can be further simplified and studied as a two-enzyme metabolic pathway.

In addition, the number of distinct strongly connected interactions graphs can be exactly computed, as shown in the next Proposition.

Proposition 5 [Number of distinct interaction graphs] System (13) has at most four (for n = 2) and nine (for n ≥ 3) distinct strongly connected interaction graphs with n + 1 nested loops, up to re-ordering the x i nodes:

(a) two distinct graphs, if w = z 1 or w = z 2 ; (b) two distinct graphs, if y = z 1 or y = z 2 ; (c) five distinct graphs if n ≥ 3 and all w, y, z 1 , z 2 are different.

Proof. (a) In this case, the first two loops are L 1 = {w} and L 2 = {w, y} (i.e., v 1 = y, in the notation of Prop. 3). Chosing the sequence of nodes v i as in Prop. 3, recall that v j = z 2 implies v j+1 = y, by definition of y.

To guarantee that the graph has a single strongly connected component, it must be that v n = z 2 . Otherwise, if v j = z 2 for j < n, then v j+1 = y , which implies that L j+1 ≡ L j+1 = {w, v 1 , . . . , v j } is a strongly connected component and so the full graph cannot be strongly connected.

(b) In this case, assume y = z 1 without loss of generality. Node y has a self arrow (which counts as both an outgoing and an incoming arrow), an incoming arrow from node w and an outgoing arrow towards z 2 . To guarantee that the graph has a single strongly connected component, it must be that v n = y and v n-1 = z 2 . Otherwise, by the same argument as in case (b), if v j = y for some j < n, the set of nodes {w, v 1 , . . . , v j } forms a strongly connected component and so the full graph cannot be strongly connected.

(c) In this case, the second incoming arrow for node w can originate either from one of the variables z i or from one of the n -3 remaining variables x i , up to reordering. Suppose first that z 1 → w (or z 2 → w). Then, as in case (b), we need v n = y and v n-1 = z 2 to garantee a strongly connected interaction graph (count 2 distinct graphs). Suppose next that x a → w and z 1 → y (or z 2 → y). Then, we again need v n = y and v n-1 = z 2 (count 2 distinct graphs). Suppose finally that x a → w and x b → y. The remaining nodes can be assigned in any order, according to rules (R1)-(R4) (count only 1 distinct graph) .

Proposition 5 is illustrated in Fig. 2. The interaction graphs associated with [START_REF] Gupta | Dynamic regulation of metabolic flux in engineered bacteria using a pathwayindependent quorum-sensing circuit[END_REF] give an indication of the possible dynamics of the feedback system (for instance, multistability is expected whenever all loops are positive [START_REF] Thomas | Multistationarity, the basis of cell differentiation and memory. I. structural conditions of multistationarity and other nontrivial behavior[END_REF]).

Global dynamics

As a counterpart to the interaction graph, next we develop a state transition graph that allows for analysis of the global dynamics. To this end, the system in equation ( 12) will be represented by a discrete system formed by a set of vertices V and a set of directed edges E: each vertex represents one region of the state space and there exists an edge between two vertices if there is at least one trajectory connecting the corresponding regions. The resulting state transition graph provides a complete qualitative view of the global dynamics of ( 12), for a given range of parameters.

The system in equation ( 12) can be written as the piecewise linear system (see also [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: equilibria and their stability[END_REF])

ėi = γ(E σi i -e i ), e 1 > 1 βi e i+1 γ(E -σi i -e i ), e 1 < 1 βi e i+1 , (14)
where σ = (σ 1 , . . . , σ n ) represents the regulatory motif, and (E -1 i , E 1 i )= (E off i , E on i ) are the minimal and maximal enzyme concentrations in equation [START_REF] Irene | Dynamic metabolic engineering: New strategies for developing responsive cell factories[END_REF]. In this new description, the state space is partitioned into several conic regions defined by the two planes z j = 0 and the n -1 planes x j = 0 (see [START_REF] Gupta | Dynamic regulation of metabolic flux in engineered bacteria using a pathwayindependent quorum-sensing circuit[END_REF]). Each region has a focal point φ = (φ 1 , . . . , φ n ) with φ i = E σi i if ν i (β i , e) > 0 and φ i = E -σi i otherwise. If a trajectory crosses the boundary between two regions, a switch in the vector field -and hence in the focal point -occurs, leading to a new direction of the trajectory. The conic regions and its focal points are the main elements in the algorithm to construct the state transition graph.

Vertices of the transition graph

As shown in the next definition, each vertex of the graph represents one region R a , characterized in terms of the cone membership functions, and the nonnegative sign function σ + (r, 0). Definition 6 [Polytopic regions] Let p ∈ {1, . . . , n} be the index defined in [START_REF] Del | Control theory meets synthetic biology[END_REF]. A region R a , with a = (a 1 , . . . , a n ) with a p ∈ {0, 1, 2} and a i ∈ {0, 1} for i ∈ {1, . . . , n} \ {p}, is defined by: R a = { e ∈ R n+1 ≥0 :

a i = σ + (ν i (β qi , e), 0), a p = σ + (ν p (β p1 , e), 0) + σ + (ν p (β p2 , e), 0) } .
The set of vertices of the graph is thus given by

V = {(a 1 , . . . , a n ) : a p ∈ {0, 1, 2}, a i ∈ {0, 1} for i = p}.

Edges of the transition graph

Each edge represents a crossing between two adjacent conic regions: if there is a trajectory starting in some region R a and going through a neighbouring region R b , then there is a directed edge from vertex a to vertex b. The path of a trajectory is determined by the location of focal points. In general, the focal point φ a may be contained in some other region, say R b with a = b, so trajectories starting in R a will eventually cross to a neighbouring region.

Here, an asynchronous strategy will be considered in the sense that if an edge joins two vertices, say a → b, then these vertices differ on exactly one coordinate j (a j = b j and a i = b i for all i = j). To determine the edges, the next steps are the calculation of the focal points, φ a = (φ a 1 , . . . , φ a n ), for each region R a and their location R b , i.e. φ a ∈ R b , with a = (a 1 , . . . , a n ) and b = (b 1 , . . . , b n ).

Focal points

The coordinates of the focal point φ a for region R a can be explicitly calculated recalling (14) and using Def. 6:

φ a qi = a i E σq i qi + (1 -a i )E -σq i qi , if i = p, (15)
and in the case i = p:

φ a p1 = max(a p -1, 0)E σp 1 p1 + (1 -max(a p -1, 0))E -σp 1 p1 , φ a p2 = min(a p , 1)E σp 2 p2 + (1 -min(a p , 1))E -σp 2 p2 . (16)

Location of the focal points

Computation of the region R b is straigthforward, as it suffices to check the sign of the functions ν i when computed at the (already known) focal points:

b i = σ + (ν i (β qi , φ a ), 0) b p = σ + (ν p (β p1 , φ a ), 0) + σ + (ν p (β p2 , φ a ), 0). ( 17
)
Observe that the focal points and their location in the state space depend on the parameters of the system (κ 0,1 i , γ, θ i ). Thus, one single regulatory motif ( 9) may lead to different dynamical behavior depending on its parameters. However, as seen below in Section 4.4, once the state transition graph is constructed, it is easy to obtain a large range of parameters for which the graph holds. For simplicity, we assume that no focal point belongs to the boundaries of a cone, that is A4. ν i (β qi , φ a ) = 0 for all i and all a ∈ V .

Construction of the graph edges

There are two cases to consider, depending on whether the given region contains its own focal point. First, assume that φ a / ∈ R a . Then it can be shown that there is at least one boundary hyper-plane e i+1 = β qi e 1 through which trajectories can escape R a .

Proposition 7 Consider a region R a of system [START_REF] He | Synthetic biology and regulatory networks: where metabolic systems biology meets control engineering[END_REF] with focal point φ a ∈ R b , with a, b ∈ V and a = b. For each i ∈ {1, . . . , n} such that a i = b i , there exist trajectories crossing from R a to R c where the coordinates of c are given by c i = a i + sign(b i -a i ) and c k = a k for k = i.

Proof. Assume that i = p (the case i = p is similar).

The assumption a = b implies that, for at least one coordinate i:

ν i (β qi , e) • ν i (β qi , φ a ) < 0 for all e ∈ R a . ( 18
)
Next, equation ( 14) implies w 1 (t) = φ 1 -e 1 (t) → 0 and w i+1 (t) = φ i+1 -e i+1 (t) → 0, as t → ∞, which leads to

(∀ε > 0)(∃t ε > 0) : |ν i (β qi , e(t)) -ν i (β qi , φ a )| < ε, ( 19
)
for all t > t ε . Applying [START_REF] Keasling | Synthetic biology and the development of tools for metabolic engineering[END_REF] yields

|ν i (β qi , e(t))| + |ν i (β qi , φ a )| < ε.
Since ν i (β qi , φ a ) = νi is a constant, for ε = νi condition [START_REF] Khalil | Nonlinear systems[END_REF] implies |ν i (β qi , e(t b ))| = 0 for some finite time instant t b , for which the trajectory exactly reaches the boundary. In conclusion, trajectories may indeed cross the boundary from R a to a neighboring region R c , characterized by ν i (β qi , e) • ν i (β qi , φ a ) > 0 (the opposite of ( 18)).

Next, assume that region R a contains its focal point, φ a ∈ R a . In the more "classical" piecewise linear systems with hyper-rectangular regions [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: equilibria and their stability[END_REF], it is well known that regions that contain their own focal points are forward invariant. However, this is not necessarily true for conic regions (as already discussed in [START_REF] Oyarzún | Multistability and oscillations in genetic control of metabolism[END_REF]), unless all degradation rates are equal. In our case, we have the following result.

Proposition 8 Consider a region R a of system ( 14) containing its own focal point φ a , a ∈ V . Then, R a is forward-invariant.

Proof. Assumption φ a ∈ R a implies:

ν i (β qi , e) • ν i (β qi , φ a ) > 0 for all i and all e ∈ R a . ( 20
)
To check whether trajectories may leave region R a through some coordinate i, consider the vector field along that boundary. Note that, for all e ∈ R a , the enzyme dynamics may be written in terms of the focal point ( 14): ėi = γ(φ a i -e i ). Now define a new variable

r = ν i (β qi , e) ν i (β qi , φ a ) ⇒ ṙ = γ(1 -r),
by reordering the terms d/dt(ν i (β qi , e)) = γ(φ a 1 -e 1 ) -

1 βq i γ(φ a i+1 -e i+1 )
. By [START_REF] Kuntz | Model reduction of genetic-metabolic systems using timescale separation[END_REF], r > 0 for all e ∈ R a . At the boundary ν i (β qi , e) = 0, it holds that r = 0 and ṙ = γ > 0. Since ν i (β qi , φ a ) is constant, we can conclude that r remains positive and never crosses the boundary, that is, solutions e(t) of ( 14) remain in R a .

The set of edges can now be constructed:

E = {a → c : c j = a j + G j , some j and c k = a k , k = j} with G j = sign(b j -a j )
, where R b is the region that contains φ a . The procedure to compute the set E is summarized in Algorithm 1.

The state transition graph

Propositions 7 and 8 determine all possible transitions from a given region R a and establish the basis for an exact algorithmic representation of the global qualitative dynamics of system [START_REF] Glass | The logical analysis of continuous, nonlinear biochemical control networks[END_REF]. The coefficients c j = 0 identify the boundary crossing corresponding to a threshold in variable a j . We note that step 4 in Algorithm 1 reflects an asynchronous updating rule for the evolution of the discrete trajectories, since only one coordinate is allowed to change at each transition.

Algorithm 1. State transition graph for system [START_REF] Glass | The logical analysis of continuous, nonlinear biochemical control networks[END_REF]. 0. Input: set of vertices V . 1. Initialize E.

2. Pick a ∈ V and compute φ a from ( 15) and ( 16). The state transition graph can be analyzed with graph theoretical tools to obtain properties such as the existence of one or more steady states, or oscillatory behavior. In Algorithm 1, a node labeled â represents a region R a which contains its own focal point (step 3). In this case, by Prop. 8, φ a is a regular equilibrium point and a (locally stable) steady state of the system.

Parameter range for the state transition graph

Steps 2 and 3 of Algorithm 1 compute the focal points and their location using [START_REF] Heinrich | The regulation of cellular systems[END_REF] and [START_REF] Gouzé | Positive and negative circuits in dynamical systems[END_REF] which, in principle, require knowledge of the set of parameters. However, the computation of the edges in Step 4 uses only the signs:

δ i,a = δ i,a (κ 0 i , κ 1 i , β qi , γ) = sign(ν i (β qi , φ a )).
So, in fact, the state transition graph computed for a particular set of parameters remains valid for all sets of parameters which yield exactly the same family of signs.

We can formally write

Q = π ∈ R 4(n+1) >0 : δ i,a (κ 0 i , κ 1 i , β qi , γ) ≡ δi,a , ∀ i = 1, . . . , n + 1 and a ∈ V } (21) where π = (κ 0 1 , κ 1 1 , β 1 , . . . , κ 0 n+1 , κ 1 n+1 , β n+1 , γ).
Therefore, given a first set of realistic parameters, the algorithm leads to a full characterization of Q.

5 Case study: a two-regulator pathway

Model definition

To illustrate the application of our results to a metabolic pathway with a complex regulatory architecture, we consider the system shown in Fig. 3 with n = 2 metabo-lites and three enzymes. In this system, both metabolites control enzyme synthesis, and thus we have k = 2 regulators. The regulatory architecture resembles the "metabolator", a metabolic oscillator implemented in the Escherichia coli bacterium by Fung et al [START_REF] Fung | A synthetic genemetabolic oscillator[END_REF]. The metabolator is based on the metabolic reactions that link glycolysis and the respiration cycle, which include two metabolites and three enzymes, and contains positive and negative feedback regulation of two pathway enzymes.

s 0 s 1 s 2 e 1 e 2 e 3
Fig. 3. Example pathway with two regulatory metabolites. We consider the cases in which metabolite s1 activates or represses the production of enzyme e3, denoted by a dotarrow.

The full 5-dimensional model for the system in Fig. 3 reads:

ṡ1 = g 1 (s 0 )e 1 -g 2 (s 1 )e 2 , ṡ2 = g 2 (s 1 )e 2 -g 3 (s 2 )e 3 , ė1 = κ 0 1 + κ 1 1 σ+ (s 2 , θ 1 ) -γe 1 , ė2 = κ 0 2 + κ 1 2 σ-(s 2 , θ 2 ) -γe 2 , ė3 = κ 0 3 + κ 1 3 σ3 (s 1 , θ 3 ) -γe 3 . (22)
In the notation of Section 2, the model in [START_REF] Mannan | Fundamental Design Principles for Transcription-Factor-Based Metabolite Biosensors[END_REF] has the subindices p = 2, ( 1 , 2 , 3 ) = (2, 2, 1), and (q 1 , p) = (3, {1, 2}). We assume that the metabolic reactions follow standard Michaelis-Menten kinetics of the form

g 1 (s 0 ) = k cat 1 s 0 K M 1 + s 0 , g 2 (s 1 ) = k cat 2 s 1 K M 2 + s 1 , g 3 (s 2 ) = k cat 3 s 2 K M 3 + s 2 . (23)
We fix the kinetic parameters to realistic values k cat 1 = 60 min -1 , k cat 2 = 180 min -1 , k cat 3 = 480 min -1 , with a Michaelis-Menten constant K M i = 1 µM for all three enzymatic reactions (i = 1, 2, 3). The degradation rate γ = 0.0116 min -1 corresponds to a doubling time of ∼ 60 min, typical of the bacterium Escherichia coli. We set the expression rates κ 0 i and κ 1 i so that the minimal and maximal enzyme concentrations are E off = [ 1, 4, 1.5 ] and E on = [START_REF] Heinrich | The regulation of cellular systems[END_REF][START_REF] Kuntz | Model reduction of genetic-metabolic systems using timescale separation[END_REF][START_REF] Segel | On the validity of the steady state assumption of enzyme kinetics[END_REF], both in units of nM concentrations. The regulatory functions in [START_REF] Mannan | Fundamental Design Principles for Transcription-Factor-Based Metabolite Biosensors[END_REF] are assumed to be Hill functions of the form

σ+ (s, θ) = s h θ h + s h , σ-(s, θ) = θ h θ h + s h , (24)
for a Hill coefficient h 1 and regulatory thresholds θ 1 = 0.047 µM, θ 2 = 0.055 µM, and θ 3 = 0.22 µM. Note that in the model [START_REF] Mannan | Fundamental Design Principles for Transcription-Factor-Based Metabolite Biosensors[END_REF] we have not specified the type of regulation for σ 3 because we will use our approach for both cases σ 3 = σ + and σ 3 = σ -.

After the separation of timescales, ṡi ≈ 0 for all t ≥ 0 and using step regulatory functions, the model in [START_REF] Mannan | Fundamental Design Principles for Transcription-Factor-Based Metabolite Biosensors[END_REF] reduces to

ė1 = κ 0 1 + κ 1 1 σ + (e 1 /e 3 , 1/β 1 ) -γe 1 , ė2 = κ 0 2 + κ 1 2 σ -(e 1 /e 3 , 1/β 2 ) -γe 2 , ė3 = κ 0 3 + κ 1 3 σ 3 (e 1 /e 2 , 1/β 3 ) -γe 3 , (25)
where the β i parameters are

β 1 = g 1 (s 0 ) g 3 (θ 1 ) , β 2 = g 1 (s 0 ) g 3 (θ 2 ) , β 3 = g 1 (s 0 ) g 2 (θ 3 ) . (26)

State transition graphs

Assuming that the regulatory functions in ( 25) can be well approximated by step functions, the reduced model is equivalent to a 3-dimensional piecewise affine system defined on 3 × 2 1 = 6 disjoint polytopes. Following the notation of Section 4, we associate each polytope to a pair a = (a 1 , a 2 ) with a 1 ∈ {0, 1} and a 2 ∈ {0, 1, 2}. Since s 2 regulates the production of two enzymes, we have that p = 2. Moreover, under the chosen parameter sets, we have that β 2 < β 1 and thus the subindices are p 1 = 2 and p 2 = 1. The regulatory motif vector is σ = [ 1, -1, 1] or σ = [ 1, -1, -1], depending on whether σ 3 = σ + or σ 3 = σ -, respectively. Using equations ( 15)-( 16) we get the focal points for each of the six polytopes in both cases, shown in Fig. 4A.

Fig. 4B shows the six polytopes for a fixed nutrient concentration (s 0 ). To illustrate the impact of model parameters on the graph structure, we used Algorithm 1 to compute the state transition graphs for increasing values of s 0 . The results, shown in Fig. 4C-D, reveal salient differences between the resulting graphs. Some graphs have one or two sink nodes, while others display various cycles. We observe changes in the graphs for a fixed regulatory architecture and increasing various values of s 0 , as well for a fixed s 0 and different regulatory architectures.

One particular case, highlighted in yellow in Fig. 4C-D, reveals that the regulatory sign of σ 3 switches the graph from having two sink nodes to having two cycles.

We investigated these two cases further to examine the bistable behaviour of the model and the possibility of a periodic orbit along the graph cycles. As described in Section 4.4, we can determine the set Q in the parameter space that leads to both these graphs

E off 1 < 1 β 3 E on 2 , E on 1 < 1 β 3 E on 2 , E on 1 > 1 β 3 E off 2 , (27)
E off 1 > 1 β 1 E off 3 , E on 1 > 1 β 2 E off 3 , E on 1 < 1 β 1 E on 3 , (28)
E off 1 < E on 1 (29)
Note that the above conditions link together the enzyme expression parameters (E off i and E on i ) with the enzyme kinetic parameters (k cat i and K M i appearing in the β i parameters).

In the case of repression, i. e. σ 3 = σ -, the interaction graph in Fig. 4C contains only positive loops and therefore the system cannot exhibit oscillatory behavior [START_REF] Thomas | Multistationarity, the basis of cell differentiation and memory. I. structural conditions of multistationarity and other nontrivial behavior[END_REF]. Accordingly, the state transition graph (yellow) has two sinks at nodes 1 and 6, and thus we expect two distinct locally stable steady states in polytopes R 00 and R 12 . As shown in Fig. 5A, simulations of the full model in [START_REF] Mannan | Fundamental Design Principles for Transcription-Factor-Based Metabolite Biosensors[END_REF] with a Hill coefficient h = 4 verify the existence of two stable state states located at the predicted polytopes. Note that because of our piecewise affine approximation, the exact location of the steady states of the continuous model differs slightly from the focal points φ 00 and φ 12 in Fig. 4A.

In the case of positive regulation of enzyme e 3 , i. e. σ 3 = σ + , the interaction graph in Fig. 4D contains three negative loops and thus the system may exhibit periodic oscillations [START_REF] Thomas | Multistationarity, the basis of cell differentiation and memory. I. structural conditions of multistationarity and other nontrivial behavior[END_REF]. Indeed, the state transition graph (Fig. 4D, yellow) has a short cycle through nodes 2-3-6-5, and long cycle through nodes 1-2-3-6-5-4. Upon a close examination of the piecewise affine model along these cycles, we can prove that [START_REF] Chaves | Analysis of a genetic-metabolic oscillator with piecewise linear models[END_REF]: a) the short cycle corresponds to a Filippov-type equilibrium at the intersection of the four polytopes R 01 ∩ R 02 ∩ R 12 ∩ R 11 ; b) the long cycle corresponds to a stable periodic orbit along the polytopes in the cycle. We verified the latter prediction through simulations of the full continuous model in [START_REF] Mannan | Fundamental Design Principles for Transcription-Factor-Based Metabolite Biosensors[END_REF] for increasing values of the Hill coefficient. As observed in Fig. 5B, we found a periodic orbit for h ≥ 10, which qualitatively matches the predicted orbit from the piecewise model for increasing values of the Hill coefficient. From the time course simulations for the metabolites Fig. 4B (bottom), we also observe that the shape and period of the oscillations corresponds well with the metabolic oscillator in [START_REF] Fung | A synthetic genemetabolic oscillator[END_REF] (see, for instance, Figs. 1,2 therein), thus suggesting the periodic orbit is within physiological ranges and potentially observable in experiments.

Conclusions

Cellular metabolism is controlled by networks of transcriptional regulatory loops, often masked by complex feedback architectures between metabolic intermediates and enzymatic genes. Determining the closed-loop dynamics of a specific architecture is challenging because of the strong nonlinearities and the large number of chemical species involved. In this paper we have proposed a formalism for the qualitative analysis of closed-loop dynamics in a large class of regulation systems for unbranched pathways. Our analysis relies on time scale separation and piecewise affine models to accomplish: (i) a complete classification of all possible feedback interaction networks associated to a metabolic pathway; and (ii) an algorithmic construction of the state transition graph that describes the global system dynamics. In a synthetic biology and metabolic engineering context, this feedback network catalogue provides an indication of the components (metabolites, enzymes, promoters) required to assemble a metabolic pathway exhibiting a desired circuitry and dynamics [START_REF] Mannan | Fundamental Design Principles for Transcription-Factor-Based Metabolite Biosensors[END_REF].

Our analytical results are obtained for a new class of piecewise linear systems in conic regions, derived by approximating Hill functions by step functions. Among other benefits, this framework allows establishing existence of an asymptotically stable periodic orbit for a system whose regulatory motif is reminiscent of the metabolator [START_REF] Fung | A synthetic genemetabolic oscillator[END_REF]. In fact, theoretical proofs of oscillatory behavior are difficult to obtain in general, and the result for the metabolator remains unclear [START_REF] Reznik | The dynamics of hybrid metabolic-genetic oscillators[END_REF].

Despite their usefulness, piecewise linear systems provide a possibly less realistic description then their continuous counterparts. Large Hill coefficients are typically associated with strong cooperativity effects [START_REF] Heinrich | The regulation of cellular systems[END_REF], while gene expression profiles can often be described by steep sigmoidal curves [START_REF] Wang | Engineering modular and orthogonal genetic logic gates for robust digitallike synthetic biology[END_REF]. To verify the suitability of the piecewise linear approximation in our case, we performed simulations for Hill coefficients in a wide interval. The differences observed in the dynamical behavior reside essentially in the form and amplitude of the time solutions, both in the case of multiple equilibria and in the presence of periodic behavior (see Fig. 5).

One of our goals in this work was to devise an analysis technique applicable to a wide range of regulatory architectures. As a result, for tractability we made a number of simplifying assumptions. In particular, we have not considered regulatory mechanisms such as product inhibition, allostery and other post-translational effects that are common in metabolic pathways. A key step in our method is that the solution of quasi steady state equation in (5) depends only on enzyme concentrations. In such case, the complete system (1)-( 2) can be projected onto the state space of enzymes, partitioned into suitably defined polytopes. We note that this does require the explicit solution of the quasi steady state equation: it is sufficient to guarantee the existence of a unique solution for the regulatory metabolites in terms of enzyme concentrations. Specific types of post-translational mechanisms may lead to quasi steady state equations that are solvable and respect these conditions. In such cases our approach could be extended with minor modifications, taking into account the geometry of the resulting polytopic partition. A similar reasoning applies to other potential extensions of our method, including the analysis of branched pathways.

In conclusion, our results draw a link between various architectures of transcriptional regulation and classes of prototypical networks, with all the positive and negative feedback loops explicitly identified. This allowed us to single out candidate global dynamics using graph analyses, and thus uncover the enormous dynamical diver-sity that can emerge in metabolic pathways under gene regulation. 22). (A) For the case of negative regulation of e3 (Fig. 4C), we observe bistable dynamics as predicted by the state transition graph highlighted in Fig. 4C. Plots show the state trajectories for 500 randomized initial conditions and Hill coefficient h = 4. The enzyme initial conditions are coloured according to the polytopic domain to which they belong, defined in Fig. 4B. (B) For the case of positive regulation of e3 (Fig. 4D), we found a stable periodic orbit as predicted by the state transition graph highlighted in Fig. 4D and our theoretical analysis. Top plots show the limit cycle for high and increasing values of the Hill coefficient, starting from h = 10. For large values of h, the orbit approaches the one predicted by the piecewise affine model under timescale separation; we have coloured the limiting orbit according to the polytopic domains defined in Fig. 4B). Bottom plots shows the time courses for h = 10 and h = 20.

Fig. 1 .

 1 Fig. 1. Regulatory architectures in unbranched metabolic pathways. We consider a general class of nonlinear feedback systems with up to n regulatory metabolites (si) with any combination of positive of negative regulation on the enzymatic genes (ei). The diagram illustrates assumptions A1, A2, and A2'. Regular arrows (resp., bar-headed arrows) represent activation (resp., repression). Dot-headed arrows represent either of these.

Fig. 2 .

 2 Fig. 2. Interaction graphs of the transformed system in (13). (A) All possible interaction graphs (up to re-ordering) for pathways of length n = 3. Left to right: w = z1; y = z1; w = z1 and y = z2. (B) All possible interaction graphs (up to re-ordering) for pathways of length n = 4. Left to right and top to bottom: w, y, z1, z2 all distinct; w = z1; y = z1; w = z1 and y = z2; alternative to w = z1; alternative to y = z1.

3 .

 3 Compute b such that φ a ∈ R b from (17). If b ≡ a label node a by adding a symbol, â. 4. For j = 1 until n do: 4.1. Compute cj = aj + sign(bj -aj), 4.2. If cj = aj add an edge a → c to E, with c k = a k for k = j. 5. Repeat steps 2 to 4 for each vertex in V . 6. Output: set of edges E.

Fig. 4 .

 4 Fig. 4. State transition graphs for the example in Fig. 3 and equation (22). (A) Polytopic domains and their focal points as defined in equations (15) and (16). (B) Positive orthant partitioned into the six polytopic domains in Definition 6. (C-D) Pathway diagram and state transition graphs for the system with negative (σ = [1, -1, -1]) (C) and positive (σ = [1, -1, 1]) (D) regulation of enzyme e3. The interaction graphs are given for the transformed variables defined in Fig. 2a(left), i.e. w = z1 = e1 -e3/β1, y = e1 -e2/β3, and z2 = e1 -e3/β2. The state transition graphs were computed with Algorithm 1 for increasing nutrient concentrations s0 = {0.5, 1.5, 4, 10} µM; arrows denote sink nodes or graph cycles, respectively. The two highlighted graphs are for s0 = 1.5, µM and suggest the bistable and oscillatory dynamics.

Fig. 5 .

 5 Fig.5. Simulations of the example system in Fig.3and equation[START_REF] Mannan | Fundamental Design Principles for Transcription-Factor-Based Metabolite Biosensors[END_REF]. (A) For the case of negative regulation of e3 (Fig.4C), we observe bistable dynamics as predicted by the state transition graph highlighted in Fig.4C. Plots show the state trajectories for 500 randomized initial conditions and Hill coefficient h = 4. The enzyme initial conditions are coloured according to the polytopic domain to which they belong, defined in Fig.4B. (B) For the case of positive regulation of e3 (Fig.4D), we found a stable periodic orbit as predicted by the state transition graph highlighted in Fig.4Dand our theoretical analysis. Top plots show the limit cycle for high and increasing values of the Hill coefficient, starting from h = 10. For large values of h, the orbit approaches the one predicted by the piecewise affine model under timescale separation; we have coloured the limiting orbit according to the polytopic domains defined in Fig.4B). Bottom plots shows the time courses for h = 10 and h = 20.

Acknowledgements M.C. was partly supported by the French agency for research through project ICycle ANR-16-CE33-0016-01. D.O. was partly supported by the Human Frontier Science Program through a Young Investigator Grant (RGY0076-2015).