# Analytical Reduction of Nonlinear Metabolic Networks Accounting for Dynamics in Enzymatic Reactions

Claudia López Zazueta, Olivier Bernard, Jean-Luc Gouzé

► To cite this version:

Claudia López Zazueta, Olivier Bernard, Jean-Luc Gouzé. Analytical Reduction of Nonlinear Metabolic Networks Accounting for Dynamics in Enzymatic Reactions. Complexity, 2018, 2018, pp.1-22. 10.1155/2018/2342650. hal-01872615

HAL Id: hal-01872615

https://inria.hal.science/hal-01872615

Submitted on 12 Sep 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Analytical Reduction of Nonlinear Metabolic
Networks Accounting for Dynamics in
Enzymatic Reactions

Claudia López Zazueta, Olivier Bernard and Jean-Luc Gouzé

Université Côte d'Azur, Inria, INRA, CNRS, UPMC Univ Paris 06, Biocore
Team, Sophia Antipolis, France

Abstract

Metabolic modeling has been particularly efficient to understand
the conditions affecting the metabolism of an organism. But so far,
metabolic models have mainly considered static situations, assuming
balanced growth. Some organisms are always far from equilibrium
and metabolic modeling must account for their dynamics. This leads
to high dimensional models were metabolic fluxes are no more con-
stant but vary depending on the intracellular concentrations. Such
metabolic models must be reduced and simplified so that they can be
calibrated and analyzed. Reducing these models of large dimension
down to a model of smaller dimension is very challenging, specially,
when dealing with non linear metabolic rates. Here, we propose a
rigorous approach to reduce metabolic models using Quasi Steady
State Reduction based on Tikhonov's Theorem, with characterized
and bounded reduction error. We assume that the metabolic net-
work can be represented with Michaelis-Menten enzymatic reactions,
with two time scales in the reactions. In this simplest approach, some
metabolites can accumulate. We consider the case with a continuous

Correspondence should be addressed to Claudia López Zazueta at claudia.lopez-zazueta@inria.fr
(slowly) varying input in the model, such as light for microalgae, so
that the system is never at steady state. Furthermore, our analysis
proves that the metabolites which can accumulate reach higher con-
centrations (by one order of magnitude) than the fast metabolites. A
simple example illustrates our approach and the resulting accuracy of
the reduction method.

**Introduction**

Metabolic models have considerably helped in understanding the metabolism
of an organism, and enhancing its production capability. These models are
based on simplified metabolic networks, and generally include several hun-
dreds of reactions associated to many metabolic compounds. For example,
metabolic models to better understand the production of triacylglycerols
and carbohydrates from microalgae (both compounds can then be turned
into biofuel) [1] use between 56 and 2190 reactions and between 46 and 1862
metabolites, depending on authors and studies. In order to manage the large
dimension of these models, some simplifying assumptions are generally nec-
essary.

The most classical hypothesis is balanced growth, i.e. global Steady State
Assumption (SSA). This means that the derivatives, with respect to time, of
all variables are put to zero. For instance, Flux Balance Analysis (FBA) [2]
or Macroscopic Bioreaction Models (MBM) [3] 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.

Yet, metabolisms of microalgae and cyanobacteria are directly related
to solar light providing the energy for incorporating CO₂ through Calvin
cycle. Periodic fluctuation of light induces unstationarity, and permanent
accumulation and reuse of metabolites (specially lipids and carbohydrates).
Therefore, such metabolisms are never at steady state, and the classical ap-
proaches based on balanced growth hypothesis cannot be used.

Here we propose a rigorous mathematical approach to reduce the dimen-
sion of a dynamical metabolic system, in order to analyze its behavior and
calibrate it. The reduction that we propose allows to characterize the approx-
imation error and it is appropriate to model based control strategies. The
idea is to keep some dynamical components of the model, that are necessary
specially when dealing with microalgae and cyanobacteria.
A first attempt in this direction was carried out with the Dynamic Reduction of Unbalanced Metabolism (DRUM) method [4]. DRUM considers subnetworks in Quasi Steady State (QSS), which are interconnected by metabolites that can accumulate. Then, Elementary Flux Modes (EFM) are computed in each subnetwork to reduce them using Quasi Steady State Assumption (QSSA). As a result, the dynamics of accumulative metabolites form a reduced system of Ordinary Differential Equations (ODE). It provided sound results, specially to describe accumulation of lipids and carbohydrates in microalgae. However, as almost all the methods, it also relies on a series of assumptions whose mathematical bases have not been rigorously established [5]. Beyond QSSA which is not rigorously defined from a mathematical viewpoint, these approaches also neglect intracellular dilution due to growth.

Models of metabolic networks are non-linear and high dimensional systems, which makes difficult to determine their dynamical behavior and calibrate them. The main objective of our work is to provide mathematical foundations for the reduction of metabolic networks down to low dimensional dynamical models.

Here we study a class of metabolic models of dimension $n$, where the enzymatic reactions rates are represented by Michaelis-Menten reactions. This class of models is the simplest nonlinear one to get accumulation of some intermediate compounds. The objective is to reduce this model accounting for a permanently fluctuating input and rigorously including dilution of the metabolic compounds due to the growth rate. The system is not closed and never reaches a steady state. At the end, we can express a slow dynamical system of small dimension and a fast system as a function of the variables of the slow system. The error in this reduction is then assessed and kept bounded (minimal) controlled.

In Section 1, we introduce the class of models we consider, which is composed of two (general) subnetworks of fast reactions connected by metabolites with slow dynamics. In Section 2, we develop a mathematical model for these metabolic systems.

In Section 3, with proper mathematical hypotheses, after a change of variables for the metabolites with fast dynamics, the system becomes a slow-fast system. The conditions for applying Tikhonov Theorem for singularly perturbed systems are verified and we end-up with a reduced dynamical model and a bound of the approximation error.

In Section 4, we prove that metabolites in QSS have concentration one order of magnitude lower than slow metabolites. Additionally, in Section 5,
we propose an identification algorithm to estimate the parameters of the reduced system from available data.

Finally, we apply our method to a toy metabolic model in Section 6. This simple model is forced by a periodic input and includes standard bricks in metabolic networks: combination of reversible and non-reversible reactions, with chains and cycles.

# 1 Network of Enzymatic Reactions

In this section we present the class of metabolic networks studied all over the paper, which are illustrated in Figure 1. These networks are composed of two subnetworks of fast reactions, which are interconnected by several metabolites with slow rates of consumption. The subnetworks have an arbitrary finite number of metabolites and reactions between them.

These subnetworks are not assumed to have a specific topology. Therefore, they represent a generic case of metabolic networks. The only condition on them is that their metabolites $X_2, \dots, X_{m-1}, X_{m+1}, \dots, X_{n-1}$ are consumed by fast reactions.

The class of systems addressed in this paper can be considered as a simplification or one part of a larger network. However, the results presented through this paper can be extended, allowing the study of more complex systems on the bases of this approach.

<figure><img src="image_2.png" /><figcaption>**Figure 1:** System of enzymatic reactions. An arrow from $X_i$ to $X_j$ represents a Michaelis-Menten reaction catalyzed by an enzyme $e_{ji}$, with substrate $X_i$, product $X_j$ and product formation rate $k_{ji}$ or $\frac{k_{ji}}{\varepsilon}$. Fast reactions are within two subnetworks, which are interconnected by the metabolites $X_1, X_m$ and $\bar{X}_n$. The connector metabolites are consumed by reactions with low rates, while the metabolites in the subnetworks are consumed by fast reactions and in Quasi Steady State.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 2:** Enzymatic reactions between metabolites in QSS depicted in Figure 1. The metabolites inside the subnetworks are substrates or products of fast reactions catalyzed by an enzyme.</figcaption></figure>

## 1.1 Summary of the methodology for reducing Slow-Fast Dynamical Metabolic models

We consider a general class of metabolic models allowing internal accumulation, represented with the network in Figure 1. In order to rigorously reduce this large dimensional model, our objective is to take benefit of the two-time scales and finally rewrite it in the canonical form of singularly perturbed systems. Then the Theorem of Tikhonov [6] can be applied and a reduced system is derived with an accurate bound of the error.

The first challenge is to find the appropriate change of variable for the metabolites with fast dynamics, to end up with a slow-fast system

$$ \frac{dX}{dt} = F(t, X, Y, \eta) \qquad X(0) = x^0 \qquad (1) $$

$$ \frac{dY}{dt} = \frac{1}{\eta}G(t, X, Y, \eta) \qquad Y(0) = y^0, $$

where X is the vector of metabolites with slow dynamics, Y is the vector of metabolites with fast dynamics and η is a very small parameter. Actually, Y results from a rescaling of fast dynamics metabolites $X_{fast}$ in the model: $Y := X_{fast}/\eta$.

When the system is under this general form, we prove some conditions necessary to apply Tikhonov's Theorem [6] and finally we obtain a Quasi Steady State Reduction of system (1):

$$ \frac{dX}{dt} = F(t, X, \bar{Y}, 0) \qquad X(0) = x^0, \qquad (2) $$
where $\bar{Y}$ is a root of the equation

$$0 = G(t, X, \bar{Y}, 0).$$

If $\bar{X}$ is a solution for (2), the Quasi Steady State Approximation $(\bar{X}, \bar{Y})$ to the solution of (1) satisfies

$$X = \bar{X} + \mathcal{O}(\eta),$$

$$Y = \bar{Y} + \mathcal{O}(\eta),$$

after an initial fast transient. In other words, the error of the Quasi Steady State Approximation has order of magnitude $\eta$, which is supposed to be a small positive number. In the manuscript, we show that the reduced system differs from existing approaches, mainly because we do not neglect the metabolite dilution due to cell growth.

The mathematical validity of the Quasi Steady State Reduction (QSSR) for the class of systems considered in this paper (Figure 1) is showed from Section 2 to Section 3.

As a new striking result, this approach allows to prove that the concentration of the metabolites in Quasi Steady State is one order of magnitude lower than the metabolites with slow dynamics, i.e.,

$$\eta \cdot Y = X_{\text{fast}} \leq \mathcal{O}(\eta \cdot X).$$

The conditions under which this assertion holds are given in Section 4.

## 2 Considered Class of Networks

In this section, we describe the systems of the network class considered in this work. Then, in Section 3 and Section 4, we deduce a QSSR and prove some conclusions about the magnitude of metabolite concentrations (see Theorem 1) for these systems.

The results obtained in the following sections can be extended to more complex networks. For instance, considering additional slow reactions or more subnetworks of fast reactions connected by metabolites with slow dynamics.
## 2.1 Notation

Consider the network of $n$ enzymatic reactions depicted in Figure 1 and Figure 2, where an arrow from $X_i$ to $X_j$ represents an enzymatic reaction catalyzed by $e_{ji}$, with substrate $X_i$, product $X_j$ and product formation rate $k_{ji}$ or $k_{ji}/\varepsilon$. Then, every enzymatic reaction can be described with the Michaelis-Menten model (see Appendix A).

However, it is necessary to justify the Quasi Steady State Approximation for the Michaelis-Menten model. For this purpose, many solutions have been presented. For example, this holds if the initial substrate concentration $x_i^0$ is sufficiently large compared with the initial enzyme concentration $e_{ji}^0$ [7], or if the product formation rate $k_{ji}$ is small enough [8].

We suppose that among the product formation rates there are two scales of magnitude. Reactions with large rate are within two subnetworks, which are interconnected by the metabolites $X_1, X_m$ and $X_n$. We suppose that the metabolites connecting the subnetworks are consumed by reactions with low rates.

In this context, we say that a reaction is fast if its rate is large, while a reaction is slow if its rate is low. Moreover, we assume the rates of fast reactions sufficiently larger than those of the slow reactions. Then, we denote fast reactions rates by

$$ \frac{k_{ji}}{\varepsilon} $$

and slow reactions rates by

$$ k_{ji}, $$

where $\varepsilon$ is a small positive number.

Additionally, a continuously varying nonnegative input $I(t)$ (e.g. the CO$_2$ uptake in a microalgae submitted to light/dark cycles) and a growth rate $\mu > 0$, which acts as a dilution factor, are taken into account for the models.

## 2.2 Dynamical Model

According to the standard Quasi Steady State Reduction for Michaelis-Menten enzymatic reactions described in Appendix A, we write the ODE system for the model in Figure 1 as

$$ \frac{dX_i}{dt} = F_i(t, X_1, \dots, X_n, \varepsilon, \mu) \qquad X_i(0) = x_i^0, \qquad (3) $$
where

$$F_1 := I(t) - e_0^{21} k_{21} \frac{X_1}{X_1 + K_{21}} - \mu X_1,$$

$$F_i := \sum_{j \in \{1, m, n\}} e_{ij}^0 k_{ij} \frac{X_j}{X_j + K_{ij}} + \sum_{\substack{j=2 \\ j \neq m}}^{n-1} e_{ij}^0 \frac{k_{ij}}{\varepsilon} \frac{X_j}{X_j + K_{ij}} - \sum_{j=1}^{n} e_{ji}^0 k_{ji} \frac{X_i}{X_i + K_{ji}} - \mu X_i,$$

for $i = m, n$, and

$$F_i := \sum_{j \in \{1, m, n\}} e_{ij}^0 k_{ij} \frac{X_j}{X_j + K_{ij}} + \sum_{\substack{j=2 \\ j \neq m}}^{n-1} e_{ij}^0 \frac{k_{ij}}{\varepsilon} \frac{X_j}{X_j + K_{ij}} - \sum_{j=1}^{n} e_{ji}^0 \frac{k_{ji}}{\varepsilon} \frac{X_i}{X_i + K_{ji}} - \mu X_i,$$

for every $i = 2, \dots, n-1, i \neq m$.

The variable $X_i$ describes the $i$-th metabolite cell concentration, $I(t)$ is a nonnegative continuous function, $\varepsilon$ is a small positive number, $e_{ji}^0, k_{ji}$ and $K_{ji}$ are nonnegative parameters, and $\mu > 0$ is the growth rate. When there is no reaction with substrate $X_i$ and product $X_j$, we define $k_{ji} = 0$, and also $k_{ii} = 0$ for every $i = 1, \dots, n$.

*Note.* In our model we can include first order (linear) reactions. In this case, instead of writing

$$e_{ji}^0 k_{ji} \cdot \frac{X_i}{X_i + K_{ji}} \quad \text{or} \quad e_{ji}^0 \frac{k_{ji}}{\varepsilon} \cdot \frac{X_i}{X_i + K_{ji}}$$

as for enzymatic reactions, we have to write

$$e_{ji}^0 k_{ji} X_i \quad \text{or} \quad e_{ji}^0 \frac{k_{ji}}{\varepsilon} X_i,$$

respectively, in Equation (3). For the sake of simplicity, in this paper we only consider the more general case with Michaelis-Menten reactions.

In line with the QSSR of Michaelis-Menten system, we recall that $e_{ji}^0 k_{ji}$ (or $e_{ji}^0 k_{ji}/\varepsilon$ for the fast reactions) and $K_{ji}$ are parameters related to the enzyme reaction with substrate $X_i$ and product $X_j$. Indeed, $e_{ji}^0$ is the initial enzyme concentration, $k_{ji}$ (or $k_{ji}/\varepsilon$) is the product formation rate and $K_{ji} > 0$ is the specific Michaelis-Menten constant defined as

$$K_{ji} := \frac{k_{-1}^{ji} + k_{ji}}{k_1^{ji}} \quad (4)$$
(see Appendix A).

An important preliminary property that the dynamical system (3) has to obey is that, the concentration $X_i(t)$ has to remain nonnegative over the time if the initial conditions are nonnegative. In our model, this depends on the input $I(t)$. This is stated in the following Property:

**Property 1.** If the initial condition $x_i^0$ is nonnegative for every $i = 1, \dots, n$ and $I(t) \ge 0$ for every $t \in [0, T_1]$, then system (3) is positively invariant in $\overline{\mathbb{R}}_+^n$.

*Proof.* To verify this, we show that system (3) is positively invariant in $\overline{\mathbb{R}}_+^n$ if $I(t)$ is nonnegative over any interval $[0, T_1]$.

Recall that all $K_{ji}$ is supposed to be positive and every parameter, $e_{ji}^0, k_{ji}$, $\mu$, is nonnegative. Then, we have for any $i = 1, \dots, n$,

$$F_i(X_1, \dots, \underset{\text{i-th entry}}{0}, \dots, X_n, t, \varepsilon, \mu) \ge 0$$

if $X_j \ge 0$ for every $j = 1, \dots, n, j \ne i$. Therefore, system (3) is positively invariant in $\overline{\mathbb{R}}_+^n$. $\square$

## 2.3 Parameter Order of magnitude

With our notations, to represent different time-scales in the reactions, we fix $\varepsilon$ a small positive number highlighting the difference between the parameter scale orders. We suppose that the parameters $e_{ji}^0, k_{ji}$ are of standard range, i.e.

$$e_{ji}^0 k_{ji} = \mathcal{O}(1) \quad \text{as} \quad \varepsilon \to 0 \qquad \forall i, j \in \{1, \dots, n\}, \quad (5)$$

where $\mathcal{O}$ denotes the Big-O or Landau symbol. For the definition and some properties of $\mathcal{O}$, we refer to [9].

Also, we suppose that the input $I(t)$ has a magnitude not larger than the slow reactions. In other words,

$$I(t) = \mathcal{O}(1)$$

The rate of growth $\mu$ is considered as a parameter smaller than any reaction rate (a standard hypothesis [10]). Here, we assume

$$\varepsilon\mu = \mathcal{O}(\varepsilon). \qquad (6)$$
## 3 Quasi Steady State Reduction

In this section we propose a rigorous Quasi Steady State Reduction (QSSR) of (3). Its mathematical validity is proved thanks to the Theorem of Tikhonov [6]. In other words, this theorem states that the error of this Quasi Steady State Approximation is bounded by the small parameter $\epsilon$.

We formally define the QSSR after Tikhonov's Theorem, of the metabolic network in Figure 1 and system (3), as the following system of dimension three:

$$ \frac{d\bar{X}_1}{dt} = I(t) - e_{21}^0 k_{21} \left( \frac{\bar{X}_1}{\bar{X}_1 + K_{21}} \right) - \mu \bar{X}_1 \quad (7) $$

$$ \frac{d\bar{X}_m}{dt} = \mathbf{b}_{m-1} \cdot e_{21}^0 k_{21} \frac{e_{m,m-1}^0 k_{m,m-1}}{K_{m,m-1}} \cdot \left( \frac{\bar{X}_1}{\bar{X}_1 + K_{21}} \right) - e_{m+1,m}^0 k_{m+1,m} \left( \frac{\bar{X}_m}{\bar{X}_m + K_{m+1,m}} \right) - \mu \bar{X}_m $$

$$ \frac{d\bar{X}_n}{dt} = \mathbf{b}_{n-1} \cdot e_{m+1,m}^0 k_{m+1,m} \frac{e_{n,n-1}^0 k_{n,n-1}}{K_{n,n-1}} \cdot \left( \frac{\bar{X}_m}{\bar{X}_m + K_{m+1,m}} \right) - \mu \cdot \bar{X}_n, $$

with initial conditions $\bar{X}_1(0) = x_1^0$, $\bar{X}_m(0) = x_m^0$ and $\bar{X}_n(0) = x_n^0$, and for the metabolites in QSS

$$ \bar{X}_i(t) = \varepsilon \cdot \mathbf{b}_i \cdot e_{21}^0 k_{21} \left( \frac{\bar{X}_1(t)}{\bar{X}_1(t) + K_{21}} \right) \quad i = 2, \dots, m-1, \quad (8) $$

$$ \bar{X}_i(t) = \varepsilon \cdot \mathbf{b}_i \cdot e_{m+1,m}^0 k_{m+1,m} \left( \frac{\bar{X}_m(t)}{\bar{X}_m(t) + K_{m+1,m}} \right) \quad i = m+1, \dots, n-1, $$

for every $t \in [0, T_1]$. The definition of the parameters $\mathbf{b}_i$ is given later in this section (see Proposition 1 and its proof).

### 3.1 Slow-Fast System

In order to write system (3) in the canonical form of singularly perturbed systems, we define a change of variable for the fast metabolites by

$$ Y_i := \frac{X_i}{\varepsilon} \quad \forall i = 2, 3, \dots, n-1, i \neq m. \quad (9) $$

Let us set the initial conditions for these new variables as

$$ y_i^0 := x_i^0 / \varepsilon, $$
and growth rate as

$$ \tilde{\mu} = \varepsilon\mu. $$

Therefore, after the change of variables (9), Equation (3) can be rewritten as follows:

$$ \frac{dX_i}{dt} = F_i(t, X_1, \varepsilon Y_2, \dots, \varepsilon Y_{m-1}, X_m, \varepsilon Y_{m+1}, \dots, \varepsilon Y_{n-1}, X_n, \varepsilon, \mu) \quad i = 1, m, n, \quad (10) $$

$$ \frac{dY_i}{dt} = \frac{1}{\varepsilon} F_i(t, X_1, \varepsilon Y_2, \dots, \varepsilon Y_{m-1}, X_m, \varepsilon Y_{m+1}, \dots, \varepsilon Y_{n-1}, X_n, \varepsilon, \tilde{\mu}) \quad i = 2, \dots, n-1, i \neq m. $$

Since $\varepsilon$ is a very small positive number, the dynamics of $Y_i$ are faster than those of $X_i$. Hence, the equations of $X_i$ form the slow part of system (10), while the equations of $Y_i$ constitute its fast part.

The previous Equation (10) is written with further details in the next subsection. The goal is to expose how the Quasi Steady State Reduction is obtained and validated using Tikhonov's Theorem.

## 3.2 Canonical form of singularly perturbed systems

The slow-fast system (10) is in the class of singularly perturbed systems of the exact form:

$$ \frac{dX_1}{dt} = I(t) - e_{21}^0 k_{21} \frac{X_1}{X_1 + K_{21}} - \mu X_1 \quad X_1(0) = x_1^0 \quad (11) $$

$$ \frac{dX_m}{dt} = e_{m,m-1}^0 k_{m,m-1} \frac{Y_{m-1}}{\eta Y_{m-1} + K_{m,m-1}} - e_{m+1,m}^0 k_{m+1,m} \frac{X_m}{X_m + K_{m+1,m}} - \mu X_m \quad X_m(0) = x_m^0 $$

$$ \frac{dX_n}{dt} = e_{n,n-1}^0 k_{n,n-1} \frac{Y_{n-1}}{\eta Y_{n-1} + K_{n,n-1}} - \mu X_n \quad X_n(0) = x_n^0, $$
$$ \frac{dY_2}{dt} = \frac{1}{\eta} \left[ e_{21}^0 k_{21} \frac{X_1}{X_1 + K_{21}} + \left( \sum_{j=3}^{m-1} e_{2j}^0 k_{2j} \frac{Y_j}{\eta Y_j + K_{2j}} \right) - \left( \sum_{i=3}^{m-1} e_{i2}^0 k_{i2} \frac{Y_2}{\eta Y_2 + K_{i2}} \right) - \tilde{\mu} Y_2 \right] \quad (12) $$

∴

$$ \frac{dY_{m-1}}{dt} = \frac{1}{\eta} \left[ \left( \sum_{j=2}^{m-2} e_{m-1,j}^0 k_{m-1,j} \frac{Y_j}{\eta Y_j + K_{m-1,j}} \right) - \left( \sum_{\substack{i=2 \\ i \neq m-1}}^{m} e_{i,m-1}^0 k_{i,m-1} \frac{Y_{m-1}}{\eta Y_{m-1} + K_{i,m-1}} \right) - \tilde{\mu} Y_{m-1} \right] $$

$$ \frac{dY_{m+1}}{dt} = \frac{1}{\eta} \left[ e_{m+1,m}^0 k_{m+1,m} \frac{X_m}{X_m + K_{m+1,m}} + \left( \sum_{j=m+2}^{n-1} e_{m+1,j}^0 k_{m+1,j} \frac{Y_j}{\eta Y_j + K_{m+1,j}} \right) - \left( \sum_{i=m+2}^{n-1} e_{i,m+1}^0 k_{i,m+1} \frac{Y_{m+1}}{\eta Y_{m+1} + K_{i,m+1}} \right) - \tilde{\mu} Y_{m+1} \right] $$

∴

$$ \frac{dY_{n-1}}{dt} = \frac{1}{\eta} \left[ \left( \sum_{j=m+1}^{n-2} e_{n-1,j}^0 k_{n-1,j} \frac{Y_j}{\eta Y_j + K_{n-1,j}} \right) - \left( \sum_{\substack{i=m+1 \\ i \neq n-1}}^{n} e_{i,n-1}^0 k_{i,n-1} \frac{Y_{n-1}}{\eta Y_{n-1} + K_{i,n-1}} \right) - \tilde{\mu} Y_{n-1} \right] $$

with initial conditions $Y_i(0) = y_i^0$ for every $i = 2, \dots, n-1, i \neq m$.

*Note.* The Equation (11)-(12) above is a more detailed expression of (10). Indeed, we obtain system (10) when $\eta$ is substituted for $\varepsilon$ in the equations Equation (11)-(12).

An approximation to the solution of system (11)-(12) can be obtained considering the limit when $\eta \to 0$. Then, dynamics in Equation (12) are considered as fast and the QSSA is applied to the metabolites $Y_i$ for every $i = 2, \dots, n-1, i \neq m$.

Hereafter, we say that Equation (11) is the slow part and Equation (12) the fast part of system (3).
### 3.3 Hypotheses necessary for Quasi Steady State

In the following two subsections, we check the assumptions of Tikhonov's Theorem [6]. First we demonstrate that the system has a single steady state (which is not straightforward for non-linear systems). Then we demonstrate that this steady state is asymptotically stable. Eventually, once all the conditions have been established, in Section 3.5 we present the result of Tikhonov's Theorem.

Consider the following algebraic system of equations, obtained from equating to 0 the terms in square brackets in (12) and substituting $\eta = 0$:

$$0 = e_{ji}^0 k_{21} \frac{X_1}{X_1 + K_{21}} + \left( \sum_{j=3}^{m-1} e_{2j}^0 k_{2j} \frac{Y_j}{K_{2j}} \right) - \left( \sum_{i=3}^{m-1} e_{i2}^0 k_{i2} \frac{Y_2}{K_{i2}} \right) - \tilde{\mu} Y_2 \quad (13)$$

∴

$$0 = \left( \sum_{j=2}^{m-2} e_{m-1,j}^0 k_{m-1,j} \frac{Y_j}{K_{m-1,j}} \right) - \left( \sum_{\substack{i=2 \\ i \neq m-1}}^{m} e_{i,m-1}^0 k_{i,m-1} \frac{Y_{m-1}}{K_{i,m-1}} \right) - \tilde{\mu} Y_{m-1}$$

$$0 = e_{m+1,m}^0 k_{m+1,m} \frac{X_m}{X_m + K_{m+1,m}} + \left( \sum_{j=m+2}^{n-1} e_{m+1,j}^0 k_{m+1,j} \frac{Y_j}{K_{m+1,j}} \right) - \left( \sum_{i=m+2}^{n-1} e_{i,m+1}^0 k_{i,m+1} \frac{Y_{m+1}}{K_{i,m+1}} \right) - \tilde{\mu} Y_{m+1}$$

∴

$$0 = \left( \sum_{j=m+1}^{n-2} e_{n-1,j}^0 k_{n-1,j} \frac{Y_j}{K_{n-1,j}} \right) - \left( \sum_{\substack{i=m+1 \\ i \neq n-1}}^{n} e_{i,n-1}^0 k_{i,n-1} \frac{Y_{n-1}}{K_{i,n-1}} \right) - \tilde{\mu} Y_{n-1}$$

In order to apply Tikhonov's Theorem [6], we have to prove that Equation (13) has an isolated root for any nonnegative constant values $X_1$ and
$X_m$, and that this root is asymptotically stable for the following system:

$$ \frac{dY_2}{dt} = e_{21}^0 k_{21} \frac{X_1}{X_1 + K_{21}} + \left( \sum_{j=3}^{m-1} \frac{e_{2j}^0 k_{2j}}{K_{2j}} Y_j \right) - \left( \sum_{i=3}^{m-1} \frac{e_{i2}^0 k_{i2}}{K_{i2}} + \tilde{\mu} \right) Y_2 \quad (14) $$

$\vdots$

$$ \frac{dY_{m-1}}{dt} = \left( \sum_{j=2}^{m-2} \frac{e_{m-1,j}^0 k_{m-1,j}}{K_{m-1,j}} Y_j \right) - \left( \sum_{\substack{i=2 \\ i \neq m-1}}^{m} \frac{e_{i,m-1}^0 k_{i,m-1}}{K_{i,m-1}} + \tilde{\mu} \right) Y_{m-1} $$

$$ \frac{dY_{m+1}}{dt} = e_{m+1,m}^0 k_{m+1,m} \frac{X_m}{X_m + K_{m+1,m}} + \left( \sum_{j=m+2}^{n-1} \frac{e_{m+1,j}^0 k_{m+1,j}}{K_{m+1,j}} Y_j \right) - \left( \sum_{i=m+2}^{n-1} \frac{e_{i,m+1}^0 k_{i,m+1}}{K_{i,m+1}} + \tilde{\mu} \right) Y_{m+1} $$

$\vdots$

$$ \frac{dY_{n-1}}{dt} = \left( \sum_{j=m+1}^{n-2} \frac{e_{n-1,j}^0 k_{n-1,j}}{K_{n-1,j}} Y_j \right) - \left( \sum_{\substack{i=m+1 \\ i \neq n-1}}^{n} \frac{e_{i,n-1}^0 k_{i,n-1}}{K_{i,n-1}} + \tilde{\mu} \right) Y_{n-1}. $$

The purpose of finding a root of Equation (13) is to write the fast variables $Y_i$ in terms of the slow variables $X_i$. In this case it is possible to find an analytic solution of this algebraic system, because it is a linear equation for the variables $Y_i$. Similarly, the asymptotical stability of this root for system (14) can be verified with the theory of linear systems of ODE.

**Proposition 1.** Consider $X_1$ and $X_m$ as nonnegative constants values. Then system (14) has a single equilibrium point

$$ (\bar{Y}_i)_{i=2,...,n-1, i \neq m} $$

which is globally asymptotically stable. Moreover,

$$ \bar{Y}_i = \mathbf{b}_i \cdot e_{21}^0 k_{21} \left( \frac{X_1}{X_1 + K_{21}} \right) \quad i = 2, \dots, m-1 \quad (15) $$

$$ \bar{Y}_i = \mathbf{b}_i \cdot e_{m+1,m}^0 k_{m+1,m} \left( \frac{X_m}{X_m + K_{m+1,m}} \right) \quad i = m+1, \dots, n-1, $$

where $\mathbf{b}_i \in \mathbb{R}_+$ are nonnegative coefficients.
Proof. First notice that system (14) is a linear system for $Y_i$ under the hypotheses of Proposition 1. Then, we just have to show that its Jacobian matrix is stable, i.e. that all its eigenvalues have negative real part [11].

The Jacobian matrix of (14) is

$$J = \begin{pmatrix} K_1 & 0 \\ 0 & K_2 \end{pmatrix}, \qquad (16)$$

where

$$K_1 = \begin{pmatrix} -\sum_{i=3}^{m-1} \frac{e_{i2}^0 k_{i2}}{K_{i2}} - \tilde{\mu} & \cdots & \frac{e_{2,m-1}^0 k_{2,m-1}}{K_{2,m-1}} \\ \frac{e_{32}^0 k_{32}}{K_{32}} & \cdots & \frac{e_{3,m-1}^0 k_{3,m-1}}{K_{3,m-1}} \\ \vdots & & \vdots \\ \frac{e_{m-1,2}^0 k_{m-1,2}}{K_{m-1,2}} & \cdots & -\sum_{i \neq m-1}^{m} \frac{e_{i,m-1}^0 k_{i,m-1}}{K_{i,m-1}} - \tilde{\mu} \end{pmatrix}, \qquad (17)$$

$$K_2 = \begin{pmatrix} -\sum_{i=m+2}^{n-1} \frac{e_{i,m+1}^0 k_{i,m+1}}{K_{i,m+1}} - \tilde{\mu} & \cdots & \frac{e_{m+1,n-1}^0 k_{m+1,n-1}}{K_{m+1,n-1}} \\ \frac{e_{m+2,m+1}^0 k_{m+2,m+1}}{K_{m+2,m+1}} & \cdots & \frac{e_{m+2,n-1}^0 k_{m+2,n-1}}{K_{m+2,n-1}} \\ \vdots & & \vdots \\ \frac{e_{n-1,m+1}^0 k_{n-1,m+1}}{K_{n-1,m+1}} & \cdots & -\sum_{i \neq n-1}^{n} \frac{e_{i,n-1}^0 k_{i,n-1}}{K_{i,n-1}} - \tilde{\mu} \end{pmatrix}.$$

But $J$ is a strictly column diagonally dominant matrix, because $\tilde{\mu} > 0$. In other words, for every column of the matrix $J$, the sum of the entries out of the diagonal is strictly less than the absolute value of the entry in the diagonal. Hence, by the Theorem of Gershgorin, $J$ is a stable matrix [12].

The matrix form of Equation (13) is

$$\begin{pmatrix} K_1 & 0 \\ 0 & K_2 \end{pmatrix} \cdot \begin{pmatrix} Y_2 \\ \vdots \\ Y_{m-1} \\ Y_{m+1} \\ \vdots \\ Y_{n-1} \end{pmatrix} = - \begin{pmatrix} e_{21}^0 k_{21} \left( \frac{X_1}{X_1+K_{21}} \right) \\ \vdots \\ 0 \\ e_{m+1,m}^0 k_{m+1,m} \left( \frac{X_m}{X_m+K_{m+1,m}} \right) \\ \vdots \\ 0 \end{pmatrix},$$
Then, the solution of the algebraic system (13) is

$$ 
\begin{pmatrix} \bar{Y}_2 \\ \vdots \\ \bar{Y}_{m-1} \end{pmatrix} = (-e_{21}^0 k_{21}) \cdot (K_1)^{-1} \cdot \begin{pmatrix} \frac{X_1}{X_1+K_{21}} \\ \vdots \\ 0 \end{pmatrix}, \qquad (18) \\
\begin{pmatrix} \bar{Y}_{m+1} \\ \vdots \\ \bar{Y}_{n-1} \end{pmatrix} = (-e_{m+1,m}^0 k_{m+1,m}) \cdot (K_2)^{-1} \cdot \begin{pmatrix} \frac{X_m}{X_m+K_{m+1,m}} \\ \vdots \\ 0 \end{pmatrix}.
$$

Therefore, the variables of the solution can be written as

$$ \bar{Y}_i = \mathbf{b}_i \cdot e_{21}^0 k_{21} \left( \frac{X_1}{X_1 + K_{21}} \right) \qquad i = 2, \dots, m-1, $$

$$ \bar{Y}_i = \mathbf{b}_i \cdot e_{m+1,m}^0 k_{m+1,m} \left( \frac{X_m}{X_m + K_{m+1,m}} \right) \qquad i = m+1, \dots, n-1, $$

with $\mathbf{b}_i \in \mathbb{R}$. Moreover, since $K_i$ is strictly column diagonally dominant, by the Theorem of Gershgorin, $K_i$ is a stable matrix [12]. Then, its inverse matrix is nonpositive [13] (i.e. each entry of $(K_i)^{-1}$ is nonpositive). Therefore, all entries in

$$ (-e_{21}^0 k_{21}) \cdot (K_1)^{-1} \qquad \text{and} $$

$$ (-e_{m+1,m}^0 k_{m+1,m}) \cdot (K_2)^{-1} $$

are nonnegative. We conclude that coefficients $\mathbf{b}_i$ in (15) are nonnegative. $\square$

<figure><img src="image_4.png" /></figure>

*Note.* Although Proposition 1 is proved for nonnegative constant values $X_1$ and $X_m$, we consider $\bar{Y}_i$ in (15) also as functions of $t \in [0, T_1]$. Then we have the functions in (8), defined for the QSSR.

## 3.4 Study of the slowly varying system

The dynamics of the slow system (metabolites which do accumulate) are obtained by setting $\eta = 0$ in (11) and substituting the fast variables $Y_i$ for the expression given by (15):
$$ 
\begin{aligned}
\frac{dX_1}{dt} &= I(t) - e_{21}^0 k_{21} \frac{X_1}{X_1 + K_{21}} - \mu X_1, \\
\frac{dX_m}{dt} &= e_{m,m-1}^0 k_{m,m-1} \frac{Y_{m-1}}{K_{m,m-1}} - \left( e_{m+1,m}^0 k_{m+1,m} \frac{X_m}{X_m + K_{m+1,m}} \right) - \mu X_m, \\
\frac{dX_n}{dt} &= e_{n,n-1}^0 k_{n,n-1} \frac{Y_{n-1}}{K_{n,n-1}} - \mu X_n.
\end{aligned}
$$

Then we obtain the remaining dynamical system (7), which provides the dynamics to the overall network.

The other variables of the metabolic network, which are the fast variables $Y_i$ (indeed, most of the variables are fast) can then be reconstructed after the solution of the algebraic Equation (13), given in (8). Finally, these fast variables rely on system (7). This system is also referred as the quasi steady state system [6].

**Proposition 2.** If system (7) has nonnegative initial conditions, then it has a unique nonnegative solution ($\overline{X}_1, \overline{X}_m, \overline{X}_n$) defined on the interval $[0, T_1]$.

For the proof of Proposition 2, see Appendix C.

## 3.5 Tikhonov's Theorem

Propositions 1 and 2 prove that the class of systems with the form (11)-(12) satisfy the hypothesis of Tikhonov's Theorem [6]. Then, we can apply this theorem to system (10).

The following proposition is a consequence of Tikhonov's Theorem [6]. The proposition states that the approximation given by the QSSR (7)-(8) has an error with order $O(\varepsilon)$, after a fast initial transient for the fast variables.

**Proposition 3.** [Deduction of Tikhonov's Theorem] If $I(t)$ is a nonnegative continuous function over $[0, T_1]$, then

$$ X_j(t) = \overline{X}_j(t) + O(\varepsilon) \qquad j = 1, m, n, \forall t \in [0, T_1]. $$

Moreover, there exists $T_0 > 0$ such that for every $t \in [T_0, T_1]$,

$$ X_i(t) = \varepsilon [\overline{Y}_i(t) + O(\varepsilon)] \qquad \forall i = 2, \dots, n-1, i \neq m. $$

where $X_i$ are the solutions of the original system (3), $\overline{X}_i$ are the solutions of (7) and $\overline{Y}_i$ are the functions defined in (8) after the algebraic Equation (13).

*Note.* The solution of the boundary layer problem for system (11)-(12) is similar to that of Equation (14). We include its demonstration in Appendix B.
# 4 Magnitude of Concentrations throughout the metabolic network

In this section we study the magnitude of metabolite concentrations, depending on if they are associated to slow or fast reactions. They are deduced from the reduced system after Tikhonov's Theorem (7)-(8). We now show that the concentration of metabolites in QSS (that do not trap the input flux) is one order of magnitude $\varepsilon$ lower than metabolites with slow dynamics.

In order to prove this assertion, we define the conditions under which

$$ \mathbf{b}_i \cdot e_{21}^0 k_{21} = \mathcal{O}(1), $$

$$ \mathbf{b}_i \cdot e_{m+1,m}^0 k_{m+1,m} = \mathcal{O}(1), $$

to obtain

$$ \bar{X}_i(t) = \mathcal{O}\left(\varepsilon \cdot \frac{\bar{X}_1(t)}{\bar{X}_1(t) + K_{21}}\right), $$

$$ \bar{X}_i(t) = \mathcal{O}\left(\varepsilon \cdot \frac{\bar{X}_m(t)}{\bar{X}_m(t) + K_{m+1,m}}\right), $$

for every $t \in [T_0, T_1]$.

## 4.1 Parameter Orders

We show that all off-diagonal entries of the Jacobian matrix $K_i$ have the same order of magnitude, for both matrices defined in (17).

**Lemma 1.** Suppose that the parameters of each Michaelis-Menten enzymatic reaction (see Appendix A) satisfy

$$ \mathcal{O}(k_{\gamma}^{ji}) = \mathcal{O}(k_{ji}) \qquad \forall i, j = 1, \dots, n, \gamma \in \{-1, 1\}. \qquad (19) $$

Then,

$$ \mathcal{O}\left(\frac{e_{ji}^0 k_{ji}}{K_{ji}}\right) = \mathcal{O}\left(\frac{e_{j'i'}^0 k_{j'i'}}{K_{j'i'}}\right) = \mathcal{O}(e_{j'i'}^0 k_{j'i'}) \quad \forall i, j, i', j' \in \{1, \dots, n\}. $$

*Proof.* As a consequence of (5),

$$ \mathcal{O}(e_{ji}^0 k_{ji}) = \mathcal{O}(e_{j'i'}^0 k_{j'i'}) \quad \forall i, j, i', j' \in \{1, \dots, n\}. \qquad (20) $$
Moreover, by the definition of the Michaelis-Menten constant (4) and (19), we have that $O(K_{ji}) = 1$ and then

$$O\left(\frac{e_{ji}^0 k_{ji}}{K_{ji}}\right) = O(e_{ji}^0 k_{ji}) \quad \forall i, j. \qquad (21)$$

Hence, combining (20) and (21), we have

$$O\left(\frac{e_{ji}^0 k_{ji}}{K_{ji}}\right) = O\left(\frac{e_{j'i'}^0 k_{j'i'}}{K_{j'i'}}\right) = O(e_{j'i'}^0 k_{j'i'}) \quad \forall i, j, i', j' \in \{1, \dots, n\}. \quad \square$$

<figure><img src="image_5.png" /></figure>

Actually, all the entries of the matrix $K_i$ have the same order of magnitude, as asserts the following corollary.

**Proposition 4.** Consider the matrices defined in (17). All the entries of $K_1$ and $K_2$ have the same order.

*Proof.* According to Lemma 1, for the sums in the diagonal of the matrices we have

$$\sum_{\substack{i=2 \\ i \neq j}}^{n} \frac{e_{ij}^0 k_{ij}}{K_{ij}} = O(e_{ji}^0 k_{ji}).$$

Moreover, $\tilde{\mu} \ll e_{ji}^0 k_{ji}$. Then,

$$\sum_{\substack{i=2 \\ i \neq j}}^{n} \frac{e_{ij}^0 k_{ij}}{K_{ij}} + \tilde{\mu} = O(e_{ji}^0 k_{ji}).$$

For the off-diagonal entries consider (21). Therefore, all the entries of $K_1$ and $K_2$ have order $O(e_{ji}^0 k_{ji})$. $\square$

<figure><img src="image_6.png" /></figure>

## 4.2 A Theorem for Magnitude of Concentrations

In order to prove that a metabolite in QSS does not reach high concentrations, we have to suppose that it is not in a trap for the input flux. The definition of trap was introduced in [14], and we formally adapt it to the class of models considered in this article (see Appendix D.1). Then, we define a *flux trap*, which is a trap reached by the flux.
Assumption 1. There exists F a flux from $X_1$ to $X_n$ in the system of enzymatic reactions (3) (depicted in Figure 1). Moreover, we define $\mathcal{I}_{\mathcal{T}_F}$ as the set of indices such that $X_i$ is in a flux trap, for every $i \in \mathcal{I}_{\mathcal{T}_F}$, and $X_j$ is not in a flux trap for every $j \in \{1, \dots, n\} \setminus \mathcal{I}_{\mathcal{T}_F}$.

Notice that the presence of the flux F from $X_1$ to $X_n$ implies

$$ \{1, m, n\} \cap \mathcal{I}_{\mathcal{T}_F} = \emptyset. $$

Then, flux traps are only possible within the subnetworks in QSS. Also, $\mathcal{I}_{\mathcal{T}_F} = \emptyset$ if there is no flux trap.

The following lemma sets down the order of magnitude of the parameters in (8), for the metabolites which are not in a flux trap. These parameters are used for writing the expression of fast metabolites in the QSSR.

**Lemma 2.** Suppose the system of enzymatic reactions (3) (Figure 1) under Assumption 1. Consider the parameters $\mathbf{b}_i$ of Equation (8), obtained in Section 3. Then, if $\mathbf{b}_i \neq 0$, it holds

$$ \mathbf{b}_i \cdot e_{21}^0 k_{21} = \mathcal{O}(1) \qquad \text{if } i \in \{2, \dots, m-1\} \setminus \mathcal{I}_{\mathcal{T}_F}, \qquad (22) $$

$$ \mathbf{b}_i \cdot e_{m+1,m}^0 k_{m+1,m} = \mathcal{O}(1) \qquad \text{if } i \in \{m+1, \dots, n-1\} \setminus \mathcal{I}_{\mathcal{T}_F}. $$

*Proof.* From the results stated in Appendix D, particularly Theorem 2, we have for $\mathbf{b}_i \neq 0$,

$$ \mathbf{b}_i \cdot \frac{e_{21}^0 k_{21}}{K_{21}} = \mathcal{O}(1) \qquad \text{if } i \in \{2, \dots, m-1\} \setminus \mathcal{I}_{\mathcal{T}_F}, $$

$$ \mathbf{b}_i \cdot \frac{e_{m+1,m}^0 k_{m+1,m}}{K_{m+1,m}} = \mathcal{O}(1) \qquad \text{if } i \in \{m+1, \dots, n-1\} \setminus \mathcal{I}_{\mathcal{T}_F}. $$

Using Equation (21) we conclude that, for $\mathbf{b}_i \neq 0$,

$$ \mathbf{b}_i \cdot e_{21}^0 k_{21} = \mathcal{O}(1) \qquad \text{if } i \in \{2, \dots, m-1\} \setminus \mathcal{I}_{\mathcal{T}_F}, $$

$$ \mathbf{b}_i \cdot e_{m+1,m}^0 k_{m+1,m} = \mathcal{O}(1) \qquad \text{if } i \in \{m+1, \dots, n-1\} \setminus \mathcal{I}_{\mathcal{T}_F}. $$

□

The next theorem is a powerful conclusion obtained after the QSSR (7)-(8). Theorem 1 states that the concentration of a metabolite in QSS, which
is not in a flux trap, is one order of magnitude $\epsilon$ lower than the concentration of a metabolite with slow dynamics. This result holds even if there is a trap or a flux trap in the system.

**Theorem 1.** [Magnitude of Concentration Theorem] Consider the system of enzymatic reactions (3) (Figure 1). Under Assumption 1, the following inequalities hold:

$$ 
\begin{aligned}
\bar{X}_i &\le \mathcal{O}(\epsilon \cdot \bar{X}_1) & \forall i \in \{2, \dots, m-1\} \setminus \mathcal{I}_{\mathcal{T}_F}, \\
\bar{X}_i &\le \mathcal{O}(\epsilon \cdot \bar{X}_m) & \forall i \in \{m+1, \dots, n-1\} \setminus \mathcal{I}_{\mathcal{T}_F}.
\end{aligned}
$$

*Proof of Theorem 1.* Since the reduction from Tikhonov's Theorem, we have Equation (8), i.e.,

$$ 
\begin{aligned}
\bar{X}_i &= \epsilon \cdot \mathbf{b}_i \cdot e_{21}^0 k_{21} \left( \frac{\bar{X}_1}{\bar{X}_1 + K_{21}} \right) & i = 2, \dots, m-1, \\
\bar{X}_i &= \epsilon \cdot \mathbf{b}_i \cdot e_{m+1,m}^0 k_{m+1,m} \left( \frac{\bar{X}_m}{\bar{X}_m + K_{m+1,m}} \right) & i = m+1, \dots, n-1.
\end{aligned}
$$

Then, as stated in Lemma 2, for $i$ such that $\mathbf{b}_i \neq 0$,

$$ 
\begin{aligned}
\bar{X}_i &= \mathcal{O}\left(\epsilon \frac{\bar{X}_1}{\bar{X}_1 + K_{21}}\right) & \text{if } i \in \{2, \dots, m-1\} \setminus \mathcal{I}_{\mathcal{T}_F}, \\
\bar{X}_i &= \mathcal{O}\left(\epsilon \frac{\bar{X}_m}{\bar{X}_m + K_{m+1,m}}\right) & \text{if } i \in \{m+1, \dots, n-1\} \setminus \mathcal{I}_{\mathcal{T}_F}.
\end{aligned}
$$

But $1 \le \mathcal{O}(\bar{X}_i + K_{i+1,i})$, because system (7) is positively invariant and $\mathcal{O}(K_{i+1,i}) = 1$. Hence,

$$ 
\begin{aligned}
\mathcal{O}\left(\frac{\bar{X}_1}{\bar{X}_1 + K_{21}}\right) &\le \mathcal{O}(\bar{X}_1), \\
\mathcal{O}\left(\frac{\bar{X}_m}{\bar{X}_m + K_{m+1,m}}\right) &\le \mathcal{O}(\bar{X}_m).
\end{aligned}
$$

We conclude that

$$ 
\begin{aligned}
\bar{X}_i &\le \mathcal{O}(\epsilon \cdot \bar{X}_1) & \forall i \in \{2, \dots, m-1\} \setminus \mathcal{I}_{\mathcal{T}_F}, \\
\bar{X}_i &\le \mathcal{O}(\epsilon \cdot \bar{X}_m) & \forall i \in \{m+1, \dots, n-1\} \setminus \mathcal{I}_{\mathcal{T}_F}.
\end{aligned}
$$

□
Note. The approach presented in this work can be used to reduce a metabolic network which has flux traps, obtaining an error characterization (as established in Proposition 3) and the conclusion of Theorem 1. But, in agreement with Theorem 1, the magnitude of concentrations of the metabolites in the flux traps cannot be bounded by the concentration of the metabolites in the slow part of the system. This fact can be inferred from the proof of Theorem 1 (see Appendix D.2.)

The presence of a flux trap leads to accumulation, without reuse, of compounds in the metabolic network. However, accumulation of some compounds to large concentrations often results to cell death. For example, the accumulation of lactate has been recognized as one cause of cell death [15, 16].

# 5 Reduced Model Calibration

Now that we have described the way to synthesize the initial model of the metabolic network into a small dynamical system (for accumulating metabolites) and a set of algebraic equations, we will explain how to calibrate this reduced model from experimental data. Of course, we assume that the initial stoichiometric coefficients are known, but the parameters associated to reaction rates are unknown.

Here we propose a method to estimate the parameters of the reduced system. In a first stage we identify the parameters of the reduced dynamical system representing the accumulating metabolites (7). The identification method is based on the minimization of a cost function, computing the error model with respect to experimental data.

Furthermore, if data of any metabolite in QSS is available, we can also estimate the respective parameters in (8), to write its concentration as a linear combination of the slow metabolites.

## 5.1 Calibration of the slow dynamics

We suppose experimental data of the metabolites in the slow part of the system (11), denoted by

$$Z_i(t_j) = X_i(t_j) + \beta_i(t_j) \qquad i = 1, m, n, j = 1, \dots, r, \qquad (23)$$

where $X_i$ is the solution of the original system (3) and $\beta_i$ represents an error of measurements. In order to estimate the parameters of the reduced system
(7), we rewrite it as

$$
\begin{align}
\frac{d\overline{X}_1}{dt} &= I(t) - \frac{\theta_1 \overline{X}_1}{\overline{X}_1 + \theta_2} - \theta_3 \cdot \overline{X}_1 & \overline{X}_1(0) &= Z_1(t_1) \tag{24} \\
\frac{d\overline{X}_m}{dt} &= \frac{\theta_4 \overline{X}_1}{\overline{X}_1 + \theta_2} - \frac{\theta_5 \overline{X}_m}{\overline{X}_m + \theta_6} - \theta_3 \overline{X}_m & \overline{X}_m(0) &= Z_m(t_1) \\
\frac{d\overline{X}_n}{dt} &= \frac{\theta_7 \overline{X}_m}{\overline{X}_m + \theta_6} - \theta_3 \cdot \overline{X}_n & \overline{X}_n(0) &= Z_n(t_1),
\end{align}
$$

Let $\theta = (\theta_1, \theta_2, \theta_3, \theta_4, \theta_5)$ and define a cost function $\mathcal{F}(\theta)$. This cost function has to measure the error between the solution of (24) and the data $Z_1, Z_m, Z_n$, for every value $\theta$ in a domain $D \subset \mathbb{R}^7$. For example, we can define $\mathcal{F}$ as

$$
\mathcal{F}(\theta) = \sum_{i \in \{1, m, n\}} \sum_{j=1}^{r} \left( \bar{X}_i(\theta, t_j) - Z_i(t_j) \right)^2 . \quad (25)
$$

Then, we have to find $\hat{\theta} = (\hat{\theta}_1, \hat{\theta}_2, \hat{\theta}_3, \hat{\theta}_4, \hat{\theta}_5)$ such that

$$
\mathcal{F}(\hat{\theta}) = \min\{F(\theta) : \theta \in D\}.
$$

Note. For obtaining the vector of parameters $\hat{\theta}$ to calibrate (24), it is not
necessary to have data of any metabolite in QSS, $X_i$ with $i = 2, \dots, n-1, i \neq m$. Only the data (23) of the metabolites in the slow part, $X_1, X_m, X_n$, is
used.

5.2 Fast dynamics parameters

In some (rare) cases, measurements of some fast metabolites can be available.
Generally, these data are only obtained at quasi steady state after the initial
transient and for a subset of the metabolic compounds.

Supposing that we have experimental data of the metabolites in QSS after
the initial fast transient,

$$
\begin{align}
Z_i(t'_j) = X_i(t'_j) + \mathcal{N}(t'_j) \tag{26} \\
\qquad i = 2, \dots, n-1, i \neq m, \notag \\
\qquad T_0 \le t'_1 < \dots < t'_{r'}, \notag
\end{align}
$$

and that we have obtained $\hat{\theta}$ after calibrating (24), we can estimate the
parameters in (8). As a matter of fact, in line with the reduced system
(7)-(8) and the calibrated system (24), for the metabolites in QSS we have

$$ \bar{X}_i = \alpha_i \cdot \frac{\bar{X}_1}{\bar{X}_1 + \hat{\theta}_2} \qquad i = 2, \dots, m-1 \qquad (27) $$

$$ \bar{X}_i = \alpha_i \cdot \frac{\bar{X}_m}{\bar{X}_m + \hat{\theta}_6} \qquad i = m + 1, \dots, n-1, $$

where $\alpha_i$ are the parameters to be estimated.

Here, we can explicitly resolve the linear least square problem. The least squares solution that minimize the difference between the data $Z_i$ and the expressions in (27) is the following [17]:

$$ \hat{\alpha}_i = \frac{\sum_{j=1}^{r'} \frac{Z_i(t'_j) \bar{X}_1(\hat{\theta}, t'_j)}{\bar{X}_1(\hat{\theta}, t'_j) + \hat{\theta}_2}}{\sum_{j=1}^{r'} \left( \frac{\bar{X}_1(\hat{\theta}, t'_j)}{\bar{X}_1(\hat{\theta}, t'_j) + \hat{\theta}_2} \right)^2} \quad \forall i = 2, \dots, m-1 $$

$$ \hat{\alpha}_i = \frac{\sum_{j=1}^{r'} \frac{Z_i(t'_j) \bar{X}_m(\hat{\theta}, t'_j)}{\bar{X}_m(\hat{\theta}, t'_j) + \hat{\theta}_6}}{\sum_{j=1}^{r'} \left( \frac{\bar{X}_m(\hat{\theta}, t'_j)}{\bar{X}_m(\hat{\theta}, t'_j) + \hat{\theta}_6} \right)^2} \quad \forall i = m + 1, \dots, n-1. $$

Indeed, we look for values of $\hat{\alpha}_i$ that minimize the differences

$$ L_i(\alpha) = \sum_{j=1}^{r'} \left( \alpha \frac{\bar{X}_1(t'_j, \hat{\theta})}{\bar{X}_1(t'_j, \hat{\theta}) + \hat{\theta}_2} - Z_i(t'_j) \right)^2 \quad i = 2, \dots, m-1 $$

$$ L_i(\alpha) = \sum_{j=1}^{r'} \left( \alpha \frac{\bar{X}_m(t'_j, \hat{\theta})}{\bar{X}_m(t'_j, \hat{\theta}) + \hat{\theta}_6} - Z_i(t'_j) \right)^2 \quad i = m + 1, \dots, n-1. $$

*Note.* To obtain the parameter $\hat{\alpha}_i$, we only need the data $Z_i$ (of the corresponding metabolites in QSS, $X_i$) and the calibrated system (24) with $\hat{\theta}$.

# 6 Illustrative example with a Toy Enzymatic Network

In this section we apply the method developed in this paper to the toy network represented in Figure 3. This toy network accounts for one reversible enzymatic reaction and a cycle of enzymatic reactions. Moreover, the toy
network contains two subnetworks in QSS (in blue in Figure 3), which are interconnected by metabolites with slow rates of consumption (in black in Figure 3).

<figure><img src="image_7.png" /><figcaption>**Figure 3:** We consider that reactions represented by black arrows are slow, while reactions represented by blue arrows are fast. Metabolites in black are accumulative, whereas metabolites in blue are non accumulative and they are supposed to be in Quasi Steady State. Every reaction is catalyzed by and enzyme $e_{ji}$.</figcaption></figure>

First we consider the ODE of the toy enzymatic network, as in Section 1. Then, using the time-scale separation hypothesis, we reduce this ODE with the method described in Section 3. Finally, we estimate the parameters of the reduced system as it is suggested in Section 5

All the parameters in the toy network are supposed to satisfy the conditions established in (5) and Section 4. The periodic and continuous input considered is given by

$$I(t) = k[\cos(t \cdot \omega + \pi) + 1],$$

where k is a parameter with the same order of magnitude as the slow reactions rates.

## 6.1 Reduction

We apply to the toy network our reduction scheme, as described in Section 3. First, to simplify the notation, we define the following parameters:

$$a_{ji} := \frac{e_{ji}^0 k_{ji}}{K_{ji}} \quad \forall i, j = 1, \dots, n,$$

$$d_1 := a_{32}(a_{43} + \tilde{\mu}) + \tilde{\mu}(a_{23} + a_{43} + \tilde{\mu}),$$

$$d_2 := (a_{65} + \tilde{\mu})(a_{86} + \tilde{\mu})(a_{57} + \tilde{\mu})(a_{78} + a_{98} + \tilde{\mu}) - (a_{65})(a_{57})(a_{86})(a_{78}).$$
Then, we obtain the following reduced system for the toy network,

$$ \frac{dX_1}{dt} = I(t) - \frac{e_{21}^0 k_{21} \cdot X_1}{X_1 + K_{21}} - \mu X_1 \quad (28) $$

$$ \frac{dX_4}{dt} = \frac{(a_{43})(a_{32})}{d_1} \cdot \frac{e_{21}^0 k_{21} \cdot X_1}{X_1 + K_{21}} - \frac{e_{54}^0 k_{54} \cdot X_4}{X_4 + K_{54}} - \mu X_4 $$

$$ \frac{dX_9}{dt} = \frac{(a_{98})(a_{65})(a_{86})(a_{57} + \tilde{\mu})}{d_2} \cdot \frac{e_{54}^0 k_{54} \cdot X_4}{X_4 + K_{54}} - \mu X_9, $$

and the expressions for the metabolites in QSS,

$$ X_2 = \frac{\varepsilon(a_{23} + a_{43} + \tilde{\mu})}{d_1} \cdot \frac{e_{21}^0 k_{21} X_1}{X_1 + K_{21}} \quad (29) $$

$$ X_3 = \frac{\varepsilon(a_{32})}{d_1} \cdot \frac{e_{21}^0 k_{21} X_1}{X_1 + K_{21}} $$

$$ X_5 = \frac{\varepsilon(a_{86} + \tilde{\mu})(a_{57} + \tilde{\mu})(a_{78} + a_{98} + \tilde{\mu})}{d_2} \cdot \frac{e_{54}^0 k_{54} X_4}{X_4 + K_{54}} $$

$$ X_6 = \frac{\varepsilon(a_{65})(a_{57} + \tilde{\mu})(a_{78} + a_{98} + \tilde{\mu})}{d_2} \cdot \frac{e_{54}^0 k_{54} X_4}{X_4 + K_{54}} $$

$$ X_7 = \frac{\varepsilon(a_{65})(a_{86})(a_{78})}{d_2} \cdot \frac{e_{54}^0 k_{54} X_4}{X_4 + K_{54}} $$

$$ X_8 = \frac{\varepsilon(a_{65})(a_{86})(a_{57} + \tilde{\mu})}{d_2} \cdot \frac{e_{54}^0 k_{54} X_4}{X_4 + K_{54}}. $$

## 6.2 Calibration of the Reduced Toy Network

We follow the procedure in Section 5. For simplicity we suppose that the data are measured at the same time instants $t_1, \dots, t_r$ (we assume that 48 measurement instants are available) for the slow and the fast parts of the system.

The measurements are the variables (units g/L) of the original system (3) (for the toy network in Figure 3) plus a white noise:

$$ Z_i(t_j) := X_i(t_j) + \beta(t_j) \qquad j = 1, 2, \dots, r \quad (30) $$

where $\beta \sim N(\sigma_i)$ and $\sigma_i = 10^{-1} \cdot \text{median}(X_i)$ for every $i = 1, \dots, n$.

As in Section 5, to estimate the parameters of (28), we use the reduced system (24), with $m = 4$ and $n = 9$. The cost function considered is $F$, defined in (25), with $m = 4$ and $n = 9$.
The function `fminsearch` in Scilab was used for minimizing $\mathcal{F}$. This function is based on the Nelder-Maid algorithm to compute the unconstrained minimum of a given function. For the simulations in Figure 4, the value $\hat{\theta}$ obtained is in Table 1 and $\mathcal{F}(\hat{\theta}) = 0.097$.

Here, for illustration purpose, we suppose that the metabolites in QSS are also measured, we calculate the parameters $\hat{\alpha}$ to estimate their concentrations as explained in Section 5.2. Then, their concentrations are obtained according to (27).

We computed the numerical solution of the systems describing the dynamics in the toy network of Figure 3. The results are represented in Figure 4 and Figure 5. As expected, the concentrations of the metabolites in QSS are one order of magnitude $\varepsilon$ lower than the metabolites in the slow part.

Note that the parameters $\theta_2$ and $\theta_6$ are affinity constant in Michaelis-Menten functions, whose sensitivity is low [18]. Here we have used 48 samples for parameter identification.

It is worth noting that the identification process results in a satisfying agreement between simulations of the calibrated system (24)-(27) and recorded data, as represented in Figure 4 and Figure 5.

<table><caption>**Table 1:** Parameter estimation for system (28), rewritten as (24). The estimation of this parameters only requires the slow dynamics of the toy network in Figure 3.</caption><thead><tr><th>i</th><th>Theoretical value θ<sub>i</sub></th><th>Initial guess</th><th>Estimated value θ̂<sub>i</sub></th><th>Units</th><th>Error percent</th></tr></thead><tbody><tr><td>1</td><td>0.110</td><td>0.010</td><td>0.072</td><td>g(L.min)<sup>-1</sup></td><td>34.54</td></tr><tr><td>2</td><td>2.000</td><td>1.000</td><td>1.298</td><td>gL<sup>-1</sup></td><td>35.10</td></tr><tr><td>3</td><td>0.010</td><td>0.010</td><td>0.011</td><td>min<sup>-1</sup></td><td>10.00</td></tr><tr><td>4</td><td>0.110</td><td>0.010</td><td>0.073</td><td>g(L.min)<sup>-1</sup></td><td>33.63</td></tr><tr><td>5</td><td>0.013</td><td>0.010</td><td>0.006</td><td>g(L.min)<sup>-1</sup></td><td>53.85</td></tr><tr><td>6</td><td>2.000</td><td>1.000</td><td>2.143</td><td>gL<sup>-1</sup></td><td>7.15</td></tr><tr><td>7</td><td>0.013</td><td>0.010</td><td>0.016</td><td>g(L.min)<sup>-1</sup></td><td>23.08</td></tr></tbody></table>
<table><caption>**Table 2:** Estimation of the parameters in (29), corresponding to the equalities in (27).</caption><thead><tr><th>i</th><th>Theoretical value α<sub>i</sub></th><th>Estimated value α̂<sub>i</sub></th><th>Error percent</th></tr></thead><tbody><tr><td>2</td><td>0.00124</td><td>0.00092</td><td>25.81</td></tr><tr><td>3</td><td>0.00122</td><td>0.00089</td><td>27.05</td></tr><tr><td>5</td><td>0.00168</td><td>0.00185</td><td>10.12</td></tr><tr><td>6</td><td>0.00018</td><td>0.00019</td><td>5.56</td></tr><tr><td>7</td><td>0.00002</td><td>0.00002</td><td>0</td></tr><tr><td>8</td><td>0.00017</td><td>0.00019</td><td>11.76</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 4:** Dynamics of the toy network represented in Figure 3. The functions $X_i$ represent the metabolite concentrations in units $(gL^{-1})$. The numerical solution of the original system (3) is depicted by the green line; the reduced system obtained by the method exposed in this work (28), by the blue dashed line; the supposed data with white noise (30), by green points; and the calibrated system (24)-(27) with the estimated parameters in Table 1 and Table 2, by the red line. The parameters considered are in Table 3 and Table 4. As expected, the concentrations of the metabolites in QSS are one order of magnitude $\epsilon$ lower than the metabolites in the slow part.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5:** Zoom on the concentration of metabolites in QSS in Figure 4.</figcaption></figure>
<table><caption>**Table 3:** Parameters considered for the simulations in Figure 4. The symbol $\gamma \in \{-1, 1\}$ denotes a rate in an enzymatic reaction (see the Michaelis-Menten Equation (31)). The initial conditions for all the enzymes are the same, as well as the initial conditions of all the metabolites are identical, i.e. $j, i \in \{1, \dots, n\}$ in this table.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td>ε</td><td>0.001</td><td>–</td></tr><tr><td>k</td><td>0.01</td><td>g(L.min)<sup>-1</sup></td></tr><tr><td>μ</td><td>0.01</td><td>min<sup>-1</sup></td></tr><tr><td>ω</td><td>0.004</td><td>–</td></tr><tr><td>k<sub>γ</sub><sup>21</sup>, k<sub>21</sub></td><td>1.10</td><td>min<sup>-1</sup></td></tr><tr><td>k<sub>γ</sub><sup>54</sup>, k<sub>54</sub></td><td>0.13</td><td>min<sup>-1</sup></td></tr><tr><td>k<sub>γ</sub><sup>32</sup>, k<sub>32</sub></td><td>1.90</td><td>min<sup>-1</sup></td></tr><tr><td>k<sub>γ</sub><sup>23</sup>, k<sub>23</sub></td><td>0.12</td><td>min<sup>-1</sup></td></tr><tr><td>k<sub>γ</sub><sup>43</sup>, k<sub>43</sub></td><td>1.80</td><td>min<sup>-1</sup></td></tr><tr><td>k<sub>γ</sub><sup>65</sup>, k<sub>65</sub></td><td>0.17</td><td>min<sup>-1</sup></td></tr><tr><td>k<sub>γ</sub><sup>57</sup>, k<sub>57</sub></td><td>1.40</td><td>min<sup>-1</sup></td></tr><tr><td>k<sub>γ</sub><sup>86</sup>, k<sub>86</sub></td><td>1.60</td><td>min<sup>-1</sup></td></tr><tr><td>k<sub>γ</sub><sup>78</sup>, k<sub>78</sub></td><td>0.15</td><td>min<sup>-1</sup></td></tr><tr><td>k<sub>γ</sub><sup>98</sup>, k<sub>98</sub></td><td>1.50</td><td>min<sup>-1</sup></td></tr><tr><td>e<sub>ji</sub><sup>0</sup></td><td>0.10</td><td>gL<sup>-1</sup></td></tr><tr><td>x<sub>i</sub><sup>0</sup></td><td>0.001</td><td>gL<sup>-1</sup></td></tr></tbody></table>

<table><caption>**Table 4:** Slow and fast reaction rates considered for the toy network in Figure 3 and the simulations in Figure 4. Fast reaction rates are characterized by the factor $1/\epsilon$. All reaction rates are in units of $g(L.min)^{-1}$.</caption><thead><tr><th>Slow rates</th><th>Value</th></tr></thead><tbody><tr><td>e<sub>21</sub><sup>0</sup>k<sub>21</sub></td><td>1.1 &times; 10<sup>-1</sup></td></tr><tr><td>e<sub>54</sub><sup>0</sup>k<sub>54</sub></td><td>1.3 &times; 10<sup>-2</sup></td></tr><tr><th>Fast rates</th><th>Value</th></tr><tr><td>e<sub>32</sub><sup>0</sup>k<sub>32</sub>/ε</td><td>1.3 &times; 10<sup>2</sup></td></tr><tr><td>e<sub>23</sub><sup>0</sup>k<sub>23</sub>/ε</td><td>1.2 &times; 10<sup>1</sup></td></tr><tr><td>e<sub>43</sub><sup>0</sup>k<sub>43</sub>/ε</td><td>1.8 &times; 10<sup>2</sup></td></tr><tr><td>e<sub>65</sub><sup>0</sup>k<sub>65</sub>/ε</td><td>1.7 &times; 10<sup>1</sup></td></tr><tr><td>e<sub>57</sub><sup>0</sup>k<sub>57</sub>/ε</td><td>1.4 &times; 10<sup>2</sup></td></tr><tr><td>e<sub>86</sub><sup>0</sup>k<sub>86</sub>/ε</td><td>1.6 &times; 10<sup>2</sup></td></tr><tr><td>e<sub>78</sub><sup>0</sup>k<sub>78</sub>/ε</td><td>1.5 &times; 10<sup>1</sup></td></tr><tr><td>e<sub>98</sub><sup>0</sup>k<sub>98</sub>/ε</td><td>1.5 &times; 10<sup>2</sup></td></tr></tbody></table>
# 7 Discussion

## 7.1 Time-Scales Hypotheses

Metabolic networks can involve much more than two different time scales. Actually, our method considers the division of these in two groups of reaction rates. The kinetics slower than a certain threshold, and the kinetics faster than this threshold. Our approach eventually preserve the dynamics of the slower kinetics (keeping the different time scales), while the fastest dynamics are lumped and approximated.

Also, to better illustrate this important aspect, we have considered several reaction rate orders in the toy network. The reaction rates are divided into slow and fast, and each group of reactions has different scales (see Table 3 and Table 4). The simulation results illustrate Theorem 1 (see Figure 4) and the reduced system accurately represents all different time scales (see Figure 4 and Figure 5).

Finally, note that it would be possible to set up a finer approximation considering several time-scales for Tikhonov's Theorem, but at the risk of higher mathematical complexity. Indeed, extended versions of Tikhonov's Theorem exist for several time-scales, using powers of $\varepsilon$ [6, 19, 20] or even different epsilons [21]. But computations with this method highly complicates the reduction.

## 7.2 Comparison with Experimental Data

To the best knowledge of the authors, there are to date no example of metabolome measured at high frequency, at least for a large number of metabolites to assess the kinetics. In general, only a very limited number macromolecules (typically proteins, carbohydrates, lipids, chlorophyll,...) are recorded, specially for microalgae. However, to show that our findings are in agreement with experimental studies we considered the results from [4] for an autotrophic microalgae metabolic network.

The authors in [4] fitted parameters of a metabolic model to the set of available experimental data. We examined the reaction rates which ranged from $10^2$ to $10^{-1}$ ($h^{-1}.mM.B^{-1}$) and compared them with the level of concentrations in the cell. Indeed (see Table 5 and Table 6), the concentration of carbohydrates has magnitude $10^2$ (mM) times higher than those of the intermediate metabolites (GAP, PEP and G6P). Moreover, GAP, PEP and
G6P are consumed by reactions with rate order $10^1$ or $10^2$ ($h^{-1}.mM.B^{-1}$), while carbohydrates are consumed at rate of order $10^0$ ($h^{-1}.mM.B^{-1}$). Additionally, carbohydrates are produced by a single reaction with rate of order $10^1$ ($h^{-1}.mM.B^{-1}$), as well as GAP, and G6P and PEP by reactions of order $10^2$ ($h^{-1}.mM.B^{-1}$). This evidences that the concentration is related to the rate of consumption, in the way predicted by Theorem 1 in our paper.

Nevertheless, we emphasize that, the reduction method proposed in this paper can be used even if only some metabolites with large concentration have been measured. Indeed, such data will support the calibration of the reduced model, i.e. describing the dynamics of the slow metabolites (see Section 5.1).

<table><caption>**Table 5:** Experimental measures (M) and estimated (E) values obtained from [4], for an autotrophic microalga metabolic network [4]. Carbon quotas of the different compounds are considered within a period of 24 hours. Light intensity values are on an interval from 0 to 1400 $uE.m^{-2}.s^{-1}$. Two different magnitudes of concentration can be distinguished among these compounds.</caption><thead><tr><th>Compound</th><th>Value</th><th>Mean value</th></tr></thead><tbody><tr><td>Carbohydrates</td><td>M</td><td>8.436 × 10<sup>-1</sup> mM</td></tr><tr><td>G6P</td><td>E</td><td>5.208 × 10<sup>-3</sup> mM</td></tr><tr><td>PEP</td><td>E</td><td>4.167 × 10<sup>-3</sup> mM</td></tr><tr><td>GAP</td><td>E</td><td>1.389 × 10<sup>-3</sup> mM</td></tr></tbody></table>

<table><caption>**Table 6:** Rates are in $h^{-1}.mM.B^{-1}$. Typical concentrations in Table 5 were used to estimate the consumption rates for GAP and PEP in the lipid synthesis reaction.</caption><thead><tr><th>Compound</th><th>Production rate</th><th>Consumption rate</th><th>Sub-network</th></tr></thead><tbody><tr><td>Carbohydrates</td><td>7.00 × 10<sup>1</sup></td><td>6.50 × 10<sup>0</sup></td><td>Carbohydrate synthesis</td></tr><tr><td rowspan="2">G6P</td><td>2.24 × 10<sup>2</sup></td><td>1.03 × 10<sup>1</sup></td><td>Upper glycolysis</td></tr><tr><td>6.50 × 10<sup>0</sup></td><td>7.00 × 10<sup>1</sup></td><td>Lipid synthesis</td></tr><tr><td rowspan="2">PEP</td><td>4.37 × 10<sup>2</sup></td><td>5.00 × 10<sup>0</sup></td><td>Lower glycolysis</td></tr><tr><td>9.97 × 10<sup>0</sup></td><td>1.04 × 10<sup>2</sup></td><td>Lipid Synthesis</td></tr><tr><td rowspan="3">GAP</td><td>2.06 × 10<sup>1</sup></td><td>4.47 × 10<sup>2</sup></td><td>Upper glycolysis</td></tr><tr><td>5.00 × 10<sup>0</sup></td><td>4.37 × 10<sup>2</sup></td><td>Upper glycolysis</td></tr><tr><td>6.00 × 10<sup>-1</sup></td><td>1.88 × 10<sup>1</sup></td><td>Lipid synthesis</td></tr></tbody></table>

## 7.3 Extensions of Results

In order to obtain reduced metabolic systems by a rigorous procedure, many extensions of the results can be obtained. Particularly, considering more reactions between the metabolites with slow dynamics is possible (as long as these reactions are slow), without changing the equations of the fast part.
Hence, modifications in the reactions between slow metabolites do not alter
the equations of the metabolites in QSS, and the slow dynamics remain in
the reduced system. Moreover, the result obtained in Theorem 1 still holds
if the equations of the fast part are not changed.

Furthermore, effects such as inhibition can be considered in the slow
part of the system. For example, using the model of Haldane or feedback
inhibition in enzyme-catalyzed subnetworks.

In addition, models with more subnetworks of fast reactions, connected
by metabolites with slow dynamics, can be reduced and analyzed using the
present approach.

8 Conclusions

Quasi Steady State Assumption without verifying mathematical conditions
can lead to erroneous conclusion and strongly biased reduced systems [22, 23].
The aim of our work was to define the mathematical foundations of Quasi
Steady State Reduction for metabolic networks.

We reduced a general class of dynamical metabolic systems using time
scale separation and Tikhonov's Theorem. The considered models, include
Michaelis-Menten reaction rates and the possibility for some compounds to
accumulate. The reduction leads to a simpler model given by a small system
of differential equations: regardless the initial dimension of the network, we
end up with a low dimensional dynamical system, representing the dynamics
of the slow variables. The dilution due to growth plays an important role
and must not be neglected. It is worth noting that keeping the growth rate
in the equations strongly improves approximation precision and preserves
qualitative (stability) features of the original system.

We show that a metabolite in QSS has a concentration one order of magni-
tude lower than a metabolite in the slow part of the system. This is indirectly
a way to validate the hypotheses on the magnitude of the reaction kinetics.

Eventually, the calibration algorithm is very simple. It is remarkable
that the reduced model can predict all the fast compounds which have been
measured, regardless of the other compounds whose concentrations cannot
be recorded.

This approach covers a large class of metabolic enzymatic networks. But
more work remains to be done to treat further metabolic systems. For exam-
ple, networks with more reactions between fast and slow metabolites can be
studied in detail. Moreover, to obtain models that rigorously describe several
hierarchies in metabolic networks, systems with more than two time-scales
can be analyzed on the basis of the present paper.

Acknowledgments

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). The authors acknowledge the support of the Investissements d'Avenir Bio-informatique program for the project RESET (ANR-11-BINF-0005), the ANR project Maximic (ANR-17-CE40-0024-01), and the Inria Project Lab Algae in silico.

References

1. C. Baroukh, R. Muñoz-Tamayo, O. Bernard, and J. P. Steyer. Mathematical modeling of unicellular microalgae and cyanobacteria metabolism for biofuel production. *Current opinion in biotechnology*, 2015.

2. J. D. Orth, I. Thiele, and B. Ø. Palsson. What is flux balance analysis? *Nature biotechnology*, 28(3):245-248, 2010.

3. A. Provost, G. Bastin, S. N. Agathos, and Y. J. Schneider. Metabolic design of macroscopic bioreaction models: application to chinese hamster ovary cells. *Bioprocess and biosystems engineering*, 29(5-6):349-366, 2006.

4. C. Baroukh, R. Muñoz Tamayo, J. P. Steyer, and O. Bernard. DRUM: a new framework for metabolic modeling under non-balanced growth. application to the carbon metabolism of unicellular microalgae. *PloS one*, 9(8):e104499, 2014.

5. C. Baroukh, R. Muñoz-Tamayo, O. Bernard, and J. P. Steyer. Reply to the comment on “mathematical modeling of unicellular microalgae and cyanobacteria metabolism for biofuel production” by baroukh et al. [curr. opin. biotechnol. 2015, 33: 198-205]. *Current opinion in biotechnology*, 38:200, 2016.
6. P. Kokotović, H. K. Khalil, and J. O'Reilly. *Singular perturbation methods in control: analysis and design*. Society for Industrial and Applied Mathematics, Philadelphia, 1999.

7. J.D. Murray. *Mathematical Biology*, volume I. An Introduction. Springer, 2001.

8. L. A. Segel. *Modeling dynamic phenomena in molecular and cellular biology*. Cambridge University Press, 1984.

9. A. Narang-Siddarth and J. Valasek. *Nonlinear Time Scale Systems in Standard and Nonstandard Forms: Analysis and Control*. Society for Industrial and Applied Mathematics, Philadelphia, 2014.

10. G. Stephanopoulos, A. A. Aristidou, and J. Nielsen. *Metabolic engineering: principles and methodologies*. Academic press, San Diego, California, 1998.

11. L. Perko. *Differential equations and dynamical systems*. Texts in Applied Mathematics 7. Springer-Verlag, 1991.

12. R. Horn and C. R. Johnson. *Matrix Analysis*. Cambridge University Press, New York, 1985.

13. A. Berman and R. Plemmons. *Nonnegative matrix in the mathematical science*. Classics in Applied Mathematics. Society for Industrial and Applied Mathematics, Philadelphia, 1994.

14. J. A. Jacquez and C. P. Simon. Qualitative theory of compartmental systems. *Siam Review*, 35(1):43–79, 1993.

15. D. Rotin, B. Robinson, and I. F. Tannock. Influence of hypoxia and an acidic environment on the metabolism and viability of cultured cells: potential implications for cell death in tumors. *Cancer research*, 46(6):2821–2826, 1986.

16. F. Menegus, L. Cattaruzza, A. Chersi, and G. Fronza. Differences in the anaerobic lactate-succinate production and in the changes of cell sap pH for plants with high and low resistance to anoxia. *Plant Physiology*, 90(1):29–32, 1989.
17. C. D. Meyer. *Matrix analysis and applied linear algebra*. Society for Industrial and Applied Mathematics, Philadelphia, 2000.

18. E. Walter and L. Pronzato. *Identification of parametric models from experimental data*. Springer Verlag, 1997.

19. H. K. Khalil. *Nonlinear systems*. Prentice Hall, New Jersey, third edition, 2002.

20. F. Verhulst. Singular perturbation methods for slow-fast dynamics. *Non-linear Dynamics*, 50(4):747-753, 2007.

21. F. Hoppensteadt. On systems of ordinary differential equations with several parameters multiplying the derivatives. *Journal of Differential Equations*, 5(1):106-116, 1969.

22. S. Boie, V. Kirk, J. Sneyd, and M. Wechselberger. Effects of quasi-steady-state reduction on biophysical models with oscillations. *Journal of theoretical biology*, 393:16-31, 2016.

23. A. Goeke, C. Schilli, S. Walcher, and E. Zerz. Computing quasi-steady state reductions. *Journal of Mathematical Chemistry*, 50(6):1495-1513, 2012.

24. P. Lancaster and M. Tismenetsky. *The theory of matrices: with applications*. Academic Press, San Diego, California, second edition, 1985.

25. D. S. Bernstein. *Matrix mathematics: Theory, facts, and formulas with application to linear systems theory*. Princeton University Press, New Jersey, 2005.

A Michaelis-Menten Reaction

In this paper we present a metabolic network which contains enzymatic re-
actions. Therefore, we present the Michaelis-Menten enzymatic reaction to
set the notation that we use throughout the text.

The Michaelis-Menten model considers a substrate $X_i$ which reacts with
an enzyme $e_{ji}$ to produce a complex $C_{ji}$. Then, this complex is transformed
into a product $X_j$ and the enzyme $e_{ji}$. This enzymatic reaction is abstracted as follows:

$$X_i + e_{ji} \underset{k_{-1}^{ji}}{\overset{k_1^{ji}}{\rightleftharpoons}} C_{ji} \xrightarrow{k_{ji}} X_j + e_{ji} \quad (31)$$

$$\frac{dX_i}{dt} = -k_1^{ji} \cdot e_{ji} X_i + k_{-1}^{ji} \cdot C_{ji} \qquad X_i(0) = x_i^0$$

$$\frac{de_{ji}}{dt} = -k_1^{ji} \cdot e_{ji} X_i + (k_{-1}^{ji} + k_{ji}) C_{ji} \qquad e_{ji}(0) = e_{ji}^0$$

$$\frac{dC_{ji}}{dt} = k_1 ji \cdot e_{ji} X_i - (k_{-1}^{ji} + k_{ji}) C_{ji} \qquad C_{ji}(0) = 0$$

$$\frac{dX_j}{dt} = k_{ji} \cdot C_{ji} \qquad X_j(0) = x_j^0.$$

It is necessary to justify the Quasi Steady State Approximation for the Michaelis-Menten model. For example, this hold if the initial substrate concentration $x_i^0$ is sufficiently large compared with the initial enzyme concentration $e_{ji}^0$ [7], or if the product formation rate $k_{ji}$ is enough small [8]. A widely used Quasi Steady State Reduction of system (31) is the following [7, 8]:

$$\frac{dX_i}{dt} = -e_{ji}^0 k_{ji} \frac{X_i}{X_i + K_{ji}} \quad (32)$$

$$\frac{dX_j}{dt} = e_{ji}^0 k_{ji} \frac{X_i}{X_i + K_{ji}}$$

$$C_{ji} = \frac{e_{ji}^0 \cdot X_i}{X_i + K_{ji}}$$

$$e_{ji} = e_{ji}^0 - C_{ji},$$

where

$$K_{ji} := \frac{k_{-1}^{ji} + k_{ji}}{k_1^{ji}}$$

is the Michaelis-Menten constant.

## B Boundary layer

A second condition related to the uniform convergence of approximations when $\eta \to 0$ has to be verified with the boundary layer of Equation (14)
[6]. For this we define the boundary layer correction $\hat{Y}(\tau) = Y(t) - \bar{Y}(t)$, $\tau = t/\eta$, and the boundary layer problem:

$$ \frac{d\hat{Y}_2}{d\tau} = e_{21}^0 k_{21} \frac{x_1^0}{x_1^0 + K_{21}} + \left[ \sum_{j=3}^{m-1} \frac{e_{2j}^0 k_{2j}}{K_{2j}} (\hat{Y}_j + \bar{Y}_j(0)) \right] - \left[ \sum_{i=3}^{m-1} \frac{e_{i2}^0 k_{i2}}{K_{i2}} + \tilde{\mu} \right] (\hat{Y}_2 + \bar{Y}_2(0)) \quad (33) $$

$$ \frac{d\hat{Y}_{m-1}}{d\tau} = \left[ \sum_{j=2}^{m-2} \frac{e_{m-1,j}^0 k_{m-1,j}}{K_{m-1,j}} (\hat{Y}_j + \bar{Y}_j(0)) \right] - \left[ \sum_{\substack{i=2 \\ i \neq m-1}}^{m} \frac{e_{i,m-1}^0 k_{i,m-1}}{K_{i,m-1}} + \tilde{\mu} \right] (\hat{Y}_{m-1} + \bar{Y}_{m-1}(0)) $$

$$ \frac{d\hat{Y}_{m+1}}{d\tau} = e_{m+1,m}^0 k_{m+1,m} \frac{x_m^0}{x_m^0 + K_{m+1,m}} + \left[ \sum_{j=m+2}^{n-1} \frac{e_{m+1,j}^0 k_{m+1,j}}{K_{m+1,j}} (\hat{Y}_j + \bar{Y}_j(0)) \right] - \left[ \sum_{i=m+2}^{n-1} \frac{e_{i,m+1}^0 k_{i,m+1}}{K_{i,m+1}} + \tilde{\mu} \right] (\hat{Y}_{m+1} + \bar{Y}_{m+1}(0)) $$

$$ \frac{d\hat{Y}_{n-1}}{d\tau} = \left[ \sum_{j=m+1}^{n-2} \frac{e_{n-1,j}^0 k_{n-1,j}}{K_{n-1,j}} (\hat{Y}_j + \bar{Y}_j(0)) \right] - \left[ \sum_{\substack{i=m+1 \\ i \neq n-1}}^{n} \frac{e_{i,n-1}^0 k_{i,n-1}}{K_{i,n-1}} + \tilde{\mu} \right] (\hat{Y}_{n-1} + \bar{Y}_{n-1}(0)) $$

with initial conditions $\hat{y}_i^0(0) = y_i^0 - \bar{Y}_i(0)$ for every $i = 2, \dots, n-1, i \neq m$.

**Proposition 5.** The equilibrium point $\hat{Y} = \bar{0}$ of system (33) is asymptotically stable.

*Proof.* First notice that system (33) is linear, since (14) is linear. That $\hat{Y} = \bar{0}$ is an equilibrium point of system (33) is a consequence of Equation (18). Moreover, the Jacobian matrix of system (33) is the same that (14). Therefore, as in the proof of Proposition 1, we conclude that the origin is asymptotically stable for system (33). □
On the other hand, the boundary layer correction $\hat{Y}$ allows to correct the error of the approximation (8) at the initial fast transition. Indeed, notice that the initial condition $y_i^0$ in (12) can be different from $\bar{Y}_i(0)$ in (8). But

$$ \bar{Y}_i(0) + \hat{Y}_i(0) = \bar{Y}_i(0) + Y_i(0) - \bar{Y}_i(0) = y_i^0 $$

Moreover, the boundary layer correction $\hat{Y}$ vanishes quickly [6] since

$$ \lim_{\tau \to \infty} \hat{Y}(\tau) = \lim_{\eta \to 0} (Y(t) - \bar{Y}(t)) = 0. $$

# C Solution of the Slow System

*Proof of Proposition 2.* As in the proof of Proposition 1, we use the fact that $I(t)$ is a nonnegative continuous function on $[0, T_1]$ and all the parameters in (7) are nonnegative real numbers. Hence, system (7) is positively invariant in $\mathbb{R}_+^3$.

Let us denote $F(t, X)$ the right hand of Equation (7). Then, we have that $F$ and $\frac{\partial F}{\partial X}$ are continuous on $[0, T_1] \times \mathbb{R}_+^3$. Moreover, $\frac{\partial F}{\partial X}$ is uniformly bounded on $[0, T_1] \times \mathbb{R}_+^3$.

As a consequence, we can deduce from the Global Existence and Uniqueness Theorem [19] that (7) has a unique solution $\bar{X}(t)$ over $[0, T_1]$. □

# D Supplement for the Proof of Theorem 1

## D.1 Fluxes, Traps and Flux Traps

In order to see when the metabolites in QSS do not accumulate, we have to introduce the following definitions.

*Definition 1.* We define a directed graph $\Gamma$ related to the network in Figure 1, equivalent to system (3), as follows: the substrates and products $X_i$, $i = 1, \dots, n$, are the nodes of the $\Gamma$. Then, if $e_{ji}^0 k_{ji} \neq 0$ (i.e. if there is a reaction with substrate $X_i$ and product $X_j$) there is an edge with initial node $X_i$ and final node $X_j$. In a similar way, we define the graph associated to a subsystem of (3), with metabolites $\{X_{i_1}, \dots, X_{i_l}\} \subset \{X_1, \dots, X_n\}$.

The concept of graph allows the following definitions:
**Definition 2. (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$.

**Definition 3. (Trap).** Consider a graph with set of vertices $V$ and a subset of this $T = \{X_{i_1}, \dots, X_{i_l}\} \subset V$, $n > l \ge 1$. We say that $T$ is a *trap* if

* for every vertex $X_{i_k} \in T$ there is no flux from $X_{i_k}$ to any metabolite of $V \setminus T$ and

* no $X_{i_k}$ is the initial vertex of an edge with final vertex $X_* \notin V$.

In this case, we also say that $X_{i_k}$ is in a *trap* for every $X_{i_k} \in T$.

**Definition 4. (Flux trap).** Consider a flux $F$ with initial vertex $X_1$ and final vertex $X_n$ in a graph with vertices $V = \{X_1, \dots, X_n\}$. We say that the graph has a *trap* for the flux $F$ if there is a subset $T_F = \{X_{i_1}, \dots, X_{i_l}\} \subset V \setminus \{X_1, X_n\}$, such that

* $T_F$ is a *trap* (hence, there is no flux from $X_{i_k}$ to $X_n$ for every vertex $X_{i_k} \in T_F$) and

* there is a flux from $X_1$ to $X_{i_k}$ for every vertex $X_{i_k} \in T_F$.

When it is clear which is the flux $F$ taken into account, we only say that the graph has a *flux trap*. We also say that $X_{i_k}$ is in a *flux trap* for every vertex $X_{i_k} \in T_F$.

If the graph associated to a network has a *flux*, *trap* or *flux trap*, we also say that the network has a *flux*, *trap* or *flux trap*, respectively.

## D.2 Matrix analysis

Consider the Jacobian matrix defined in (16). For the sake of simplicity we denote

$$l_j := \sum_{\substack{i=2 \\ i \neq j}}^{n} \frac{e_{ij}^0 k_{ij}}{K_{ij}} + \tilde{\mu}$$

$$l_{ij} := \frac{e_{ij}^0 k_{ij}}{K_{ij}},$$
where $k_{ij} = 0$ if there is no reaction from $X_j$ to $X_i$. Then,

$$J = K' := \begin{pmatrix} K'_1 & 0 \\ 0 & K'_2 \end{pmatrix},$$

where

$$K_1 = K'_1 := \begin{pmatrix} -l_2 & l_{23} & \dots & l_{2,m-1} \\ l_{32} & -l_3 & \dots & l_{3,m-1} \\ \vdots & \vdots & \ddots & \vdots \\ l_{m-1,2} & l_{m-1,3} & \dots & -l_{m-1} \end{pmatrix},$$

$$K_2 = K'_2 := \begin{pmatrix} -l_{m+1} & l_{m+1,m+2} & \dots & l_{m+1,n-1} \\ l_{m+2,m+1} & -l_{m+2} & \dots & l_{m+2,n-1} \\ \vdots & \vdots & \ddots & \vdots \\ l_{n-1,m+1} & l_{n-1,m+2} & \dots & -l_{n-1} \end{pmatrix}.$$

**Theorem 2.** Suppose that the graph associated to (3) satisfies Assumption 1. Consider the expression of the metabolites in QSS (8) and define

$$c_i := \mathbf{b}_i \cdot l_{21} \qquad i = 2, \dots, m-1$$

$$c_i := \mathbf{b}_i \cdot l_{m+1,m} \qquad i = m+1, \dots, m-1.$$

Then for every $i \in \{2, \dots, m-1, m+1, \dots, n-1\} \setminus \mathcal{I}_{F_T}$,

$$c_i = \mathcal{O}(1) \quad (\text{if } c_i \neq 0).$$

We recall that $i \in \{2, \dots, m-1, m+1, \dots, n-1\} \setminus \mathcal{I}_{F_T}$ means the metabolite $X_i$ is not in a flux trap.

Before proving Theorem 2, we demonstrate several propositions. The proof of Theorem 2 is in Appendix D.3. For this, we have to analyze the order of the parameters

$$c_i = \frac{1}{\det(K'_1)} C_{1,i-1} \cdot (-l_{21}) \qquad \forall i = 2, \dots, m-1,$$

$$c_i = \frac{1}{\det(K'_2)} C'_{1,i-m} \cdot (-l_{m+1,m}) \qquad \forall i = m+1, \dots, n-1,$$
where $C_{1,i-1}$ and $C'_{1,i-m}$ are the cofactors of $K'_1$ and $K'_2$, respectively.

**Lemma 3.** Consider a singular matrix $A$ of dimension $n \times n$ and $\varepsilon\mu > 0$. Suppose $a_{ij} = O(1)$ when $\varepsilon \to 0$, for every entry of $A$. Then

$$ \det(A - \varepsilon\mu \cdot I) \le O(\varepsilon\mu). $$

*Proof.* Define $f$ as the function

$$ f(\varepsilon) = \det(A - \varepsilon\mu \cdot I). $$

Since $a_{ij} = O(1)$ when $\varepsilon \to 0$ for every entry of $A$, $f$ is infinitely differentiable at 0. Then, considering its Taylor series around zero, it follows

$$ f(\varepsilon\mu) = f(0) + f^{(1)}(0) \cdot \varepsilon\mu + \frac{f^{(2)}(0)}{2} \cdot (\varepsilon\mu)^2 + \dots $$

But $f(0) = \det(A) = 0$ and $f^{(n)}(0) = O(1)$ when $\varepsilon \to 0$, for every $n \in \mathbb{N}$, as a consequence of the hypothesis on the orders of $A$ entries. We conclude that

$$ \begin{aligned} f(\varepsilon\mu) &= \varepsilon\mu \cdot (f^{(1)}(0) + \frac{f^{(2)}(0)}{2} \cdot (\varepsilon\mu) + \frac{f^{(3)}(0)}{3} \cdot (\varepsilon\mu)^2 + \dots) \\ &\le O(\varepsilon\mu) \quad \text{when } \varepsilon \to 0. \end{aligned} $$

□

**Lemma 4.** Suppose that $M$ is a column diagonally dominant matrix of size $n \times n$, such that $\det(M) \ne 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 $\operatorname{sgn}(\det(M)) = (-1)^n$.

*Proof.* Since $-M$ is nonsingular and column diagonally dominant, by the Theorem of Gershgorin, $-M$ is a positive stable matrix [12]. Then its inverse matrix is nonnegative [13] (i.e. each entry of $(-M)^{-1}$ is nonnegative). But

$$ -((-M)^{-1}) = (M)^{-1} = \frac{1}{\det(M)} \cdot C \le 0, $$
where

$$ C = \begin{pmatrix} C_{11} & C_{12} & \dots & C_{1n} \\ C_{21} & C_{22} & \dots & C_{2n} \\ \vdots & \vdots & & \vdots \\ C_{n1} & C_{n2} & \dots & C_{nn} \end{pmatrix}^T, $$

is the transpose matrix of cofactors of $M$ [24]. Then

$$ \frac{C_{ij}}{\det(M)} \le 0 \qquad \forall i, j = 1, \dots, 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 [24], have the same sign. Moreover, since all the principal minors of $-M$ are positive [13], then $\det(-M) > 0$. We conclude that

$$ \operatorname{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. $\square$

**Proposition 6.** Let

$$ M_n = \begin{pmatrix} -\sum_{i=2}^{n} l_{i1} - l_{*1} - \varepsilon\mu & \cdots & l_{1n} \\ \vdots & \ddots & \vdots \\ l_{n1} & \cdots & -\sum_{i=1}^{n-1} l_{in} - l_{n+1,n} - \varepsilon\mu \end{pmatrix}. \quad (34) $$

where $l_{*i} \ge 0$. Consider the directed graph $\Gamma(M_n)$ associated to $M_n$ as a graph with $n$ nodes $X_1, \dots, X_n$ and an edge with origin $X_i$ and final $X_j$ if $l_{ji} > 0$. Suppose that $\Gamma(M_n)$ has no traps and that $l_{n+1,n} > 0$. Then,

$$ \det(M_n) = (-1)^n \cdot O(l_{ij}^n). $$

*Proof.* Notice that an output from the $i$-th metabolite is equivalent to $l_{*i} > 0$. Here, without loss of generality, we begin by supposing that the $n$-th metabolite has an output. Then $l_{n+1,n} > 0$.

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

$$ M_2 = \begin{pmatrix} -l_{21} - \varepsilon\mu & l_{12} \\ l_{21} & -l_{12} - l_{32} - \varepsilon\mu \end{pmatrix} \quad (35) $$
of a system with two metabolites and one output. The determinant of $M_2$ is

$$ \det(M_2) = l_{21}(l_{32} + \varepsilon\mu) + \varepsilon\mu(l_{12} + l_{32} + \varepsilon\mu). $$

If $l_{21} \cdot l_{32} \neq 0$, then $\det(M_2) = \mathcal{O}(l_{ij}^2)$. We examine in which cases $l_{21} \cdot l_{32} = 0$. If $l_{32} = 0$, the system has no output, contrary to our hypothesis. On the other hand, $l_{21} = 0$ implies that $X_1$ is in a trap (see Figure 6). We conclude that $\det(M_2) = \mathcal{O}(l_{ij}^2)$. The case in dimension $n = 2$ with more than one output is verified immediately.

<figure><img src="image_10.png" /><figcaption>**Figure 6:** Possible scenarios where $l_{21} = 0$ in a system with two metabolites and one output. Both cases represent a trap in $X_1$.</figcaption></figure>

We make the following induction hypothesis: consider a graph $\Gamma(M_{n-1})$ of $n-1$ metabolites with no traps and one output at least. If $M_{n-1}$ is the matrix of size $(n-1) \times (n-1)$ associated to $\Gamma(M_{n-1})$, then $\det(M_{n-1}) = (-1)^{n-1} \cdot \mathcal{O}(l_{ij}^{n-1})$.

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 4. It holds

$$ \det(M_n) = -l_{n+1,n}C_{nn} + \left[ \sum_{j=1}^{n-1} l_{jn}C_{jn} - \left( \sum_{i=1}^{n-1} l_{in} + \varepsilon\mu \right)C_{nn} \right], \quad (36) $$

where $C_{jn} = (-1)^{j+n}(M_n)_{jn}$ are cofactors of $M_n$ [24].

Suppose that $l_{ni} = 0$ and $l_{*i} = 0$ for every $i \in \{1, \dots, n-1\}$. Then $X_n$ is isolated and the rest of metabolites $\{X_1, \dots, X_{n-1}\}$ form a trap. Hence, $l_{ni} > 0$ or $l_{*i} > 0$ for some $i \in \{1, \dots, n-1\}$ and we can apply the hypothesis of induction to deduce that

$$ C_{nn} = (-1)^{n-1} \cdot \mathcal{O}(l_{ji}^{n-1}). $$
On the other hand, the term in the squared brackets in (36) is the determinant of the matrix $(M_n + l_{n+1,n} \cdot \delta_{nn})$, where $\delta_{nn}$ is a matrix of size $n \times n$ with zero at every entry, except for in the entry $nn$ which is equal to 1.

If $l_{*,i} = 0$ for every $i = 1, \dots, n-1$, then

$$ \det(M_n + l_{n+1,n} \cdot \delta_{nn}) \le \mathcal{O}(\varepsilon\mu) $$

according to Lemma 3 and the statement of Proposition 6 is proved. In other case, suppose $l_{*,n-1} > 0$ without loss of generality. Hence, if we develop the determinant of $(M_n + l_{n+1,n} \cdot \delta_{nn})$ by the $n-1$-th column and we substitute in (36), we have

$$ \begin{aligned} \det(M_n) = & -l_{n+1,n}\mathcal{C}_{nn} - l_{*,n-1}(M_n + l_{n+1,n} \cdot \delta_{nn})_{n-1,n-1} \quad (37) \\ & + \left[ \sum_{\substack{j=1 \\ j \ne n-1}}^{n} l_{j,n-1}(M_n + l_{n+1,n} \cdot \delta_{nn})_{j,n-1} \right. \\ & \qquad \left. - (\sum_{i=1}^{n-2} l_{i,n-1} + \varepsilon\mu)(M_n + l_{n+1,n} \cdot \delta_{nn})_{n-1,n-1} \right], \end{aligned} $$

where $(M_n + l_{n+1,n} \cdot \delta_{nn})_{j,n-1}$ are minors of $(M_n + l_{n+1,n} \cdot \delta_{nn})$. Moreover, the matrix $(M_n + l_{n+1,n} \cdot \delta_{nn})$ satisfies the conditions of Lemma 4. Then, all its cofactors have the same sign. Particularly, $sgn((M_n + l_{n+1,n} \cdot \delta_{nn})_{n-1,n-1}) = (-1)^{n-1}$, and then

$$ \operatorname{sgn}(-l_{n+1,n}\mathcal{C}_{nn}) = \operatorname{sgn}(-l_{*,n-1}(M_n + l_{n+1,n} \cdot \delta_{nn})). $$

Once again, the term in square brackets in (37) is equal to $\det(M_n + l_{n+1,n} \cdot \delta_{nn} + l_{*,n-1} \cdot \delta_{n-1,n-1})$. We proceed as for $\det(M_n + l_{n+1,n} \cdot \delta_{nn})$ to extract the following term

$$ -l_{*,n-2}(M_n + l_{n+1,n} \cdot \delta_{nn} + l_{*,n-1} \cdot \delta_{n-1,n-1})_{n-2,n-2} $$

which has the same sign as $-l_{n+1,n}\mathcal{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 to $-\varepsilon\mu$. That is to say, if we define

$$ \tilde{M}_i := (M_n + l_{n+1,n} \cdot \delta_{nn} + \sum_{j=1}^{n-i} l_{*,n-j} \delta_{n-j,n-j}), $$
for every $i = 2, \dots, n$, where we define $\sum_{j=1}^{0} l_{*,n-j} \delta_{n-j,n-j} = 0$. Then

$$ \det(M_n) = -l_{n+1,n}C_{nn} - \sum_{i=2}^{n} l_{*,i-1}(\tilde{M}_i)_{i-1,i-1} + [\mathcal{O}((\varepsilon\mu)^k)], $$

with $(\tilde{M}_i)_{i-1,i-1}$ a principal minor of $\tilde{M}_i$ and the term in square brackets represents

$$ \det \begin{pmatrix} -\sum_{i=2}^{n} l_{i1} - \varepsilon\mu & l_{12} & \dots & l_{1n} \\ l_{21} & -\sum_{i \neq 2}^{n} l_{i2} - \varepsilon\mu & \dots & l_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ l_{n1} & l_{n2} & \dots & -\sum_{i=1}^{n-1} l_{in} - \varepsilon\mu \end{pmatrix}, $$

for some $0 < k$, according to Lemma 3. Moreover,

$$ \operatorname{sgn}(-l_{n+1,n}C_{nn}) = \operatorname{sgn}(-l_{*,i-1}(\tilde{M}_i)_{i-1,i-1}) = (-1)^n, $$

for every $i = 2, \dots, n$, as a consequence of Lemma 4. Therefore, we conclude

$$ \det(M_n) = (-1)^n \cdot \mathcal{O}(l_{ij}^n). $$

□

The goal of the following proposition is to define the order of some $M'$ minors, as required for the Proof of Theorem 2.

**Proposition 7.** Let us suppose that $M_n$ represents a graph 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 \begin{pmatrix} l_{21} & -(\sum_{i \neq 2}^{n} l_{i2} + \varepsilon\mu) & \dots & l_{2,n-1} \\ l_{31} & l_{32} & \dots & l_{3,n-1} \\ \vdots & \vdots & \ddots & \vdots \\ l_{n-1,1} & l_{n-1,2} & \dots & -(\sum_{i \neq n-1}^{n} l_{i,n-1} + \varepsilon\mu) \\ l_{n1} & l_{n2} & \dots & l_{n,n-1} \end{pmatrix} $$

Then

$$ 0 < (M_n)_{1n} = \mathcal{O}(l_{ij}^{n-1}). $$
*Proof.* The demonstration is by induction over the squared matrix size. For the case of a minor with dimension two we have:

$$ 
\begin{aligned}
\det \begin{pmatrix} l_{21} & -(\sum_{i \neq 2}^{3} l_{i2} + \varepsilon \mu) \\ l_{31} & l_{32} \end{pmatrix} &= l_{21}l_{32} + l_{31}(\sum_{\substack{i=1 \\ i \neq 2}}^{3} l_{i2} + \varepsilon \mu) \\
&= \mathcal{O}(l_{ij}^2),
\end{aligned}
$$

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 satisfies the hypothesis of this lemma after $x-1$ changes of columns, for $x = 1, 2, \dots, n-1$. Hence, applying the *induction hypothesis* to these minors, we obtain that they are quantities equal to $(-1)^{x-1} \cdot \mathcal{O}(l_{ij}^{n-2})$, where $x$ is the struck row index.

Since there is no *traps* by hypothesis, the minor obtained after omitting the first line and column and the last line and row of $M_n$ has a column which is strictly diagonally dominant. We can then apply Proposition 6 and conclude that it has order $(-1)^{n-2} \cdot \mathcal{O}(l_{ij}^{n-2})$.

Therefore, we conclude that the determinant $(M_n)_{1n}$ is the sum of positive quantities of order $\mathcal{O}(l_{ij}^{n-1})$:

$$ 
\begin{aligned}
0 < (M_n)_{1n} &= l_{21} \cdot \mathcal{O}(l_{ij}^{n-2}) + \dots \\
&\quad + (-1)^{x+1} (-1)^{x-1} l_{x1} \cdot \mathcal{O}(l_{ij}^{n-2}) + \dots \\
&\quad + (-1)^n (-1)^{n-2} l_{n1} \cdot \mathcal{O}(l_{ij}^{n-2}) \\
&= \mathcal{O}(l_{ij}^{n-1}).
\end{aligned}
$$

□

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 1.* When the graph $\Gamma(M_n)$ related to $M_n$ has no traps, the minor $(M_n)_{1x}$ has order $(-1)^{n-x} \cdot \mathcal{O}(l_{ij}^{n-1})$, for every $x = 1, \dots, n$.

Recall that in Assumption 1 we only take into consideration *flux traps*. For this reason, we also analyze the determinant of the matrix associated to
a system with *traps*. For instance, with the matrix $M_2$ defined in (35), if $\Gamma(M_2)$ has a *trap*, $l_{21} = 0$ and its determinant has order $O(\varepsilon\mu)$.

In general, we can expect that a graph $\Gamma(M_n)$ with a *trap* has a determinant with order $\varepsilon\mu$. As a consequence, the matrix $M_n$ is ill-conditioned. This happens because a *trap* implies a block of zeros in the matrix. Indeed, the j-th column of the matrix system represents the edges whose origin is the metabolite $X_j$. Then, if $X_j$ is in a *trap*, $l_{ij} = 0$ for every $i$ with $X_i$ out of the *trap*.

**Proposition 8.** Let $M_n$ a matrix defined as in (34). If $M_n$ has a trap, then

$$ \det(M_n) = \det(M') \cdot \det(T), $$

where $T$ and $M'$ are a square submatrices of $M_n$, which correspond to metabolites in a trap and to metabolites not in a trap, respectively.

*Proof.* If there is a *trap* in $\Gamma(M_n)$, the matrix $M_n$ is reducible [24]. Then, after the same number of interchanges of rows than columns, $M_n$ can be transformed in a square block triangular matrix (keeping the dominant diagonal structure):

$$ M_n = \begin{pmatrix} M' & 0 \\ * & T \end{pmatrix}, \qquad (38) $$

where $M'$ and $T$ are square submatrices that correspond to the metabolites which are not in a *trap* and the metabolites which are in a *trap*, respectively. Since the matrix in (38) is square block triangular, its determinant is the product of the determinants of the diagonal blocks [25]. $\square$

**Corollary 2.** If $\Gamma(M_n)$ has a trap, then

$$ \det(M_n) \le O(\varepsilon\mu). $$

*Proof.* The square block $T$ is equal to a singular matrix minus $\varepsilon\mu \cdot I$. Then, by Lemma 3, its determinant has order less or equal to $O(\varepsilon\mu)$. $\square$

If $C_{ij}$ is a cofactor of $M_n$ and $\det(T)$ has order $\varepsilon\mu$, then the coefficients

$$ \frac{C_{1i}}{\det(M') \cdot \det(T)} \cdot (-l_{21}) $$

can be affected by a factor of order $(\varepsilon\mu)^{-1}$. However, in the following proposition we prove that when there is a trap $T$, $\det(T)$ is also a factor of the cofactor $C_{1i}$ if $X_i$ is not in the trap.
Proposition 9. Let $M_n$ be a matrix defined as in (34) and $\mathbf{F}$ a flux from $X_1$ to $X_n$. If $M_n$ has traps (not reached by $\mathbf{F}$) or flux traps for $\mathbf{F}$, then $M_n$ has the form

$$ M_n = \begin{pmatrix} [M']_{r \times r} & [C_1]_{r \times s} & 0_{r \times p} & 0_{r \times q} \\ 0_{s \times r} & [C_2]_{s \times s} & 0_{s \times p} & 0_{s \times q} \\ 0_{p \times r} & [C_3]_{p \times s} & [\mathbf{T}]_{p \times p} & 0_{p \times q} \\ [*]_{q \times r} & [C_4]_{q \times s} & 0_{q \times p} & [\mathbf{T}_F]_{q \times q} \end{pmatrix}, \quad r + s + p + q = n, \quad (39) $$

where $M'$ is a matrix with no traps, $\mathbf{T}$ is the square block corresponding to metabolites in traps not reached by $\mathbf{F}$, $\mathbf{T}_F$ to metabolites which are in flux traps and $C_2$ to metabolites that connect the traps to the rest of the network but which do not have a flux from the input. Then,

$$ \det(M_n) = \det(M') \cdot \det(C_2) \cdot \det(\mathbf{T}) \cdot \det(\mathbf{T}_F). $$

Furthermore, its minors satisfy

$$ (M_n)_{1j} = (M')_{1j} \cdot \det(C_2) \cdot \det(\mathbf{T}) \cdot \det(\mathbf{T}_F) \quad \forall j = 1, \dots, r, $$

with $(M')_{1j}$ a minor of $M'$, and

$$ (M_n)_{1j} = 0 \quad \forall j = r + 1, \dots, n - q. $$

*Note.* Notice that the block $[*]_{q \times r}$ is different from zero if there is a flux from $X_1$ to the flux trap $(\mathbf{T}_F)$.

*Proof.* Since $M_n$ defined in (39) is a square block triangular matrix, its determinant is the product of the determinants of the diagonal blocks [25]. Then,

$$ \begin{aligned} \det(M_n) &= \begin{vmatrix} M' & C_1 \\ 0 & C_2 \end{vmatrix} \cdot \det(\mathbf{T}) \cdot \det(\mathbf{T}_F) \\ &= \det(M') \cdot \det(C_2) \cdot \det(\mathbf{T}) \cdot \det(\mathbf{T}_F). \end{aligned} $$

For $j = 1, \dots, r$, the submatrix obtained from deleting the first row and the $j$-th column of $M_n$ is also a square block triangular matrix. Then, its determinant is

$$ \begin{aligned} (M_n)_{1j} &= \begin{vmatrix} (M')_{1j} & [C'_1] \\ 0 & [C_2] \end{vmatrix} \cdot \det(\mathbf{T}) \cdot \det(\mathbf{T}_F) \\ &= (M')_{1j} \cdot \det(C_2) \cdot \det(\mathbf{T}) \cdot \det(\mathbf{T}_F), \end{aligned} $$
where $(M')_{1j}$ is a minor of $M'$ and $[C'_1]$ is the matrix $C_1$ without its first row. On the other hand, for $j = r + 1, \dots, (r + s + p)$, the minor $(M_n)_{1j}$ is also the determinant of a square block triangular matrix, i.e.

$$ (M_n)_{1j} = \begin{pmatrix} M' & C_1 & 0 \\ 0 & C_2 & 0 \\ 0 & C_3 & T \end{pmatrix}_{1j} \cdot \det(T_F). $$

On the other hand, we have that the minor

$$ \begin{pmatrix} M' & C_1 & 0 \\ 0 & C_2 & 0 \\ 0 & C_3 & T \end{pmatrix}_{1j} = 0 \qquad \forall j = r+1, \dots, (r+s+p), $$

as a consequence of the block of zeros below $M'$. We conclude

$$ (M_n)_{ij} = 0 \qquad \forall j = r + 1, \dots, (r + s + p). $$

□

Finally, to analyze the minors of $M'$, the block of $M_n$ corresponding to the subgraph with no traps, we refer to Proposition 7 and Corollary 1.

## D.3 Proof of Theorem 2

*Proof of Theorem 2.* Since $K'_1$ is a nonsingular matrix,

$$ (K'_1)^{-1} = \frac{1}{\det(K'_1)} \cdot C, $$

where $C$ is the transpose matrix of cofactors of $K'_1$ [24] (i.e. $C_{ji} = (-1)^{j+i}(K'_1)_{ji}$). We then have according to Equation (18)

$$ \bar{Y}_i = \frac{1}{\det(K'_1)} C_{1,i-1} \cdot (-e_{21}^0 k_{21} \frac{\bar{X}_1}{\bar{X}_1 + K_{21}}) $$

Then, by definition of $c_i$,

$$ c_i = \frac{1}{\det(K'_1)} C_{1,i-1} \cdot (-l_{21}). $$
If $K'_1$ has no traps (i.e. the subnetwork with metabolites $X_2, \dots, X_{m-1}$ has no traps), then

$$ \det(K'_1) = (-1)^{m-2} \cdot \mathcal{O}(l_{ij}^{m-2}), $$

as stated by Proposition 6. Moreover, Corollary 1 implies that the cofactors $C_{1,i-1}$ have order

$$ C_{1,i-1} = (-1)^{m-1} \cdot \mathcal{O}(l_{ij}^{m-3}). $$

On the other hand, if $K'_1$ has a trap T not reached by the flux or a flux trap T_F, as a consequence of Corollary 1, Propositions 6 and 9,

$$ \frac{C_{1,i-1}}{\det(K'_1)} = (-1) \cdot \mathcal{O}(l_{ij}^{-1}) \quad \text{if } X_i \notin T_F, C_{1,i-1} \neq 0, $$

$$ \frac{C_{1,i-1}}{\det(K'_1)} = 0 \quad \text{if } X_i \in T. $$

We conclude that

$$ \frac{-l_{21} \cdot C_{1,i-1}}{\det(K'_1)} = \mathcal{O}(1) $$

if $C_{1,i-1} \neq 0$, for $i \in \{2, \dots, m-1\} \setminus \mathcal{I}_{T_F}$. The same reasoning applies for $K'_2$ and the variables of the second subnetwork $X_{m+1}, \dots, X_{n-1}$. $\square$