Dynamics of complex feedback architectures in metabolic pathways

Automatica 99 (2019) 323–332

Contents lists available at ScienceDirect

Automatica journal homepage: www.elsevier.com/locate/automatica

Dynamics of complex feedback architectures in metabolic pathways✩ Madalena Chaves a, *, Diego A. Oyarzún b,c a b c

Université Côte d’Azur, Inria, INRA, CNRS, Sorbonne Université, Biocore team, Sophia Antipolis, France School of Informatics, University of Edinburgh, Edinburgh EH8 9AB, United Kingdom School of Biological Sciences, University of Edinburgh, Edinburgh EH9 3JH, United Kingdom

article

info

Article history: Received 12 December 2017 Received in revised form 6 July 2018 Accepted 5 October 2018

Keywords: Metabolic pathways Synthetic biology Stability analysis Piecewise linear models Gene regulatory networks

a b s t r a c t Cellular metabolism contains intricate arrays of feedback control loops. These are a key mechanism by which cells adapt and 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. © 2018 Elsevier Ltd. All rights reserved.

1. Introduction Cellular metabolism adapts to environmental fluctuations through networks of regulatory feedback loops (Chubukov, Gerosa, Kochanowski, & Sauer, 2014). 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 (Farmer & Liao, 2000; Keasling, 2012; Xu, Li, Zhang, Stephanopoulos, & Koffas, 2014) and the generation of entirely new metabolic dynamics such as oscillations (Fung et al., 2005) or cellto-cell communication (Gupta, Reizman, Reisch, & Prather, 2017; Soma & Hanai, 2015). ✩ 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). The material in this paper was not presented at any conference. This paper was recommended for publication in revised form by Associate Editor Daniel Liberzon under the direction of Editor Nathan van de Wouw. Corresponding author. E-mail addresses: [email protected] (M. Chaves), [email protected] (D.A. Oyarzún).

*

https://doi.org/10.1016/j.automatica.2018.10.046 0005-1098/© 2018 Elsevier Ltd. All rights reserved.

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 (Oyarzún & Chaves, 2015; Stevens & Carothers, 2015). 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 (Brockman & Prather, 2015; He, Murabito, & Westerhoff, 2016; Liu, Mannan, Han, Oyarzún, & Zhang, 2018), while the use of control theoretical principles has gained substantial traction in the field (Del Vecchio, Dy, & Qian, 2016; Oyarzún & Stan, 2013). 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

324

M. Chaves, D.A. Oyarzún / Automatica 99 (2019) 323–332

corresponding dynamical behavior. 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 Oyarzún, Chaves, and Hoffmeyer-Zlotnik (2012) 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 Fung et al. (2005) that illustrates the dynamical diversity within each class of feedback architectures (Section 5). 2. 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 si into si+1 is catalyzed by an enzyme ei+1 . The substrate s0 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, s0 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.

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.

(De Jong, 2002; Mannan, Liu, Zhang, & Oyarzún, 2017; Yagil & Yagil, 1971). The positive parameters (κi0 , κi1 ) 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 Eqs. (1)–(2). To this end, we will use two techniques for simplifying the models (Oyarzún et al., 2012): 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. 2.1. Piecewise linear models for gene regulation

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 (Oyarzún et al., 2012) that considered the case of only one regulatory metabolite (k = 1). In particular, Assumption A2 allows for any number of regulatory metabolites to control enzyme expression. In most natural instances, e.g. the lactose operon or amino acid metabolism (Chubukov et al., 2014), 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 (Liu et al., 2018). If si and ei denote the concentration of metabolites and enzymes respectively, the system in Fig. 1 can be described by the nonlinear model:

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 (Thomas, 1973) and Glass and Kauffman (Glass & Kauffman, 1973):

s˙i = gi (si−1 )ei − gi+1 (si )ei+1 , i = 1, . . . , n,

(1)

Eioff =

e˙ i = κi0 + κi1 σi (sℓi , θi ) − γ ei , i = 1, . . . , n + 1,

(2)

where sℓi is the metabolite that regulates ei , gi : R+ ≥0 → [0, Mi ] are continuous and strictly increasing functions representing the kinetics of each enzyme (Heinrich & Schuster, 1996), and σi ∈ {σ˜ + , σ˜ − } are increasing or decreasing sigmoidal functions describing the nonlinear regulation from metabolites to enzyme synthesis

σ+ (r , θ ) = 1 − σ− (r , θ ) =

1, 0,

{

r > θ, 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, Eioff ≤ ei (t) ≤ Eion for all t ≥ t0 > 0, with (κ 0 + κi1 ) κi0 , Eion = 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 (Oyarzún et al., 2012; Ropers, Baldazzi, & de Jong, 2011) and genetic transcription and translation networks (Edwards, Machina, McGregor, & van den Driessche,

M. Chaves, D.A. Oyarzún / Automatica 99 (2019) 323–332

2015; Hudson & Edwards, 2016). 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 (Batt, Belta, & Weiss, 2008; Casey, de Jong, & Gouzé, 2006); this graph provides a global view of the dynamics of the system and will be constructed in Section 4. 2.2. Timescale separation

εi = max

γ KM , i Eion kcat i

,

γ KM ,i+1 Eion kcat i

}

≪ 1,

(4)

for i = 2, . . . , n. These conditions can be obtained by applying an appropriate change of variables: e˜ i = ei /Eion , s˜i = si /KM ,i , and t˜ = γ t, similar to those used in Kuntz et al. (2014) and Segel (1988). The conditions typically hold for metabolic systems, since dilution by growth is much slower than enzyme turnover rates, and hence γ ≪ kcat i . The quasi-steady state assumption is: A3. We assume that s˙i ≈ 0 for each t ≥ 0, leading to g1 (s0 )e1 = gi (si−1 )ei , for i = 2, . . . , n + 1.

(5)

Using assumption A3 we can substitute sℓi =

1 gℓ−+ i 1

( g1 (s0 )

e1

)

eℓi +1

(6)

into Eq. (2) for enzyme ei . To guarantee the existence of solutions of (6), a sufficient condition is: Eioff ≥

g1 (s0 ) Mℓi +1

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 sp ) 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}

Metabolic timescales are much faster than those of enzyme synthesis (Heinrich & Schuster, 1996; Ropers et al., 2011), so the system (1)–(2) can be simplified by assuming the metabolites to be in quasi-steady state. The variables si 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 si in the ODEs for the enzymes, the slow dynamics in (2) can be simplified to depend only on the enzyme concentrations ei . Such timescale separation ideas have been used for analysis of various enzymatic mechanisms (Baldazzi, Ropers, Geiselmann, Kahn, & de Jong, 2012; Kuntz, Oyarzún, & Stan, 2014), genetic transcription and translation networks (Edwards et al., 2015), and justify the well known Michaelis–Menten kinetic mechanism (Segel, 1988). Briefly, the goal is to rewrite (1)–(2) in the form: ε˙s = G(s, e, ε), e˙ = F (s, e, ε), where ε ≪ 1. Under suitable conditions (see Khalil, 2002 for more details), G(s, e, 0) = 0 defines ˜ the slow manifold s = G(e), and Tikhonov’s Theorem guarantees ˜ , e, 0) remain that solutions of the slower dynamics e˙ = F (G(e) close to those of the original system. In the particular case where all enzymes follow Michaelis– Menten kinetics, that is, gi (si ) = kcat i si−1 /(KM ,i + si−1 ) (see example in Section 5), the timescale separation can be applied to (1)–(2) provided that

{

325

Eion , for i = 2, . . . , n + 1.

(7)

Now, let Ireg = {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 Inrg = {1, 2, . . . , n}\ Ireg . From (6), it becomes clear that the subset of k + 1 equations corresponding to {ei : i ∈ Ireg }, does not depend on the variables {ej , i ∈ Inrg }, so it can be studied separately. Therefore, without loss of generality we will thus consider a modified version of Assumption A2:

(8)

with I (i) = ℓi and I (p1 ) = I (p2 ) = p. Its inverse is defined as: I −1 (j) = qj and I −1 (p) = {p1 , p2 }. 2.3. 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, −1,

σi (r , θ ) = σ+ (r , θ ) σi (r , θ ) = σ− (r , θ ).

{

(9)

In addition, the quasi steady state assumption (6) and the fact that gi are non-decreasing imply sℓi > θi ⇔ g1 (s0 )e1 > gℓi +1 (θi )eℓi +1 ⇔ e1 >

1

βi

eℓi +1 .

(10)

where we have defined the parameter βi = g1 (s0 )/gℓi +1 (θi ). Under the relations in (10), the step functions satisfy

σi (sℓi , θi ) = σi (sℓi − θi , 0) = σi (e1 −

1

βi

eℓi +1 , 0),

(11)

and therefore the enzyme dynamics in (2) can be rewritten as e˙ i = κi0 + κi1 σi

( e1 −

1

βi

)

eℓi +1 , 0

− γ ei .

(12)

Eq. (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. 3. General structure of the interaction graph Given a metabolic pathway such as the one in Fig. 1, one would like to 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 (Gouzé, 1998; Remy, Ruet, & Thieffry, 2008; Richard & Comet, 2007; Soulé, 2003). An interaction graph describes the dependencies between model variables according to the sign pattern of the Jacobian matrix De Jong (2002). Although the Jacobian of (12) has a nonconstant 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 e1 − 1 e , and thus we introduce the following definition: β ℓi +1 i

Definition 1 (Cone Membership Function). Let β > 0, r ∈ Rn>+01 , and ri the ith component of r. Define the cone membership function by

νi (β, r) : R>0 × Rn>+01 −→ R (β, r) −→ r1 −

1

β

ri+1 .

326

M. Chaves, D.A. Oyarzún / Automatica 99 (2019) 323–332

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:

xj = νj−1 (βqj−1 , e), j ∈ {2, . . . , n + 1} \ {p + 1},

w ; w ⇆ y; w → z2 → y → w,

z1 = νp (βp1 , e), z2 = νp (βp2 , e),

and for the first case in Fig. 2B

From the model in (12), it follows that:

w ; w ⇆ z1 ; w → y → z1 → w; w → z2 → y → z1 → w.

x˙ j = κ˜ 10 + κ11 σ1 (xℓ1 +1 , 0) −

κj1 βqj−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 halfline with slope β . The new variables can now be defined as:

In fact, this result can be generalized as follows:

σj (xℓj +1 , 0) − γ xj ,

Proposition 3 (Nested Loops). Consider the n + 1 dimensional system (13). If the interaction graph is formed of a single strongly connected component, and if either w ̸ = zj or y ̸ = zj , for all j ∈ {1, 2}, +1 then it contains n + 1 nested loops, {Lℓ }nℓ= 1 , defined by:

j ∈ {2, . . . , n + 1} \ {p + 1} z˙j = κ˜ p0j + κ11 σ1 (xℓ1 +1 , 0) −

κp1+1 βpj

σp+1 (xℓp+1 +1 , 0) − γ zj ,

j = 1, 2.

(13)

where κ˜ j0 = κ10 − κj0 /βqj−1 and κ˜ p0j = κ10 − κp0+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 xj to xi if ∂ (x˙ i )/∂ xj ̸ = 0, and (for simplicity) omitting all arrows related to the natural degradation terms −γ xj or −γ zj . 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 z1 and z2 ; R4 all other nodes xℓj +1 have exactly one outgoing arrow, signed −σj , towards xj . In particular, node zj has an outgoing arrow, signed −σpj , towards xpj . A self arrow counts as both an incoming and an outgoing arrow for the node. If a node has an outgoing arrow to xp1 , then it has necessarily an outgoing arrow to xp2 , 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 Eqs. (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 z1 and z2 (which can coincide with either w or y, when pj = ℓ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 xa , xb , etc. To obtain all the distinct interaction graphs, consider all combinations of z1 and z2 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, s1 and s2 , so it follows from (6) that e4 will not influence the dynamics of ei , i = 1, 2, 3. Since s2 regulates two enzymes, it follows that p = 2 and, from the definition in (8), the index correspondence is (ℓ1 , ℓ2 , ℓ3 ) = (1, 2, 2), (q1 , p) = (1, {2, 3}). The change of variables is as follows: 1

1

1

e2 , z1 = e1 − e3 , z2 = e1 − e3 , β1 β3 β2 with w = x2 and y = z1 , yielding an interaction graph as in

x2 = e1 −

Fig. 2A(middle).

L1 = {w}, L2 = {w, v1 : w → v1 , v1 → w}, Lℓ+1 = {w, v1 , . . . , vℓ : w → vℓ , vj+1 → vj , v1 → w}, where j = 1, . . . , ℓ − 1 in Lℓ+1 . If the interaction graph contains more than one strongly connected component, then one of these components has j + 1 nodes and contains nested loops L1 , . . . , Lj+1 . Proof. This structure follows directly from the properties (R1)– (R4) and (i)–(iv). First, note that w has an arrow towards every node, including itself, so w forms a trivial first loop, L1 . Second, choose the node v1 ̸ = 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 v1 . This gives a 2-node circuit L2 . Now, to use an induction argument, suppose we have chosen j ≥ 1 distinct nodes, v1 , . . . , vj , defining loops L1 up to Lj+1 . Note that, by construction, each of the nodes vi , i < j has exactly two incoming arrows, w and vi+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 vj+1 ̸ = w which has an edge towards vj : again, by point (R1) there can be only one such node. If vj+1 ̸ = vj , then we have found loop Lj+2 and can continue. Otherwise, vj+1 = vj , which means that vj = y since it has two outgoing arrows. But it also has two incoming arrows, w and vj . This means that the set of nodes w, v1 , . . . , vj (with j < n) cannot have any further incoming arrows and thus forms a strongly connected component, which contains the nested loops L1 , . . . , Lj+1 . Finally, the remaining nodes (vj+1 , . . . , vn ) 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 = z1 and y = z2 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 (13) 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 L1 , . . . , Lj+1 . Indeed, the dynamics of this sub-system w, v1 , . . . , vj can be studied as an autonomous system, whose behavior is an input to the system formed by the remaining variables vj+1 , . . . , vn . Examples are shown in Fig. 2B (bottom middle and right) which can both be decomposed into a 3-dimensional system (w, y, z2 ) which controls xa . 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.

M. Chaves, D.A. Oyarzún / Automatica 99 (2019) 323–332

327

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 .

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 xi nodes: (a) two distinct graphs, if w = z1 or w = z2 ; (b) two distinct graphs, if y = z1 or y = z2 ; (c) five distinct graphs if n ≥ 3 and all w, y, z1 , z2 are different. Proof. (a) In this case, the first two loops are L1 = {w} and L2 = {w, y} (i.e., v1 = y, in the notation of Proposition 3). Choosing the sequence of nodes vi as in Proposition 3, recall that vj = z2 implies vj+1 = y, by definition of y. To guarantee that the graph has a single strongly connected component, it must be that vn = z2 . Otherwise, if vj = z2 for j < n, then vj+1 = y , which implies that Lj+1 ≡ Lj+1 = {w, v1 , . . . , vj } is a strongly connected component and so the full graph cannot be strongly connected. (b) In this case, assume y = z1 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 z2 . To guarantee that the graph has a single strongly connected component, it must be that vn = y and vn−1 = z2 . Otherwise, by the same argument as in case (a), if vj = y for some j < n, the set of nodes {w, v1 , . . . , vj } 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 zi or from one of the n − 3 remaining variables xi , up to re-ordering. Suppose first that z1 → w (or z2 → w ). Then, as in case (b), we need vn = y and vn−1 = z2 to guarantee a strongly connected interaction graph (count 2 distinct graphs). Suppose next that xa → w and z1 → y (or z2 → y). Then, we again need vn = y and vn−1 = z2 (count 2 distinct graphs). Suppose finally that xa → w and xb → 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 (13) give an indication of the possible dynamics of the feedback system (for instance, multistability is expected whenever all loops are positive (Thomas & Kaufman, 2001)). 4. 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 Eq. (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 Eq. (12) can be written as the piecewise linear system (see also (Casey et al., 2006))

e˙ i =

⎧ σ ⎪ ⎨γ (Ei i − ei ),

e1 >

⎪ ⎩γ (Ei−σi − ei ),

e1 <

1

βi

eℓi +1

βi

eℓi +1

1

,

(14)

where σ = (σ1 , . . . , σn ) represents the regulatory motif, and (Ei−1 , Ei1 ) = (Eioff , Eion ) are the minimal and maximal enzyme concentrations in Eq. (3). In this new description, the state space is partitioned into several conic regions defined by the two planes zj = 0 and the n − 1 planes xj = 0 (see (13)). Each region has a σ focal point φ = (φ1 , . . . , φn )′ with φi = Ei i if νℓi (βi , e) > 0 and −σi φi = Ei 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.

328

M. Chaves, D.A. Oyarzún / Automatica 99 (2019) 323–332

4.1. Vertices of the transition graph As shown in the next definition, each vertex of the graph represents one region Ra , 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 (8). A region Ra , with a = (a1 , . . . , an )′ with ap ∈ {0, 1, 2} and ai ∈ {0, 1} for i ∈ {1, . . . , n} \ {p}, is defined by: Ra = { e ∈ Rn≥+01 : ai = σ+ (νi (βqi , e), 0), ap = σ+ (νp (βp1 , e), 0) + σ+ (νp (βp2 , e), 0) .

}

4.2.3. 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 ̸ ∈ Ra . Then it can be shown that there is at least one boundary hyperplane ei+1 = βqi e1 through which trajectories can escape Ra . Proposition 7. Consider a region Ra of system (14) with focal point

φ a ∈ Rb , with a, b ∈ V and a ̸= b. For each i ∈ {1, . . . , n} such that ai ̸ = bi , there exist trajectories crossing from Ra to Rc where the coordinates of c are given by ci = ai + sign(bi − ai ) and ck = ak 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:

The set of vertices of the graph is thus given by V = {(a1 , . . . , an ) : ap ∈ {0, 1, 2}, ai ∈ {0, 1} for i ̸ = p}.

νi (βqi , e) · νi (βqi , φ a ) < 0 for all e ∈ Ra .

4.2. Edges of the transition graph

Next, Eq. (14) implies w1 (t) = φ1 − e1 (t) → 0 and wi+1 (t) = φi+1 − ei+1 (t) → 0, as t → ∞, which leads to

Each edge represents a crossing between two adjacent conic regions: if there is a trajectory starting in some region Ra and going through a neighboring region Rb , 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 Rb with a ̸ = b, so trajectories starting in Ra will eventually cross to a neighboring 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 (aj ̸ = bj and ai = bi for all i ̸ = j). To determine the edges, the next steps are the calculation of the focal points, φ a = (φ1a , . . . , φna ), for each region Ra and their location Rb , i.e. φ a ∈ Rb , with a = (a1 , . . . , an ) and b = (b1 , . . . , bn ). 4.2.1. Focal points The coordinates of the focal point φ a for region Ra can be explicitly calculated recalling (14) and using Definition 6: σq

−σqi

φqai = ai Eqi i + (1 − ai )Eqi

,

if i ̸ = p,

(15)

and in the case i = p: σp

−σp1

φpa1 = max(ap − 1, 0)Ep1 1 + (1 − max(ap − 1, 0))Ep1 φ

a p2

σp2

−σp2

= min(ap , 1)Ep2 + (1 − min(ap , 1))Ep2

.

,

(18)

(∀ε > 0)(∃tε > 0) : ⏐νi (βqi , e(t)) − νi (βqi , φ a )⏐ < ε,

⏐

⏐

(19)

for all t > tε . Applying (18) yields

⏐ ⏐ ⏐ ⏐ ⏐νi (βq , e(t))⏐ + ⏐νi (βq , φ a )⏐ < ε. i i Since νi (βqi , φ⏐ a ) = νˆ i is a constant, for ε = νˆ i condition (19) implies ⏐ ⏐νi (βq , e(tb ))⏐ = 0 for some finite time instant tb , for which the i trajectory exactly reaches the boundary. In conclusion, trajectories may indeed cross the boundary from Ra to a neighboring region Rc , characterized by νi (βqi , e) · νi (βqi , φ a ) > 0 (the opposite of (18)). Next, assume that region Ra contains its focal point, φ a ∈ Ra . In the more ‘‘classical’’ piecewise linear systems with hyperrectangular regions (Casey et al., 2006), 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 Oyarzún et al., 2012), unless all degradation rates are equal. In our case, we have the following result. Proposition 8. Consider a region Ra of system (14) containing its own focal point φ a , a ∈ V . Then, Ra is forward-invariant. Proof. Assumption φ a ∈ Ra implies:

(16)

νi (βqi , e) · νi (βqi , φ a ) > 0 for all i and all e ∈ Ra .

(20) a

4.2.2. Location of the focal points Computation of the region Rb is straightforward, as it suffices to check the sign of the functions νi when computed at the (already known) focal points:

To check whether trajectories may leave region R through some coordinate i, consider the vector field along that boundary. Note that, for all e ∈ Ra , the enzyme dynamics may be written in terms of the focal point (14): e˙ i = γ (φia − ei ). Now define a new variable

bi = σ+ (νi (βqi , φ a ), 0)

r =

bp = σ+ (νp (βp1 , φ a ), 0) + σ+ (νp (βp2 , φ a ), 0).

(17)

Observe that the focal points and their location in the state space 0,1 depend on the parameters of the system (κ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 , φ ) ̸ = 0 for all i and all a ∈ V . a

νi (βqi , e) ⇒ r˙ = γ (1 − r), νi (βqi , φ a )

by reordering the terms d/dt(νi (βqi , e)) = γ (φ1a − e1 ) − β1 γ (φia+1 − qi

ei+1 ). By (20), r > 0 for all e ∈ Ra . At the boundary νi (βqi , e) = 0, it holds that r = 0 and r˙ = γ > 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 Ra . The set of edges can now be constructed: E = {a → c : cj = aj + Gj , for some j and ck = ak , k ̸ = j}

with Gj = sign(bj − aj ), where Rb is the region that contains φ a . The procedure to compute the set E is summarized in Algorithm 1.

M. Chaves, D.A. Oyarzún / Automatica 99 (2019) 323–332

329

4.3. The state transition graph Propositions 7 and 8 determine all possible transitions from a given region Ra and establish the basis for an exact algorithmic representation of the global qualitative dynamics of system (12). The coefficients cj ̸ = 0 identify the boundary crossing corresponding to a threshold in variable aj . 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 (12).

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

0. Input: set of vertices V . 1. Initialize E . 2. Pick a ∈ V and compute φ a from (15) and (16). 3. Compute b such that φ a ∈ Rb from (17). If b ≡ a label node a by adding a symbol, aˆ . 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 ck = ak for k ̸ = j. 5. Repeat steps 2 to 4 for each vertex in V . 6. Output: set of edges E .

s˙1 = g1 (s0 )e1 − g2 (s1 )e2 , s˙2 = g2 (s1 )e2 − g3 (s2 )e3 , e˙ 1 = κ10 + κ11 σ˜ + (s2 , θ1 ) − γ e1 ,

e˙ 3 = κ30 + κ31 σ˜ 3 (s1 , θ3 ) − γ e3 .

Steps 2 and 3 of Algorithm 1 compute the focal points and their location using (15) and (17) 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 (κi0 , κi1 , β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 3(n+1)+1

: δi,a (κi0 , κi1 , βqi , γ ) ≡ δ¯i,a , ∀i = 1, . . . , n + 1 and a ∈ V

In the notation of Section 2, the model in (22) has the subindices p = 2, (ℓ1 , ℓ2 , ℓ3 ) = (2, 2, 1), and (q1 , p) = (3, {1, 2}). We assume that the metabolic reactions follow standard Michaelis–Menten kinetics of the form kcat 1 s0 g1 (s0 ) = , KM 1 + s 0 g2 (s1 ) =

4.4. Parameter range for the state transition graph

{

(22)

e˙ 2 = κ20 + κ21 σ˜ − (s2 , θ2 ) − γ e2 ,

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 aˆ represents a region Ra which contains its own focal point (step 3). In this case, by Proposition 8, φ a is a regular equilibrium point and a (locally stable) steady state of the system.

Q = π ∈ R>0

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 dot-headed arrow.

}

(21)

where π = (κ10 , κ11 , β1 , . . . , κn0+1 , κn1+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 5.1. 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 metabolites 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. (2005). The metabolator is based on the metabolic reactions that link glycolysis and the respiration cycle, which include two metabolites and three enzymes, and contain positive and negative feedback regulation of two pathway enzymes.

g3 (s2 ) =

kcat 2 s1 KM 2 + s 1 kcat 3 s2 KM 3 + s 2

,

(23)

.

We fix the kinetic parameters to realistic values kcat 1 = 60 min−1 , kcat 2 = 180 min−1 , kcat 3 = 480 min−1 , with a Michaelis–Menten constant KM 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 κi0 and κi1 so that the minimal and maximal enzyme concentrations are E off = [ 1, 4, 1.5 ] and E on = [ 15, 20, 30 ], both in units of nM concentrations. The regulatory functions in (22) are assumed to be Hill functions of the form

σ˜ + (s, θ ) =

sh

θh + s

, h

σ˜ − (s, θ ) =

θh , θ h + sh

(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 (22) 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, s˙i ≈ 0 for all t ≥ 0 and using step regulatory functions, the model in (22) reduces to e˙ 1 = κ10 + κ11 σ+ (e1 /e3 , 1/β1 ) − γ e1 , e˙ 2 = κ20 + κ21 σ− (e1 /e3 , 1/β2 ) − γ e2 ,

(25)

e˙ 3 = κ30 + κ31 σ3 (e1 /e2 , 1/β3 ) − γ e3 , where the βi parameters are

β1 =

g1 (s0 ) g 3 (θ 1 )

,

β2 =

g1 (s0 ) g3 (θ2 )

,

β3 =

g1 (s0 ) g2 (θ3 )

.

(26)

5.2. State transition graphs Assuming that the regulatory functions in (25) can be well approximated by step functions, the reduced model is equivalent

330

M. Chaves, D.A. Oyarzún / Automatica 99 (2019) 323–332

to a 3-dimensional piecewise affine system defined on 3 × 21 = 6 disjoint polytopes. Following the notation of Section 4, we associate each polytope to a pair a = (a1 , a2 ) with a1 ∈ {0, 1} and a2 ∈ {0, 1, 2}. Since s2 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 p1 = 2 and p2 = 1. The regulatory motif vector is σ = [ 1, −1, 1] or σ = [ 1, −1, −1], depending on whether σ3 = σ+ or σ3 = σ− , respectively. Using Eqs. (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 (s0 ). 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 s0 . 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 s0 , as well for a fixed s0 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 behavior 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 E1off < E1off >

1

β3 1

β1

E2on ,

E1on <

E3off ,

E1on >

E1off < E1on

1

β3 1

β2

E2on ,

E1on >

E3off ,

E1on <

1

β3 1

β1

E2off ,

(27)

E3on ,

(28) (29)

Note that the above conditions link together the enzyme expression parameters (Eioff and Eion ) with the enzyme kinetic parameters (kcat i and KM 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 (Thomas & Kaufman, 2001). 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 R00 and R12 . As shown in Fig. 5A, simulations of the full model in (22) with a Hill coefficient h = 4 verify the existence of two stable 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 e3 , i. e. σ3 = σ+ , the interaction graph in Fig. 4D contains three negative loops and thus the system may exhibit periodic oscillations (Thomas & Kaufman, 2001). 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 (Chaves, Oyarzún, and Gouzé, 2019): (a) the short cycle corresponds to a Filippov-type equilibrium at the intersection of the four polytopes R01 ∩R02 ∩R12 ∩ R11 ; (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 (22) 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 Fung et al. (2005) (see, for instance, Figs. 1, 2 therein), thus suggesting the periodic orbit is within physiological ranges and potentially observable in experiments.

6. 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 (Mannan et al., 2017). 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 (Fung et al., 2005). In fact, theoretical proofs of oscillatory behavior are difficult to obtain in general, and the result for the metabolator remains unclear (Reznik, Kaper, & Segrè, 2013). 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 (Heinrich & Schuster, 1996), while gene expression profiles can often be described by steep sigmoidal curves (Wang, Kitney, Joly, & Buck, 2011). 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 diversity that can emerge in metabolic pathways under gene regulation.

M. Chaves, D.A. Oyarzún / Automatica 99 (2019) 323–332

331

Fig. 4. State transition graphs for the example in Fig. 3 and Eq. (22). (A) Polytopic domains and their focal points as defined in Eqs. (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 existence of bistable and oscillatory dynamics.

Fig. 5. Simulations of the example system in Fig. 3 and Eq. (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 colored 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 colored the limiting orbit according to the polytopic domains defined in (Fig. 4B). Bottom plots shows the time courses for h = 10 and h = 20.

332

M. Chaves, D.A. Oyarzún / Automatica 99 (2019) 323–332

References Baldazzi, V., Ropers, D., Geiselmann, J., Kahn, D., & de Jong, H. (2012). Importance of metabolic coupling for the dynamics of gene expression following a diauxic shift in Escherichia coli. Journal of Theoretical Biology, 295, 100–115. Batt, G., Belta, C., & Weiss, R. (2008). Temporal logic analysis of gene networks under parameter uncertainty. IEEE Transactions on Automatic Control, 53, 215–229. Brockman, Irene M., & Prather, Kristala L. J. (2015). Dynamic metabolic engineering: New strategies for developing responsive cell factories. Biotechnology Journal, 10(9), 1360–1369. Casey, R., de Jong, H., & Gouzé, J.-L. (2006). Piecewise-linear models of genetic regulatory networks: equilibria and their stability. Journal of Mathematical Biology, 52, 27–56. Chaves, M., Oyarzún, D. A., & Gouzé, Jean-Luc (2019). Analysis of a genetic-metabolic oscillator with piecewise linear models. Journal of Theoretical Biology (in press). Chubukov, Victor, Gerosa, Luca, Kochanowski, Karl, & Sauer, Uwe (2014). Coordination of microbial metabolism. Nature Reviews. Microbiology, 12(5), 327–340. De Jong, H. (2002). Modeling and simulation of genetic regulatory systems: A literature review. Journal of Computational Biology, 9, 67–103. Del Vecchio, D., Dy, A. J., & Qian, Y. (2016). Control theory meets synthetic biology. Journal of the Royal Society Interface, 13(120). Edwards, R., Machina, A., McGregor, G., & van den Driessche, P. (2015). A modelling framework for gene regulatory networks involving transcription and translation. Bulletin of Mathematical Biology, 77, 953–983. Farmer, W., & Liao, J. (2000). Improving lycopene production in Escherichia coli by engineering metabolic control. Nature biotechnology, 18, 533–537. Fung, Eileen, Wong, Wilson W., Suen, Jason K., Bulter, Thomas, gu Lee, Sun, & Liao, James C. (2005). A synthetic gene-metabolic oscillator. Nature, 435(7038), 118–122. Glass, Leon, & Kauffman, Stuart A. (1973). The logical analysis of continuous, nonlinear biochemical control networks. Journal of Theoretical Biology, 39, 103–129. Gupta, Apoorv, Reizman, Irene M. Brockman, Reisch, Christopher R., & Prather, Kristal L. J. (2017). Dynamic regulation of metabolic flux in engineered bacteria using a pathway-independent quorum-sensing circuit. Nature biotechnology, 35(3), 273–279. He, F., Murabito, E., & Westerhoff, H. V. (2016). Synthetic biology and regulatory networks: where metabolic systems biology meets control engineering. Journal of the Royal Society Interface, 13(117). Heinrich, R., & Schuster, S. (1996). The regulation of cellular systems. Chapman & Hall. Hudson, D., & Edwards, R. (2016). Dynamics of transcription-translation networks. Physica D, 331, 102–111. Gouzé, J. -L. (1998). Positive and negative circuits in dynamical systems. Journal of Biological Systems, 6(1), 11–15. Keasling, J. D. (2012). Synthetic biology and the development of tools for metabolic engineering. Metabolic Engineering, 14, 189–195. Khalil, H. K. (2002). Nonlinear systems. Prentice Hall, New Jersey. Kuntz, J., Oyarzún, D. A., & Stan, G.-B. (2014). Model reduction of genetic-metabolic systems using timescale separation. In System theoretic approaches to systems and synthetic biology. Springer-Verlag. Liu, Di, Mannan, Ahmad A., Han, Yichao, Oyarzún, Diego A., & Zhang, Fuzhong (2018). Dynamic metabolic control: towards precision engineering of metabolism. Journal of Industrial Microbiology & Biotechnology, 45(7), 1–9. Mannan, A. A., Liu, D., Zhang, F., & Oyarzún, D. A. (2017). Fundamental design principles for transcription-factor-based metabolite biosensors. ACS Synthetic Biology, 6(10), 1851–1859. Oyarzún, D. A., & Chaves, M. (2015). Design of a bistable switch to control cellular uptake. Journal of the Royal Society Interface, 12(20150618). Oyarzún, D. A., Chaves, M., & Hoffmeyer-Zlotnik, M. (2012). Multistability and oscillations in genetic control of metabolism. Journal of Theoretical Biology, 295, 139–153. Oyarzún, D. A., & Stan, G.-B. V. (2013). Synthetic gene circuits for metabolic control: design trade-offs and constraints. Journal of the Royal Society Interface, 10(78). Remy, E., Ruet, P., & Thieffry, D. (2008). Graphics requirement for multistability and attractive cycles in a Boolean dynamical framework. Advances in Applied Mathematics, 41, 335–350.

Reznik, E., Kaper, T. J., & Segrè, D. (2013). The dynamics of hybrid metabolic-genetic oscillators. Chaos, 23, 013132. Richard, A., & Comet, J.-P. (2007). Necessary conditions for multistationarity in discrete dynamical systems. Discrete Applied Mathematics, 155, 2403–2413. Ropers, D., Baldazzi, V., & de Jong, H. (2011). Model reduction using piecewiselinear approximations preserves dynamic properties of the carbon starvation response in Escherichia coli. IEEE/ACM Transactions on Computational Biology and Bioinformatics, 8, 166–181. Segel, L. A. (1988). On the validity of the steady state assumption of enzyme kinetics. Bulletin of Mathematical Biology, 50, 579–593. Soma, Y., & Hanai, T. (2015). Self-induced metabolic state switching by a tunable cell density sensor for microbial isopropanol production.. Metabolic Engineering, 30, 7–15. Soulé, C. (2003). Graphical requirements for multistationarity. ComPlexUs, 1, 123–133. Stevens, Jason T., & Carothers, James M. (2015). Designing RNA-based genetic control systems for efficient production from engineered metabolic pathways. ACS Synthetic Biology, 4(2), 107–115. Thomas, R. (1973). Boolean formalization of genetic control circuits. Journal of Theoretical Biology, 42, 563–585. Thomas, R., & Kaufman, M. (2001). Multistationarity, the basis of cell differentiation and memory. I. structural conditions of multistationarity and other nontrivial behavior. Chaos, 11, 165–179. Wang, B., Kitney, R. I., Joly, N., & Buck, M. (2011). Engineering modular and orthogonal genetic logic gates for robust digital-like synthetic biology. Nature Communications, 2, 580. Xu, P., Li, L., Zhang, F., Stephanopoulos, G., & Koffas, M. (2014). Improving fatty acids production by engineering dynamic pathway regulation and metabolic control. Proceedings of the National Academy of Sciences, 111(31), 11299–11304. Yagil, G., & Yagil, E. (1971). On the relation between effector concentration and the rate of induced enzyme synthesis. Biophysical Journal, 11, 11–27.

Madalena Chaves is a Senior Researcher in the team Biocore at Inria Sophia Antipolis (France). She obtained a PhD in Mathematics from Rutgers University in 2003 and an HDR – Habilitation à diriger des recherches from the University of Nice – Sophia Antipolis in 2013. She was a visiting scientist at Sanofi-Aventis (New Jersey, USA) between 2003 and 2005, and held a post-doctoral fellowship at the Institute for Systems Theory and Automatic Control at the University of Stuttgart in 2005/2006. She is a researcher at Inria since 2007. In 2015 she joined the editorial board of SIAM Journal on Applied Dynamical Systems. Her research is in the area of dynamical systems, in particular problems related to stability, robustness and control of nonlinear systems, with applications to the modeling and analysis of biological networks.

Diego A. Oyarzún obtained his BSc and MSc in Electronic Engineering from U Santa María, Chile (2006), and his PhD from the Hamilton Institute, Maynooth University, Ireland (2010). He was a Research Fellow in Biomathematics at the Department of Mathematics, Imperial College London (2013–2018) and currently holds a joint faculty appointment between the School of Informatics and the School of Biological Sciences at the University of Edinburgh. His group develops mathematical methods for the analysis and design of biomolecular systems with applications in Systems & Synthetic Biology. Oyarzún holds advisory roles at the World Economic Forum and is associate editor for several journals in biotechnology and synthetic biology.