Equilibria and stability analysis of a branched metabolic network with feedback inhibition

: This paper deals with the analysis of a metabolic network with feedback inhibition. The considered system is an acyclic network of mono-molecular enzymatic reactions in which metabolites can act as feedback regulators on “preceding” enzymes of their own pathway, and in which one metabolite is the root of the whole network. We show, under mild assumptions, the uniqueness of the equilibrium. In the simpliﬁed case where inhibition only acts on the reactions having the root of the network as substrate, we show that the equilibrium is globally attractive. This requires that we impose conditions on the kinetic parameters of the metabolic reactions.


INTRODUCTION
The cellular metabolism is defined as the (huge) set of biochemical reactions that occur inside a living cell for growth and reproduction. It is usually represented by an intricate network connecting the involved biochemical species (called "metabolites"). The pathways of the network are called "metabolic pathways". In the metabolic engineering literature, it is widely accepted that "despite their immense complexity, metabolic systems are characterized by their ability to reach stable steady states" (quoted from Stephanopoulos, Aristidou & Nielsen (1997), Chapter 4). It should however be fair to recognize that a mathematical analysis of this fundamental stability property is a difficult question which was not much investigated. We shall limit ourselves to simple metabolic pathways which are made up of sequences of mono-molecular enzyme-catalysed reactions in the form X s → X p .
Those reactions can be inhibited by the presence of other metabolites in the network. Without this inhibition, the stability analysis is straightforward; we will then concentrate on networks with inhibition, like the one given by the aspartate amino-acid pathways (Umbarger, 1978), see Figure 1. In this network, each produced amino-acid inhibits an enzyme of its own pathway. This action can be seen as a negative feedback, that regulates the behavior of the network. Indeed, if we, for example, consider a large excess of isoleucine (X 20 ), the reaction X 16 −→ X 17 is shut down, so that the concentration of isoleucine is progressively reduced.
In Section 2, a model of metabolic networks such as the one of Figure 1 will be presented. The equilibria of these models will then be studied in Section 3, followed by a stability analysis in Section 4, where global attractivity of a single equilibrium is shown.
x 5 x 6 x 7 x 8 x 9 x 10 x 11 x 12 x 13 x 14 x 15 x 16 x 17 x 18 x 19 x 20 Fig. 1. Metabolic network representing the aspartate amino-acid pathways: the solid lines represent the reactions and the dash-dotted lines the inhibition produced by the state at the start of the arrow onto the reaction that lies at the end of the arrow. The root of the metabolic pathway is x 1 (aspartate), and the products are the corresponding amino acids: lysine (x 9 ), methionine (x 14 ), threonine(x 16 ), and isoleucine (x 20 ).
The non-genericity of the stability of the equilibrium is then illustrated in Section 5.

MODEL OF A METABOLIC NETWORK
In our model of a metabolic network made of enzymecatalysed reactions in the form X s → X p which are inhibited by other metabolites of the network, we will denote the inhibition factor by x [p] , a n p dimensional vector containing the molar fractions of all the metabolites inhibiting the reaction X s → X p . This reaction is characterized by a velocity ϕ sp (x s , x [p] ).
On the other hand, some metabolites of the considered network are used in stages of the metabolism that are not modeled in the considered network. Therefore, we must include consumption terms in the model for those reactions in the form X s → ... We will generically denote those terms ϕ s0 (x s ).
We impose that these reaction velocities satisfy the following assumption: • For all s, p such that the reaction X s → X p belongs to the metabolic network, the function j . Note that this must also be valid when the value of p is 0.
In order to define the class of metabolic networks that we consider, we need the following definition from graph theory: Definition 1. A directed graph is called an arborescence if, from a given node x, known as the root node, there is exactly one elementary path from this node to any other node y.
The metabolic networks that we consider in this paper then satisfy the following assumption Assumption 2.
• the involved species are denoted X 1 , X 2 , · · · , X n • the graphic representation of the network (with the different metabolites as nodes and the different reactions as oriented edges) is an arborescence with X 1 as root • the inhibition acting on a reaction X s → X p only results from the action of metabolites from the (sub)-arborescence rooted in X p Stemming from the definition and known properties of arborescences, Assumption 2 has the following consequence on the class of metabolic network that we consider: (i) Each metabolite is produced by a single other metabolite; (ii) There is no cycle of reactions; With these definitions and notations, we shall now define a mass-balance dynamical model in the forṁ where x = (x 1 , · · · , x n ) T ∈ IR n , and x i denotes the molar fraction of the metabolite X i inside the cell. The factor µ ≥ 0 represents the specific growth rate of the cell: we assume that the cell metabolism is analyzed during a period of exponential cell growth with a constant specific growth rate µ. The vector e 1 = (1, 0, · · · , 0) T and the scalar c denote the constant supply rate of the metabolite X 1 at the root of the network. The function Φ includes all the reaction velocities, whether they correspond to reactions inside the network or reactions consuming the metabolites of the network for use in subsequent stages of the metabolism.
In order to specify Φ(x), we introduce the following notations: Notation 1.
• P(j) = {k| the reaction X j → X k belongs to the network }. P(j) defines the set of all metabolites that are produced by reactions having X j as substrate. If there is a consumption term in the form It can easily be seen that P(j) \ {0} is a subsets of A(j) under Assumption 2.
From the arborescence structure, it is clear that we can separate the metabolites into three different families: • the root X 1 : we suppose that there is a constant supply rate c of X 1 so that the corresponding mass-balance equation is the following: Under Assumption 2, we can only have A particular case of this network is the metabolic chain with sequential feedback inhibition (cf. Chitour, Grognard & Bastin (2003)) with the last metabolite X n acting as an inhibitor of the first reaction X 1 −→ X 2 . In this case, all x [j] are empty, except x [2] = x n .

EQUILIBRIUM OF A METABOLIC NETWORK
Based on equations (1)-(2)-(3), we can now compute the mass-balance of the whole arborescence: and of the arborescence that has its root in X j Those expressions will be critical in the proof of the following proposition: Proposition 1. If Assumptions 1 and 2 are satisfied, then: (C) if µ > 0, then system (1)-(2)-(3) has a unique equilibriumx = (x 1 , · · · ,x n ) in IR n + . Moreover, the solutions of (1)-(2)-(3) are bounded for any initial condition in IR n + .
Proof: (A) is easily seen by considering the system on the boundaries of the positive orthant.
The proofs of (B) and (C) are very similar. We write the proof for (B) and highlight the differences that arise for the proof of (C) .
We will first consider system (2)-(3) with x 1 =x 1 as constant input. For any value ofx 1 , we will denote by (x 2 , · · · ,x n ) the equilibrium of (2)-(3); this equilibrium is a function ofx 1 , so that we will state that it is (x 2 , · · · ,x n )(x 1 ). We will now show, by induction, that every elementx i is an increasing function ofx 1 (resp. non-decreasing in case (C)).
The initial step of the proof considers the equilibrium of (3) with x i =x i as constant input is the only solution. Also, the left-hand side of this equation is an increasing function ofx i (resp. nondecreasing in case (C)) and a decreasing function of x j . It is then easily seen that, if we increasex i ,x j needs also to be increased (resp.increased or kept constant) to keep this equality satisfied. We then have that, in this case,x j (x i ) is an increasing function such thatx j (0) = 0 (resp. non-decreasing function such thatx j (0) = 0). When µ = 0 (which can only happen in case (B)), the definition ofx j could be limited to an interval [0,x m i ). Let us now make the following induction hypothesis: for a given j, the functionsx k (x j ) are increasing (resp. non-decreasing) functions for all k ∈ A(j) with x k (0) = 0. We then study the equilibrium of the massbalance of the arborescence that has its root in X j . From (5): With a similar argument to that of the initial step, we see thatx j is an increasing (resp. non-decreasing) function ofx i and thatx j (0) = 0. The same can be said for allx k with k ∈ A j because they are already increasing (resp. non-decreasing) functions of x j .
By induction, we then have that everyx k (x 1 ) is an increasing function ofx 1 defined on the interval [0,x m 1 ) (resp. non decreasing function defined for allx 1 ). An equilibrium of the whole system then has to satisfy the equilibrium of the total mass-balance. From (4), this comes to: The system admits as many equilibria as this equation has roots. In case (B), the left-hand side is increasing from 0 whenx 1 increases from 0 tox m 1 (because of the second term). Therefore, if there exists an equilibrium, it is unique. In case (C), the left-hand side is strictly increasing from 0 to +∞ whenx 1 increases from 0 to +∞ (because of the µx 1 term), so that the equilibrium exists and it is unique.
The final point of (C) is a direct consequence of (4); this implies x l which clearly implies boundedness of the solutions when µ > 0. ✷ Uniqueness of the equilibrium, especially when it is coupled with boundedness of solutions, gives hope of the possibility of having some general result about the structural global asymptotic stability of the equilibrium. In the next section, we study a special case of feedback inhibition, in which we will be able to prove global attractivity.

STABILITY OF A CLASS OF METABOLIC NETWORKS
In this section, we will study the stability of the equilibrium of a class of metabolic networks belonging to the family that was described in Section 2. In order to do that, we will first present a technical lemma.

Technical lemma
In order to analyze the stability of the considered metabolic networks, we will use the small-gain theorem of De Leenheer, Angeli & Sontag (2003) for the interconnection of monotone systems. Due to space limitation, the details of this theory are omitted here, and we only give two lemmas that we will use in the following sections (and which are proven in Grognard, Chitour & Bastin (2003). Let us consider the stability of a systemẋ which is not cooperative, but which satisfies the following assumption: Assumption 3. Each off-diagonal element of the Jacobian J of f (x) is sign-definite (independent of x).
We then build a new systemẋ = F (x, −v) according to the following construction: for all i, j with j = i, if J ij < 0, then replace x j in f i (x) with the constant −v j . Define u 1 as the vector that contains all the constants v j that are necessary for the preceding construction (not necessarily all j have been required); we denote those constants by v s1 , · · · , v s l . From this construction, it directly appears that the systeṁ is cooperative (in the sense defined in De Leenheer et al. (2003)) We then define the output of (7) as y 1 = (x s1 , · · · , x s l ) T . The interconnection of system (7) with the trivial system through the feedback connections u 1 = −y 2 and u 2 = y 1 results in the original system (6). We can then show the following lemma Lemma 1. Suppose that the solutions of (6) are bounded and that system (7) has a unique globally attractive equilibrium for any constant input u 1 (this equilibrium defines an input output characteristicȳ 1 = k y1 (u 1 ) for system (7)). Suppose that the expression of y 1 = k y1 (u 1 ) is unknown, but that it is known that it satisfies equation then system (6)  and if the unique equilibrium of system (7) with u 1 constant is globally attractive for any u 1 .

.2 Feedback inhibition
The metabolic networks that we will consider for our stability study are made up of an arborescence of mono-molecular enzyme-catalysed reactions which satisfy the assumption that only one type of inhibition is present: the inhibition of the reactions that use X 1 as substrate. Also, a metabolite X m can only inhibit the reaction X 1 −→ X p , which is the first reaction of the unique elementary path linking X 1 to X m (this is a consequence of the third point of Assumption 2). This translates into Assumption 4.
A mass-balance model for such a network fits in the structure that was described in the previous section, so that we know that there is a single equilibrium. The model can be written as (1)-(2)-(3), but is particularized due to the presence of Assumption 4: where x ki is a product of x 1 and x k of x l = x 1 . We impose the boundedness of the partial derivatives of ϕ ij in the following assumption: We now define a new notation Notation 2. From the arborescence structure, we know that there exists a unique path from X 1 to any metabolite X s . This path takes the form X 1 → X kj → · · · X k → X l → · · · → X w → X s ; if X s is an arbitrary metabolite, we will store the indices of this path (without 1 and s) in C s ; alternatively, if x s corresponds to some x [kj ] b , we will also denote this path C [kj ] b . Similarly, we denote g s (k) or g [kj ] b (k) the index of the metabolite that follows X k in the path that connects X 1 to X s . This allows for the following theorem Theorem 2. If Assumptions 1, 2, 4 and 5 are satisfied and and ϕ 1ki is bounded for all i ∈ {1, · · · , r} (0 ≤ ϕ 1ki ≤ B ki ), then the equilibrium of system (10)-(11)-(12) is globally attractive in the positive orthant.
Proof: The proof of this Theorem is an application of Lemma 1 (see Grognard et al. (2003)). The decomposition of (10)-(11)-(12) as in (7) is done by replacing all the inhibiting terms x [p] with the constant vector −u 1 .
We see that, despite the fact that the inhibition is classically presented as a negative feedback that regulates the system, we have only been able to prove global attractivity under the restrictive condition (13), while stability is easily seen in the absence of inhibition. This condition is very strong, especially if the specific growth rate is small. We will analyze this condition further in the next section.

LIMIT CYCLES IN METABOLIC NETWORKS
Having obtained the sufficient result for global attractivity of Theorem 2, it is relevant to ask two questions: is condition (13) necessary and sufficient, on the one hand, and is Theorem 2 still valid without condition (13) on the other hand? The answer to the first question is easily seen to be "no": the stability of the equilibrium is retained even if condition (13) is slightly violated (a few simulations of simple systems is already convincing). In this section, we will show that the answer to the second question is also "no": without condition (13), the stability can be lost, so that we see that the stability of the metabolic networks is not a simple consequence of the structure of the models. Indeed, in this section, we shall exhibit an example where the equilibrium becomes unstable with a limit cycle (Hopf bifurcation) when condition (13) is not satisfied. We will concentrate on the stability of the equilibrium of a simple sequential pathway of four metabolites without branching and with sequential feedback inhibition (that was presented at the end of Section 2). Each metabolite produces a single other metabolite, and X 1 −→ X 2 is inhibited by the last metabolite, X 4 . We can directly apply Theorem 2 to this system: where we take ϕ 1 (x 1 , x 4 ) = which is Hurwitz for p < 56.4519, and is not Hurwitz for p larger than that value. This transition from a stable to an unstable equilibrium is illustrated on Figure  2, where the time responses of the four states is illustrated for the value of p = 0 (no inhibition), p = 10 (weak inhibition), and p = 60 (strong inhibition). In the latter case, oscillations appear. This corresponds to a limit cycle in the state-space. A Hopf bifurcation has taken place in p = 56.4519. Despite this oscillation, the reaction rates for the three reactions X 2 −→ X 3 , X 3 −→ X 4 , and X 4 −→ ... are close to their maximum after the transient. The only limiting reaction is X 1 −→ X 2 which, due to the inhibition is far from its maximum reaction rate.

CONCLUSION
In this paper, we have shown that a large class of models of metabolic system only has a single equilibrium. We then have proven that, under a small gain condition, this equilibrium is globally attractive in the particular case where inhibition only acts on reaction having the root as substrate. Finally, we have shown that stability of this equilibrium is not a generic property of the metabolic systems: a condition needs to be imposed on the parameters to have stability (similar to the small gain condition that we found).