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

Authors: Claudia López Zazueta; Olivier Bernard; Jean-Luc Gouzé
Journal: Complexity
 (2018)

---

## 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 in which metabolic fluxes are no more constant 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 nonlinear metabolic rates. Here, we propose a rigorous approach to reduce metabolic models using quasi-steady-state reduction based on Tikhonov’s theorem, with a characterized and bounded reduction error. We assume that the metabolic network can be represented with Michaelis-Menten enzymatic reactions that evolve at different time scales. In this simplest approach, some metabolites can accumulate. We consider the case with a continuous varying input in the model, such as light for microalgae, so that the system is never at a steady state. Furthermore, our analysis proves that metabolites in the slow part of the metabolic system reach higher concentrations (by one order of magnitude) than metabolites in the fast part under some flux conditions. A simple example illustrates our approach and the resulting accuracy of the reduction method.


---

## Body

## 1. 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 hundreds 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 necessary.The most classical hypothesis is balanced growth, that is, 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 CO2 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 a steady state, and the classical approaches based on balanced growth hypothesis cannot be used to describe their metabolisms.Here, we propose a rigorous mathematical approach to reduce the dimension of a dynamical metabolic system, in order to analyze its behavior and calibrate it. The reduction that we propose allows to characterize the approximation 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 nonlinear and high-dimensional systems, which make their dynamical behavior difficult to determine and calibrate. 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 dimensionn, where the enzymatic reaction 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 characterized and bounded.In Section2, 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 3, we develop a mathematical model for these metabolic systems.In Section4, 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’s theorem for singularly perturbed systems are verified and we end up with a reduced dynamical model and a bound of the approximation error.In Section5, we prove that metabolites in QSS have a concentration one order of magnitude lower than slow metabolites. Additionally, in Section 6, 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 Section7. This simple model is forced by a periodic input and includes standard bricks in metabolic networks: combination of reversible and nonreversible reactions, with chains and cycles.


## 2. Network of Enzymatic Reactions

In this section, we present the class of metabolic networks studied all over the paper, which are illustrated in Figure1. 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.Figure 1
System of enzymatic reactions. An arrow fromXi to Xj represents a Michaelis-Menten reaction catalyzed by an enzyme eji, with substrate Xi, product Xj, and product formation rate kji or kji/ε. Fast reactions are within two subnetworks, which are interconnected by the metabolites X1, Xm, and Xn. 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.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 metabolitesX2,…,Xm−1,Xm+1,…,Xn−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.


### 2.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 Figure1. 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:(1)dXdt=Ft,X,Y,η,X0=x0,dYdt=1ηGt,X,Y,η,Y0=y0,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 dynamic metabolites Xfast in the following model:
(2)Y≔Xfastη.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):
(3)dXdt=Ft,X,Y¯,0,X0=x0,where Y¯ is a root of the equation
(4)0=Gt,X,Y¯,0.IfX¯ is a solution for (3), the quasi-steady-state approximation X¯,Y¯ to the solution of (1) satisfies
(5)X=X¯+Oη,Y=Y¯+Oη,after an initial fast transient. In other words, the error of the quasi-steady-state approximation has order of magnitude η, 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 (Figure1) is showed from Section 3 to Section 4.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 that of the metabolites with slow dynamics, that is,(6)η⋅Y=Xfast≤Oη⋅X.The conditions under which this assertion holds are given in Section5.


## 2.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 Figure1. 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:(1)dXdt=Ft,X,Y,η,X0=x0,dYdt=1ηGt,X,Y,η,Y0=y0,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 dynamic metabolites Xfast in the following model:
(2)Y≔Xfastη.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):
(3)dXdt=Ft,X,Y¯,0,X0=x0,where Y¯ is a root of the equation
(4)0=Gt,X,Y¯,0.IfX¯ is a solution for (3), the quasi-steady-state approximation X¯,Y¯ to the solution of (1) satisfies
(5)X=X¯+Oη,Y=Y¯+Oη,after an initial fast transient. In other words, the error of the quasi-steady-state approximation has order of magnitude η, 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 (Figure1) is showed from Section 3 to Section 4.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 that of the metabolites with slow dynamics, that is,(6)η⋅Y=Xfast≤Oη⋅X.The conditions under which this assertion holds are given in Section5.


## 3. Considered Class of Networks

In this section, we describe the systems of the network class considered in this work. Then, in Section4 and Section 5, 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.


### 3.1. Notations

Consider the network ofn enzymatic reactions depicted in Figures 1 and 2, where an arrow from Xi to Xj represents an enzymatic reaction catalyzed by eji, with substrate Xi, product Xj, and product formation rate kji or kji/ε. Then, every enzymatic reaction can be described with the Michaelis-Menten model (see Appendix A).Figure 2
Enzymatic reactions between metabolites in QSS depicted in Figure1. The metabolites inside the subnetworks are substrates or products of fast reactions catalyzed by an enzyme.However, it is necessary 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 concentrationxi0 is sufficiently large compared with the initial enzyme concentration eji0 [7] or if the product formation rate kji 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 metabolitesX1, Xm, and Xn. 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 reaction rates by(7)kjiεand slow reaction rates by
(8)kji,where ε is a small positive number.Additionally, a continuously varying nonnegative inputIt (e.g., the CO2 uptake in a microalgae submitted to light/dark cycles) and a growth rate μ>0, which acts as a dilution factor, are taken into account for the models.


### 3.2. Dynamical Model

According to the standard quasi-steady-state reduction for Michaelis-Menten enzymatic reactions described in AppendixA, we write the ODE system for the model in Figure 1 as
(9)dXidt=Fit,X1,…,Xn,ε,μ,Xi0=xi0,where
(10)F1≔It−e021k21X1X1+K21−μX1,Fi≔∑j∈1,m,neij0kijXjXj+Kij+∑j=2j≠mn−1eij0kijεXjXj+Kij−∑j=1neji0kjiXiXi+Kji−μXi,for i=m,n and
(11)Fi≔∑j∈1,m,neij0kijXjXj+Kij+∑j=2j≠mn−1eij0kijεXjXj+Kij−∑j=1neji0kjiεXiXi+Kji−μXi,for every i=2,…,n−1, i≠m.The variableXi describes the ith metabolite cell concentration; It is a nonnegative continuous function; ε is a small positive number; eji0, kji, and Kji are nonnegative parameters; and μ>0 is the growth rate. When there is no reaction with substrate Xi and product Xj, we define kji=0 and also kii=0 for every i=1,…,n.Note.
In our model, we can include first-order (linear) reactions. In this case, instead of writing(12)eji0kji⋅XiXi+Kjioreji0kjiε⋅XiXi+Kjias for enzymatic reactions, we have to write
(13)eji0kjiXioreji0kjiεXi,respectively, in the algebraic equation (9). 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 thateji0kji (or eji0kji/ε for the fast reactions) and Kji are parameters related to the enzyme reaction with substrate Xi and product Xj. Indeed, eji0 is the initial enzyme concentration, kji (or kji/ε) is the product formation rate, and Kji>0 is the specific Michaelis-Menten constant defined as
(14)Kji≔k−1ji+kjik1ji(see Appendix A).An important preliminary property that the dynamical system (9) has to obey is that the concentration Xit has to remain nonnegative over the time if the initial conditions are nonnegative. In our model, this depends on the input It. This is stated in the following property:Property 1.
If the initial conditionxi0 is nonnegative for every i=1,…,n and It≥0 for every t∈0,T1, then system (9) is positively invariant in ℝ+n¯.Proof.
To verify this, we show that system (9) is positively invariant in ℝ+n¯ if It is nonnegative over any interval 0,T1.Recall that allKji is supposed to be positive and every parameter, eji0, kji, and μ, is nonnegative. Then, we have for any i=1,…,n
(15)FiX1,…,0i−thentry,…,Xn,t,ε,μ≥0,if Xj≥0 for every j=1,…,n, j≠i. Therefore, system (9) is positively invariant in ℝ+n¯.


### 3.3. Parameter Order of Magnitude

With our notations, to represent different time scales in the reactions, we fixε a small positive number highlighting the difference between the parameter scale orders. We suppose that the parameters eji0kji are of standard range, that is,
(16)eji0kji=O1,asε→0∀i,j∈1,…,n,where O denotes the Big O or Landau symbol. For the definition and some properties of O, we refer to [9].Also, we suppose that the inputIt has a magnitude not larger than the slow reactions. In other words,
(17)It=O1.The rate of growthμ is considered as a parameter smaller than any reaction rate (a standard hypothesis [10]). Here, we assume
(18)εµ=Oε.


## 3.1. Notations

Consider the network ofn enzymatic reactions depicted in Figures 1 and 2, where an arrow from Xi to Xj represents an enzymatic reaction catalyzed by eji, with substrate Xi, product Xj, and product formation rate kji or kji/ε. Then, every enzymatic reaction can be described with the Michaelis-Menten model (see Appendix A).Figure 2
Enzymatic reactions between metabolites in QSS depicted in Figure1. The metabolites inside the subnetworks are substrates or products of fast reactions catalyzed by an enzyme.However, it is necessary 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 concentrationxi0 is sufficiently large compared with the initial enzyme concentration eji0 [7] or if the product formation rate kji 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 metabolitesX1, Xm, and Xn. 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 reaction rates by(7)kjiεand slow reaction rates by
(8)kji,where ε is a small positive number.Additionally, a continuously varying nonnegative inputIt (e.g., the CO2 uptake in a microalgae submitted to light/dark cycles) and a growth rate μ>0, which acts as a dilution factor, are taken into account for the models.


## 3.2. Dynamical Model

According to the standard quasi-steady-state reduction for Michaelis-Menten enzymatic reactions described in AppendixA, we write the ODE system for the model in Figure 1 as
(9)dXidt=Fit,X1,…,Xn,ε,μ,Xi0=xi0,where
(10)F1≔It−e021k21X1X1+K21−μX1,Fi≔∑j∈1,m,neij0kijXjXj+Kij+∑j=2j≠mn−1eij0kijεXjXj+Kij−∑j=1neji0kjiXiXi+Kji−μXi,for i=m,n and
(11)Fi≔∑j∈1,m,neij0kijXjXj+Kij+∑j=2j≠mn−1eij0kijεXjXj+Kij−∑j=1neji0kjiεXiXi+Kji−μXi,for every i=2,…,n−1, i≠m.The variableXi describes the ith metabolite cell concentration; It is a nonnegative continuous function; ε is a small positive number; eji0, kji, and Kji are nonnegative parameters; and μ>0 is the growth rate. When there is no reaction with substrate Xi and product Xj, we define kji=0 and also kii=0 for every i=1,…,n.Note.
In our model, we can include first-order (linear) reactions. In this case, instead of writing(12)eji0kji⋅XiXi+Kjioreji0kjiε⋅XiXi+Kjias for enzymatic reactions, we have to write
(13)eji0kjiXioreji0kjiεXi,respectively, in the algebraic equation (9). 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 thateji0kji (or eji0kji/ε for the fast reactions) and Kji are parameters related to the enzyme reaction with substrate Xi and product Xj. Indeed, eji0 is the initial enzyme concentration, kji (or kji/ε) is the product formation rate, and Kji>0 is the specific Michaelis-Menten constant defined as
(14)Kji≔k−1ji+kjik1ji(see Appendix A).An important preliminary property that the dynamical system (9) has to obey is that the concentration Xit has to remain nonnegative over the time if the initial conditions are nonnegative. In our model, this depends on the input It. This is stated in the following property:Property 1.
If the initial conditionxi0 is nonnegative for every i=1,…,n and It≥0 for every t∈0,T1, then system (9) is positively invariant in ℝ+n¯.Proof.
To verify this, we show that system (9) is positively invariant in ℝ+n¯ if It is nonnegative over any interval 0,T1.Recall that allKji is supposed to be positive and every parameter, eji0, kji, and μ, is nonnegative. Then, we have for any i=1,…,n
(15)FiX1,…,0i−thentry,…,Xn,t,ε,μ≥0,if Xj≥0 for every j=1,…,n, j≠i. Therefore, system (9) is positively invariant in ℝ+n¯.


## 3.3. Parameter Order of Magnitude

With our notations, to represent different time scales in the reactions, we fixε a small positive number highlighting the difference between the parameter scale orders. We suppose that the parameters eji0kji are of standard range, that is,
(16)eji0kji=O1,asε→0∀i,j∈1,…,n,where O denotes the Big O or Landau symbol. For the definition and some properties of O, we refer to [9].Also, we suppose that the inputIt has a magnitude not larger than the slow reactions. In other words,
(17)It=O1.The rate of growthμ is considered as a parameter smaller than any reaction rate (a standard hypothesis [10]). Here, we assume
(18)εµ=Oε.


## 4. Quasi-Steady-State Reduction

In this section, we propose a rigorous quasi-steady-state reduction (QSSR) of (9). 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 ε.We formally define the QSSR after Tikhonov’s theorem, of the metabolic network in Figure1 and system (9), as the following system of dimension three:
(19)dX¯1dt=It−e210k21X¯1X¯1+K21−μX¯1,(20)dX¯mdt=bm‐1⋅e210k21em,m−10km,m−1Km,m−1.X¯1X¯1+K21−em+1,m0km+1,mX¯mX¯m+Km+1,m−μX¯m,(21)dX¯ndt=bn‐1⋅em+1,m0km+1,men,n−10kn,n−1Kn,n−1⋅X¯mX¯m+Km+1,m−μ⋅X¯n,with initial conditions X¯10=x10, X¯m0=xm0, and X¯n0=xn0, and for the metabolites in QSS,
(22)X¯it=ε⋅bi⋅e210k21X¯1tX¯1t+K21,i=2,…,m−1,X¯it=ε⋅bi⋅em+1,m0km+1,mX¯mtX¯mt+Km+1,m,i=m+1,…,n−1,for every t∈0,T1. The definition of the parameter bi is given later in this section (see Proposition 1 and its proof).


### 4.1. Slow-Fast System

In order to write system (9) in the canonical form of singularly perturbed systems, we define a change of variable for the fast metabolites by
(23)Yi≔Xiε∀i=2,3,…,n−1,i≠m.Let us set the initial conditions for these new variables as(24)yi0≔xi0εand growth rate as
(25)µ~=εµ.Therefore, after the change of variables (23), (9) can be rewritten as follows:
(26)dXidt=Fit,X1,εY2,…,εYm−1,Xm,εYm+1,…,εYn−1,Xn,ε,μ,i=1,m,n,dYidt=1εFit,X1,εY2,…,εYm−1,Xm,εYm+1,…,εYn−1,Xn,ε,μ~,i=2,…,n−1,i≠m.Sinceε is a very small positive number, the dynamics of Yi are faster than those of Xi. Hence, the equations of Xi form the slow part of system (26), while the equations of Yi constitute its fast part.The previous (26) 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.


### 4.2. Canonical Form of Singularly Perturbed Systems

The slow-fast system (26) is in the class of singularly perturbed systems of the exact form:
(27)dX1dt=It−e210k21X1X1+K21−μX1,X10=x10,dXmdt=em,m−10km,m−1Ym−1ηYm−1+Km,m−1−em+1,m0km+1,mXmXm+Km+1,m−μXm,Xm0=xm0,dXndt=en,n−10kn,n−1Yn−1ηYn−1+Kn,n−1−μXn,Xn0=xn0,(28)dY2dt=1ηe210k21X1X1+K21+∑j=3m−1e2j0k2jYjηYj+K2j−∑i=3m−1ei20ki2Y2ηY2+Ki2−μ~Y2,⋮dYm−1dt=1η∑j=2m−2em−1,j0km−1,jYjηYj+Km−1,j−∑i=2i≠m−1mei,m−10ki,m−1Ym−1ηYm−1+Ki,m−1−μ~Ym−1,dYm+1dt=1ηem+1,m0km+1,mXmXm+Km+1,m+∑j=m+2n−1em+1,j0km+1,jYjηYj+Km+1,j−∑i=m+2n−1ei,m+10ki,m+1Ym+1ηYm+1+Ki,m+1−μ~Ym+1,⋮dYn−1dt=1η∑j=m+1n−2en−1,j0kn−1,jYjηYj+Kn−1,j−∑i=m+1i≠n−1nei,n−10ki,n−1Yn−1ηYn−1+Ki,n−1−μ~Yn−1,with initial conditions Yi0=yi0 for every i=2,…,n−1, i≠m.Note.
Equations (27)-(28) above are a more detailed expression of (26). Indeed, we obtain system (26) when η is substituted for ε in (27)-(28).An approximation to the solution of system (27)-(28) can be obtained considering the limit when η→0. Then, the dynamics in (28) are considered as fast, and the QSSA is applied to the metabolites Yi for every i=2,…,n−1, i≠m.Hereafter, we say that (27) is the slow part and (28) the fast part of system (9).


### 4.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 nonlinear systems). Then, we demonstrate that this steady state is asymptotically stable. Eventually, once all the conditions have been established, in Section 4.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 (28) and substituting η=0:
(29)0=eji0k21X1X1+K21+∑j=3m−1e2j0k2jYjK2j−∑i=3m−1ei20ki2Y2Ki2−μ~Y2⋮0=∑j=2m−2em−1,j0km−1,jYjKm−1,j−∑i=2i≠m−1mei,m−10ki,m−1Ym−1Ki,m−1−μ~Ym−10=em+1,m0km+1,mXmXm+Km+1,m+∑j=m+1n−1em+1,j0km+1,jYjKm+1,j−∑i=m+2n−1ei,m+10ki,m+1Ym+1Ki,m+1−μ~Ym+1⋮0=∑j=m+1n−2en−1,j0kn−1,jYjKn−1,j−∑i=m+1i≠n−1nei,n−10ki,n−1Yn−1Ki,n−1−μ~Yn−1.In order to apply Tikhonov’s theorem [6], we have to prove that (29) has an isolated root for any nonnegative constant values X1 and Xm and that this root is asymptotically stable for the following system:
(30)dY2dt=e210k21X1X1+K21+∑j=3m−1e2j0k2jK2jYj−∑i=3m−1ei20ki2Ki2+μ~Y2,⋮dYm−1dt=∑j=2m−2em−1,j0km−1,jKm−1,jYj−∑i=2i≠m−1mei,m−10ki,m−1Ki,m−1+μ~Ym−1,dYm+1dt=em+1,m0km+1,mXmXm+Km+1,m+∑j=m+2n−1em+1,j0km+1,jKm+1,jYj−∑i=m+2n−1ei,m+10ki,m+1Ki,m+1+μ~Ym+1,⋮dYn−1dt=∑j=m+1n−2en−1,j0kn−1,jKn−1,jYj−∑i=m+1i≠n−1nei,n−10ki,n−1Ki,n−1+μ~Yn−1.The purpose of finding a root of (29) is to write the fast variables Yi in terms of the slow variables Xi. In this case, it is possible to find an analytic solution of this algebraic system, because it is a linear equation for the variables Yi. Similarly, the asymptotical stability of this root for system (30) can be verified with the theory of linear systems of ODE.Proposition 1.
ConsiderX1 and Xm as nonnegative constant values. Then, system (30) has a single equilibrium point
(31)Y¯ii=2,…,n−1,i≠m,which is globally asymptotically stable. Moreover,
(32)Y¯i=bi⋅e210k21X1X1+K21,i=2,…,m−1,Y¯i=bi⋅em+1,m0km+1,mXmXm+Km+1,m,i=m+1,…,n−1,where bi∈ℝ+¯ are nonnegative coefficients.Proof.
First, notice that system (30) is a linear system for Yi under the hypotheses of Proposition 1. Then, we just have to show that its Jacobian matrix is stable, that is, that all its eigenvalues have negative real part [11].The Jacobian matrix of (30) is
(33)J=K100K2,where
(34)K1=−∑i=3m−1ei20ki2Ki2−μ~…e2,m−10k2,m−1K2,m−1e320k32K32…e3,m−10k3,m−1K3,m−1⋮⋮em−1,20km−1,2Km−1,2…−∑i=2i≠m−1mei,m−10ki,m−1Ki,m−1−μ~,(35)K2=−∑i=m+2n−1ei,m+10ki,m+1Ki,m+1−μ~⋯em+1,n−10km+1,n−1Km+1,n−1em+2,m+10km+2,m+1Km+2,m+1⋯em+2,n−10km+2,n−1Km+2,n−1⋮⋮en−1,m+10kn−1,m+1Kn−1,m+1⋯−∑i=m+1i≠n−1nei,n−10ki,n−1Ki,n−1−μ~.ButJ is a strictly column diagonally dominant matrix, because µ~>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 (29) is
(36)K100K2Y2⋮Ym−1Ym+1⋮Yn−1=−e210k21X1X1+K21⋮0em+1,m0km+1,mXmXm+Km+1,m⋮0.Then, the solution of the algebraic system (29) is
(37)Y¯2⋮Y¯m−1=−e210k21⋅K1−1⋅X1X1+K21⋮0,Y¯m+1⋮Y¯n−1=−em+1,m0km+1,m⋅K2−1XmXm+Km+1,m⋮0.Therefore, the variables of the solution can be written as(38)Y¯i=bi⋅e210k21X1X1+K21,i=2,…,m−1,Y¯i=bi⋅em+1,m0km+1,mXmXm+Km+1,m,i=m+1,…,n−1,with bi∈ℝ. Moreover, since Ki is strictly column diagonally dominant, by the theorem of Gershgorin, Ki is a stable matrix [12]. Then, its inverse matrix is nonpositive [13] (i.e., each entry of Ki−1 is nonpositive). Therefore, all entries in
(39)−e210k21⋅K1−1,−em+1,m0km+1,m⋅K2−1,are nonnegative. We conclude that coefficients bi in (32) are nonnegative.Note.
Although Proposition1 is proved for nonnegative constant values X1 and Xm, we consider Y¯i in (32) also as functions of t∈0,T1. Then, we have the functions in (22), defined for the QSSR.


### 4.4. Study of the Slowly Varying System

The dynamics of the slow system (metabolites which do accumulate) are obtained by settingη=0 in (27) and substituting the fast variables Yi for the expression given by (32):
(40)dX1dt=It−e210k21X1X1+K21−μX1,dXmdt=em,m−10km,m−1Ym−1Km,m−1−em+1,m0km+1,mXmXm+Km+1,m−μXm,dXndt=en,n−10kn,n−1Yn−1Kn,n−1−μXn.Then, we obtain the remaining dynamical system (19)–(21), which provides the dynamics to the overall network.The other variables of the metabolic network, which are the fast variablesYi (indeed, most of the variables are fast) can then be reconstructed after the solution of the algebraic system (29), given in (22). Finally, these fast variables rely on system (19)–(21). This system is also referred as the quasi-steady-state system [6].Proposition 2.
If system (19)–(21) has nonnegative initial conditions, then it has a unique nonnegative solution (X¯1,X¯m,X¯n) defined on the interval 0,T1.For the proof of Proposition2, see Appendix C.


### 4.5. Tikhonov’s Theorem

Propositions1 and 2 prove that the class of systems with the form (27)-(28) satisfies the hypothesis of Tikhonov’s theorem [6]. Then, we can apply this theorem to system (26).The following proposition is a consequence of Tikhonov’s theorem [6]. The proposition states that the approximation given by the QSSR (19)–(22) has an error with order Oε, after a fast initial transient for the fast variables.Proposition 3 (deduction of Tikhonov’s theorem).
IfIt is a nonnegative continuous function over 0,T1, then
(41)Xjt=X¯jt+Oε,j=1,m,n∀t∈0,T1.Moreover, there existsT0>0 such that for every t∈T0,T1,
(42)Xit=εY¯it+Oε∀i=2…,n−1,i≠m,where Xi are the solutions of the original system (9), X¯i are the solutions of (19)–(21), and Y¯i are the functions defined in (32).Note.
The solution of the boundary layer problem for system (27)-(28) is similar to that of system (30). We include its demonstration in Appendix B.


## 4.1. Slow-Fast System

In order to write system (9) in the canonical form of singularly perturbed systems, we define a change of variable for the fast metabolites by
(23)Yi≔Xiε∀i=2,3,…,n−1,i≠m.Let us set the initial conditions for these new variables as(24)yi0≔xi0εand growth rate as
(25)µ~=εµ.Therefore, after the change of variables (23), (9) can be rewritten as follows:
(26)dXidt=Fit,X1,εY2,…,εYm−1,Xm,εYm+1,…,εYn−1,Xn,ε,μ,i=1,m,n,dYidt=1εFit,X1,εY2,…,εYm−1,Xm,εYm+1,…,εYn−1,Xn,ε,μ~,i=2,…,n−1,i≠m.Sinceε is a very small positive number, the dynamics of Yi are faster than those of Xi. Hence, the equations of Xi form the slow part of system (26), while the equations of Yi constitute its fast part.The previous (26) 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.


## 4.2. Canonical Form of Singularly Perturbed Systems

The slow-fast system (26) is in the class of singularly perturbed systems of the exact form:
(27)dX1dt=It−e210k21X1X1+K21−μX1,X10=x10,dXmdt=em,m−10km,m−1Ym−1ηYm−1+Km,m−1−em+1,m0km+1,mXmXm+Km+1,m−μXm,Xm0=xm0,dXndt=en,n−10kn,n−1Yn−1ηYn−1+Kn,n−1−μXn,Xn0=xn0,(28)dY2dt=1ηe210k21X1X1+K21+∑j=3m−1e2j0k2jYjηYj+K2j−∑i=3m−1ei20ki2Y2ηY2+Ki2−μ~Y2,⋮dYm−1dt=1η∑j=2m−2em−1,j0km−1,jYjηYj+Km−1,j−∑i=2i≠m−1mei,m−10ki,m−1Ym−1ηYm−1+Ki,m−1−μ~Ym−1,dYm+1dt=1ηem+1,m0km+1,mXmXm+Km+1,m+∑j=m+2n−1em+1,j0km+1,jYjηYj+Km+1,j−∑i=m+2n−1ei,m+10ki,m+1Ym+1ηYm+1+Ki,m+1−μ~Ym+1,⋮dYn−1dt=1η∑j=m+1n−2en−1,j0kn−1,jYjηYj+Kn−1,j−∑i=m+1i≠n−1nei,n−10ki,n−1Yn−1ηYn−1+Ki,n−1−μ~Yn−1,with initial conditions Yi0=yi0 for every i=2,…,n−1, i≠m.Note.
Equations (27)-(28) above are a more detailed expression of (26). Indeed, we obtain system (26) when η is substituted for ε in (27)-(28).An approximation to the solution of system (27)-(28) can be obtained considering the limit when η→0. Then, the dynamics in (28) are considered as fast, and the QSSA is applied to the metabolites Yi for every i=2,…,n−1, i≠m.Hereafter, we say that (27) is the slow part and (28) the fast part of system (9).


## 4.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 nonlinear systems). Then, we demonstrate that this steady state is asymptotically stable. Eventually, once all the conditions have been established, in Section 4.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 (28) and substituting η=0:
(29)0=eji0k21X1X1+K21+∑j=3m−1e2j0k2jYjK2j−∑i=3m−1ei20ki2Y2Ki2−μ~Y2⋮0=∑j=2m−2em−1,j0km−1,jYjKm−1,j−∑i=2i≠m−1mei,m−10ki,m−1Ym−1Ki,m−1−μ~Ym−10=em+1,m0km+1,mXmXm+Km+1,m+∑j=m+1n−1em+1,j0km+1,jYjKm+1,j−∑i=m+2n−1ei,m+10ki,m+1Ym+1Ki,m+1−μ~Ym+1⋮0=∑j=m+1n−2en−1,j0kn−1,jYjKn−1,j−∑i=m+1i≠n−1nei,n−10ki,n−1Yn−1Ki,n−1−μ~Yn−1.In order to apply Tikhonov’s theorem [6], we have to prove that (29) has an isolated root for any nonnegative constant values X1 and Xm and that this root is asymptotically stable for the following system:
(30)dY2dt=e210k21X1X1+K21+∑j=3m−1e2j0k2jK2jYj−∑i=3m−1ei20ki2Ki2+μ~Y2,⋮dYm−1dt=∑j=2m−2em−1,j0km−1,jKm−1,jYj−∑i=2i≠m−1mei,m−10ki,m−1Ki,m−1+μ~Ym−1,dYm+1dt=em+1,m0km+1,mXmXm+Km+1,m+∑j=m+2n−1em+1,j0km+1,jKm+1,jYj−∑i=m+2n−1ei,m+10ki,m+1Ki,m+1+μ~Ym+1,⋮dYn−1dt=∑j=m+1n−2en−1,j0kn−1,jKn−1,jYj−∑i=m+1i≠n−1nei,n−10ki,n−1Ki,n−1+μ~Yn−1.The purpose of finding a root of (29) is to write the fast variables Yi in terms of the slow variables Xi. In this case, it is possible to find an analytic solution of this algebraic system, because it is a linear equation for the variables Yi. Similarly, the asymptotical stability of this root for system (30) can be verified with the theory of linear systems of ODE.Proposition 1.
ConsiderX1 and Xm as nonnegative constant values. Then, system (30) has a single equilibrium point
(31)Y¯ii=2,…,n−1,i≠m,which is globally asymptotically stable. Moreover,
(32)Y¯i=bi⋅e210k21X1X1+K21,i=2,…,m−1,Y¯i=bi⋅em+1,m0km+1,mXmXm+Km+1,m,i=m+1,…,n−1,where bi∈ℝ+¯ are nonnegative coefficients.Proof.
First, notice that system (30) is a linear system for Yi under the hypotheses of Proposition 1. Then, we just have to show that its Jacobian matrix is stable, that is, that all its eigenvalues have negative real part [11].The Jacobian matrix of (30) is
(33)J=K100K2,where
(34)K1=−∑i=3m−1ei20ki2Ki2−μ~…e2,m−10k2,m−1K2,m−1e320k32K32…e3,m−10k3,m−1K3,m−1⋮⋮em−1,20km−1,2Km−1,2…−∑i=2i≠m−1mei,m−10ki,m−1Ki,m−1−μ~,(35)K2=−∑i=m+2n−1ei,m+10ki,m+1Ki,m+1−μ~⋯em+1,n−10km+1,n−1Km+1,n−1em+2,m+10km+2,m+1Km+2,m+1⋯em+2,n−10km+2,n−1Km+2,n−1⋮⋮en−1,m+10kn−1,m+1Kn−1,m+1⋯−∑i=m+1i≠n−1nei,n−10ki,n−1Ki,n−1−μ~.ButJ is a strictly column diagonally dominant matrix, because µ~>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 (29) is
(36)K100K2Y2⋮Ym−1Ym+1⋮Yn−1=−e210k21X1X1+K21⋮0em+1,m0km+1,mXmXm+Km+1,m⋮0.Then, the solution of the algebraic system (29) is
(37)Y¯2⋮Y¯m−1=−e210k21⋅K1−1⋅X1X1+K21⋮0,Y¯m+1⋮Y¯n−1=−em+1,m0km+1,m⋅K2−1XmXm+Km+1,m⋮0.Therefore, the variables of the solution can be written as(38)Y¯i=bi⋅e210k21X1X1+K21,i=2,…,m−1,Y¯i=bi⋅em+1,m0km+1,mXmXm+Km+1,m,i=m+1,…,n−1,with bi∈ℝ. Moreover, since Ki is strictly column diagonally dominant, by the theorem of Gershgorin, Ki is a stable matrix [12]. Then, its inverse matrix is nonpositive [13] (i.e., each entry of Ki−1 is nonpositive). Therefore, all entries in
(39)−e210k21⋅K1−1,−em+1,m0km+1,m⋅K2−1,are nonnegative. We conclude that coefficients bi in (32) are nonnegative.Note.
Although Proposition1 is proved for nonnegative constant values X1 and Xm, we consider Y¯i in (32) also as functions of t∈0,T1. Then, we have the functions in (22), defined for the QSSR.


## 4.4. Study of the Slowly Varying System

The dynamics of the slow system (metabolites which do accumulate) are obtained by settingη=0 in (27) and substituting the fast variables Yi for the expression given by (32):
(40)dX1dt=It−e210k21X1X1+K21−μX1,dXmdt=em,m−10km,m−1Ym−1Km,m−1−em+1,m0km+1,mXmXm+Km+1,m−μXm,dXndt=en,n−10kn,n−1Yn−1Kn,n−1−μXn.Then, we obtain the remaining dynamical system (19)–(21), which provides the dynamics to the overall network.The other variables of the metabolic network, which are the fast variablesYi (indeed, most of the variables are fast) can then be reconstructed after the solution of the algebraic system (29), given in (22). Finally, these fast variables rely on system (19)–(21). This system is also referred as the quasi-steady-state system [6].Proposition 2.
If system (19)–(21) has nonnegative initial conditions, then it has a unique nonnegative solution (X¯1,X¯m,X¯n) defined on the interval 0,T1.For the proof of Proposition2, see Appendix C.


## 4.5. Tikhonov’s Theorem

Propositions1 and 2 prove that the class of systems with the form (27)-(28) satisfies the hypothesis of Tikhonov’s theorem [6]. Then, we can apply this theorem to system (26).The following proposition is a consequence of Tikhonov’s theorem [6]. The proposition states that the approximation given by the QSSR (19)–(22) has an error with order Oε, after a fast initial transient for the fast variables.Proposition 3 (deduction of Tikhonov’s theorem).
IfIt is a nonnegative continuous function over 0,T1, then
(41)Xjt=X¯jt+Oε,j=1,m,n∀t∈0,T1.Moreover, there existsT0>0 such that for every t∈T0,T1,
(42)Xit=εY¯it+Oε∀i=2…,n−1,i≠m,where Xi are the solutions of the original system (9), X¯i are the solutions of (19)–(21), and Y¯i are the functions defined in (32).Note.
The solution of the boundary layer problem for system (27)-(28) is similar to that of system (30). We include its demonstration in Appendix B.


## 5. Magnitude of Concentrations throughout the Metabolic Network

In this section, we study the magnitude of metabolite concentrations, depending if they are associated to slow or fast reactions. They are deduced from the reduced system after Tikhonov’s theorem (19)–(22). We now show that the concentration of metabolites in QSS (that do not trap the input flux) is one order of magnitude ε lower than that of metabolites with slow dynamics. In order to prove this assertion, we define the conditions under which
(43)bi⋅e210k21=O1,bi⋅em+1,m0km+1,m=O1,to obtain
(44)X¯it=Oε⋅X¯1tX¯1t+K21,X¯it=Oε⋅X¯mtX¯mt+Km+1,m,for every t∈T0,T1.


### 5.1. Parameter Orders

We show that all off-diagonal entries of the Jacobian matrixKi have the same order of magnitude, for both matrices defined in (34)-(35).Lemma 1.
Suppose that the parameters of each Michaelis-Menten enzymatic reaction (see AppendixA) satisfy
(45)Okγji=Okji∀i,j=1,…,n,γ∈−1,1,
then(46)Oeji0kjiKji=Oej′i′0kj′i′Kj′i′=Oej′i′0kj′i′∀i,j,i′,j′∈1,…,n.Proof.
As a consequence of (16),
(47)Oeji0kji=Oej′i′0kj′i′∀i,j,i′,j′∈1,…,n.Moreover, by the definition of the Michaelis-Menten constant (14) and (45), we have OKji=1 and then
(48)Oeji0kjiKji=Oeji0kji∀i,j.Hence, combining (47) and (48), we have
(49)Oeji0kjiKji=Oej′i′0kj′i′Kj′i′=Oej′i′0kj′i′∀i,j,i′,j′∈1,…,n.Actually, all the entries of the matrixKi have the same order of magnitude, as asserts the following corollary.Proposition 4.
Consider the matrices defined in (34)-(35). All the entries of K1 and K2 have the same order.Proof.
According to Lemma1, for the sums in the diagonal of the matrices, we have
(50)∑i=2i≠jneij0kijKij=Oeji0kji.Moreover,μ~≪eji0kji. Then,
(51)∑i=2i≠jneij0kijKij+μ~=Oeji0kji.For the off-diagonal entries, consider (48). Therefore, all the entries of K1 and K2 have the order Oeji0kji.


### 5.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 existsF a flux from X1 to Xn in the system of enzymatic reaction (9) (depicted in Figure 1). Moreover, we define ℐTF as the set of indices such that Xi is in a flux trap, for every i∈ℐTF, and Xj is not in a flux trap for everyj∈1,…,n\ℐTF.Notice that the presence of the fluxF from X1 to Xn implies
(52)1,m,n∩ℐTF=∅.Then, flux traps are only possible within the subnetworks in QSS. Also,ℐTF=∅ if there is no flux trap.The following lemma sets down the order of magnitude of the parameters in (22), 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 reaction (9) (Figure 1) under Assumption 1. Consider the parameters bi of (22), obtained in Section 4. Then, if bi≠0, it holds
(53)bi⋅e210k21=O1ifi∈2,…,m−1\ℐTF,bi⋅em+1,m0km+1,m=O1ifi∈m+1,…,n−1\ℐTF.Proof.
From the results stated in AppendixD, particularly Theorem 2, we have for bi≠0,
(54)bi⋅e210k21K21=O1ifi∈2,…,m−1\ℐTF,bi⋅em+1,m0km+1,mKm+1,m=O1ifi∈m+1,…,n−1\ℐTF.Using equality (48), we conclude that, for bi≠0,
(55)bi⋅e210k21=O1ifi∈2,…,m−1\ℐTF,bi⋅em+1,m0km+1,m=O1ifi∈m+1,…,n−1\ℐTF.The next theorem is a powerful conclusion obtained after the QSSR (19)–(22). Theorem 1 states that the concentration of a metabolite in QSS, which is not in a flux trap, is one order of magnitude ε lower than that of 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 (9) (Figure 1). Under Assumption 1, the following inequalities hold:
(56)X¯i≤Oε⋅X¯1∀i∈2,…,m−1\ℐTF,X¯i≤Oε⋅X¯m∀i∈m+1,…,n−1\ℐTF.Proof of Theorem 1.
Since the reduction from Tikhonov’s theorem, we have (22), that is,
(57)X¯i=ε⋅bi⋅e210k21X¯1X¯1+K21,i=2,…,m−1X¯i=ε⋅bi⋅em+1,m0km+1,mX¯mX¯m+Km+1,m,i=m+1,…,n−1.Then, as stated in Lemma2, for i such that bi≠0,
(58)X¯i=OεX¯1X¯1+K21,i∈2,…,m−1\ℐTF,X¯i=OεX¯mX¯m+Km+1,mifi∈m+1,…,n−1\ℐTF.But1≤OX¯i+Ki+1,i, because system (19)–(21) is positively invariant and OKi+1,i=1. Hence,
(59)OX¯1X¯1+K21≤OX¯1,OX¯mX¯m+Km+1,m≤OX¯m.We conclude that(60)X¯i≤Oε⋅X¯1∀i∈2,…,m−1\ℐTF,X¯i≤Oε⋅X¯m∀i∈m+1,…,n−1\ℐTF.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 Proposition3) 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.1. Parameter Orders

We show that all off-diagonal entries of the Jacobian matrixKi have the same order of magnitude, for both matrices defined in (34)-(35).Lemma 1.
Suppose that the parameters of each Michaelis-Menten enzymatic reaction (see AppendixA) satisfy
(45)Okγji=Okji∀i,j=1,…,n,γ∈−1,1,
then(46)Oeji0kjiKji=Oej′i′0kj′i′Kj′i′=Oej′i′0kj′i′∀i,j,i′,j′∈1,…,n.Proof.
As a consequence of (16),
(47)Oeji0kji=Oej′i′0kj′i′∀i,j,i′,j′∈1,…,n.Moreover, by the definition of the Michaelis-Menten constant (14) and (45), we have OKji=1 and then
(48)Oeji0kjiKji=Oeji0kji∀i,j.Hence, combining (47) and (48), we have
(49)Oeji0kjiKji=Oej′i′0kj′i′Kj′i′=Oej′i′0kj′i′∀i,j,i′,j′∈1,…,n.Actually, all the entries of the matrixKi have the same order of magnitude, as asserts the following corollary.Proposition 4.
Consider the matrices defined in (34)-(35). All the entries of K1 and K2 have the same order.Proof.
According to Lemma1, for the sums in the diagonal of the matrices, we have
(50)∑i=2i≠jneij0kijKij=Oeji0kji.Moreover,μ~≪eji0kji. Then,
(51)∑i=2i≠jneij0kijKij+μ~=Oeji0kji.For the off-diagonal entries, consider (48). Therefore, all the entries of K1 and K2 have the order Oeji0kji.


## 5.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 existsF a flux from X1 to Xn in the system of enzymatic reaction (9) (depicted in Figure 1). Moreover, we define ℐTF as the set of indices such that Xi is in a flux trap, for every i∈ℐTF, and Xj is not in a flux trap for everyj∈1,…,n\ℐTF.Notice that the presence of the fluxF from X1 to Xn implies
(52)1,m,n∩ℐTF=∅.Then, flux traps are only possible within the subnetworks in QSS. Also,ℐTF=∅ if there is no flux trap.The following lemma sets down the order of magnitude of the parameters in (22), 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 reaction (9) (Figure 1) under Assumption 1. Consider the parameters bi of (22), obtained in Section 4. Then, if bi≠0, it holds
(53)bi⋅e210k21=O1ifi∈2,…,m−1\ℐTF,bi⋅em+1,m0km+1,m=O1ifi∈m+1,…,n−1\ℐTF.Proof.
From the results stated in AppendixD, particularly Theorem 2, we have for bi≠0,
(54)bi⋅e210k21K21=O1ifi∈2,…,m−1\ℐTF,bi⋅em+1,m0km+1,mKm+1,m=O1ifi∈m+1,…,n−1\ℐTF.Using equality (48), we conclude that, for bi≠0,
(55)bi⋅e210k21=O1ifi∈2,…,m−1\ℐTF,bi⋅em+1,m0km+1,m=O1ifi∈m+1,…,n−1\ℐTF.The next theorem is a powerful conclusion obtained after the QSSR (19)–(22). Theorem 1 states that the concentration of a metabolite in QSS, which is not in a flux trap, is one order of magnitude ε lower than that of 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 (9) (Figure 1). Under Assumption 1, the following inequalities hold:
(56)X¯i≤Oε⋅X¯1∀i∈2,…,m−1\ℐTF,X¯i≤Oε⋅X¯m∀i∈m+1,…,n−1\ℐTF.Proof of Theorem 1.
Since the reduction from Tikhonov’s theorem, we have (22), that is,
(57)X¯i=ε⋅bi⋅e210k21X¯1X¯1+K21,i=2,…,m−1X¯i=ε⋅bi⋅em+1,m0km+1,mX¯mX¯m+Km+1,m,i=m+1,…,n−1.Then, as stated in Lemma2, for i such that bi≠0,
(58)X¯i=OεX¯1X¯1+K21,i∈2,…,m−1\ℐTF,X¯i=OεX¯mX¯m+Km+1,mifi∈m+1,…,n−1\ℐTF.But1≤OX¯i+Ki+1,i, because system (19)–(21) is positively invariant and OKi+1,i=1. Hence,
(59)OX¯1X¯1+K21≤OX¯1,OX¯mX¯m+Km+1,m≤OX¯m.We conclude that(60)X¯i≤Oε⋅X¯1∀i∈2,…,m−1\ℐTF,X¯i≤Oε⋅X¯m∀i∈m+1,…,n−1\ℐTF.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 Proposition3) 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].


## 6. 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 (19)–(21). 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 (22), to write its concentration as a linear combination of the slow metabolites.


### 6.1. Calibration of the Slow Dynamics

We suppose experimental data of the metabolites in the slow part of the system (27), denoted by
(61)Zitj=Xitj+βitj,i=1,m,n,j=1,…,r,where Xi is the solution of the original system (9) and βi represents an error of measurements. In order to estimate the parameters of the reduced system (19)–(21), we rewrite it as
(62)dX¯1dt=It−θ1X¯1X¯1+θ2−θ3⋅X¯1,X¯10=Z1t1,dX¯mdt=θ4X¯1X¯1+θ2−θ5X¯mX¯m+θ6−θ3X¯m,X¯m0=Zmt1,dX¯ndt=θ7X¯mX¯m+θ6−θ3⋅X¯n,X¯n0=Znt1.Letθ=θ1,θ2,θ3,θ4,θ5 and define a cost function ℱθ. This cost function has to measure the error between the solution of (62) and the data Z1,Zm,Zn, for every value θ in a domain D⊂ℝ7. For example, we can define ℱ as
(63)ℱθ=∑i∈1,m,n∑j=1rX¯iθ,tj−Zitj2.Then, we have to findθ^=θ^1,θ^2,θ^3,θ^4,θ^5 such that
(64)ℱθ^=minFθ:θ∈D.Note.
For obtaining the vector of parametersθ^ to calibrate (62), it is not necessary to have data of any metabolite in QSS, Xi with i=2,…,n−1, i≠m. Only the data (61) of the metabolites in the slow part, X1,Xm,Xn, is used.


### 6.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,(65)Zitj′=Xitj′+Ntj′,i=2,…,n−1,i≠m,T0≤t1′<⋯<tr′′,and that we have obtained θ^ after calibrating (62), we can estimate the parameters in (22), as a matter of fact, in line with the reduced system.In (19)–(22) and the calibrated system (62), for the metabolites in QSS, we have
(66)X¯i=αi⋅X¯1X¯1+θ^2,i=2,…,m−1,X¯i=αi⋅X¯mX¯m+θ^6,i=m+1,…,n−1,where αi are the parameters to be estimated.Here, we can explicitly resolve the linear least square problem. The least squares solution that minimizes the difference between the dataZi and the expressions in (66) is the following [17]:
(67)α^i=∑j=1r′Zitj′X¯1θ^,tj′/X¯1θ^,tj′+θ^2∑j=1r′X¯1θ^,tj′/X¯1θ^,tj′+θ^22∀i=2,…,m−1,α^i=∑j=1r′Zitj′X¯mθ^,tj′/X¯mθ^,tj′+θ^6∑j=1r′X¯mθ^,tj′/X¯mθ^,tj′+θ^62∀i=m+1,…,n−1.Indeed, we look for values ofα^i that minimize the differences
(68)Liα=∑j=1r′αX¯1tj′,θ^X¯1tj′,θ^+θ^2−Zitj′2,i=2,…,m−1,Liα=∑j=1r′αX¯mtj′,θ^X¯mtj′,θ^+θ^6−Zitj′2,i=m+1,…,n−1.Note.
To obtain the parameterα^i, we only need the data Zi (of the corresponding metabolites in QSS, Xi) and the calibrated system (62) with θ^.


## 6.1. Calibration of the Slow Dynamics

We suppose experimental data of the metabolites in the slow part of the system (27), denoted by
(61)Zitj=Xitj+βitj,i=1,m,n,j=1,…,r,where Xi is the solution of the original system (9) and βi represents an error of measurements. In order to estimate the parameters of the reduced system (19)–(21), we rewrite it as
(62)dX¯1dt=It−θ1X¯1X¯1+θ2−θ3⋅X¯1,X¯10=Z1t1,dX¯mdt=θ4X¯1X¯1+θ2−θ5X¯mX¯m+θ6−θ3X¯m,X¯m0=Zmt1,dX¯ndt=θ7X¯mX¯m+θ6−θ3⋅X¯n,X¯n0=Znt1.Letθ=θ1,θ2,θ3,θ4,θ5 and define a cost function ℱθ. This cost function has to measure the error between the solution of (62) and the data Z1,Zm,Zn, for every value θ in a domain D⊂ℝ7. For example, we can define ℱ as
(63)ℱθ=∑i∈1,m,n∑j=1rX¯iθ,tj−Zitj2.Then, we have to findθ^=θ^1,θ^2,θ^3,θ^4,θ^5 such that
(64)ℱθ^=minFθ:θ∈D.Note.
For obtaining the vector of parametersθ^ to calibrate (62), it is not necessary to have data of any metabolite in QSS, Xi with i=2,…,n−1, i≠m. Only the data (61) of the metabolites in the slow part, X1,Xm,Xn, is used.


## 6.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,(65)Zitj′=Xitj′+Ntj′,i=2,…,n−1,i≠m,T0≤t1′<⋯<tr′′,and that we have obtained θ^ after calibrating (62), we can estimate the parameters in (22), as a matter of fact, in line with the reduced system.In (19)–(22) and the calibrated system (62), for the metabolites in QSS, we have
(66)X¯i=αi⋅X¯1X¯1+θ^2,i=2,…,m−1,X¯i=αi⋅X¯mX¯m+θ^6,i=m+1,…,n−1,where αi are the parameters to be estimated.Here, we can explicitly resolve the linear least square problem. The least squares solution that minimizes the difference between the dataZi and the expressions in (66) is the following [17]:
(67)α^i=∑j=1r′Zitj′X¯1θ^,tj′/X¯1θ^,tj′+θ^2∑j=1r′X¯1θ^,tj′/X¯1θ^,tj′+θ^22∀i=2,…,m−1,α^i=∑j=1r′Zitj′X¯mθ^,tj′/X¯mθ^,tj′+θ^6∑j=1r′X¯mθ^,tj′/X¯mθ^,tj′+θ^62∀i=m+1,…,n−1.Indeed, we look for values ofα^i that minimize the differences
(68)Liα=∑j=1r′αX¯1tj′,θ^X¯1tj′,θ^+θ^2−Zitj′2,i=2,…,m−1,Liα=∑j=1r′αX¯mtj′,θ^X¯mtj′,θ^+θ^6−Zitj′2,i=m+1,…,n−1.Note.
To obtain the parameterα^i, we only need the data Zi (of the corresponding metabolites in QSS, Xi) and the calibrated system (62) with θ^.


## 7. Illustrative Example with a Toy Enzymatic Network

In this section, we apply the method developed in this paper to the toy network represented in Figure3. 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 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 nonaccumulative and they are supposed to be in quasi-steady state. Every reaction is catalyzed by an enzymeeji.First, we consider the ODE of the toy enzymatic network, as in Section2. Then, using the time-scale separation hypothesis, we reduce this ODE with the method described in Section 4. Finally, we estimate the parameters of the reduced system as it is suggested in Section 6.All the parameters in the toy network are supposed to satisfy the conditions established in (16) and Section 5. The periodic and continuous input considered is given by
(69)It=kcost⋅ω+π+1,where k is a parameter with the same order of magnitude as the slow reactions rates.


### 7.1. Reduction

We apply to the toy network our reduction scheme, as described in Section4.First, to simplify the notation, we define the following parameters:(70)aji≔eji0kjiKji∀i,j=1,…,n,d1≔a32a43+μ~+μ~a23+a43+μ~,d2≔a65+μ~a86+μ~a57+μ~a78+a98+μ~−a65a57a86a78.Then, we obtain the following reduced system for the toy network,(71)dX1dt=It−e210k21⋅X1X1+K21−μX1,dX4dt=a43a32d1⋅e210k21⋅X1X1+K21−e540k54⋅X4X4+K54−μX4,dX9dt=a98a65a86a57+μ~d2⋅e540k54⋅X4X4+K54−μX9,and the expressions for the metabolites in QSS,
(72)X2=εa23+a43+μ~d1⋅e210k21X1X1+K21,X3=εa32d1⋅e210k21X1X1+K21,X5=εa86+μ~a57+μ~a78+a98+μ~d2⋅e540k54X4X4+K54,X6=εa65a57+μ~a78+a98+μ~d2⋅e540k54X4X4+K54,X7=εa65a86a78d2⋅e540k54X4X4+K54,X8=εa65a86a57+μ~d2⋅e540k54X4X4+K54.


### 7.2. Calibration of the Reduced Toy Network

We follow the procedure in Section6. For simplicity, we suppose that the data are measured at the same time instants t1,…,tr (we assume that 48 measurement instants are available) for the slow and the fast parts of the system.The measurements are the variables (unit g/L) of the original system (9) (for the toy network in Figure 3) plus a white noise:
(73)Zitj≔Xitj+βtj,j=1,2,…,r,where β∼Nσi and σi=10−1·medianXi for every i=1,…,n.As in Section6, to estimate the parameters of (71), we use the reduced system (62), with m=4 and n=9. The cost function considered is ℱ, defined in (63), with m=4 and n=9.The functionfminsearch in Scilab (www.scilab.org) was used for minimizing ℱ. This function is based on the Nelder-Mead algorithm to compute the unconstrained minimum of a given function. For the simulations in Figure 4, the value θ^ obtained is in Table 1 and ℱθ^=0.097.Figure 4
Dynamics of the toy network represented in Figure3. The functions Xi represent the metabolite concentrations in units (g·L−1). The numerical solution of the original system (9) is depicted by the green line; the reduced system obtained by the method exposed in this work (71), by the blue dashed line; the supposed data with white noise (73), by green points; and the calibrated system (62)–(66) 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 ε lower than those of the metabolites in the slow part.Table 1
Parameter estimation for system (71), rewritten as (62). The estimation of these parameters only requires the slow dynamics of the toy network in Figure 3.
i
Theoretical valueθi
Initial guess
Estimated valueθ^i
Units
Error percent
1
0.110
0.010
0.072
g(L·min)−1
34.54
2
2.000
1.000
1.298
g·L−1
35.10
3
0.010
0.010
0.011
min− 1
10.00
4
0.110
0.010
0.073
g(L·min)−1
33.63
5
0.013
0.010
0.006
g(L·min)−1
53.85
6
2.000
1.000
2.143
g·L− 1
7.15
7
0.013
0.010
0.016
g(L·min)−1
23.08Table 2
Estimation of the parameters in (72), corresponding to the equalities in (66).
i
Theoretical valueαi
Estimated valueα^i
Error percent
2
0.00124
0.00092
25.81
3
0.00122
0.00089
27.05
5
0.00168
0.00185
10.12
6
0.00018
0.00019
5.56
7
0.00002
0.00002
0
8
0.00017
0.00019
11.76Here, for illustration purposes, we suppose that the metabolites in QSS are also measured; we calculate the parametersα^ to estimate their concentrations as explained in Section 6.2. Then, their concentrations are obtained according to (66).We computed the numerical solution of the systems describing the dynamics in the toy network of Figure3. The results are represented in Figures 4 and 5. As expected, the concentrations of the metabolites in QSS are one order of magnitude ε lower than those of the metabolites in the slow part.Figure 5
Zoom on the concentration of metabolites in QSS in Figure4.Note that the parametersθ2 and θ6 are affinity constants 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 (62)–(66) and recorded data, as represented in Figures 4 and 5.


## 7.1. Reduction

We apply to the toy network our reduction scheme, as described in Section4.First, to simplify the notation, we define the following parameters:(70)aji≔eji0kjiKji∀i,j=1,…,n,d1≔a32a43+μ~+μ~a23+a43+μ~,d2≔a65+μ~a86+μ~a57+μ~a78+a98+μ~−a65a57a86a78.Then, we obtain the following reduced system for the toy network,(71)dX1dt=It−e210k21⋅X1X1+K21−μX1,dX4dt=a43a32d1⋅e210k21⋅X1X1+K21−e540k54⋅X4X4+K54−μX4,dX9dt=a98a65a86a57+μ~d2⋅e540k54⋅X4X4+K54−μX9,and the expressions for the metabolites in QSS,
(72)X2=εa23+a43+μ~d1⋅e210k21X1X1+K21,X3=εa32d1⋅e210k21X1X1+K21,X5=εa86+μ~a57+μ~a78+a98+μ~d2⋅e540k54X4X4+K54,X6=εa65a57+μ~a78+a98+μ~d2⋅e540k54X4X4+K54,X7=εa65a86a78d2⋅e540k54X4X4+K54,X8=εa65a86a57+μ~d2⋅e540k54X4X4+K54.


## 7.2. Calibration of the Reduced Toy Network

We follow the procedure in Section6. For simplicity, we suppose that the data are measured at the same time instants t1,…,tr (we assume that 48 measurement instants are available) for the slow and the fast parts of the system.The measurements are the variables (unit g/L) of the original system (9) (for the toy network in Figure 3) plus a white noise:
(73)Zitj≔Xitj+βtj,j=1,2,…,r,where β∼Nσi and σi=10−1·medianXi for every i=1,…,n.As in Section6, to estimate the parameters of (71), we use the reduced system (62), with m=4 and n=9. The cost function considered is ℱ, defined in (63), with m=4 and n=9.The functionfminsearch in Scilab (www.scilab.org) was used for minimizing ℱ. This function is based on the Nelder-Mead algorithm to compute the unconstrained minimum of a given function. For the simulations in Figure 4, the value θ^ obtained is in Table 1 and ℱθ^=0.097.Figure 4
Dynamics of the toy network represented in Figure3. The functions Xi represent the metabolite concentrations in units (g·L−1). The numerical solution of the original system (9) is depicted by the green line; the reduced system obtained by the method exposed in this work (71), by the blue dashed line; the supposed data with white noise (73), by green points; and the calibrated system (62)–(66) 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 ε lower than those of the metabolites in the slow part.Table 1
Parameter estimation for system (71), rewritten as (62). The estimation of these parameters only requires the slow dynamics of the toy network in Figure 3.
i
Theoretical valueθi
Initial guess
Estimated valueθ^i
Units
Error percent
1
0.110
0.010
0.072
g(L·min)−1
34.54
2
2.000
1.000
1.298
g·L−1
35.10
3
0.010
0.010
0.011
min− 1
10.00
4
0.110
0.010
0.073
g(L·min)−1
33.63
5
0.013
0.010
0.006
g(L·min)−1
53.85
6
2.000
1.000
2.143
g·L− 1
7.15
7
0.013
0.010
0.016
g(L·min)−1
23.08Table 2
Estimation of the parameters in (72), corresponding to the equalities in (66).
i
Theoretical valueαi
Estimated valueα^i
Error percent
2
0.00124
0.00092
25.81
3
0.00122
0.00089
27.05
5
0.00168
0.00185
10.12
6
0.00018
0.00019
5.56
7
0.00002
0.00002
0
8
0.00017
0.00019
11.76Here, for illustration purposes, we suppose that the metabolites in QSS are also measured; we calculate the parametersα^ to estimate their concentrations as explained in Section 6.2. Then, their concentrations are obtained according to (66).We computed the numerical solution of the systems describing the dynamics in the toy network of Figure3. The results are represented in Figures 4 and 5. As expected, the concentrations of the metabolites in QSS are one order of magnitude ε lower than those of the metabolites in the slow part.Figure 5
Zoom on the concentration of metabolites in QSS in Figure4.Note that the parametersθ2 and θ6 are affinity constants 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 (62)–(66) and recorded data, as represented in Figures 4 and 5.


## 8. Discussion

### 8.1. Time-Scale 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 preserves 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 Table3 and Table 4). The simulation results illustrate Theorem 1 (see Figure 4), and the reduced system accurately represents all different time scales (see Figures 4 and 5).Table 3
Parameters considered for the simulations in Figure4. The symbol γ∈−1,1 denotes a rate in an enzymatic reaction (see the Michaelis-Menten equation (A.1)). The initial conditions for all the enzymes are the same, as well as the initial conditions of all the metabolites are identical, that is, j,i∈1,…,n in this table.
Parameter
Value
Units
ε
0.001
—
k
0.01
g(L·min)−1
μ
0.01
min−1
ω
0.004
—
kγ21,k21
1.10
min−1
kγ54,k54
0.13
min−1
kγ32,k32
1.90
min−1
kγ23,k23
0.12
min−1
kγ43,k43
1.80
min−1
kγ65,k65
0.17
min−1
kγ57,k57
1.40
min−1
kγ86,k86
1.60
min−1
kγ78,k78
0.15
min−1
kγ98,k98
1.50
min−1
eji0
0.10
g·L−1
xi0
0.001
g·L−1Table 4
Slow and fast reaction rates considered for the toy network in Figure3 and the simulations in Figure 4. Fast reaction rates are characterized by the factor 1/ε. All reactions rates are in units of g(L·min)−1.
Slow rates
Value
e210k21
1.1 × 10−1
e540k54
1.3 × 10−2
Fast rates
Value
e320k32/ε
1.3 × 102
e230k23/ε
1.2 × 101
e430k43/ε
1.8 × 102
e650k65/ε
1.7 × 101
e570k57/ε
1.4 × 102
e860k86/ε
1.6 × 102
e780k78/ε
1.5 × 101
e980k98/ε
1.5 × 102Finally, 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ε [6, 19, 20] or even different epsilons [21]. But computations with this method highly complicate the reduction.


### 8.2. Comparison with Experimental Data

To the best knowledge of the authors, there are to date no examples 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 of macromolecules (typically proteins, carbohydrates, lipids, chlorophyll, etc.) 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 the parameters of a metabolic model to the set of available experimental data. We examined the reaction rates which ranged from 102 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 a magnitude 102 (mM) times higher than that of the intermediate metabolites (GAP, PEP, and G6P). Moreover, GAP, PEP, and G6P are consumed by reactions with a rate order 101 or 102 (h−1·mM·B−1), while carbohydrates are consumed at a rate of order 100 (h−1·mM·B−1). Additionally, carbohydrates are produced by a single reaction with a rate of order 101 (h−1·mM·B−1), as well as GAP, and G6P, and PEP by reactions of order 102 (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.Table 5
Experimental measures (M) and estimated (E) values obtained from [4], for an autotrophic microalgae 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.
Compound
Value
Mean value
Carbohydrates
M
8.436 × 10−1 mM
G6P
E
5.208 × 10−3 mM
PEP
E
4.167 × 10−3 mM
GAP
E
1.389 × 10−3 mMTable 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.
Compound
Production rate
Consumption rate
Subnetwork
Carbohydrates
7.00 × 101
6.50 × 100
Carbohydrate synthesis
G6P
2.24 × 102
1.03 × 101
Upper glycolysis
6.50 × 100
7.00 × 101
Lipid synthesis
PEP
4.37 × 102
5.00 × 100
Lower glycolysis
9.97 × 100
1.04 × 102
Lipid synthesis
GAP
2.06 × 101
4.47 × 102
Upper glycolysis
5.00 × 100
4.37 × 102
Upper glycolysis
6.00 × 10−1
1.88 × 101
Lipid synthesisNevertheless, we emphasize that the reduction method proposed in this paper can be used even if only some metabolites with large concentrations have been measured. Indeed, such data will support the calibration of the reduced model, that is, describing the dynamics of the slow metabolites (see Section6.1).


### 8.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 Theorem1 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.1. Time-Scale 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 preserves 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 Table3 and Table 4). The simulation results illustrate Theorem 1 (see Figure 4), and the reduced system accurately represents all different time scales (see Figures 4 and 5).Table 3
Parameters considered for the simulations in Figure4. The symbol γ∈−1,1 denotes a rate in an enzymatic reaction (see the Michaelis-Menten equation (A.1)). The initial conditions for all the enzymes are the same, as well as the initial conditions of all the metabolites are identical, that is, j,i∈1,…,n in this table.
Parameter
Value
Units
ε
0.001
—
k
0.01
g(L·min)−1
μ
0.01
min−1
ω
0.004
—
kγ21,k21
1.10
min−1
kγ54,k54
0.13
min−1
kγ32,k32
1.90
min−1
kγ23,k23
0.12
min−1
kγ43,k43
1.80
min−1
kγ65,k65
0.17
min−1
kγ57,k57
1.40
min−1
kγ86,k86
1.60
min−1
kγ78,k78
0.15
min−1
kγ98,k98
1.50
min−1
eji0
0.10
g·L−1
xi0
0.001
g·L−1Table 4
Slow and fast reaction rates considered for the toy network in Figure3 and the simulations in Figure 4. Fast reaction rates are characterized by the factor 1/ε. All reactions rates are in units of g(L·min)−1.
Slow rates
Value
e210k21
1.1 × 10−1
e540k54
1.3 × 10−2
Fast rates
Value
e320k32/ε
1.3 × 102
e230k23/ε
1.2 × 101
e430k43/ε
1.8 × 102
e650k65/ε
1.7 × 101
e570k57/ε
1.4 × 102
e860k86/ε
1.6 × 102
e780k78/ε
1.5 × 101
e980k98/ε
1.5 × 102Finally, 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ε [6, 19, 20] or even different epsilons [21]. But computations with this method highly complicate the reduction.


## 8.2. Comparison with Experimental Data

To the best knowledge of the authors, there are to date no examples 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 of macromolecules (typically proteins, carbohydrates, lipids, chlorophyll, etc.) 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 the parameters of a metabolic model to the set of available experimental data. We examined the reaction rates which ranged from 102 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 a magnitude 102 (mM) times higher than that of the intermediate metabolites (GAP, PEP, and G6P). Moreover, GAP, PEP, and G6P are consumed by reactions with a rate order 101 or 102 (h−1·mM·B−1), while carbohydrates are consumed at a rate of order 100 (h−1·mM·B−1). Additionally, carbohydrates are produced by a single reaction with a rate of order 101 (h−1·mM·B−1), as well as GAP, and G6P, and PEP by reactions of order 102 (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.Table 5
Experimental measures (M) and estimated (E) values obtained from [4], for an autotrophic microalgae 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.
Compound
Value
Mean value
Carbohydrates
M
8.436 × 10−1 mM
G6P
E
5.208 × 10−3 mM
PEP
E
4.167 × 10−3 mM
GAP
E
1.389 × 10−3 mMTable 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.
Compound
Production rate
Consumption rate
Subnetwork
Carbohydrates
7.00 × 101
6.50 × 100
Carbohydrate synthesis
G6P
2.24 × 102
1.03 × 101
Upper glycolysis
6.50 × 100
7.00 × 101
Lipid synthesis
PEP
4.37 × 102
5.00 × 100
Lower glycolysis
9.97 × 100
1.04 × 102
Lipid synthesis
GAP
2.06 × 101
4.47 × 102
Upper glycolysis
5.00 × 100
4.37 × 102
Upper glycolysis
6.00 × 10−1
1.88 × 101
Lipid synthesisNevertheless, we emphasize that the reduction method proposed in this paper can be used even if only some metabolites with large concentrations have been measured. Indeed, such data will support the calibration of the reduced model, that is, describing the dynamics of the slow metabolites (see Section6.1).


## 8.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 Theorem1 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.


## 9. 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 the 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 of 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 magnitude 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 example, 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.



---

*Source: 2342650-2018-08-12.xml*