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

HAL
open science

# Flux balance analysis-based ranking for model order reduction of biochemical networks

Carlos Eduardo Robles-Rodriguez, Erik Steur

► To cite this version:

Carlos Eduardo Robles-Rodriguez, Erik Steur. Flux balance analysis-based ranking for model order reduction of biochemical networks. IFAC-PapersOnLine, 2021, 54 (15), pp.556-561. 10.1016/j.ifacol.2021.10.315. hal-03450478

HAL Id: hal-03450478

https://hal.inrae.fr/hal-03450478

Submitted on 26 Nov 2021

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Flux balance analysis-based ranking for model order reduction of biochemical networks

Carlos Eduardo Robles-Rodriguez* Erik Steur**

* TBI, Université de Toulouse, CNRS, INRAE, INSA, Toulouse, France. (e-mail: roblesro@insa-toulouse.fr).

** Institute for Complex Molecular Systems and the Department of Mechanical Engineering, Eindhoven University of Technology, Eindhoven, The Netherlands. (e-mail: e.steur@tue.nl)

**Abstract:** This paper presents a parameter free approach to identify the importance of reactions (and involved species) in biochemical reaction networks for the purpose of model order reduction. The new methodology is based on the structure of the network assuming that at steady state the flux distribution is preserved. Our method employs Flux Balance Analysis (FBA) to calculate these flux distributions. The ranking is based on the comparison of the FBA results of the original and reduced networks. The main purpose of this identification step is to guide the selection of species that could be deleted by model order reduction methods. Our ranking method is illustrated in a model for Glycolysis, where model reduction is performed via the Kron reduction method considering the elimination of the less sensitive species proposed by the ranking method. The reduced model shows that the global dynamics could be accurately reproduced with a smaller number of species.

Copyright © 2021 The Authors. This is an open access article under the CC BY-NC-ND license (https://creativecommons.org/licenses/by-nc-nd/4.0/)

**Keywords:** Biochemical Reactions, Flux Balance Analysis, Sensitivity analysis, Model reduction

## 1. INTRODUCTION

Dynamic models of biochemical reaction networks typically consist of large sets of nonlinear ordinary differential equations and many parameters, which complicates their analysis. Model reduction techniques can alleviate this complexity by decreasing the number of species, reactions and parameters involved (Radulescu et al., 2012). Furthermore, these methods help to gain a better understanding of the role of the constituting components of the network and to analyze better their contribution to the response of the system. There exist a wide variety of reduction methods for biochemical systems (Snowden et al., 2017), where no single method is superior to other since most of them depend on the complexity of the system. The main challenge in applying these methods is the selection of the species and/or reactions to be removed in the reduction. Such a selection, however, often requires the kinetic rate constants and/or species concentrations to be known, while in practice one often only knows the ranges in which the parameters take values.

This paper introduces a ranking methodology that is solely based on the structure of the network, and is independent of the kinetics. In order to rank the reactions and species, we assume that the system reaches a steady state with a specific flux distribution. Under this assumption, we can identify the contribution of the reactions towards the desired steady state and assume that the low-contribution-reactions and the species participating in them are consumed by fast reactions. This is similar to the quasi-steady state assumption (QSSA) used in different approaches of

model order reduction (e.g. time-scale separation, layered decomposition, cf. (Snowden et al., 2017)). The main difference, however, is that our identification is performed under different steady state conditions to mimic the behavior of the network. The ranking is computed as a sensitivity index based on the comparison of the flux distribution of the original and reduced networks, which is computed via Flux Balance Analysis (FBA) (Varma and Palsson, 1994). FBA is a well known tool in metabolic engineering to reconstruct networks and to determine their thermodynamically feasible, stationary flux distribution based on maximization of objective functions (Kauffman et al., 2003). Therefore, the proposed method present the following advantages: 1. parameter-free, 2. preservation of steady-states, 3. independent of species concentrations.

The less sensitive species from the ranking methodology are considered as candidates to be removed in the model order reduction procedure. In this paper, we use the Kron reduction method presented by Rao et al. (2014) to perform the reduction. This method simplifies a given reaction network by deleting complexes from the complex graph. To illustrate the idea, let us consider the following set of reactions:

$$A + B \rightleftharpoons AB \rightarrow P$$

in which $A + B, AB, P$ are the complexes, formed by the species $A, B, AB, P$. Kron reduction deletes complexes (nodes) from this complex graph and replaces the reactions (edges) by new ones. By removing $AB$, we get $A + B \rightarrow P$. These new reactions (the dashed arrow in the example) in the reduced graph are computed by taking the Schur
complement of the (weighted) Laplacian matrix that encodes the interactions in the original complex graph. The main attractiveness of this method is the simplicity of its implementation and the possibility to deal with various reversible and irreversible enzyme kinetics, e.g. mass-action and Michaelis-Menten kinetics.

In (Rao et al., 2014), the Kron model reduction is automated using an iterative scheme. In this scheme, the effect on the dynamics for removing one complex from the model via Kron reduction is measured using an integral error. The best performing reduced model (i.e., the one resulting in the lowest error) is then used as the (new) model in the next iteration step. The procedure is repeated until a predefined criterion "the maximal tolerable error" is reached. Our proposed FBA-based ranking method provides an alternative for this iterative scheme. Our approach characterizes algebraically (and independent of the kinetic parameters) the importance of species and reactions in the generation of the output(s) of biochemical reaction networks. Thus the reduced model can be obtained without numerical integration of the dynamics and evaluation of an error. The proposed FBA based ranking method is applied, together with the Kron reduction method, to a case study concerning yeast Glycolysis.

The remainder of the paper is as follows. Section 2 presents the FBA ranking methodology. Section 3 described the Kron reduction method. The case study and results are reported in section 4. Finally section 5 concludes this work.

## 2. FBA-BASED RANKING METHOD

The dynamics of biochemical systems with $m$ species are usually described (Snowden et al., 2017; Radulescu et al., 2012; Klamt et al., 2018) by equations of the form

$$ \dot{x} = Sv(x), \quad (1) $$

where $x \in \mathbb{R}_+^m$ is the vector of non-negative species concentrations, $v: \mathbb{R}_+^m \to \mathbb{R}^r$ is the reaction rate flux vector (kinetics) and $S \in \mathbb{R}^{m \times r}$ the stoichiometric matrix. Typically, the identification of the importance of the species and reactions to be deleted requires detailed knowledge about the kinetics of the system $v(x)$. Our ranking methodology is built on steady-state analysis, hence, it only depends on the stoichiometry of the network $S$. This is based on Flux Balance Analysis (FBA) (Varma and Palsson, 1994), which helps to determine how much a reaction contributes to reaching a specific objective. We make use of this characteristic to rank the reactions assuming that the flux distribution towards the objective is preserved. In this sense, the low contribution reactions can be assumed to be fast reactions and the species involved are assumed to be at QSS.

The application of the FBA-based ranking method follow several steps. Step 1: All the reactions of the systems must have a reactant and a product (input/output representation) and be expressed as irreversible reactions. Step 2: Several objectives (the maximization of one or several products) are introduced and evaluated via FBA to determine the flux distribution of each objective. Step 3: The flux distribution resulting from FBA is normalized as yields relating the flux of the output over the flux of an input, therewith providing a measure of input-output

sensitivity of single reactions. Step 4: Reduced networks are built by adding up one reaction to another reaction with which it is connected to obtain reduced networks with one reaction less. As the number of reactions is large, we only consider the reactions which, combined to another, eliminate at least one species. Step 5: The fluxes for the reduced network are computed and compared with the complete one as yields. The errors on the yields provide us a sensitivity measurement of the importance of each reaction. Step 6: After all the reduced networks have been evaluated, a new normalization is performed to obtain the ranking of the reactions. By using the stoichiometric matrix, these results can be translated into the ranking of the deleted species. In the next subsections, we describe in more detail how the ranking is performed.

### 2.1 Flux Balance Analysis

Flux Balance Analysis (FBA) is used to analyze the flux distribution through a metabolic network towards a specific objective in steady-state (Varma and Palsson, 1994). The main assumption of FBA is that the networks evolve to achieve an optimal metabolic objective. The heart of FBA relies on the constraints imposed by the stoichiometry to the flow of species through the network (Kauffman et al., 2003). This approach identifies particular flux distributions through a metabolic network by optimizing a linear objective function (e.g. maximize growth) by means of linear programming (Orth et al., 2010). The FBA problem is formulated as

$$ \begin{array}{ll} \max_{v} & J = w^T v \\ \text{s.t.} & Sv = 0 \\ & v^L \le v \le v^u \end{array} \quad (2) $$

where the solution $v$ is at steady state (as $Sv = 0$) and constrained between a lower and an upper bound, $v^L$ and $v^u$, respectively. The vector $w$ specifies the weights indicating how much each reaction contributes to the objective function. In practice, when only one reaction is desired for optimization, $w$ is a vector of zeros with a one in the position of the reaction of interest.

The main advantage of the FBA is that it does not require any knowledge of the species concentrations, and more importantly, the kinetics of the system (Orth et al., 2010). Nonetheless, defining a single objective function may not be sufficient, as metabolism seems to operate at an optimal surface in multi-dimensional space with competing objectives.

### 2.2 Selection of reactions and species reduction candidates

Let us consider the representation of the system as a species-reaction graph (Feinberg, 2019) as,

$$ s_0 \to R_k \to s_1 \to R_\ell \to s_2 \to \dots \quad (3) $$

in which $R_i$ are reactions and $s_j$ are species. As species can take part in several reactions, let us define interacting reactions as the reactions that share at least one species. For instance, in (3) $R_k$ and $R_\ell$ are interacting reactions. However, the interacting reactions are not always linearly connected, as sometimes a species participates in several
reactions (e.g., Isocitrate in Krebs cycle). Then, the reduction would be to delete the shared species ($s_1$) and produce the new path: $s_0 \rightarrow R_k + R_\ell \rightarrow s_2 \rightarrow \dots$. Here the + indicate that the new reaction depends on both $R_k$ and $R_\ell$. This procedure is performed with one interacting reaction at a time considering the total number of interacting reactions $Q_k$ for a given reaction $R_i$. This sum-deletion approach implies the construction of a new network that contains $r-1$ reactions.

In this process, we must ensure that neither new species nor pathways are created that can alter the description of the network. Two important properties are defined and included in the calculation of sensitivity in order to prevent these concerns. The properties are:

* **Feasibility ($\lambda$)**. The reduced system does not contain new pathways and the species in the reactions that are not inputs and outputs are either consumed or produced. For example, let us consider the reaction network of Fig. 1a. Species A and C are inputs and outputs, respectively, and B and D are internal species. Summing up $R_1$ and $R_2$, the reduced reaction system is as in Fig. 1b, where the species B is never produced, but it is consumed to form D. In this case, the reduction is infeasible. On the contrary, adding up $R_3$ and $R_4$ eliminates D and the network is still connected showing a feasible reduction (Fig. 1c).

* **Reduction capability ($\omega$)**. Only the sum of reactions that delete at least one species are considered. However, if all species are eliminated (sum of reversible reaction), that combination is not accounted.

<figure><img src="image_3.png" /><figcaption>Fig. 1. Example of Network (a) with infeasible reduction (b) and feasible reduction (c).</figcaption></figure>

These two properties are included in the sensitivity analysis as Boolean variables regarding reaction $k$ and its interacting reaction $\ell$ as,

$$ \lambda_{k,\ell} = \begin{cases} 1, & \text{if feasible} \\ 0, & \text{otherwise} \end{cases}, \quad \omega_{k,\ell} = \begin{cases} 1, & \text{if species deleted} \\ 0, & \text{otherwise} \end{cases} $$

## 2.3 FBA-based sensitivity analysis

In order to assess the reduced networks, the method takes advantage of FBA to measure the contribution of reactions to attain a specific objective. A sensitivity index of the reactions is calculated based on the contribution of each reaction towards the specific objective by quantifying the impact of its deletion on the objective. Maximizing only one objective, however, is usually not sufficient to represent the whole system, specially if the interest is on the dynamic behavior (Hjersted and Henson, 2009).

To this end, it is necessary to set different objectives, here denoted as $J_j$ with $j$ being the number of objective functions. These objectives should represent the largest possible amount of objectives that the reaction network can attain. As the system is in an input/output representation, the maximization of different outputs (products) and/or combinations of them are good targets.

**Set the target objectives** The first step is to calculate the flux distribution of the complete reaction network described by the stoichiometric matrix $S$ with $m$ species and $r$ reactions. This imply solving problem (2) for each objective $J_j$.

Let $F^j$ be the optimal flux distribution for the objective $J_j$. It is important to note that the values of $F_j$ are fluxes and not concentrations. Therefore, following literature, it is preferred to expressed those results as yields. In order to express the results as yields and to respect an input/output representation, the flux of the objective $F_O^j$ is normalized with respect to the flux of an input to the network $F_I^j$. For the sake of simplicity, let us set the target objective as the maximum yield $Y^j$ which is calculated as follows:

$$ Y^j = F_O^j/F_I^j. \quad (5) $$

This maximum yield (target) is used to make a comparison between the flux distribution (yield distribution) of the complete system and the reduced system. The idea is to rank the importance of each reaction $k$ to attain the desired objective $J_j$ by the feasible elimination of one reaction at a time.

**Yields of networks with $r-1$ reactions** FBA is performed for each objective $J_j$ considering each reaction $k$ and its addition to each $\ell$ interacting reaction. This results in computing FBA $k \times \ell$ times for an objective $j$. The flux distribution is expressed as $\bar{F}_{k,\ell}^j$. The results are normalized taking the flux of the objective reaction $\bar{F}_{O,k,\ell}^j$ over the flux of and input reaction $\bar{F}_{I,k,\ell}^j$ to get the yield of the reduced system as,

$$ \bar{Y}_{k,p}^j = \bar{F}_{O,k,\ell}^j/\bar{F}_{I,k,\ell}^j \quad (6) $$

These yields $\bar{Y}_{k,\ell}^j$ are then compared with the objective yield $Y^j$ to build a sensitivity function.

## 2.4 Ranking with sensitivity

The sensitivity measurement of each reaction $k$ is expressed as $\sigma_{rxn,k}^j$ for each $j$ objective. The ranking is computed by means of the Root Mean Squared Error as,

$$ \sigma_{rxn,k}^j = \sqrt{\sum_{\ell=1}^{Q_k} s_{k,\ell}^{j,2} / \max \left( 1, \sum_{\ell=1}^{Q_k} \lambda_{k,\ell} \omega_{k,\ell} \right)} \quad (7) $$

where

$$ s_{k,\ell}^j = \frac{Y^j - \bar{Y}_{k,\ell}^j}{Y^j} \lambda_{k,\ell} \omega_{k,\ell} \quad (8) $$

This measurement takes into account the difference of the reduced systems with respect to the complete one considering a feasible reduction. Furthermore, this measurement
can be extrapolated to the species via the stoichiometric matrix S as,

$$ \sigma_{sp}^{j} = S\sigma_{rxn}^{j} \quad (9) $$

From the results of sensitivity, the ranking method identifies which reactions and species are more sensitive to the different objective functions. If the amount of species to be deleted is high, a cut-off value for $\sigma_{sp}^j$ could be set according to the magnitude of the results in order to classify the species to be deleted.

### 3. KRON REDUCTION METHOD

In this section we show how the ranking method presented in Section 2 can be combined with Kron reduction to generate reduced order models of biochemical reaction networks. The Kron reduction method of (Rao et al., 2014) follows the modeling approach of (van der Schaft et al., 2015). For this approach, let c be complexes in the biochemical reaction network associated to model (1). These complexes relate to the m species by the complexes stoichiometric matrix $Z \in \mathbb{R}^{m \times c}$, and the link between the complexes and the reactions is defined by the incidence matrix $B \in \mathbb{R}^{c \times r}$ (van der Schaft et al., 2015) (One can show that $S = ZB$).

With these conventions, model (1) can be written in the form

$$ \dot{x} = ZBF(x)KExp(Z^T Ln(x)) \quad (10) $$

where the elements of the matrix $K \in \mathbb{R}^{r \times c}$ are the mass-action rate constants of the k reaction of the network, and $F: \mathbb{R}_+^m \to \mathbb{R}^r$ is a matrix of rational functions comprising the denominators of the kinetics of reaction k, which allows the utilization of different kinetics, such as Michaelis-Menten and Hill (Rao et al., 2014). The functions $Exp(\cdot)$ and $Ln(\cdot)$ are the component-wise natural exponential and logarithmic functions, respectively. Define $L = L(x) := -BF(x)K$. It can be easily verified that L has non-negative diagonal elements and non-positive off-diagonal elements. Moreover, the sum of the elements in each column of L is equal to zero. Hence, hereafter $L \in \mathbb{R}^{c \times c}$ is considered as the weighted Laplacian matrix. Substituting the definition of L in (10), the system is expressed as,

$$ \dot{x} = -Z L \exp(Z^T Ln(x)) \quad (11) $$

The Kron reduction method performs a reduction by using a partition of the system into the combination of the remaining complexes R and the deleted ones D. Therefore, the system (11) is rewritten as

$$ \begin{bmatrix} \dot{x}_R \\ \dot{x}_D \end{bmatrix} = - \begin{bmatrix} Z_R \\ Z_D \end{bmatrix} \begin{bmatrix} L_{R,R} & L_{R,D} \\ L_{D,R} & L_{D,D} \end{bmatrix} \begin{bmatrix} \exp(Z_R^T Ln(x)) \\ \exp(Z_D^T Ln(x)) \end{bmatrix}. \quad (12) $$

In order to perform the reduction, the deleted species $x_D$ are assumed to be at their steady state (i.e. $\dot{x}_D = 0$). The reduced system is then obtained by taking the Schur complement of the Laplacian matrix and takes the form:

$$ \dot{x}_R = -Z_R L_R \exp(Z_R^T Ln(x)), \quad (13) $$

where $L_R$ is defined as

$$ L_R = L_{R,R} - L_{R,D} L_{D,D}^{-1} L_{D,R}. \quad (14) $$

As mentioned before, in (Rao et al., 2014), the reduction is obtained by measuring the effect of the deletion of one

complex at a time using the integral error in an iterative way. We use the ranking method to define the partition of the vector of species concentrations x into $x_R$ and $x_D$, and then perform the reduction.

To ensure that the steady-state of the remaining species is preserved in the reduction process, the use of the method is restricted to complex balanced networks. Hence, its use with intermediate complexes having more than two species requires rewriting the system with either conservation laws (Gasparyan et al., 2020) or clustering the species in a complex (Rao et al., 2014).

### 4. CASE STUDY: YEAST GLYCOLYSIS

#### 4.1 Description of Yeast Glycolysis Model

The implementation of the ranking method is illustrated in a detailed kinetic model for yeast Glycolysis from literature (van Eunen et al., 2012). A schematic representation of the metabolic network is displayed in Fig. 2. This model considers that Glucose (Glci) is taken up by the cells and can be metabolized via glycolysis to produce five products: Ethanol (EtOH), Acetate, Succinate, Glycerol and Trehalose. Ethanol and Acetate are produced from acetaldehyde (AcAld), Succinate from pyruvate (PYR), and Glycerol from TRIO which is a pool including dihydroxyacetone phosphate and glyceraldehyde 3-phosphate. Trehalose is produced from glucose 6-phosphate and it is, at the same time, consumed to produce intracellular glucose. The model comprises 12 ordinary differential equations to represent the dynamics of the internal metabolites. The detailed mathematical model can be found in van Eunen et al. (2012).

<figure><img src="image_4.png" /><figcaption>Fig. 2. Metabolic network of yeast Glycolysis.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 3. Sensitivity of the reactions involved in yeast glycolysis</figcaption></figure>

## 4.2 Ranking results

The first step in the sensitivity analysis is to write the system with irreversible reactions as it is displayed in Fig. 2. Then, the objective functions that will be used for FBA need to be defined. Based on an input/output representation, we consider that the metabolism will prioritize the production of the products. Hence, the objective functions were defined as the maximization of the production of Ethanol (EtOH), Glycerol (GLY), Succinate (SUC), Acetate (ACT) and Trehalose (TREA). Additionally, we assumed that the cofactors NAD, NADH, ADP and ATP cannot be eliminated in order to avoid energetic imbalances in the network. The computation of the sensitivity of reactions and species was performed with equations (7) and (9). The ranking of the reactions is displayed in Fig. 3 at a scale ranging from 0 *not important*, to 1 *very important/necessary*. As it can be observed, there is a large difference in the reactions that are important (values close to 1) and the ones least important (values close to 0). Reactions R5, R6, R12, R13, R14, R15, R16, R17, R19 and R21 have values close to 0, which means that they do not contribute to the maximization of the defined objectives. Some of these reactions are the backwards reactions of the production of a species specified as objective (i.e., R21). Interestingly reactions R10 and R11, corresponding to the transformation of TRIO into glucose 1,6-biphosphate (BPG) are not important for maximizing Succinate.

The sensitivity of the species is presented in Fig. 4. It is observed that the most important variables are the cofactors NAD and NADH followed by TRIO which is the branching point to produce Glycerol or the rest of the products. A cut-off value was set to 0.4 to distinguish the less important species (in this case also equivalent to complexes). The species that were below that threshold corresponded to P3G, P2G, F6P, PEP and G6P, which were the strongest candidates for the reduction. It is worth noting that although the initial glucose (Glco) and the products have low sensitivity, they cannot be considered for elimination to keep the input/output representation.

<figure><img src="image_6.png" /><figcaption>Fig. 4. Sensitivity of the species involved in yeast Glycolysis</figcaption></figure>

<table><caption>Table 1. Initial Conditions (IC) and Steady state (SS) values for Glycolysis model. Units in mM.</caption><thead><tr><th>Species</th><th>IC</th><th>SS</th><th>Species</th><th>IC</th><th>SS</th></tr></thead><tbody><tr><td>Glci</td><td>0.057</td><td>0.623</td><td>P2G</td><td>0.161</td><td>0.041</td></tr><tr><td>G6P</td><td>0.121</td><td>1.384</td><td>PEP</td><td>0.982</td><td>0.054</td></tr><tr><td>F6P</td><td>0.026</td><td>0.304</td><td>PYR</td><td>0.503</td><td>1.571</td></tr><tr><td>F16P</td><td>0.093</td><td>11.95</td><td>AcAld</td><td>13.0</td><td>12.99</td></tr><tr><td>TRIO</td><td>0.336</td><td>3.816</td><td>NADH</td><td>0.0003</td><td>0.0015</td></tr><tr><td>BPG</td><td>0.0009</td><td>0.003</td><td>NAD</td><td>1.589</td><td>1.588</td></tr><tr><td>P3G</td><td>0.946</td><td>0.611</td><td></td><td></td><td></td></tr></tbody></table>

## 4.3 Model order reduction results

The ODE model and the parameters described in (van Eunen et al., 2012) were used for the complete model which was written in the form of (11). The Kron reduction method has been implemented to build reduced models with combinations of candidate species in order to test the impact in the dynamics when those species are eliminated. From the ranking method, the candidate species to be eliminated species were P3G, P2G, F6P, PEP and G6P. Following (van Eunen et al., 2012), the reduction was applied to a glucose up-shift in which the concentrations of Glucose are 0.2 mM and 0.5 mM for $t < 0$ and $t \ge 0$, respectively. The concentrations of ATP corresponded to 5 mM and 2.5 mM for the same intervals. Furthermore, it is assumed that the model is at steady state at $t < 0$.

The initial conditions and the steady state values of the species are reported in Table 1. In order to assess the reduction, two criteria have been considered: the Integral error (Rao et al., 2014) and the Root Mean Squared Relative Errors ($RMSE_r$) (Apri et al., 2012), which are calculated as in equations (15) and (16), respectively.

$$I = \sum_{h \in M_1} \frac{1}{\tau n(M_1)} \int_{0}^{\tau} \left| 1 - \frac{x_{h,R}(t)}{x_h(t)} \right| dt \quad (15)$$

$$RMSE_r = \sqrt{\frac{1}{N} \sum_{h \in M_1} \left( \frac{x_h(t) - x_{h,R}(t)}{x_h(t)} \right)^2} \quad (16)$$
<table><caption>Table 2. Results of the reduction with different combinations.</caption><thead><tr><td>Number</td><td>Deleted species</td><td>RMSE<sub>r</sub></td><td>I</td></tr></thead><tbody><tr><td>1</td><td>P3G, P2G, F6P</td><td>0.186</td><td>0.027</td></tr><tr><td>2</td><td>P3G, P2G, F6P, G6P</td><td>0.248</td><td>0.054</td></tr><tr><td>3</td><td>P3G, P2G, F6P, PEP</td><td>0.452</td><td>0.097</td></tr><tr><td>4</td><td>P3G, P2G, F6P, G6P, PEP</td><td>0.482</td><td>0.127</td></tr></tbody></table>

$M_1$ correspond to the measured species or the species that are significant for the reduction under the desired interval $[0,τ]$. The terms $x_{h,R}(t)$ and $x_h(t)$ are the concentrations of the $M_1$ species at time $t$ in the reduced and complete model, respectively.

Four different combinations have been tested based on the ranking of the species, where the least sensitive species (P3G, P2G, F6P) have been considered to have a negligible impact in the reduction. As the model only represents the dynamics of the internal metabolites, the set of significant species $M_1$ has been set to the metabolites closer to the inputs/outputs and the most sensitive species: Glci, F16P, TRIO, PYR, AcAld and NADH. The combinations and the results of the reduction performance are presented in Table 2. All the reductions present small values for $RMSE_r$ and $I$. Therefore, the reduced model 4 which deletes the largest amount of species or complexes is selected as the best reduction. The results of reduced model 4 compared against the complete model are displayed in Fig. 5. It is worth noting that the deleted species take a constant value corresponding to their steady state value. A good agreement is observed between the transient behavior of most of the variables. The main difference is for AcAld where the deletion of PEP induces a change in the dynamics, but the steady state is well preserved. The same model for Glycolysis has been reduced by Rao et al. (2014) using the iterative Kron reduction method. Their final reduced model is the same as the one obtained in this paper. In their implementation they have chosen other complexes like F16P, which after their iterative procedure they have found that it was an important complex and they could not deleted. Our ranking method, however, has identified beforehand that F16P cannot be deleted. It is important to note that when PEP is included, both error measures ($I$ and $RMSE_r$) increase. This highlights that although the steady state analysis provided by the ranking method simplifies the selection of the candidates to be deleted, the test of the dynamics is required.

## 5. CONCLUSIONS

We presented a method for ranking the reactions of a network and selecting the candidate species to be deleted by model reduction methods. Our selection method requires only knowing the structure of the network (parameter free). The selection is made based on the (relative) sensitivity of the reactions in flux distribution at pre-defined objectives. The implementation of the method is shown to be efficient in the identification of the candidate species and complexes to be deleted. Implemented together with the Kron reduction method a reduction of 5 complexes out of 12 was possible in the presented case study.

## REFERENCES

* Apri, M., de Gee, M., and Molenaar, J. (2012). Complexity reduction preserving dynamical behavior of biochemical networks. *Journal of*

<figure><img src="image_7.png" /><figcaption>Fig. 5. Comparison of the concentration profiles of yeast glycolysis metabolites. (—) Complete model, (—) Reduced model. All the concentrations are expressed in mmol/L.</figcaption></figure>

Theoretical Biology, 304, 16–26. doi:10.1016/j.jtbi.2012.03.019.

* Feinberg, M. (2019). Species-Reaction Graph Foundations. In *Applied Mathematical Sciences (Switzerland)*, volume 202, 419–440. Springer. doi:10.1007/978-3-030-03858-8 19.

* Gasparyan, M., Van Messem, A., and Rao, S. (2020). An automated model reduction method for biochemical reaction networks. *Symmetry*, 12(8), 1–24. doi:10.3390/sym12081321.

* Hjersted, J.L. and Henson, M.A. (2009). Steady-state and dynamic flux balance analysis of ethanol production by Saccharomyces cerevisiae. *IET systems biology*, 3(3), 167–79. doi:10.1049/iet-syb.2008.0103.

* Kauffman, K.J., Prakash, P., and Edwards, J.S. (2003). Advances in flux balance analysis. *Current Opinion in Biotechnology*, 14(5), 491–496. doi:http://dx.doi.org/10.1016/j.copbio.2003.08.001.

* Klamt, S., Müller, S., Regensburger, G., and Zanghellini, J. (2018). A mathematical framework for yield (vs. rate) optimization in constraint-based modeling and applications in metabolic engineering. *Metabolic Engineering*, 47, 153–169. doi: 10.1016/j.YMBEN.2018.02.001.

* Orth, J.D., Thiele, I., and Palsson, B.Ø. (2010). What is flux balance analysis? *Nature biotechnology*, 28(3), 245–8. doi:10.1038/nbt.1614.

* Radulescu, O., Gorban, A.N., Zinovyev, A., and Noel, V. (2012). Reduction of dynamical biochemical reactions networks in computational biology. *Frontiers in Genetics*, 3(JUL), 131. doi: 10.3389/fgene.2012.00131.

* Rao, S., van der Schaft, A., van Eunen, K., Bakker, B.M., and Jayawardhana, B. (2014). A model reduction method for biochemical reaction networks. *BMC Systems Biology*, 8(1), 1–17. doi:10.1186/1752-0509-8-52.

* Snowden, T.J., van der Graaf, P.H., and Tindall, M.J. (2017). Methods of Model Reduction for Large-Scale Biological Systems: A Survey of Current Methods and Trends. *Bulletin of Mathematical Biology*, 79(7), 1449–1486. doi:10.1007/s11538-017-0277-2.

* van der Schaft, A., Rao, S., and Jayawardhana, B. (2015). Complex and detailed balancing of chemical reaction networks revisited. *Journal of Mathematical Chemistry*, 53(6), 1445–1458. doi:10.1007/s10910-015-0498-2.

* van Eunen, K., Kiewiet, J.A.L., Westerhoff, H.V., and Bakker, B.M. (2012). Testing Biochemistry Revisited: How In Vivo Metabolism Can Be Understood from In Vitro Enzyme Kinetics. *PLoS Computational Biology*, 8(4), e1002483. doi:10.1371/journal.pcbi.1002483.

* Varma, A. and Palsson, B.O. (1994). Metabolic Flux Balancing: Basic Concepts, Scientific and Practical Use. *Bio/Technology*, 12(10), 994–998. doi:10.1038/nbt1094-994.