Claudia López Zazueta

Olivier Bernard

Jean-Luc Gouzé

Claudia López Zazueta
email: claudia.lopez-zazueta@inria.fr

Dynamical Reduction of Linearized Metabolic Networks Through Quasi Steady State Approximation

Keywords: Biomolecular Metabolic Modeling, Slow-Fast Analysis, Tikhonov Theorem, Microalgae, Unbalanced Growth QSSA, quasi steady state assumption, FBA, flux balance analysis, QSS, quasi steady state, DRUM, dynamic reduction of unbalanced metabolism, SSA, steady state assumptions, MBM, macroscopic bioreaction models, DFBA, dynamical flux balance analysis, EFM, elementary flux modes, ODE, ordinary differential equations

Metabolic modeling has gained accuracy in the last decades, but the resulting models are of high dimension and difficult to use for control purpose. Here we propose a mathematical approach to reduce high dimensional linearized metabolic models, which relies on time scale separation and the Quasi Steady State Assumption. Contrary to the Flux Balance Analysis assumption that the whole system reaches an equilibrium, our reduced model depends on a small system of differential equations which represents the slow variables dynamics. Moreover, we prove that the concentration of metabolites in Quasi Steady State is one order of magnitude lower than the concentration of metabolites with slow dynamics (under some flux conditions). Also, we propose a minimization strategy to estimate the reduced system parameters. The reduction of a toy network with the method presented here is compared with other approaches. Finally, our reduction technique is applied to an autotrophic microalgae metabolic network.

Introduction

Metabolic modeling has proved to be a very powerful tool to get a better insight into the metabolism of an organism, while assessing the main fluxes throughout its metabolic network. [START_REF] Baroukh | DRUM: A new framework for metabolic modeling under non-balanced growth. Application to the carbon metabolism of unicellular microalgae[END_REF][START_REF] Nazaret | An old paper revisited: "A mathematical model of carbohydrate energy metabolism. Interaction between glycolysis, the Krebs cycle and the H-transporting shuttles at varying ATPases load[END_REF] This approach has gained accuracy in the last decades and turned out to be particularly efficient to improve the production of target molecules, by understanding the biological processes that influence the metabolites. For example, metabolic modeling has clarified the production of triacylglycerols from microalgae and carbohydrates from cyanobacteria. [START_REF] Baroukh | A state of the art of metabolic networks of unicellular microalgae and cyanobacteria for biofuel production[END_REF] Both compounds can then be turned into biofuel (biodiesel and bioethanol, respectively) with expected reduced environmental impacts. [START_REF] Lardon | Life-cycle assessment of biodiesel production from microalgae[END_REF] The metabolisms of microalgae and cyanobacteria are driven by the solar flux which supports the fixation of CO 2 . As a consequence, the periodic fluctuation of light induces instationarity of their metabolisms, with accumulation of metabolites (especially lipids and carbohydrates). Such metabolisms are therefore never at steady state. However, most of the approaches dedicated to metabolism analysis assume balanced growth, i.e. Steady State Assumption (SSA), which are not totally justified both from a biological and from a mathematical point of view. For example, Flux Balance Analysis (FBA) [START_REF] Orth | What is flux balance analysis?[END_REF] or Macroscopic Bioreaction Models (MBM) [START_REF] Provost | Metabolic design of macroscopic bioreaction models: application to Chinese hamster ovary cells[END_REF] are based on linear algebra to solve the equation N • V = 0, where N is the stoichiometric matrix and V is the vector of intracellular reaction rates. These approaches assume therefore two hypotheses. First that the derivative of each intracellular compound (per biomass unit) is zero, and second that the dilution rate is neglected. Some approaches have attempted to introduce dynamics, for instance, assuming that the cell optimizes at each time instant an objective criterion using Dynamic Flux Balance Analysis (DFBA) [START_REF] Mahadevan | Dynamic flux balance analysis of diauxic growth in Escherichia coli[END_REF] or considering external con-ditions that might change continuously in a hybrid system. [START_REF] Kaplan | Modeling and simulation of metabolic networks for estimation of biomass accumulation parameters[END_REF] However, the resulting metabolic models are of high dimension, which makes their mathematical analysis and their parameter identification complex. Identifying conditions to maximize productivity by a rigorous mathematical analysis is generally not possible. Here we propose an approach to reduce the dimension of a dynamical metabolic model, so that model based control strategies can be derived. Contrary to most of the existing work, the idea is to keep the dynamical component of the model, which is crucial especially when dealing with microalgae and cyanobacteria.

A first attempt in this direction was carried out with the Dynamic Reduction of Unbalanced Metabolism (DRUM) method. [START_REF] Baroukh | DRUM: A new framework for metabolic modeling under non-balanced growth. Application to the carbon metabolism of unicellular microalgae[END_REF] This modeling approach considers subnetworks in Quasi Steady State (QSS) that are interconnected by metabolites that can accumulate. Then, Elementary Flux Modes (EFM) are computed in each subnetwork to reduce them by a Quasi Steady State Assumption (QSSA). Finally, dynamics of the connector metabolites form a reduced system of Ordinary Differential Equations (ODE). This approach has proven to provide sound results, a way to dynamical reduce the model together with an efficient representation of the accumulation of lipids and carbohydrates in microalgae submitted to light / dark cycles. However, as most of the methods developed for metabolic analysis, it relies on a series of assumptions whose mathematical bases have not been rigorously established. Beyond the Quasi Steady State Assumption (QSSA), which assumes "fast" and "slow" parts on the metabolism without delimiting them rigorously, these approaches also neglect the intracellular dilution term due to growth.

The main objective of our work is to provide mathematical foundations for the reduction of metabolic networks to dynamical models of low dimension. To achieve this objective, in a first stage we simplify the approach assuming linear kinetics for the whole metabolism. However, our method can be applied to any metabolic network whose kinetics can be locally linearized.

In this article, we compute an exact reduced model for any network in the class of systems addressed, whose metabolism is forced by a continuous input. In addition, we propose an algorithm to estimate the parameters of the reduced system.

In the method presented in this work we conserve the factor of dilution, which improves the precision of the approximation and preservers qualitative (stability) features of the original system.

We also show that Flux Balance Analysis (FBA) cannot be applied to the class of system studied in this work and and we discuss the accuracy of the Dynamic Reduction of Unbalanced Metabolism (DRUM) for linear systems, approximation which relies on a series of hypotheses including Quasi Steady State Assumption (QSSA) and omission of the dilution due to growth.

In the first section of this article, we present the class of metabolic systems under study, which consists of metabolic networks with slow and fast reactions forced by a continuous input. We also write the system of Ordinary Differential Equations corresponding to any metabolic network of this class, and we formulate it as a slow-fast system. Then, we verify the conditions for applying the Theorem of Tikhonov for singularly perturbed systems [START_REF] Hoppensteadt | Properties of solutions of ordinary differential equations with small parameters[END_REF][START_REF] Kokotović | Singular perturbation methods in control: analysis and design[END_REF][START_REF] Tikhonov | Differential Equations[END_REF] to obtain a reduced system. From these approximations, we demonstrate that the concentration of each metabolite in Quasi Steady State (QSS) is one order of magnitude lower than the concentration of any metabolite with slow dynamics, subject to some flux constraints. This result is one of the main outcome of the paper.

In the second section of this article, we introduce a method based on minimization for the calibration of the reduced system. In addition, when data of a metabolite involved in fast reactions is available, we can compute the corresponding parameters for its estimation in Quasi Steady State (QSS).

In the third section of this article we apply the approach to a toy metabolic network. The toy model includes a periodic input and reflects standard bricks in metabolic networks: combination of reversible and nonreversible reactions, with chains and cycles. We then compare with the reduced model with Flux Balance Analysis (FBA) and the Dynamic Reduction of Unbalanced Metabolism (DRUM) method.

Finally, we apply our method to an autotrophic microalgae metabolic network. We use the simplified network of Yang et al. [START_REF] Yang | Energetics and carbon metabolism during growth of microalgal cells under photoautotrophic, mixotrophic and cyclic light-autotrophic/dark-heterotrophic conditions[END_REF] The reduction is carried out considering the macromolecules as the elements with the highest concentration. Using the data from Lacour et al., [START_REF] Lacour | Diel variations of carbohydrates and neutral lipids in nitrogen-sufficient and nitrogen-starved cyclostat cultures of isochrysis sp[END_REF] we calibrated the reduced model. The simulations of the reduced model accurately fits the experimental data.

To conclude, we draw some conclusions and provide perspective for a rigorous slow-fast analysis of a larger class of metabolic networks.

The Slow-Fast Linear Metabolic Model

Metabolic system hypotheses

The class of metabolic networks studied in this paper consists of metabolic systems with fast and slow reactions, which are forced by a continuous in-put. To describe the elements of this class, we focus on a general network assuming that one internal component of this (X m ) is consumed with a low rate (see Figure 1). We consider then two subnetworks, before and after this component. The first subnetwork contains m -2 metabolites, while the second contains n -m -1 metabolites. We also suppose that the entire network is a connected graph, meaning that any pair of nodes (metabolites) are connected to each other by a non directed path (a sequence of reactions without considering their direction). The model is kept general in the sense that a priori all the fluxes are possible within the two subnetworks. We suppose that there are fast and slow reactions in the system. The fast reactions are depicted in Figure 1 by blue arrows and their reaction rates are denoted by k ij /ε, where ε is a very small positive number. A reaction with rate k ij /ε consumes only the element X j and produces just X i (single reactant-single product hypothesis). Similarly for the slow reactions, which are represented by the black arrows, the reaction with rate k ij consumes X j and produces X i . Notice that reversible reactions are allowed in the fast subnetworks (see Figure 2).

We assume that one dimensional metabolic flux enters in the metabolism from X 1 and finishes into variable X n . The input I(t) is a positive continuous function of time t in an interval [0, T ], feeding in the system at a slow rate. Moreover, we consider a constant dilution affecting every metabolite. The rate of dilution µ > 0 is considered as a parameter smaller than any reaction rate (a classical hypothesis [START_REF] Stephanopoulos | Metabolic engineering: principles and methodologies[END_REF] ).

For a metabolic network in the class of systems addressed in this work (see Figure 1), the hypothesis of single reactant -single product reactions implies that the ODE of metabolite concentrations is a non-homogeneous linear system with a continuous positive input. However, the approach developed in this paper is suitable for any metabolic network whose kinetics can be locally linearized. In Appendix E of the Supporting Information we show a metabolic network example with nonlinear enzymatic reactions, whose kinetics are linearized around a functional point.

Figure 2: Single reactant-single product hypothesis: fast reaction k ji /ε just consumes one metabolite (X i ) and produces another (X j ). This conditions leads to a linear system of equations.

Henceforth, we focus on a metabolic network with linear kinetics and a continuous positive input. In order to obtain some conclusions about the accumulation of metabolites, we make some assumptions about the fluxes in the network.

Assumption 1. We assume that a connector metabolite is not the origin of any fast reaction, nor the origin of an output with large rate.

On the other hand, for detecting when the input flux leads to the accumulation of metabolites which have fast outflows, we introduce the next definitions. For these, we consider that our metabolic network model (Figure 1) can be represented with a directed graph, where the metabolites are the vertex and the reactions are the directed edges (arrows). Only inputs and outputs are allowed to be edges with one vertex: final or initial vertex, respectively (the dilution due to growth does not appear on the graph).

Definition 1. (Flux ).

A flux from X i to X j is a directed path which has as initial vertex X i and as final vertex X j .

The following definition of trap is equivalent to the definition for compartmental linear systems, [START_REF] Jacquez | Qualitative theory of compartmental systems Siam Review[END_REF] considering the metabolites as compartments and the arrows (reactions, inputs and outputs) as flows into and out of the compartment. A compartmental system can also be represented by a directed graph. Then, a trap is defined as a compartment or a set of compartments from which there are no transfers to the environment nor to compartments that are not in the trap. [START_REF] Jacquez | Qualitative theory of compartmental systems Siam Review[END_REF] Indeed, it has been proved that an autonomous linear compartmental system has a trap if, and only if, its compartmental matrix is singular. [START_REF] Jacquez | Qualitative theory of compartmental systems Siam Review[END_REF] Definition 2. (Trap). Consider a graph with set of vertices N and a subset of this

T = {X k 1 , . . . , X k l } ⊂ N, n > l ≥ 1. We say that T is a trap if • for every vertex X k i ∈ T there is no flux from X k i to any metabolite of N \ T and
• no X k i ∈ T has an output to the exterior of the graph.

In this case, we also say that X k i ∈ T is in a trap.

Definition 3. (Flux trap)

. Consider a flux F with initial vertex X 1 and final vertex X n in a graph with vertices N = {X 1 , . . . , X n }. We say that the graph has a trap for the flux F if there is a subset

T F = {X k 1 , . . . , X k l } ⊂ N \ {X 1 , X n }, such that
• T F is a trap (as a consequence, there is no flux from any

X k i ∈ T F to X n ) and
• for every vertex X k i ∈ T F there is a flux from X 1 to X k i .

We also say that X k i ∈ T F is in a flux trap. When it is clear which is the flux F taken into account, we only say that the graph has a flux trap and we denoted by T.

Remark. Notice that in a metabolic system with linear kinetics where the dilution factor is omitted, if there is a trap, then the system is undetermined, i.e. its (Jacobian) matrix is singular. [START_REF] Jacquez | Qualitative theory of compartmental systems Siam Review[END_REF] This case results in additional complexity from a mathematical point of view and it will not be discussed in this paper, because the factor of dilution is not omitted in our approach.

On the other hand, we include the mathematical concept of flux trap, to detect the accumulation of metabolites with fast dynamics as a consequence of a input flux. However, from a biological point of view, the presence of traps or flux traps is very unlikely, unless organisms have been modified to indefinitely accumulate a metabolite.

In the system depicted in Figure 1 we assume the existence of a flux from X 1 to X n . It then has a flux trap if there exists in one of the subnetworks a metabolite, X i for some i = 2, . . . , n -1, i = m, such that there is a flux from X 1 to X i , but there is no flux from X i to X n . A flux going from X 1 to X n has to pass through X m as well. Hence, there is a flux trap if and only if there is a flux trap in one subnetwork with fast dynamics.

Assumption 2. The model represented in Figure 1 has no flux trap for the assumed flux going from X 1 to X n . We refer to this conditions saying that the system has no flux traps.

Assumption 2 implies that if X i is any metabolite in one of the subnetworks with fast dynamics and it is reached by the flux, then it has one fast outflow at least.

Dynamics of Metabolite Concentrations

We write the ordinary differential equations that describes the metabolite concentrations of the model in Figure 1 as follows:

dX 1 dt = I(t) -(k 21 + µ)X 1 (1)
dX 2 dt = k 21 X 1 + m-1 j=3 k 2j ε X j - m-1 i=3 k i2 ε + µ X 2 . . . dX m dt = k m,m-1 ε X m-1 -k m+1,m + µ X m . . . dX n-1 dt = n-2 j=m+1 k n-1,j ε X j - n i=m+1 i =n-1 k i,n-1 ε + µ X n-1 dX n dt = k n,n-1 ε X n-1 -µ • X n ,
with initial conditions X i (0) = x 0 i for every i = 1, . . . , n. System (1) is a non-homogenous linear system, continuous in a domain

[0, T ] × R n + if I(t) is continuous in the interval [0, T ]. Moreover, system (1)
is positive if and only if its (Jacobian) matrix is a Metzler matrix and I(t) ≥ 0. [START_REF] Farina | Positive Linear Systems: Theory and Applications[END_REF] For the mathematical analysis, we suppose that all the variables and the time have been non dimensionalized by appropriately normalizing the concentrations and time constants.

Parameter definition and order of magnitude

In order to distinguish the two time-scales in system (1), we make the following hypothesis about the order of the parameters. For this purpose, we use the "big O" or Landau order symbol. [START_REF] Narang-Siddarth | Nonlinear Time Scale Systems in Standard and Nonstandard Forms: Analysis and Control[END_REF] Although the Landau symbol is defined for functions, we can also use it to indicate the order of our parameters. For this, we consider a parameter as the constant function which takes its value.

Recall that fast reactions have rates k ij /ε, while slow reactions k ij . We then suppose that ε is a small positive number such that

k ij = O(1)
as ε → 0 ∀i, j.

Henceforth, we consider ε fixed. Additionally, I(t) is a slow input and µ is lower or equal to the magnitude of any slow reaction (an usual hypothesis [START_REF] Stephanopoulos | Metabolic engineering: principles and methodologies[END_REF] ). Hence, we suppose

I(t) = O(1) ∀t ∈ [0, T ], (2)
µ ≤ O(1).

For integrating the factor of dilution µ without affecting the order of the parameters in the system (1), we define

k j ε := n i=2 i =j k ij ε + µ (3)
The definition given in (3) implies

n i=2 i =j k ij < k j . (4)
Also, k j = O(1), as a consequence of condition (2) for µ.

Fast metabolites variable rescaling

In order to write the original system (1) as a singularly perturbed system, we define new variables for the metabolites with fast dynamics and their initial conditions as follow:

Y i : = X i ε Y i (0) = y 0 i := x 0 i ε .
Hence, we represent system (1) as the following slow-fast system:

dY 2 dt = 1 k 21 X 1 + m-1 j=3 k 2j Y j -k 2 Y 2 (5)
. . .

dY m-1 dt = 1 m-2 j=2 k m-1,j Y j -k m-1 Y m-1 dY m+1 dt = 1 k m+1,m X m + n-1 j=m+2 k m+1,j Y j -k m+1 Y m+1 . . . dY n-1 dt = 1 n-2 j=m+1 k n-1,j Y j -k n-1 Y n-1 dX 1 dt = I(t) -(k 21 + µ)X 1 X 1 (0) = x 0 1 ( 6
)
dX m dt = k m,m-1 Y m-1 -(k m+1,m + µ)X m X m (0) = x 0 m dX n dt = k n,n-1 Y n-1 -µ • X n X n (0) = x 0 n .
Notice that ε is a parameter that we previously fixed. From a mathematical viewpoint, we will consider → 0 to apply Tikhonov's Theorem. In the linear system (5) we can consider the initial conditions as any nonnegative values, since global and local asymptotic stability are equivalent.

Hereafter, we say that Equation ( 5) and Equation ( 6) are the fast and the slow part of the system (1), respectively. Furthermore, since the dynamics of Y i = X i /ε are fast according to (5), we consider the metabolites X 2 , . . . , X m-1 , X m+1 , X n-1 to be in Quasi Steady State.

Model properties

Let us rewrite the fast part (5) under a matrix form. We define

K := K 1 0 0 K 2 ,
where

K 1 :=      -k 2 k 23 . . . k 2,m-1 k 32 -k 3 . . . k 3,m-1 . . . . . . . . . k m-1,2 k m-1,3 . . . -k m-1      , K 2 :=      -k m+1 k m+1,m+2 . . . k m+1,n-1 k m+2,m+1 -k m+2 . . . k m+2,n-1 . . . . . . . . . k n-1,m+1 k n-1,m+2 . . . -k n-1      .
Then system (5) can be written as

dY dt = I + K • Y, (7)
with Y = (Y 2 , . . . , Y m-1 , Y m+1 , . . . , Y n-1 ) T and

I = (k 21 • X 1 , 0, . . . , 0, (m-1)-th entry k m+1,m • X m , 0, . . . , 0) T .
Property 1. For every fixed scalars X 1 , X m , X n , system (7) has a stable equilibrium point.

Proof. The matrices K 1 , K 2 are strictly column diagonally dominant as a consequence of (4). Therefore, they are stable matrices. [START_REF] Horn | Matrix Analysis[END_REF] Remark. For any constant values of X 1 , X m , X n , we obtain

Y := (Y 2 , Y 3 , . . . , Y m-1 , Y m+1 , . . . , Y n-1 ),
the equilibrium point of system (5) after computing the inverses of the matrices K 1 , K 2 . Indeed, K 1 and K 2 are nonsingular matrices, because they are strictly column diagonally dominant. [START_REF] Horn | Matrix Analysis[END_REF] Hence,

     Y 1 Y 2 . . . Y m-1      = (K 1 ) -1 •      -k 21 X 1 0 . . . 0      , (8)
     Y m+1 Y m+2 . . . Y n-1      = (K 2 ) -1 •      -k m+1,m X m 0 . . . 0      .
The equalities in (8) imply that the fast variables Y i in QSS are linear combinations of X 1 and X m .

Boundary Layer Correction

The boundary layer correction is aiming at correcting the initial transient of the fast variables, so that the approximation is also accurate at the very first time of the simulation. For this, we consider a time variable τ = t/ . The boundary layer correction is defined as the function Ŷ (τ ) = Y (t)-Y (t).

Then, the boundary layer problem of system ( 5)-( 6) is written as

d Ŷ2 dτ =k 21 x 0 1 + m-1 j=3 k 2j Ŷj + Y j (0) (9)
-k 2 Ŷ2 + Y 2 (0) . . . d Ŷm-1 dτ = m-2 j=2 k m-1,j Ŷj + Y j (0) -k m-1 Ŷm-1 + Y m-1 (0) d Ŷm+1 dτ =k m+1,m x 0 m + n-1 j=m+2 k m+1,j Ŷj + Y j (0) -k m+1 Ŷm+1 + Y m+1 (0) . . . d Ŷn-1 dτ = n-2 j=m+1 k n-1,j Ŷj + Y j (0) -k n-1 Ŷn-1 + Y n-1 (0) with initial conditions Ŷi (0) = Y i (0) -Y i (0) = x 0 i ε -Y i (0).
For obtaining Quasi Steady State Approximation (QSSA) using the Theorem of Tikhonov, the stability of the origin for the boundary layer equation ( 9) is a necessary condition. [START_REF] Kokotović | Singular perturbation methods in control: analysis and design[END_REF] Notice that the (Jacobian) matrices of system (5) and (9) are equal. Then, as in Property 1, we have the following result.

Property 2. Consider the boundary layer system (9). Then its equilibrium point Ŷ (τ ) = 0 is stable.

Reduced system

From Equation (8) we obtain the expressions of Y m-1 and Y n-1 in terms of X 1 and X m , respectively. We then substitute in the slow part (6) and we obtain a reduced system with the variables X 1 , X m , X n . Hence, we express the QSS approximation of the original system (1), deduced from the Theorem of Tikhonov as follows:

X i := ε • Y i = ε • c i • X 1 ∀i = 2, . . . , m -1, (10)
X i := ε • Y i = ε • c i • X m ∀i = m + 1, . . . , n -1, dX 1 dt = I(t) -(k 21 -µ)X 1 ( 11
)
dX m dt = k m,m-1 c m-1 • X 1 -(k m+1,m + µ)X m dX n dt = k n,n-1 c n-1 • X m -µ • X n .
where c i are parameters that satisfy

Y i = c i •X 1 or Y i = c i •X m
, respectively, deduced from Equation (8). The initial conditions for the reduced system (11) are the same that for (6), i.e., X i (0) = x 0 i for i = 1, m, n.

Property 3. Let (X i ) be the solution of the original system (1) and (X i ) the solution of system ( 10)- (11). Then

X i (t) = X i (t) + O(ε) ∀t ∈ [0, T ], i = 1, m, n,
and there exists 0 ≤ T 0 < T such that

X i (t) = X i (t) + O(ε 2 ) ∀t ∈ [T 0 , T ], i = 2, . . . , n -1, i = m.
Additionally, if ( Ŷi ) is the solution of the boundary layer system (9) and

Xi := ε Ŷi for every i = 2, . . . , n -1, i = m, X i (t) = X i (t) + Xi t ε + O(ε 2 ) ∀t ∈ [0, T ].
Proof. After Properties 1 and 2, we conclude that system (5)-( 6) satisfies the conditions of Tikhonov's Theorem. [START_REF] Hoppensteadt | Properties of solutions of ordinary differential equations with small parameters[END_REF][START_REF] Kokotović | Singular perturbation methods in control: analysis and design[END_REF][START_REF] Levin | The asymptotic behaviour of the stable initial manifolds of a system of nonlinear differential equations[END_REF] Then a simple multiplication X i = εY i leads to the desired result.

Observe that the concentration of metabolites supposed to be in QSS is expressed as a linear combination of X 1 and X m multiplied by ε. Since ε is a very small positive number, this suggests that the concentration of the metabolites in the slow part of the system is higher.

Theorem 1. Suppose that the model in Figure 1 has no flux trap. Hence, with the notation of system (1) and Property 3, for every t ∈ [T 0 , T ],

X i (t) = ε • O X 1 (t) i = 2, . . . , m -1, (12)
X i (t) = ε • O X m (t) i = m + 1, . . . , n -1.
Moreover, for every t ∈ [T 0 , T ] we have

X i (t) = ε • O X j (t) , j ∈ {1, m, n}, (13)
i = 2, . . . , n -1, i = m.
Theorem 1 states that, after the initial fast transient, any metabolite in QSS has a concentration one order of magnitude lower than any metabolite in the slow part.

Proof of Theorem 1. To demonstrate the first affirmation we consider the equalities in (10) and we will show that

c i = O(1) if c i = 0, for i = 2, . . . , m- 1. Since K 1 is a nonsingular matrix, (K 1 ) -1 = 1 det(K 1 ) • C,
where C is the transpose matrix of cofactors of K 1 . [START_REF] Lancaster | The theory of matrices: with applications[END_REF] We then have according to Equation ( 8)

Y i = 1 det(K 1 ) C 1,i-1 • (-k 21 X 1 ), then c i = 1 det(K 1 ) C 1,i-1 • (-k 21 ).
If K 1 has no traps (i.e. the subnetwork with metabolites X 2 , . . . , X m-1 has no traps), then det

(K 1 ) = (-1) m-2 • O(k m-2 ij ),
as stated by Proposition 1 of Appendix A in the Supporting Information. Moreover, Corollary 2 (Appendix A in the Supporting Information) implies that the cofactors C 1,i-1 have order

C 1,i-1 = (-1) m-1 • O(k m-3 ij
).

On the other hand, if K 1 has a trap T not reached by the flux, as a consequence of Proposition 2 and Corollary 2 of Appendix A in the Supporting Information,

C 1,i-1 det(K 1 ) = (-1) • O(k -1 ij ) X i ∈ T C 1,i-1 det(K 1 ) = 0 X i ∈ T, We conclude that -k 21 • C 1,i-1 det(K 1 ) = O(1)
if C 1,i-1 = 0, for every i = 2, . . . , m -1. The same reasoning applies for K 2 and the variables which are linear combinations of X m . Hence, in accordance with (10),

X i = ε • O(c i • X 1 ) = ε • O(X 1 ) ∀i = 2, . . . , m -1, X i = ε • O(c i • X m ) = ε • O(X m ) ∀i = m + 1, . . . , n -1,
and we obtain the equalities in (12) from Property 3.

In order to verify the second affirmation, we consider the reduced system (11). Then we deduce that the local maximum or minimum points of X 1 , X m and X n satisfy

X m = k m,m-1 • c m-1 k m+1,m + µ X 1 X n = k n,n-1 • c n-1 µ X m . But c m-1 = O(1), c n-1 = O(1) and µ/k ij ≤ O(1) for all i, j as stated by (2). Thus, O(X 1 ) = O(X m ) ≤ O(X n ).
Finally, we obtain the equalities in (13) using Property 3.

Reduced model calibration

Calibration of the slow system

We suppose that the slow metabolite concentrations are measured at r different time instants 0

≤ t 1 < t 2 < • • • < t r ≤ T : Z i (t j ) = X i (t j ) + N (t j ) i = 1, m, n
where N (t j ) is the measurement noise. From these measurements, in a first stage we want to calibrate the parameters for the reduced model (11). The objective is to find the value for the parameter vector

θ = (θ 1 , θ 2 , θ 3 , θ 4 , θ 5 )
such that the solution of the reduced system

dX 1 dt = I(t) -θ 1 • X 1 X 1 (0) = x 0 1 ( 14
)
dX m dt = θ 2 • X 1 -θ 3 • X m X m (0) = x 0 m dX n dt = θ 4 • X m -θ 5 • X n X n (0) = x 0 n .
best fits the data Z 1 , Z m and Z n . For this, we define a cost criterion that characterizes the modeling error:

F 0 (θ) = i=1,m,n j=1,...,r Z i (t j ) -X i (t j , θ) 2
The approach consists in determining the best value θ = ( θ1 , θ2 , θ3 , θ4 , θ5 ) such that

F 0 ( θ) = min{F 0 (θ) : θ ∈ D},
with D ⊂ R 5 + a domain for the parameters. This can be carried out by using a minimization algorithm, based on gradient descent, or such as the Nelder-Mead algorithm.

Calibration for a metabolite in QSS

After the calibration of the slow system (14), we can estimate the concentration dynamics of a metabolite in QSS if there are data of this. For this purpose, we assume that a fast metabolite is measured at different time instants

0 ≤ t 1 < t 2 < • • • < t r ≤ T : Z i (t j ) = X i (t j ) + N (t j )
for some i ∈ {2, . . . , n -1} \ {m}.

We first use the value of the slow metabolites

X i (t, θ) i=1,m,n
obtained from calibration of the reduced system (14), to estimate the parameter α i such that

X i := α i • X 1 if i ∈ {2, . . . , m -1} X i := α i • X m if i ∈ {m + 1, . . . , n -1}.
Indeed we look for value of αi that minimize the differences

F 1 (α i ) = r j=1 α i • X 1 (t j , θ) -Z i (t j ) 2 for i ∈ {2, . . . , m -1} and F 2 (α i ) = r j=1 α i • X m (t j , θ) -Z i (t j ) 2
for i ∈ {m + 1, . . . , n -1}. Consequently, the least squares solution is the following:

21 αi = r j=1 Z i (t j ) • X 1 (t j , θ) r j=1 X 1 (t j , θ) 2 if i ∈ {2, . . . , m -1} αi = r j=1 Z i (t j ) • X m (t j , θ) r j=1 X m (t j , θ) 2 if i ∈ {m + 1, . . . , n -1}.
In general, measurements of fast metabolites are difficult, and therefore rarely done. It is worth noting that the reduction and calibration of the model for the fast metabolites can be done for any subset of measured metabolites.

Practical remark

Let us emphasize that, even if no metabolite with fast dynamics is measured, this approach provides a reduced model (14) which can be calibrated. If not all the accumulative metabolites are measured, the reduced model may also be calibrated, but the identifiability of the parameters should be demonstrated. [START_REF] Walter | Identification of parametric models from experimental data[END_REF] Furthermore, the dynamics of all the measured metabolites are eventually predicted by means of the reduced model of equations.

This is an important point, since the calibration of fast variables (which is indeed experimentally much more tricky) can be done only on a very restricted number of metabolites.

A Toy Metabolic Network

We consider a toy system of enzymatic reactions mimicking the characteristics of a real metabolic network but of lower dimension. In this model we have included some recurrent features of the metabolic network of a photoautotrophic microalga cell, such as a cycle and several reversible reactions. For instance, in the chain that includes the citric acid cycle, lower glycolysis, upper glycolysis and carbohydrate synthesis. We thus construct our example with a chain that contains one reversible reaction and a cycle. This toy Network N1 is depicted in Figure 3. Additionally, to represent the effect of a periodic factor on the metabolism (such as sun light for microalgae) we introduce a periodic input

I(t) = k[cos(t•ω)+1].
We suppose that slow reaction rates have order O(10 -2 ), and fast reaction rates O(10 1 ). Then, we choose ε := 10 -2 /10 1 = 10 -3 .

For instance, we suppose the slow reaction rate k 21 = 0.01 • min -1 = 10 -2 • min -1 . On the other hand, we assume that the fast reaction from X 2 to X 3 has rate equal to 20 • min -1 = 2 × 10 1 • min -1 . Then we establish

k 32 ε = 20 • min -1 , which implies k 32 = 20 ε • min -1 = 0.02 • min -1 . Recall that O(k 32 ) = O(k 21 ) = 10 -2 .
For the definition of the rest of parameters, see Table 1.

Following the method proposed in the first section of this article, we write a singularly perturbed system to apply Tikhonov's Theorem (see Figure 4).

I(t) = k[cos(t • ω) + 1]. Parameter Value Units Parameter Value Units ε 0.1 × 10 -2 - k 43 /ε 20 min -1 ω 0.4 × 10 -2 - k 32 /ε 20 min -1 µ 0.2 × 10 -2 min -1 k 68 /ε 20 min -1 k 0.2 × 10 -1 min -1 k 75 /ε 10 min -1 k 21 0.1 × 10 -1 min -1 k 23 /ε 10 min -1 k 54 0.2 × 10 -1 min -1 k 56 /ε 10 min -1 x 0 i 0.1 × 10 -1 µmol/m 3 k 87 /ε 20 min -1 k 98 /ε 10 min -1
Table 1: Parameters considered for the simulation of dynamics in Network N1 (Figure 3). The initial conditions x 0 i are all the same for i = 1, . . . , 9. 5)-( 6)) that is obtained after the variable rescaling of Network N1 (Figure 3). Notice that the all the metabolite concentrations have the same order of magnitude, as a consequence of defining Y i = X i /ε for the metabolites in the fast part. The parameters considered are specified in Table 1.

Moreover, we obtain an algebraic system that let us rewrite every variable of the fast part in terms of the variables of the slow part: sions of (15) in the equations of the slow part:

dX 1 dt = k[cos(t • ω) + 1] -(k 21 + µ)X 1 X 1 (0) = x 0 1 ( 16
)
dX 4 dt = k 21 • k 32 d 1 k 43 X 1 -(k 54 + µ)X 4 X 4 (0) = x 0 4 dX 9 dt = k 54 • k 75 • k 87 • k 6 d 2 k 98 X 4 -µX 9 X 9 (0) = x 0 9 .
As expected, the concentration of metabolites belonging to the fast part is one order of magnitude lower than the concentration of metabolites in the slow part (see Theorem 1 and Figure 5).

In the first section of this article, we described the boundary layer correction that can be applied to the variables in QSS. To see the effect of the boundary layer correction in this toy network, see Appendix C in the Supporting Information. 3). Thick line: numerical solution of the original system (1); thin line: reduced system obtained using the method developed in this work (15)-( 16). The parameters considered for the simulation are specified in Table 1.

Calibration of the slow dynamics

For this example, we consider that a few measurements of the metabolite concentrations are available. We represent it as the solution of the original system (1) for Network N1 plus a white noise:

Z i (t j ) := X i (t j ) + β(t j ) j = 1, 2, . . . , r (17)
where β ∼ N (σ i ) and σ i = m(X i ) × 10 -1 for every i = 1, . . . , n. For simplicity we suppose that the data is obtained at the same time instants t 1 , . . . , t r for the slow and the fast parts of the system. Following the notation in the second section of this article, we find the vector parameter θ that minimizes the cost function F 0 . We used the function fminsearch in Scilab to compute the minimum of the square differences function F 0 , with the Nelder-Mead algorithm (see Table 2). 16), written as Equation ( 14). The parameters considered for simulation of the original dynamics in Network N1 are specified in Table 1.

Calibration of coefficients for metabolites in QSS

After estimating the parameters θi we consider the equalities in (15) for finding α i such that

X i = α i • X 1 i = 2, 3 (18)
X i = α i • X 4 i = 5, 6 , 7, 8.
We then obtain the parameters αi which resolve the linear least squares problem, as in the second section of this article (see Table 3). Lastly, we obtain a new approximation given by the following system:

dX 1 dt = k[cos(t • ω) + 1] -θ1 X 1 X 1 (0) = x 0 1 ( 19
)
dX 4 dt = θ2 X 1 -θ3 X 4 X 4 (0) = x 0 4 dX 9 dt = θ4 X 4 -θ5 X 9 X 9 (0) = x 0 i Theoretical Estimated
Error value α i value αi percent 2 0.750 × 10 -3 0.727 × 10 -3

3.0 3 0.500 × 10 -3 0.498 × 10 -3 0.4 5 0.599 × 10 -2 0.568 × 10 -2 5.1 6 0.399 × 10 -2 0.384 × 10 -2 4.0 7 0.300 × 10 -2 0.288 × 10 -2 4.0 8 0.200 × 10 -2 0.194 × 10 -2 2.9

Table 3: Original parameters (α i ) and numerical approximations ( αi ), solutions of the least square linear problems for approximating the coefficients in (18). The parameters considered for the simulation of the original dynamics in Network N1 are specified in Table 1. 3). Thick light line: original system (1); dots: supposed data with white noise (17); solid thin line: calibrated system (19) with the parameters in Table 2 and Table 3. The original parameters are specified in Table 1.

Comparison between DRUM and FBA

The reduced system for the toy Network N1 (depicted in Figure 3) obtained after the DRUM 1 approach is

dX 1 dt = k[cos(t • ω) + 1] -k 21 X 1 X 1 (0) = x 0 1 ( 20
)
dX 4 dt = k 21 X 1 -k 54 X 4 X 4 (0) = x 0 4 dX 9 dt = k 54 X 4 X 9 (0) = x 0 9 ,
(see Appendix D in the Supporting Information for details). Remark that equations ( 16) and ( 20) coincide when µ = 0. However, for nonzero growth rate, omitting µ can imply differences between both systems even in their qualitative behaviors (see Figure 9). On the other hand, the principal hypothesis of Flux Balance Analysis is that metabolic networks reach a steady state, under any external conditions. [START_REF] Orth | What is flux balance analysis?[END_REF][START_REF] Kauffman | Advances in flux balance analysis[END_REF] From the numerical solutions depicted in Figure 5 and Figure 6, we can observe that the dynamics of the metabolic Network N1 do not reach an equilibrium in the delimited period of time. Indeed, the metabolic system of N1 have no equilibrium point, because it is forced by a continuous (non constant) periodic input. This is typically a case where FBA is a rough approximation. To illustrate this fact, we applied FBA to the toy metabolic Network N1 and we compare with our approach. The purpose of FBA is to resolve the algebraic equation

N • V (X) = 0
for the variables X 1 , . . . , X 8 , where X is the vector of metabolite concentrations, N is the stoichiometric matrix and V (X) is the vector of kinetics reactions. We recall that FBA methods omit the dilution due to growth. As a consequence, the concentration of X 9 cannot be estimated from resolving the FBA algebraic equation mentioned above. The results are explained with details in Appendix D.2 of the Supporting Information. (Figure 3). Thick light solid line: graphical behavior of the original system (1); dark solid line: system obtained using the method proposed in this work ( 15)-( 16); square dotted line: system obtained after DRUM (20); thin line: Flux Balance Analysis solution (see Appendix D in the Supporting Information for details). The parameters considered are specified in Table 1. In this figure, for the metabolites in QSS there is no approximation from the framework Drum. For X 9 there is no approximation solution with FBA. The reduction method developed in this work gives the most accurate approximation to the original dynamics.

Reduction of an Autotrophic Microalgae Metabolic Network

The proposed method is applied to the metabolism of autotrophic microalgae. A metabolic network taken from Yang et al. [START_REF] Yang | Energetics and carbon metabolism during growth of microalgal cells under photoautotrophic, mixotrophic and cyclic light-autotrophic/dark-heterotrophic conditions[END_REF] was used to represent the metabolism of autotrophic microalgae. The network has 61 reactions and 59 metabolites. We assume that each enzymatic reaction can be represented by a Michaelis-Menten kinetics. Then, the metabolic system can be linearized around the working mode to end up to the structure proposed in this paper (see Appendix E in the Supporting Information for details).

The inputs in the system are a function representing CO 2 uptake rate and a function fueled by the photon flux density (PFD), which represents the light step in photosynthesis.

Macromolecules such as lipids and carbohydrates are stored during the day and reach relatively high concentration in microalgae. According to our approach, we suppose that these elements are in the slow part of the system (see Theorem 1). The rest of the metabolites are considered to have lower concentration, and then to be in QSS.

We assume there is no metabolite that accumulates as high as macromolecules between the inputs and the metabolites in the slow part, and that this is because the intermediate metabolites are consumed by fast reactions. A closer look at the considered metabolic network shows that the macromolecules (lipids, carbohydrates and chlorophyll) are all produced directly from a metabolite in the fast part and they are not consumed after. According to the generalized reduction presented in Appendix B of the Supporting Information, metabolites in the fast part (put into QSS) are linear combinations of the inputs.

Property 4. Consider a linear metabolic system of n metabolites with several inputs I 1 (t), I 2 (t), . . . , I k (t). Suppose that metabolites X 1 , X 2 , . . . , X m are not consumed by any reaction, while X m+1 , X m+2 , . . . , X n are consumed by fast reactions. Then, applying our reduction strategy, the equations of slow metabolites are equal to linear combinations of the inputs and the term of growth dilution:

dX i dt = k j=1 β i j I j (t) -µX i ∀i = 1, . . . , m,
where β i j = 0 if there is no flux from I j (t) to X i .

Proof. See Appendix B in the Supporting Information.

Derived from Property 4, we can deduce the reduced system composed of the microalgae macromolecules equations:

dX i dt = b i • L(t) + c i • C 0 + n i • N 0 -µ • X i -e i , X i (0) = x 0 i
, where X i is the concentration of a macromolecule, L(t) is the function representing the continuous evolution of light intensity, C 0 and N 0 are the controlled concentration constants of CO 2 and NO 3 , respectively, µ > 0 is the dilution due to growth, b i , c i and n i are nonnegative numbers combining the parameters of the original system (whose kinetics are unknown), as well as the real number e i that also depends on functional equilibrium points of the linearized system (see the Appendix E in the Supporting Information for details). These constants will be re-estimated by our approach.

We consider experimental data described in Lacour et al. [START_REF] Lacour | Diel variations of carbohydrates and neutral lipids in nitrogen-sufficient and nitrogen-starved cyclostat cultures of isochrysis sp[END_REF] The same light-day pattern has been considered as in the experiments. PDF at noon was 1,500 µmol.m -2 .s -1 . The data available only include the concentration of carbohydrate [CAR], neutral lipids [NL] and chlorophyll [CHLO]. From the simplified autotrophic microalgae metabolic network, [START_REF] Yang | Energetics and carbon metabolism during growth of microalgal cells under photoautotrophic, mixotrophic and cyclic light-autotrophic/dark-heterotrophic conditions[END_REF] we assume that Glucose-6-phosphate [G6P], Diacylglycerol [DG] and Glutamate [GLU] are the direct precursors of [CAR], [NL] and [CHLO], respectively. Also, we consider Acetyl-coenzyme A [AcCoA] as the direct precursor of [DG].

Therefore, we calibrated a system of three equations representing X 1 =CAR, X 2 =NL and X 3 =CHLO concentrations:

dX 1 dt = b 1 • L(t) + c 1 • C 0 -µ • X 1 -e 1 , (21)
dX 2 dt = b 2 • L(t) + c 2 • C 0 -µ • X 2 -e 2 , dX 3 dt = b 3 • L(t) + c 3 • C 0 + n 3 • N 0 -µ • X 3 -e 3 .
Notice that in the first two equations of (21) we consider n 1 = 0 and n 2 = 0, because in the simplified network there is no flux from the NO 3 input to CAR nor DG (see Proposition 33 in the Appendix B of this article). We use the minimization method described in the second section of this article and the fminsearch tool in Scilab, to estimate the parameters of system (21). The fminsearch tool in Scilab is based on the Nelder-Mead algorithm (http://www.scilab.org). The simulation results for the reduced system are presented in Figure 10 and its calibrated parameters in Table 4

s -1 c 2 0.100 × 10 -3 s -1 c 3 0.100 × 10 -3 s -1 n 3 0.100 × 10 -2 s -1 µ 0.100 × 10 -6 s -1
Table 4: Parameters obtained after the calibration of system (21) with experimental data obtained from Lacour et al. [START_REF] Lacour | Diel variations of carbohydrates and neutral lipids in nitrogen-sufficient and nitrogen-starved cyclostat cultures of isochrysis sp[END_REF] The initial conditions of the system are x 0 1 = 0.35 gC/gC, x 0 2 = 0.2 gC/gC and x 0 3 = 0.0125 gC/gC.

Discussion and Conclusion

In this work we have reduced a metabolic network to a small number of macroscopic reactions, eliminating internal metabolites under the QSSA.

We have shown that this computation is possible for a non-homogeneous linear system of n equations, accounting for a continuous input. Moreover, our approach is suitable for any metabolic network whose kinetics can be locally linearized.

The method proposed in the first section of this article incorporates the dilution factor with the coefficients of reaction rates, allowing to check the stability conditions for obtaining an accurate approximation. On the other hand, when the dilution factor is omitted for resolving a system in QSS, as most approaches do, the reduced models obtained lose accuracy, as shown for the toy network N1 in this article.

We emphasize that, writing the dynamical system into the proper canonical form of singularly perturbed systems is crucial to rigorously separate the time scales. QSSA applied to non canonical forms can lead to erroneous conclusions, [START_REF] Boie | Effects of quasi-steady-state reduction on biophysical models with oscillations[END_REF] as well as when the conditions of Tikhonov's Theorem are not fulfilled. [START_REF] Tikhonov | Differential Equations[END_REF] Additionally, for a system with no flux traps, we proved that the concentration of metabolites in QSS is very low in comparison to the accumulative metabolites (Theorem 1) as a consequence of the two time-scales. Using the reverse reasoning, we can detect slow reactions associated to accumulative metabolites and fast reactions related to metabolites with low concentration.

The Theory of Singularly Perturbed Systems is a tool which has already been used to justify the Quasi Steady State Assumption in metabolic networks. [START_REF] Radulescu | Robust simplifications of multiscale biochemical networks[END_REF][START_REF] Waldherr | Dynamic optimization of metabolic networks coupled with gene expression[END_REF] However, our approach combines this tool with others and eventually it leads to an alternative approach with a different reduction strategy. Note also that we use the concepts of trap [START_REF] Jacquez | Qualitative theory of compartmental systems Siam Review[END_REF] and flux trap, which were, to the best of our knowledge, never considered previously for Quasi Steady State Reductions of metabolic networks. As future work, more generalized systems including nonlinear kinetics can be studied with the framework developed in this article. Figure 2. Single reactant-single product hypothesis: fast reaction k ji /ε just consumes one metabolite (X i ) and produces another (X j ). This conditions leads to a linear system of equations. Figure 3. Metabolic Network N1. Metabolites X 1 , X 4 and X 9 can accumulate and have slow dynamics, while the others are in Quasi Steady State. Reactions from metabolites in QSS are faster than those from metabolites with slow dynamics. The input is

List of Figures

I(t) = k[cos(t • ω) + 1]. Figure 4.
Behavior of the slow-fast system in standard form (see Equation ( 5)-( 6)) that is obtained after the variable rescaling of Network N1 (Figure 3). Notice that the all the metabolite concentrations have the same order of magnitude, as a consequence of defining Y i = X i /ε for the metabolites in the fast part. The parameters considered are specified in Table 1.

Figure 5. Dynamics of Network N1 (Figure 3). Thick line: numerical solution of the original system (1); thin line: reduced system obtained using the method developed in this work ( 15)-( 16). The parameters considered for the simulation are specified in Table 1. Figure 6. Zoom on the dynamics of metabolites in QSS (see also Figure 5). Figure 7. Calibrated system of toy Network N1 (Figure 3). Thick light line: original system (1); dots: supposed data with white noise (17); solid thin line: calibrated system (19) with the parameters in Table 2 and Table 3. The original parameters are specified in Table 1. Figure 8. Zoom on the dynamics of metabolites in QSS (see also Figure 7). Figure 9. Comparison of different approaches to approximate the dynamics of Network N1 (Figure 3). Thick light line: graphical behavior of the original system (1); dark line: system obtained using the method proposed in this work ( 15)-( 16); square dotted line: system obtained after DRUM (20); thin line: Flux Balance Analysis solution (see Appendix D in the Supporting Information for details). In this figure, for the metabolites in QSS there is no approximation from the framework Drum. For X 9 there is no approximation solution with FBA. The parameters considered are specified in Table 1. Figure 10. Dynamical behavior of Carbohydrates, Lipids and Chlorophyll in the metabolism of autotrophic microalgae. Experimental data was obtained from Lacour et al. [START_REF] Lacour | Diel variations of carbohydrates and neutral lipids in nitrogen-sufficient and nitrogen-starved cyclostat cultures of isochrysis sp[END_REF] The reduced model (21) was calibrated using the data. The parameters obtained after the calibrations are in Table 4. Light pattern (PFD) is also represented. Figure 11. Possible scenarios where k 21 = 0 in a system with two metabolites and one output. Both cases represent a flux trap in X 1 . Figure 12. Boundary layer correction applied to the toy Network N1. The boundary layer term is only added to the approximations of the metabolites in QSS. Therefore, dynamics of the slow reduced system ( 16) are the same as in Figure 5. Thick light line: original system (1); solid line: approximation obtained after the approach proposed in this work (15); dashed line: approximation with the boundary layer correction (34). The parameters considered are in Table 1. Notice that the approximation without boundary layer correction (solid line) does not have the same initial condition as the original system (thick light line), while the system with boundary layer correction (dashed line) accurately approximates the original system in the first fast transient interval. Figure 14. Subnetworks SN1 (above) and SN2 (below) of N1. The internal metabolites X 2 , X 3 , X 5 , X 6 , X 7 and X 8 are assumed to be in Quasi Steady State. Metabolites X 1 , X 4 and X 9 are supposed to be external to the subsystems in QSS. Figure 15. Metabolic Network N2. Arrows between metabolites represent enzymatic reactions catalyzed by an enzyme e ji , with substrate X i , product X j and product formation rate k ji or k ji /ε, respectively. The input

I(t) = k[cos(t • ω) + 1] is a periodic continuous function.
Figure 16. Dynamics of the network with enzymatic reactions N2 (Figure 15). Thick light solid line: nonlinear system describing the Michaelis-Menten reactions of N2 (41); thin dark solid line: linearized system (44); dashed line: value of the functional equilibrium points (X * i ). The parameters considered for the simulation are stated in Table 5. Figure 17. Zoom on dynamics of the network with enzymatic reactions N2 (see also Figure 16).

A Complements for the Proof of Theorem 1

For the following demonstrations recall that all the terms k ij = O(1) have the same order, µ ≤ O(1), and that ε is a very small positive number (see the first section of this article).

Lemma 1. Consider a constant matrix A = (a ij ) of dimension n × n, such that O(a ij ) = 1 when ε → 0 for every i, j. Suppose that A is nonsingular and let εµ > 0. Then

det(A -εµ • I) = (-1) n • O(εµ).
Proof. If λ (1) = 0, λ [START_REF] Nazaret | An old paper revisited: "A mathematical model of carbohydrate energy metabolism. Interaction between glycolysis, the Krebs cycle and the H-transporting shuttles at varying ATPases load[END_REF] , . . . , λ (n) are the eigenvalues of A, we have that

det(A -λ • I) = (-1) n • (λ -λ (1) ) • • • (λ -λ (n) ) = (-1) n • (λ)(λ -λ (2) ) • • • (λ -λ (n) ).
Substituting λ by εµ in the formula above, we obtain the desired result.

Lemma 2. Suppose that M is a column diagonal dominant matrix of size n × n, such that det(M ) = 0. If every off-diagonal entry of M is nonnegative, then all the cofactors of M have the same sign equal to (-1) n-1 and sgn(det(M )) = (-1) n .

Proof. Since -M is nonsingular and column diagonal dominant, by the Theorem of Gershgorin, -M is a positive stable matrix. [START_REF] Horn | Matrix Analysis[END_REF] Then its inverse matrix is nonnegative 27 (i.e. each entry of (-M ) -1 is nonnegative). But

-((-M ) -1 ) = (M ) -1 = 1 det(M ) • C ≤ 0,
where

C =      C 11 C 12 . . . C 1n C 21 C 22 . . . C 2n . . . . . . . . . C n1 C n2 . . . C nn      T ,
is the transpose matrix of cofactors of M . [START_REF] Lancaster | The theory of matrices: with applications[END_REF] Then

C ij det(M ) ≤ 0 ∀i, j = 1, . . . , n,
which implies that all the cofactors C ij = (-1) i+j M ij , with M ij the minor of M obtained from removing the i-th row and the j-th column, [START_REF] Lancaster | The theory of matrices: with applications[END_REF] have the same sign. Moreover, since all the principal minors of -M are positive, [START_REF] Berman | Nonnegative matrix in the mathematical science[END_REF] then det(-M ) > 0. We conclude that

sgn(C ij ) = (-1) n-1
and that det(M ) = (-1) n det(-M ) is negative if n is odd and positive is n is even.

Proposition 1. Consider a linear system with n metabolites and no traps.

We suppose that there is an output from the n-th metabolite to the exterior of the system. If M n is the (Jacobian) matrix of the system, then

det(M n ) = (-1) n •O(k n ij ).
Proof. The matrix of the system is written as

M n =    - n i=2 k i1 -k * 1 -εµ • • • k 1n . . . . . . . . . k n1 • • • - n-1 i=1 k in -k n+1,n -εµ    ,
where k * i ≥ 0. Notice that an output from the i-th metabolite is equivalent to k * i > 0. Here, without loss of generality, we begin supposing that the n-th has an output. Then k n+1,n > 0.

We prove the proposition by induction over n. For n = 2, consider the matrix

M 2 = -k 21 -εµ k 12 k 21 -k 12 -k 32 -εµ (22)
of a system with two metabolites and one output. The determinant of

M 2 is det(M 2 ) = k 21 (k 32 + εµ) + εµ(k 12 + k 32 + εµ). If k 21 •k 32 = 0, then det(M 2 ) = O(k 2 ij ).
We examine in which cases k 21 •k 32 = 0. If k 32 = 0, the system has no output, contrary to our hypothesis. On the other hand, k 21 = 0 implies that X 1 is in a trap (see Figure 11). We conclude that det(M 2 ) = O(k 2 ij ). The case in dimension n = 2 with more than one output is verified immediately.

We make the following induction hypothesis: consider a linear system of n -1 metabolites with no traps and one output at least. If M n-1 is the matrix of this system, then det

(M n-1 ) = (-1) n-1 • O(k n-1 ij
). Now we prove the case of a network with n metabolites. We take into account that all the cofactors C ij of M n have the same sign, as claimed by Lemma 2. It holds

det(M n ) = -k n+1,n C nn (23) + n-1 j=1 k jn C jn - n-1 i=1 k in + εµ C nn ,
where C jn = (-1) j+n (M n ) jn are cofactors of M n . [START_REF] Lancaster | The theory of matrices: with applications[END_REF] Suppose that k ni = 0 and k * i = 0 for every i ∈ {1, . . . , n -1}. Then X n is isolated and the rest of metabolites {X 1 , . . . , X n-1 } form a trap. Hence, k ni > 0 or k * i > 0 for some i ∈ {1, . . . , n-1} and we can apply the hypothesis of induction to deduce that

C nn = (-1) n-1 • O(k n-1 ).
On the other hand, the term in the squared brackets in ( 23) is the determinant of the matrix (M n + k n+1,n • δ nn ), where δ nn is a matrix of size n × n with zero at every entry, except for in the entry nn which is equal to 1

If k * i = 0 for every i = 1, . . . , n -1, then

det(M n + k n+1,n • δ nn ) = (-1) n • O(εµ)
according to Lemma 1 and the statement of Proposition 2 is proved. In other case, suppose k * ,n-1 > 0 without loss of generality. Hence, if we develop the determinant of (M n + k n+1,n • δ nn ) by the n -1-th column and we substitute in (23), we have

det(M n ) = -k n+1,n C nn -k * ,n-1 (M n + k n+1,n • δ nn ) n-1,n-1 (24)
+ n j=1

j =n-1 k j,n-1 (M n + k n+1,n • δ nn ) j,n-1 - n-2 i=1 k i,n-1 + εµ (M n + k n+1,n • δ nn ) n-1,n-1 ,
where (M n + k n+1,n • δ nn ) j,n-1 are cofactors of (M n + k n+1,n • δ nn ). Moreover, the matrix (M n + k n+1,n • δ nn ) satisfies the conditions of Lemma 2. Then, all its cofactors have the same sign. Particularly, sgn((M n + k n+1,n • δ nn ) n-1,n-1 ) = (-1) n-1 , and then

sgn(-k n+1,n C nn ) = sgn(-k * ,n-1 (M n + k n+1,n • δ nn )).
Once again, the term in square brackets in ( 24) is equal to det

(M n + k n+1,n • δ nn + k * ,n-1 • δ n-1,n-1
). We proceed as for det(M n + k n+1,n • δ nn ) to extract the following term

-k * ,n-2 (M n + k n+1,n • δ nn + k * ,n-1 • δ n-1,n-1 ) n-2,n-2
which has the same sign as -k n+1,n C nn . In n steps, we arrive to an expression of the determinant where all the terms have the same sign and one term is the determinant of a matrix whose entries by column sum -εµ. That is to say, if we define

M i := (M n + k n+1,n • δ nn + n-i j=1 k * ,n-j δ n-j,n-j ),
for every i = 2, . . . , n, where we define

0 j=1 k * ,n-j δ n-j,n-j = 0. Then det(M n ) = -k n+1,n C nn - n i=2 k * ,i-1 ( M i ) i-1,i-1 + (-1) n • O(εµ) , with ( M i ) i-1,i-1 is a principal minor of M i and the term in square brackets represents det       - n i=2 k i1 -εµ k 12 . . . k 1n k 21 - n i=1 i =2 k i2 -εµ . . . k 2n . . . . . . . . . . . . k n1 k n2 . . . - n-1 i=1 k in -εµ       , according to Lemma 1. Moreover, sgn(-k n+1,n C nn ) = sgn(-k * ,i-1 ( M i ) i-1,i-1 ) = (-1) n ,
for every i = 2, . . . , n, as a consequence of Lemma 2. Therefore, we conclude

det(M n ) = (-1) n • O(k n ij ).
Recall that in our model (Figure 1) we only suppose that there is not flux trap. For this reason, we analyze the determinant of the matrix associated to a system with traps. For instance, with the matrix M 2 defined in ( 22), if the system has a trap, k 21 = 0 and its determinant has order O(εµ).

In general, we can expect that a system with a trap has determinant with order εµ. This happens because a trap implies a block of zeros in the matrix. Indeed, remember that the j-th column of the matrix system represents the reactions whose origin is the metabolite X j . Then, if X j is in a trap, k ij = 0 for every i with X i out of the trap.

In the presence of a trap, the matrix of the system is reducible. [START_REF] Lancaster | The theory of matrices: with applications[END_REF] That is to say, after the same number of interchanges of rows than columns, the matrix of a system with a trap can be transformed in a square block triangular matrix (keeping the dominant diagonal structure):

M n = M 0 * T , (25)
where M and T are square matrices that correspond to the metabolites which are not in a trap and the metabolites which are in a trap, respectively. If C ij is a cofactor of M n and det(T) has order εµ, then the coefficients

k 21 C 1i det(M ) • det(T)
can be affected by a factor of order (εµ) -1 (see the proof of Theorem 1).

However, it is possible to see that when there is a trap which is not reached by the flux, then the determinant of the block corresponding to the trap is also a factor of the cofactors C 1i , where X i is not in that trap. This is why we distinguish a flux trap of a simple trap (see Definitions 2 and 3), by determining if the input of interest reaches them or not.

A matrix with a trap that is not a flux trap has the following form:

M n =   M r×r [C 1 ] r×s 0 r×p 0 s×r [C 2 ] s×s 0 s×p 0 p×r [C 3 ] p×s T p×p   , r + s + p = n, (26)
where T represents the trap not reached by a specific flux, [C i ] represents columns of metabolites that connect the trap with the rest of the system, but which are not nourished by the flux not even by the trap. M is again the block corresponding to the rest of metabolites, including those with an input or an output. Since the matrix in ( 26) is square block triangular, its determinant is the product of the determinants of the diagonal blocks. [START_REF] Bernstein | Matrix mathematics: Theory, facts, and formulas with application to linear systems theory[END_REF] Then, this determinant has a factor det(T) = (-1) p •O(εµ), as claimed by Lemma 1.

Furthermore, the cofactors C 1j have the factor det(T) for every j = 1, . . . , r, because the sub matrix (M n ) 1n is also square block triangular. Then, when dividing by the determinant of M n , this factor is neutralized. In this way we rule out having a large coefficient of order (εµ) -1 for estimating the concentration of X i , 1 ≤ i ≤ r, for every metabolite that is not in a trap. Whereas the minors M 1j = 0 for j > r, as a consequence of the block of zeros in the lower left corner. In fact, the metabolites not reached by the flux are not related to this either in QSS. Therefore we have proved the following proposition.

Proposition 2. Consider the matrix M n in its triangular form (25), such that the square block T corresponds to metabolites in a trap and M to metabolites not in the trap. Then

det(M n ) = det(M ) • det(T),
Moreover, if the trap is not reached by a specific flux, then M n has the form (26) and its minors satisfy

(M n ) 1j = (M ) 1j • det(C 2 ) • det(T ) ∀j = 1, . . . , r,
with (M ) 1j a minor of M , and

(M n ) 1j = 0 ∀j = r + 1, . . . , n. Corollary 1. If M n has a trap, then det(M n ) = ±O(εµ).
Proof. The square block T is equal to a singular matrix minus εµ • I. Then, by Lemma 1, its determinant has order ±O(εµ).

Now we consider some minors of the matrix M n and we study the order of their determinants, as required for the proof of Theorem 1. Proposition 3. Let us suppose that M n represents a system with no traps. Moreover, assume a flux from X 1 to X n . Consider the minor of M n resulting from removing the first line and the n-th column:

(M n ) 1n = det          k 21 -( n i=1 i =2 k i2 + εµ) . . . k 2,n-1 k 31 k 32 . . . k 3,n-1 . . . . . . . . . . . . k n-1,1 k n-1,2 . . . -( n i=1 i =n-1 k i,n-1 + εµ) k n1 k n2 . . . k n,n-1          (27)
Then 0 < (M n ) 1n = O(k n-1 ij ).
Proof. The demonstration is by induction over the squared matrix size. For the case of a minor with dimension two we have:

det k 21 -( 3
i=1 i =2 k i2 + εµ) k 31 k 32 = k 21 k 32 + k 31 ( 3
i=1 i =2 k i2 + εµ) = O(k 2 ij ),
since there is a flux from X 1 to X 3 and no traps. We then suppose the validity of this lemma for a minor of dimension up to n -2 (induction hypothesis).

If we develop the determinant (M n ) 1n by the first column, we verify that the minor resulting from striking the first column and the x-th row of the matrix in (27) satisfies the hypothesis of this lemma after x -1 changes of columns, for x = 1, 2, . . . , n -2. Hence we apply the induction hypothesis to these minors and we obtain that they are quantities (-1) x-1 • O(k n-2 ij ), where x is the number of the struck row.

Since there is no traps by hypothesis, the minor obtained after omitting the first column and the last row of the matrix in (27) has a column which is strictly diagonal dominant. We can then apply Proposition 1 and conclude that it has order (-1) n-2 • O(k n-2 ij ). Therefore, we conclude that the minor (M n ) 1n is the sum of positive quantities of order O(k n-1 ij ):

For the other minors we obtain a similar result. Indeed, every minor obtained from striking the first row and the x-th column can be transformed in a matrix of the form (M n ) 1n , by n -x changes of rows. Therefore, the following assertion holds.

Corollary 2. When M n has no traps, the minor (M n ) 1x has order (-1) n-x • O(k n-1 ij ), for every x = 1, . . . , n.

B Generalization of the Quasi Steady State Reduction of Linear Metabolic Systems

In this Appendix, we generalized the Quasi Steady State Reduction obtained in the first section of this article (see Property 3). The following reduction is valid for any linear system of metabolic reactions, with any number of subnetworks of fast reactions and all possible reactions between metabolites. Also, the generic network can include a finite number of continuous inputs, entering at any metabolite.

B.1 Generic Linear System with two time-scales

As the main manuscript, we consider linear system with metabolites consumed by slow or fast reactions. Without loss of generality, we suppose that metabolites X 1 , X 2 , . . . , X m are only consumed by slow reactions and that metabolites X m+1 , X m+2 , . . . , X n can be consumed by a fast reaction. Then, we consider the following system of equations:

dX i dt = F i (t, X 1 , . . . , X n , ε) X i (0) = x 0 i , (28)
where

F i := I i (t) + m j=1 j =i k ij • X j + n j=m+1 k ij ε X j -k i • X i for i = 1, . . . , m,
where

k i := n j=1 j =i k ji + µ, and
F i := I i (t) + m j=1 k ij • X j + n j=m+1 j =i k ij ε X j - k i ε • X i for i = m + 1, . . . , n,
where

k i ε := n j=1 j =i k ji ε + µ.
Without loss of generality, we suppose I i (t) ≥ 0.

B.2 Canonical Form of Singularly Perturbed Systems.

Equation ( 28) is a slow-fast system, where the variables X 1 , X 2 , . . . , X m are in the slow part and X m+1 , X m+2 , . . . , X n are in the fast part. Indeed, making the change of variable

Y i = X i ε
for the fast variables, we obtain

dX i dt = G i (t, X 1 , . . . , X m , Y m+1 , . . , Y n , ε) X i (0) = x 0 i , (29)
ε dY i dt = G i (t, X 1 , . . . , X m , Y m+1 , . . . , Y n , ε) Y i (0) = y 0 i , (30)
where

y 0 i := x 0 i /ε, G i :=I i (t) + m j=1 j =i k ij • X j + n j=m+1 k ij Y j -k i • X i
for i = 1, . . . , m, and

G i :=I i (t) + m j=1 k ij • X j + n j=m+1 j =i k ij Y j -k i • Y i for i = m + 1, . . . , n.
Equation (30) evaluated in ε = 0 is equivalent to algebraic equation which can be rewritten in the matrix form

     -k m+1 k m+1,m+2 . . . k m+1,n k m+2,m+1 -k m+2 . . . k m+2,n . . . . . . . . . k n,m+1 k n,m+2 . . . -k n      •      Y m+1 Y m+2 . . . Y n      = -      I m+1 (t) + m j=1 k m+1,j X j I m+2 (t) m j=1 k m+2,j X j . . . I n (t) + m j=1 k nj X j      .
Notice that the matrix above is stable, since it is strictly column diagonally dominant. [START_REF] Horn | Matrix Analysis[END_REF] Define this matrix as

K :=      -k m+1 k m+1,m+2 . . . k m+1,n k m+2,m+1 -k m+2 . . . k m+2,n . . . . . . . . . k n,m+1 k n,m+2 . . . -k n     
.

Then the solution to the algebraic problem is

     Y m+1 Y m+2 . . . Y n      = K -1 •         -I m+1 (t) + k m+1,j • X j -I m+2 (t) + k m+2,j • X j . . . -I n (t) + k nj • X j         . The entries of K -1 are b ij = C ji det K ,
where C ji a cofactor of K. Thus, we write every Y i as a linear combination of the kinetics in the slow-part and the inputs on the fast part (if they exist):

X i ε = Y i = n j =m+1 |b ij | • I j (t) + m j=1 k j ,j • X j .
It follows

X i = ε • n j =m+1 |b ij | • I j (t) + m j=1 k j ,j • X j ,
for i = m + 1, . . . , n.

B.3 Reduction of the Slow-Fast Linear System

Proof. Since K defined in (33) is a square block triangular matrix, its determinant is the product of the determinants of the diagonal blocks. [START_REF] Bernstein | Matrix mathematics: Theory, facts, and formulas with application to linear systems theory[END_REF] Then,

det(K) = det(K ) • det(C 2 ).
For i = 1, . . . , r, the submatrix obtained from deleting the j -th row and the i-th column of K is also a square block triangular matrix. Then, its determinant is

(K) j i = (K ) j i [C 1 ] 0 [C 2 ] = (K ) j i • det(C 2 ) ∀i = 1, . . . , r,
where (K ) j i is a minor of K and [C 1 ] is the matrix C 1 without its first row. On the other hand, for i = r + 1, . . . , (r + s), the minor (K) j i is also the determinant of a square block triangular matrix, i.e.

(K)

j i = K C 1 0 C 2 j i = 0 ∀i = r + 1, . . . , r + s,
as a consequence of the block of zeros below K . We conclude (K) j i = 0 ∀i = r + 1, . . . , (r + s).

Proposition 5 states that the only elements in the sum

n j =m+1 |b ij | • I j (t) + m j=1 k j ,j • X j
which are different from zero are those with index j such that there is a flux from X j to X i (i.e. when |b ij | = 0). Furthermore, notice that in m j=1 k j ,j • X j the parameter k j ,j = 0 if and only if there is a reaction consuming X j and producing X j , where X j is in the slow part and X j is in the fast part of the system.

C Boundary layer correction for the toy model

In this section we illustrate the effects of the boundary layer correction, applied to the toy Network N1 as described in the first section of this article.

According to Tikhonov's Theorem, [START_REF] Kokotović | Singular perturbation methods in control: analysis and design[END_REF] the solution of the original system (1) for the toy Network N1 (in Figure 5) can be approximated with an error of order O(ε). Indeed, in this paper we prove that this approximation is valid in any time interval [0, T ], where the input I(t) is nonnegative and continuous, if we consider the boundary layer correction X.

The details about the boundary layer system ( 9) are in the first section of this article, and the approximation

X i (t) = X i (t) + Xi t ε + O(ε) ∀t ∈ [0, T ], (34)
for every i = 2, . . . , 8, i = 4, is stated in Property 3 at the first section of this article. The boundary layer correction for the toy Network N1 of the third section of this article is illustrated in the following Figure 12 and Figure 13. 16) are the same as in Figure 5. Thick light line: original system (1); solid line: approximation obtained after the approach proposed in this work (15); dashed line: approximation with the boundary layer correction (34). The parameters considered are in Table 1. 12. Notice that the approximation without boundary layer correction (solid line) does not have the same initial condition as the original system (thick light line), while the system with boundary layer correction (dashed line) accurately approximates the original system in the first fast transient interval.

D Classical approach for the stoichiometric modeling.

D.1 The DRUM methodology: recall and application

In this section we give the details for model reduction using DRUM, presented in the third section of this article, is obtained. In the DRUM approach, as in all the approaches targeting the full metabolic network, an additional reaction for biomass synthesis is required. Here, our goal is to illustrate the approach on a reduced system. To apply the Drum approach in a simplified framework, we have therefore assumed that we were studying a situation where Biomass was at a constant concentration, close to the ones typical of photolimited photobioreactors. The formulation for the stoichiometric analysis of a metabolic system with DRUM 1 considers the original system (1) for Network N1 like

dX dt = I(t) + K • V (X) X(0) = x 0 , ( 35
)
where X is the column vector for intracellular metabolites with entries X 1 ,..., X 9 , the initial condition x 0 = (x 0 1 , . . . , x 0 9 ) T , the input

I(t) = (k[cos(t • ω) + 1], 0, . . . , 0) T ,
K is the stoichiometric matrix and V (X) the vector of kinetics reactions.

To be more precise,

K =              ν 21 ν 32 ν 43 ν 54 ν 56 ν 75 ν 68 ν 87 ν 98 X 1 -1 0 0 0 0 0 0 0 0 X 2 1 -1 0 0 0 0 0 0 0 X 3 0 1 -1 0 0 0 0 0 0 X 4 0 0 1 -1 0 0 0 0 0 X 5 0 0 0 1 1 -1 0 0 0 X 6 0 0 0 0 -1 0 1 0 0 X 7 0 0 0 0 0 1 0 -1 0 X 8 0 0 0 0 0 0 -1 1 -1 X 9 0 0 0 0 0 0 0 0 1      k 21 X 1 k 32 ε X 2 - k 23 ε X 3 k 43 ε X 3 k 54 X 4 k 56 ε X 6 k 75 ε X 5 k 68 ε X 8 k 87 ε X 7 k 98 ε X 8                  .
Moreover, in this approach the factor µ is neglected. For the stoichiometric analysis it is then considered the following system

dX dt = K • V (X) X(0) = x 0 ,
and the input vector I(t) is added later.

In line with the DRUM methodology, 1 we consider that metabolites X 2 , X 3 , X 5 , X 6 , X 7 and X 8 are in QSS, while metabolites X 1 , X 4 and X 9 can accumulate and have dynamics. This leads to the division of the metabolism in two subnetworks SN1 and SN2, as shown in Figure 14.

Figure 14: Subnetworks SN1 (above) and SN2 (below) of N1. The internal metabolites X 2 , X 3 , X 5 , X 6 , X 7 and X 8 are assumed to be in Quasi Steady State. Metabolites X 1 , X 4 and X 9 are supposed to be external to the subsystems in QSS.

The next step in the DRUM methodology is to summarize each subnetwork using its EFM. We compute the EFM of subnetworks SN1 and SN2 by Gauss elimination and corroborate the result with Copasi. [START_REF] Hoops | COPASI-a complex pathway simulator[END_REF] We have that the only elementary flux mode of SN1 is

e T 1 = k21 k 32 ε k 43 ε 1 1 1 ,
and its matrix of EFM is E 1 := (e 1 ). Analogously, the EFM of SN2 and its matrix of EFM are, respectively

e T 2 = k54 k 56 ε k 75 ε k 68 ε k 87 ε k 98 ε 1 0 1 0 1 1 , e T 3 = k54 k 56 ε k 75 ε k 68 ε k 87 ε k 98 ε 0 1 1 1 1 0 , E 2 := (e 2 ).
Similar to e 1 , the mode e 2 represents an elementary flux that begins and finishes at the exterior of the system in QSS. In contrast, e 3 is a cycle in QSS, which can be omitted. Under the QSSA, a simple algebraic computation from the EFM of the two subnetworks leads to the dynamics of the main network:

K • V (X) = K 1 • E 1 K 2 • E 2 • β, (36)
where K i is the sub matrix that consists in the column vectors of K corresponding to the reactions in SN i, for i = 1, 2, and β = (β 1 β 2 ) T is a vector with entries satisfying

E 1 • (β 1 ) = ν 21 ν 32 ν 43 , E 2 • (β 2 ) =   ν 54 ν 56 ν 75 ν 68 ν 87 ν 98   ,
where ν ij is the entry of V (X) corresponding to the reaction with rate k ij or k ij /ε, respectively. Then, we compute

K 1 • (e 1 ) =              ν 21 ν 32 ν 43 X 1 -1 0 0 X 2 1 -1 0 X 3 0 1 -1 X 4 0 0 1 X 5 0 0 0 X 6 0 0 0 X 7 0 0 0 X 8 0 0 0 0 0 0 0 0 0 X 2 0 0 0 0 0 0 X 3 0 0 0 0 0 0 X 4 -1 0 0 0 0 0 X 5 1 1 -1 0 0 0 X 6 -1 0 1 0 0 X 7 0 0 1 0 -1 0 X 8 0 0 0 -1 1 -1 X 9 0 0 0 0 0 1             •        1 0 1 0 1 1        =              0 0 0 -1 0 0 0 0 1             
Hence, according to Equation (36),

K • V (X) =              -1 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 1              • β 1 β 2 =              -β 1 0 0 β 1 -β 2 0 0 0 0 β 2             
Substituting in Equation (35), we obtain

dX dt = I(t) +      -β 1 0 0 β 1 -β 2 0 0 0 0 β 2      X(0) = x 0 .
Then, the DRUM approach does not take into account the equations for the metabolites in QSS and deduces the following reduced system for the accumulative metabolites

dX 1 dt = I(t) -β 1 X 1 (0) = x 0 1 (37) dX 4 dt = β 1 -β 2 X 4 (0) = x 0 4 dX 9 dt = β 2 X 9 (0) = x 0 9 .
A crucial step in this method is to choose the entries β 1 and β 2 . This choice is arbitrary according to the DRUM method, 1 but it determines the accuracy of the reduced system. Indeed, we have that β 1 has to satisfy • (β 1 ) =

β 1 β 1 β 1 = ν 21 ν 32 ν 43 ⇒ β 1 = ν 21 = k 21 X 1 β 1 = ν 32 = k 32 ε X 2 -k 23 ε X 3 β 1 = ν 43 = k 43 ε X 3 .
Similarly, for β 2 we have

  1 0 1 0 1 1   • (β 2 ) =   β 2 0 β 2 0 β 2 β 2   =   ν 54 ν 56 ν 75 ν 68 ν 87 ν 98   ⇒ β 2 = ν 54 = k 54 X 4 β 2 = ν 75 = k 75 ε X 5 β 2 = ν 87 = k 87 ε X 7 β 2 = ν 98 = k 98 ε X 8
From the equations above, we have to choose only one definition for β 1 and one for β 2 . The DRUM method 1 does not establish any technique for selecting these elements and it assumes that this choice is arbitrary.

However, in this small example, it is straightforward to make the right choice for β 1 and β 2 . As matter of fact, looking at the reduced model (16) obtained in the third section of this article after Tikhonov's Theorem, we can deduce that choosing

β 1 := ν 21 = k 21 X 1 , (38)
β 2 := ν 54 = k 54 X 4 ,
leads us to an accurate approximation. Besides, when the terms of the metabolites in QSS appears in the reduced system after de DRUM approach, these are considered (and calibrated) as constant parameters. Therefore, the other options for β 1 and β 2 does not result in correct approximations. Finally, substituting (38) in Equation (37), we obtain the reduced model after the DRUM approach:

dX 1 dt = I(t) -k 21 X 1 X 1 (0) = x 0 1 dX 4 dt = k 21 X 1 -k 54 X 4 X 4 (0) = x 0 4 dX 9 dt = k 54 X 4 X 9 (0) = x 0 9 .

D.2 Flux Balance Analysis

In this section, we applied Flux Balance Analysis to the toy model presented in the third section of this article. Because of the periodic forcing and the strong accumulation, this approach turns out to be very inaccurate. The principal hypothesis of FBA is that all the internal metabolites of a system reach a steady state, under any external conditions. [START_REF] Orth | What is flux balance analysis?[END_REF][START_REF] Kauffman | Advances in flux balance analysis[END_REF] For our example, we consider all the metabolites X 1 , . . . , X 9 as internal and we suppose that there is an input from the exterior I(t).

For Flux Balance Analysis, the system of equations for the toy Network N1 is

dX dt = N • V (X), ( 39
)
where X is the vector of metabolite concentrations, N is the stoichiometric matrix and V (X) is the vector of kinetics reactions. We recall that FBA methods omit the dilution due to growth. Then, the factor µ • X does not appear on Equation (39) and we have written as follows: where the input u(t) = I(t) := k[cos(t•ω)+1], e 0 ji is the initial concentration for the enzyme catalyzing the reaction with substrate X i and product X j , k ji (or k ji /ε for the fast reactions) is the product formation rate and K ji is the Michaelis-Menten constant of the same enzymatic reaction. To describe how we proceed to linearized the kinetics of N2, let us consider a non linear system:

dX 1 dt = u(t
ẋ = f (x, u), (42)
where

a ji := e ji k ji K ji (X * i + K ji ) 2 ∈ [0, ∞)
∀i, j ∈ {1, . . . , 9}.

Numerical simulations (Figure 16 and Figure 17) show that the linearized system (44) is a feasible representation of the nonlinear system (41). 15). Thick light solid line: nonlinear system describing the Michaelis-Menten reactions of N2 (41); thin dark solid line: linearized system (44); dashed line: value of the functional equilibrium points (X * i ). The parameters considered for the simulation are stated in Table 5. 41) and (44), depicted in Figure 16 and Figure 17. The Michaelis-Menten constants K ji are considered to be all the same, as well as the initial conditions x 0 i for i = 1, . . . , 9.

Figure 1 :

 1 Figure 1: System of kinetic reactions with n metabolites and slow input. Subnetworks of fast reactions are connected by metabolites consumed at low rates. The metabolites within the subnetworks of fast reactions are in Quasi Steady State.

Figure 3 :

 3 Figure 3: Metabolic Network N1. Metabolites X 1 , X 4 and X 9 can accumulate and have slow dynamics, while the others are in Quasi Steady State. Reactions from metabolites in QSS are faster than those from metabolites with slow dynamics. The input is I(t) = k[cos(t • ω) + 1].

Figure 4 :

 4 Figure4: Behavior of the slow-fast system in standard form (see Equation (5)-(6)) that is obtained after the variable rescaling of Network N1 (Figure3). Notice that the all the metabolite concentrations have the same order of magnitude, as a consequence of defining Y i = X i /ε for the metabolites in the fast part. The parameters considered are specified in Table1.

Figure 5 :

 5 Figure 5: Dynamics of Network N1 (Figure3). Thick line: numerical solution of the original system (1); thin line: reduced system obtained using the method developed in this work (15)-(16). The parameters considered for the simulation are specified in Table1.

Figure 6 :

 6 Figure6: Zoom on the dynamics of metabolites in QSS (see also Figure5).

Figure 7 :

 7 Figure 7: Calibrated system of toy Network N1 (Figure3). Thick light line: original system (1); dots: supposed data with white noise (17); solid thin line: calibrated system(19) with the parameters in Table2and Table3. The original parameters are specified in Table1.

Figure 8 :

 8 Figure8: Zoom on the dynamics of metabolites in QSS (see also Figure7).

Figure 9 :

 9 Figure9: Comparison of different approaches to approximate the dynamics of Network N1 (Figure3). Thick light solid line: graphical behavior of the original system (1); dark solid line: system obtained using the method proposed in this work (15)-(16); square dotted line: system obtained after DRUM(20); thin line: Flux Balance Analysis solution (see Appendix D in the Supporting Information for details). The parameters considered are specified in Table1. In this figure, for the metabolites in QSS there is no approximation from the framework Drum. For X 9 there is no approximation solution with FBA. The reduction method developed in this work gives the most accurate approximation to the original dynamics.

Figure 1 .

 1 Figure1. System of kinetic reactions with n metabolites and slow input. Subnetworks of fast reactions are connected by metabolites consumed at low rates. The metabolites within the subnetworks of fast reactions are in Quasi Steady State. Figure2. Single reactant-single product hypothesis: fast reaction k ji /ε just consumes one metabolite (X i ) and produces another (X j ). This conditions leads to a linear system of equations. Figure3. Metabolic Network N1. Metabolites X 1 , X 4 and X 9 can accumulate and have slow dynamics, while the others are in Quasi Steady State. Reactions from metabolites in QSS are faster than those from metabolites with slow dynamics. The input isI(t) = k[cos(t • ω) + 1]. Figure 4.Behavior of the slow-fast system in standard form (see Equation (5)-(6)) that is obtained after the variable rescaling of Network N1 (Figure3). Notice that the all the metabolite concentrations have the same order of magnitude, as a consequence of defining Y i = X i /ε for the metabolites in the fast part. The parameters considered are specified in Table1.

  Figure1. System of kinetic reactions with n metabolites and slow input. Subnetworks of fast reactions are connected by metabolites consumed at low rates. The metabolites within the subnetworks of fast reactions are in Quasi Steady State. Figure2. Single reactant-single product hypothesis: fast reaction k ji /ε just consumes one metabolite (X i ) and produces another (X j ). This conditions leads to a linear system of equations. Figure3. Metabolic Network N1. Metabolites X 1 , X 4 and X 9 can accumulate and have slow dynamics, while the others are in Quasi Steady State. Reactions from metabolites in QSS are faster than those from metabolites with slow dynamics. The input isI(t) = k[cos(t • ω) + 1]. Figure 4.Behavior of the slow-fast system in standard form (see Equation (5)-(6)) that is obtained after the variable rescaling of Network N1 (Figure3). Notice that the all the metabolite concentrations have the same order of magnitude, as a consequence of defining Y i = X i /ε for the metabolites in the fast part. The parameters considered are specified in Table1.

Figure 13 .

 13 Figure 5. Dynamics of Network N1 (Figure3). Thick line: numerical solution of the original system (1); thin line: reduced system obtained using the method developed in this work(15)-(16). The parameters considered for the simulation are specified in Table1. Figure6. Zoom on the dynamics of metabolites in QSS (see also Figure5). Figure7. Calibrated system of toy Network N1 (Figure3). Thick light line: original system (1); dots: supposed data with white noise(17); solid thin line: calibrated system(19) with the parameters in Table2and Table3. The original parameters are specified in Table1. Figure8. Zoom on the dynamics of metabolites in QSS (see also Figure7). Figure9. Comparison of different approaches to approximate the dynamics of Network N1 (Figure3). Thick light line: graphical behavior of the original system (1); dark line: system obtained using the method proposed in this work (15)-(16); square dotted line: system obtained after DRUM(20); thin line: Flux Balance Analysis solution (see Appendix D in the Supporting Information for details). In this figure, for the metabolites in QSS there is no approximation from the framework Drum. For X 9 there is no approximation solution with FBA. The parameters considered are specified in Table1. Figure10. Dynamical behavior of Carbohydrates, Lipids and Chlorophyll in the metabolism of autotrophic microalgae. Experimental data was obtained from Lacour et al.[START_REF] Lacour | Diel variations of carbohydrates and neutral lipids in nitrogen-sufficient and nitrogen-starved cyclostat cultures of isochrysis sp[END_REF] The reduced model(21) was calibrated using the data. The parameters obtained after the calibrations are in Table4. Light pattern (PFD) is also represented. Figure11. Possible scenarios where k 21 = 0 in a system with two metabolites and one output. Both cases represent a flux trap in X 1 . Figure12. Boundary layer correction applied to the toy Network N1. The boundary layer term is only added to the approximations of the metabolites in QSS. Therefore, dynamics of the slow reduced system (16) are the same as in Figure5. Thick light line: original system (1); solid line: approximation obtained after the approach proposed in this work(15); dashed line: approximation with the boundary layer correction (34). The parameters considered are in Table1. Figure 13. Zoom on the initial fast transient of the simulations in Figure 12.Notice that the approximation without boundary layer correction (solid line) does not have the same initial condition as the original system (thick light line), while the system with boundary layer correction (dashed line) accurately approximates the original system in the first fast transient interval. Figure14. Subnetworks SN1 (above) and SN2 (below) of N1. The internal metabolites X 2 , X 3 , X 5 , X 6 , X 7 and X 8 are assumed to be in Quasi

Figure 11 :

 11 Figure 11: Possible scenarios where k 21 = 0 in a system with two metabolites and one output. Both cases represent a flux trap in X 1 .

Figure 12 :

 12 Figure12: Boundary layer correction applied to the toy Network N1. The boundary layer term is only added to the approximations of the metabolites in QSS. Therefore, dynamics of the slow reduced system(16) are the same as in Figure5. Thick light line: original system (1); solid line: approximation obtained after the approach proposed in this work(15); dashed line: approximation with the boundary layer correction (34). The parameters considered are in Table1.

Figure 13 :

 13 Figure13: Zoom on the initial fast transient the simulations in Figure12. Notice that the approximation without boundary layer correction (solid line) does not have the same initial condition as the original system (thick light line), while the system with boundary layer correction (dashed line) accurately approximates the original system in the first fast transient interval.

Figure 15 :

 15 Figure15: Metabolic Network N2. Arrows between metabolites represent enzymatic reactions catalyzed by an enzyme e ji , with substrate X i , product X j and product formation rate k ji or k ji /ε, respectively. The inputI(t) = k[cos(t • ω) + 1] is a periodic continuous function.

Figure 16 :

 16 Figure16: Dynamics of the network with enzymatic reactions N2 (Figure15). Thick light solid line: nonlinear system describing the Michaelis-Menten reactions of N2 (41); thin dark solid line: linearized system (44); dashed line: value of the functional equilibrium points (X * i ). The parameters considered for the simulation are stated in Table5.

Figure 17 : 8 0. 5 × 3 Table 5 :

 178535 Figure17: Zoom on dynamics of the network with enzymatic reactions N2 (see also Figure16).

Table 2 :

 2 Original parameter (θ i ) and numerical approximation ( θi ) of the parameters for Equation (

	i Theoretical Estimated	Error
		value θ i	value θi	percent
	1	0.12 × 10 -1	0.114 × 10 -1	04.7
	2	0.10 × 10 -1	0.085 × 10 -1	15.2
	3	0.22 × 10 -1	0.186 × 10 -1	15.3
	4	0.20 × 10 -1	0.174 × 10 -1	12.8
	5	0.20 × 10 -2	0.179 × 10 -2	10.5

Table 4 .

 4 Light pattern (PFD) is also represented.

	.

Figure 10: Dynamical behavior of Carbohydrates, Lipids and Chlorophyll in the metabolism of autotrophic microalgae. Experimental data was obtained from Lacour et al.

[START_REF] Lacour | Diel variations of carbohydrates and neutral lipids in nitrogen-sufficient and nitrogen-starved cyclostat cultures of isochrysis sp[END_REF]

The reduced model

(21)

was calibrated using the data. The parameters obtained after the calibrations are in

X 2 = ε • k 21 • k 3 d 1 X 1 , X 3 = ε • k 21 • k 32 d 1 X 1 ,(15)X 5 = ε • k 54 • k 6 • k 7 • k 8 d 2 X 4 , X 6 = ε • k 54 • k 68 • k 75 • k 87 d 2 X 4 , X 7 = ε • k 54 • k 75 • k 6 • k 8 d 2 X 4 , X 8 = ε • k 54 • k 75 • k 87 • k 6 d 2 X 4 ,whered 1 = k 3 • k 2 -k 23 • k 32 and d 2 = k 5 • k 6 • k 7 • k 8 -k 56 • k 68 • k 75 • k 87 .Subsequently, we achieve a reduced system by substituting the expres-

X i = αi • X 1 i = 2, 3, X i = αi • X 4 i = 5, 6, 7, 8.

< (M n )

1n =k 21 • O(k n-2 ij ) + • • • + (-1) x+1 (-1) x-1 k x1 • O(k n-2 ij ) + . . . + (-1) n (-1) n-2 k n1 • O(k n-2 ij ) =O(k n-1 ij ).

= G i (t, X

, . . . , X m , Y m+1 , . . . , Y n , 0) ∀i = m + 1, . . . , n,

Funding

Claudia López Zazueta was supported by the National Council of Science and Technology of Mexico through the program conacyt-secretaría de energía-sustentabilidad energética 2015 (scholarship number: 427010). We acknowledge the support of the Investissements d'Avenir Bio-informatique program for project RESET (ANR-11-BINF-0005) and of the Inria Project Lab Algae in silico.

To deduce the reduced system after Tikhonov's Theorem, define G i (t, X 1 , . . . , X m , 0) := G i (t, X 1 , . . . , X m , Y m+1 , . . . , Y n , 0) =

Thus, the QSSR is

for all i = 1, . . . , m. For the variables in QSS,

for all i = m + 1, . . . , n.

The following proposition is a consequence of Tikhonov's Theorem. 10 Proposition 4 (Tikhonov's Theorem). Let (X i ) be the solution of (28) and (X i ) defined by (31)-(32). If there is a solution for (31), then

and there exists 0 ≤ T 0 such that

Now we see that |b ij | = 0 if there is not a flux from X j to X i .

Proposition 5. Let K be the Jacobian matrix of the fast eq. (30)and consider metabolite X j . With out loss of generality suppose 1 ≤ j ≤ r and rewrite this matrix as

where K is the matrix representing the metabolites with a flux from X j (including X j ) and C 2 is the square block corresponding to metabolites not reached by any flux from X j . Then,

Furthermore, its minors satisfy

with (K ) j i a minor of K , and

The purpose of FBA is to resolve the algebraic system

for the variables X 1 , . . . , X 8 . In our example, the algebraic system is linear and we can compute the solution as follows:

where

E Local Linearization for an Enzymatic System

In this section we take into consideration a network similar to N1 (Figure 3), but composed of nonlinear enzymatic reactions, whose substrates and products are the metabolites X 1 , . . . , X n . Then, we show how to locally linearized its metabolic system around a functional point. The ODE system of Michaelis-Menten enzymatic reactions for Network N2 (Figure 15) is and a steady state x * associated to a constant input u * . Therefore

Now we define a linearized system associated to (42) around (x * , u * ):

In our specific example, we take

Then, we calculate the equilibrium point of Equation ( 41) when the input is equal to this constant value. In other words, we obtain the equilibrium points for the variables of Equation ( 41), when u(t) = u * , as

Hence, according to Equation (43), we obtain the following linearized system associated to Equation (41):