# An Approach for Model Reduction of Biochemical Networks

Authors: Yen-Chang Liu; Chun-Liang Lin; Chia-Hua Chuang
Journal: Computational Biology Journal
 (2013)

---

## Abstract

Biochemical networks are not only complex but also extremely large. The dynamic biological model of great complexity resulting in a large number of parameters is a main difficulty for optimization and control processes. In practice, it is highly desirable to further simplify the structure of biological models for the sake of reducing computational cost or simplification for the task of system analysis. This paper considers the S-system model used for describing the response of biochemical networks. By introducing the technique of singular value decomposition (SVD), we are able to identify the major state variables and parameters and eliminate unimportant metabolites and the corresponding signal transduction pathways. The model reduction by multiobjective analysis integrates the criteria of reactive weight, sensitivity, and flux analyses to obtain a reduced model in a systematic way. The resultant model is closed to the original model in performance but with a simpler structure. Representative numerical examples are illustrated to prove feasibility of the proposed method.


---

## Body

## 1. Introduction

Metabolic models are usually with extremely high complexity with their system model because those networks are entangled by numerous molecules and reaction pathways [1–3]. The complexity of full-scale metabolic models is commonly cognized as a major obstacle for their effective use in computational systems biology. Direct analysis of that kind of systems based on the mathematical models is usually inefficient from the viewpoint of computational cost. To cope with the difficulty, one appropriate approach is to simplify the model before conducting the analysis that is, configuration of biological systems is simplified to reduce its actual size [4–6].In the mathematics, the stoichiometric matrix is commonly applied to express the framework of biochemical networks, which is built from the stoichiometric coefficients of the reactions. TheS-system model is a generalized representation of the dynamic behaviors of genetic networks or metabolic networks. Most results commence from considering the mathematical model of biochemical networks by S-system architecture [7, 8]. The model consists of a number of particular features of biochemical systems as well as with observation from their behaviors. The biochemical reactions of the S-system models in heterogeneous media show fractal kinetic orders [9].The technique of SVD has been used in the field of engineering for model simplification and data compression. In biochemical territory, the stoichiometric matrix can be disassembled into eigen-reaction (corresponding to the left singular vector), singular value, and eigen-connectivity (corresponding to the right singular vector) [10–13]. SVD can be utilized as a useful mathematical framework for processing and modeling genome-wide expression data, in which both of the mathematical variables and operations may be assigned to specific biological meanings [12, 13]. It is thus expected that the characteristics and information of the state variables can be peeked via the decomposed network structure.Sensitivity analysis combines the branch of mathematic and engineering responses. It shows how a biochemical system responds to short-term perturbations [3]. In general, the analysis is quantified via principal component analysis (PCA) [3, 11, 14], by which the upper bound of the tolerated parameter variations related to the system matrix of the biochemical network at the steady state can be obtained. In [14], the tolerance of the parameter uncertainty was estimated under which the system stability is preserved.This paper proposes a systematic model reduction method based on the SVD which considers three criteria (SVD, flux, and sensitivity analyses) to simplify complex rational rate expressions of the biochemical networks, such as enzymatic reactions. Our objective is to obtain biochemically meaningful reduced models, which reproduce the dynamic properties of the original model. The advantage of this method is that it is algorithmic and does not significantly rely on input in the form of biochemical knowledge. It is easy to apply and that allows for user-specified reductions of individual rate expressions in a full-order model. The approach, however, is better at preserving biochemical properties, as we show through some representative case studies. Availability endorses it as an alternative to achieve the classical engineering objective of improved parameter identifiability without losing the demand for preserving biochemical interpretation.


## 2. Multiple-Objective Approach

### 2.1. S-System Model

A typicalS-system model consists of a number of specific features of biochemical systems [3]. It provides easy systemic observations while compared with other biological modelling techniques. To explain, consider the reactive relation of influx υi+∈ℜ and efflux υi-∈ℜ represented, respectively, by
(1)x˙i=υi+-υi-=αi∏j=1r+mxjgij-βi∏j=1r+mxjhij,i=1,…,r,
where xj,j=1,…,r are dependent variables such as intermediate metabolites and products, xj,j=r+1,…,r+m are independent variables such as substrates and enzymes, xi is metabolic concentration, the positive numbers αi and βi refer to rate constants, and gij and hij are kinetic order variables. Positive kinetic orders indicate activating influences while negative values mean inhibition. If the kinetic order gij  (hij) equals zero, then the network will not contain the pathway from xj to xi. This model can then be casted into a canonical matrix form regardless of the degree of its nonlinearity:
(2)dx(t)dt=Sf(x(t)),
where x(t)∈ℜr is the concentration of reactants, S∈ℜr×2r denotes the stoichiometric matrix, and f(·)∈ℜ2r is a function vector.


### 2.2. Stoichiometric Matrix and SVD Analysis

The stoichiometric matrix is commonly utilized to describe the structure of a reaction network and SVD representation of the stoichiometric matrix can be used to analyze network properties. It is particularly a useful way to obtain the basic information about four fundamental subspaces of the matrix.Ifrank(S)=r, then S can be decomposed by SVD as
(3)S=UΣVT,
where U∈ℜr×r and V∈ℜ2r×2r are unimodular matrices, Σ=[diag(σ1,σ2,…,σr)0r×r]∈ℜr×2r is the matrix consisting of singular values (σi) and null entries, the elements of U are stoichiometric coefficients, and the elements of V are reaction participation numbers in model.To further investigate the effect of the flux rate on the biochemical network, the matrixSTS is diagonalized as
(4)STS=VΣTΣVT=VΛVT,
where V=[vj,p]∈ℜ2r×2r and Λ=diag(λ1,λ2,…,λr,0,…,0)∈ℜ2r×2r with λi=σi2 are the eigenvector matrix and eigenvalue matrix, respectively. From the representation of (4), it is clear that V involves reaction participation numbers and Λ consists of the square of singular values. Relative importance of the jth reaction can thus be measured via the weighted parameter ηj defined as
(5)ηj=∑p=1rλp|vj,p|∑p=1rλp,j=1,2,…,r,
where λp and vj,p have been defined in (4).


### 2.3. Sensitivity Analysis

Sensitivity analysis is introduced to explore the relationship between concentration and reactive rate. It helps to predict how the system performance is influenced by changes in key network parameters. In (1), if the biochemical network is at steady state, then we have
(6)αi∏j=1r+mxjgij=βi∏j=1r+mxjhij,i=1,2,…,r
or
(7)∑j=1r+m(gij-hij)lnxj=lnβi-lnαi,i=1,2,…,r.When the system remains close to the steady state, sensitivity analysis provides the result that a relative change in a rate constant affects the steady-state concentration of the metabolite. The idea provides perception of the structure ofS-system and can be employed in the evaluation of typical biochemical experiments. Consider the sensitivity of rate constants as
(8)zijα=∂lnxi∂lnαj,(9)zijβ=∂lnxi∂lnβj,
where zijα=-zijβ in the steady state; i,j=1,2,…,r by [3]. Define
(10)Z=[zijα]r×r=[∂lnxi∂lnαj]r×r
which can be decomposed as
(11)ZTZ=LDLT,
where L=[lq,l]∈ℜr×r is the similarity transformation matrix and D=diag(d1,d2,…,dr).The sensitivity coefficient can be obtained via (8). A good sensitivity analysis can distinguish each path by identifying relative importance of the parameters in the model. To this aim, define the relative importance of the jth reaction as
(12)wqα=wqβ=∑l=1rdl|lq,l|∑l=1rdl,q=1,2,…,r
and let
(13)wj={wqα,q=j+12,j=oddwqβ,q=j2,j=even,
where wqβ∈ℜ equals to wqα∈ℜ by ZTZ in (12). The factor wj quantitatively weights the parametric effects on the network features in (13). Therefore, sensitivity can be a measurement for quantitatively evaluating how fast the responses of any system components change with parameter variations.


### 2.4. Flux Analysis

In biochemical networks, the reactive input and output refer to the fluxυj. Biological metabolic networks have more fluxes (reactions) in the system than metabolites. We calculate the values of fluxes to get additional constraints for the metabolic network. For each metabolite one can obtain reactive rates by using the flux analysis. The flux consists of two patterns: influx and efflux. It can be further separated into activating and inhibition influences in the S-system. The result of conversion of the compounds is affected by the magnitude of flux rate. To compare each flux, we calculate, for each average flux f-j∈ℜ, as(14)f-j=1tf-t0|∫t0tfυjdt|,j=1,2,…,2r,
where tf and t0 are, respectively, the terminal and initial times for simulation; υj represents the input flux of the ((j+1)/2)th state when j is odd; otherwise, it represents the output flux of the (j/2)th state when j is even. From (14), one can compute the average flux rate on which each metabolic accumulation and degradation are quantified. From the calculated results, the relative magnitude of flux rate for each reaction pathway is perceived.


### 2.5. Model Reduction

Relative importance of each subobjective is not only subjective but also uneasily evaluated in a fair manner. The previous indices are thus integrated into a single objective for convenience of joint evaluation. Before integration, all objective functions are normalized. Forηj, consider ηmin≤ηj≤ηmax with 0≤(ηj-ηmin)/(ηmax-ηmin)≤1, where ηmin and ηmax are, respectively, the lower and upper limits of ηj. Similar treatments are applied for wj and f-j. Let
(15)Y=[γ1ηj-ηminηmax-ηminγ2wj-wminwmax-wminγ3f-j-f-minf-max-f-min]3×1,
where γi is the weight of the ith subobjective; wj and f-j are defined, respectively, in (13) and (14). Define
(16)Y^j={(γ3f-j-f-minf-max-f-min)s(γ1ηj-ηminηmax-ηmin)s+(γ2wj-wminwmax-wmin)sz+(γ3f-j-f-minf-max-f-min)s},1/s
where j=1,2,…,2r. The exponent s∈[1,∞) is used to reflect the general vector length. If s=1, the weighting represents the relative importance among each subobjective. Next, consider combination of the three objective functions by finding the cost value Y^j∈ℜ. It can help to distinguish important and unimportant parts:
(17)y=εmaxjY^j,
where y serves as a threshold and ε is a normalized factor to normalize the value of y. If Y^j<y, Y^j is categorized as the unimportant part and others as the important part. The unimportant part implies that there is weak or no connection between two states. Referring to the system model (1), if the parameters gij and hij are treated to zero, then, the network structure is reduced accordingly. The selection of ε is case dependent which is closely related to the network characteristics. One can determine an appropriate value of ε based on the simulation results. The principle for selection of ε is that, its value is such that (i) there are identifiable unimportant parts; (ii) the responses of the reduced model are not far away from the original model.After performing model reduction, one needs to ensure that there is no significant discrepancy between the original and reduced networks. We compute the residual errorRi∈ℜ via
(18)Ri=∫t0tf(xi(t)-xir(t))2dt∫t0tfxi2(t)dt,
where xir(t) denotes concentration of the reduced network at the sampling time instant t. If the error concentration is larger than the original system concentration, then the reduced system is unfeasible and could not be reduced. Furthermore, define
(19)R-=1r∑i=1rRi.
When R-→0+, it refers to the ideal situation and there is no residual error. On the other hand, when R-→1, it refers to the worst situation. Under the situation, the model should not be reduced.


## 2.1. S-System Model

A typicalS-system model consists of a number of specific features of biochemical systems [3]. It provides easy systemic observations while compared with other biological modelling techniques. To explain, consider the reactive relation of influx υi+∈ℜ and efflux υi-∈ℜ represented, respectively, by
(1)x˙i=υi+-υi-=αi∏j=1r+mxjgij-βi∏j=1r+mxjhij,i=1,…,r,
where xj,j=1,…,r are dependent variables such as intermediate metabolites and products, xj,j=r+1,…,r+m are independent variables such as substrates and enzymes, xi is metabolic concentration, the positive numbers αi and βi refer to rate constants, and gij and hij are kinetic order variables. Positive kinetic orders indicate activating influences while negative values mean inhibition. If the kinetic order gij  (hij) equals zero, then the network will not contain the pathway from xj to xi. This model can then be casted into a canonical matrix form regardless of the degree of its nonlinearity:
(2)dx(t)dt=Sf(x(t)),
where x(t)∈ℜr is the concentration of reactants, S∈ℜr×2r denotes the stoichiometric matrix, and f(·)∈ℜ2r is a function vector.


## 2.2. Stoichiometric Matrix and SVD Analysis

The stoichiometric matrix is commonly utilized to describe the structure of a reaction network and SVD representation of the stoichiometric matrix can be used to analyze network properties. It is particularly a useful way to obtain the basic information about four fundamental subspaces of the matrix.Ifrank(S)=r, then S can be decomposed by SVD as
(3)S=UΣVT,
where U∈ℜr×r and V∈ℜ2r×2r are unimodular matrices, Σ=[diag(σ1,σ2,…,σr)0r×r]∈ℜr×2r is the matrix consisting of singular values (σi) and null entries, the elements of U are stoichiometric coefficients, and the elements of V are reaction participation numbers in model.To further investigate the effect of the flux rate on the biochemical network, the matrixSTS is diagonalized as
(4)STS=VΣTΣVT=VΛVT,
where V=[vj,p]∈ℜ2r×2r and Λ=diag(λ1,λ2,…,λr,0,…,0)∈ℜ2r×2r with λi=σi2 are the eigenvector matrix and eigenvalue matrix, respectively. From the representation of (4), it is clear that V involves reaction participation numbers and Λ consists of the square of singular values. Relative importance of the jth reaction can thus be measured via the weighted parameter ηj defined as
(5)ηj=∑p=1rλp|vj,p|∑p=1rλp,j=1,2,…,r,
where λp and vj,p have been defined in (4).


## 2.3. Sensitivity Analysis

Sensitivity analysis is introduced to explore the relationship between concentration and reactive rate. It helps to predict how the system performance is influenced by changes in key network parameters. In (1), if the biochemical network is at steady state, then we have
(6)αi∏j=1r+mxjgij=βi∏j=1r+mxjhij,i=1,2,…,r
or
(7)∑j=1r+m(gij-hij)lnxj=lnβi-lnαi,i=1,2,…,r.When the system remains close to the steady state, sensitivity analysis provides the result that a relative change in a rate constant affects the steady-state concentration of the metabolite. The idea provides perception of the structure ofS-system and can be employed in the evaluation of typical biochemical experiments. Consider the sensitivity of rate constants as
(8)zijα=∂lnxi∂lnαj,(9)zijβ=∂lnxi∂lnβj,
where zijα=-zijβ in the steady state; i,j=1,2,…,r by [3]. Define
(10)Z=[zijα]r×r=[∂lnxi∂lnαj]r×r
which can be decomposed as
(11)ZTZ=LDLT,
where L=[lq,l]∈ℜr×r is the similarity transformation matrix and D=diag(d1,d2,…,dr).The sensitivity coefficient can be obtained via (8). A good sensitivity analysis can distinguish each path by identifying relative importance of the parameters in the model. To this aim, define the relative importance of the jth reaction as
(12)wqα=wqβ=∑l=1rdl|lq,l|∑l=1rdl,q=1,2,…,r
and let
(13)wj={wqα,q=j+12,j=oddwqβ,q=j2,j=even,
where wqβ∈ℜ equals to wqα∈ℜ by ZTZ in (12). The factor wj quantitatively weights the parametric effects on the network features in (13). Therefore, sensitivity can be a measurement for quantitatively evaluating how fast the responses of any system components change with parameter variations.


## 2.4. Flux Analysis

In biochemical networks, the reactive input and output refer to the fluxυj. Biological metabolic networks have more fluxes (reactions) in the system than metabolites. We calculate the values of fluxes to get additional constraints for the metabolic network. For each metabolite one can obtain reactive rates by using the flux analysis. The flux consists of two patterns: influx and efflux. It can be further separated into activating and inhibition influences in the S-system. The result of conversion of the compounds is affected by the magnitude of flux rate. To compare each flux, we calculate, for each average flux f-j∈ℜ, as(14)f-j=1tf-t0|∫t0tfυjdt|,j=1,2,…,2r,
where tf and t0 are, respectively, the terminal and initial times for simulation; υj represents the input flux of the ((j+1)/2)th state when j is odd; otherwise, it represents the output flux of the (j/2)th state when j is even. From (14), one can compute the average flux rate on which each metabolic accumulation and degradation are quantified. From the calculated results, the relative magnitude of flux rate for each reaction pathway is perceived.


## 2.5. Model Reduction

Relative importance of each subobjective is not only subjective but also uneasily evaluated in a fair manner. The previous indices are thus integrated into a single objective for convenience of joint evaluation. Before integration, all objective functions are normalized. Forηj, consider ηmin≤ηj≤ηmax with 0≤(ηj-ηmin)/(ηmax-ηmin)≤1, where ηmin and ηmax are, respectively, the lower and upper limits of ηj. Similar treatments are applied for wj and f-j. Let
(15)Y=[γ1ηj-ηminηmax-ηminγ2wj-wminwmax-wminγ3f-j-f-minf-max-f-min]3×1,
where γi is the weight of the ith subobjective; wj and f-j are defined, respectively, in (13) and (14). Define
(16)Y^j={(γ3f-j-f-minf-max-f-min)s(γ1ηj-ηminηmax-ηmin)s+(γ2wj-wminwmax-wmin)sz+(γ3f-j-f-minf-max-f-min)s},1/s
where j=1,2,…,2r. The exponent s∈[1,∞) is used to reflect the general vector length. If s=1, the weighting represents the relative importance among each subobjective. Next, consider combination of the three objective functions by finding the cost value Y^j∈ℜ. It can help to distinguish important and unimportant parts:
(17)y=εmaxjY^j,
where y serves as a threshold and ε is a normalized factor to normalize the value of y. If Y^j<y, Y^j is categorized as the unimportant part and others as the important part. The unimportant part implies that there is weak or no connection between two states. Referring to the system model (1), if the parameters gij and hij are treated to zero, then, the network structure is reduced accordingly. The selection of ε is case dependent which is closely related to the network characteristics. One can determine an appropriate value of ε based on the simulation results. The principle for selection of ε is that, its value is such that (i) there are identifiable unimportant parts; (ii) the responses of the reduced model are not far away from the original model.After performing model reduction, one needs to ensure that there is no significant discrepancy between the original and reduced networks. We compute the residual errorRi∈ℜ via
(18)Ri=∫t0tf(xi(t)-xir(t))2dt∫t0tfxi2(t)dt,
where xir(t) denotes concentration of the reduced network at the sampling time instant t. If the error concentration is larger than the original system concentration, then the reduced system is unfeasible and could not be reduced. Furthermore, define
(19)R-=1r∑i=1rRi.
When R-→0+, it refers to the ideal situation and there is no residual error. On the other hand, when R-→1, it refers to the worst situation. Under the situation, the model should not be reduced.


## 3. Simulation Results

Data analysis of large reaction networks is usually time consuming due to higher computational cost. A variety of reduced models have been developed for the purpose [5, 14]. An efficient way to analyze a large-sized mathematical model is to decompose it into a smaller one which is easier to manage and analyze, analogous to “divide-and-conquer” method. However, the reduced module should be combined to produce the same set of inputs and outputs as those in the original model. Simply put, precision, and fidelity should not give in simplicity.


### 3.1. Case 1: Tricarboxylic Acid Cycle

In Figure1, the tricarboxylic acid cycle (TCA cycle), also called Krebs cycle and citric acid cycle, is a series of enzyme catalysed chemical reactions. TCA cycle is considered an important organism reaction, that pyruvate oxidation, carbohydrate, protein, and fatty acid are transformed into energy in the cell mitochondria. Acetyl-CoA requires Oxaloacetate (OAA) to continuously decompose and release CO2. If the carrier cannot be recycled, then the carrier will exhaust and the reaction will not be activated. Therefore, intermediate is constantly translated to carrier by the reaction cycle. The network possesses 13 metabolites (dependent variable), while the 35 enzymes are accounted as independent variables and coded as x14 through x48. The dynamic equations can be expressed as [15]
(20)x˙1=0.8282x1-0.038x6-0.0204x70.106x100.114x130.7x180.807×x280.108x310.0848x460.599x48-0.181-1.3423x1x300.915x330.0847,x˙2=1.3401x10.915x70.0848x230.0848x300.915-17.166x20.706x30.0716x220.0848x240.915x47-0.0341,x˙3=0.3231x3-0.405x50.156x210.427x350.573x460.422x470.405x48-0.418-9.6952x20.376x30.489x240.554x410.446x47-0.00206,x˙4=x9x25-0.152x40.958x26x460.0348x48-0.862,x˙5=1.875x70.0274x80.465x130.336x190.535x200.465-0.01923x3-0.717x50.413x60.306x8-0.29x10-0.0883×x210.756x290.244x460.748x470.718x48-0.741,x˙6=2.459x1-0.00921x6-0.0154x70.0162x100.086x110.276×x280.813x320.276x390.6413-1.1528x50.0963x61.01x8-0.204x10-0.062×x140.0518x270.277x290.171x450.5x460.0222x48-0.0191,x˙7=2.1167x10.129x20.129x220.129x330.129x340.741-3.4893x1-0.0187x6-0.0311x70.868x100.174x230.129×x280.165x310.129x400.577,x˙8=0.5724x50.111x60.247x80.234x100.0713x290.197x380.803-1.9369x8x200.375x440.625,x˙9=16.242x20.679x30.0782x24x47-0.0372-x9x25,x˙10=0.156x40.724x50.106x60.259x8-0.223x10-0.0679x140.0568×x260.756x290.188x460.0506x48-0.672-0.8063x1-0.0101x6-0.0168×x70.0177x100.99x11-0.879x150.911x280.0891x460.882x470.879x48-0.881,x˙11=2.0031x60.166x100.491x11-0.481x150.499x270.166×x360.335x460.483x470.481x48-0.483-2.4373x110.495x12-0.00542x160.574x320.166x420.261,x˙12=1.271x110.106x12-0.00836x160.885x370.115-9.1694x121.89x13-1.24x170.911x430.0893,x˙13=8.289x121.98x13-1.36x17-0.9387x1-0.0197x70.0196×x130.775x180.618x190.382x460.458x48-0.139.Figure 1
Tricarboxylic acid cycle network inDictyostelium discoideum.One is referred to [3] for the definitions of all state variables depicted above. Figure 2 illustrates the result of reactive weight, sensitivity analysis, flux rates, and compromise solution. The entries of the matrices Λ and D are given below:
(21)Λ=diag(296.47,264.8,94.1,85.69,69.59,16.66,9.95,ZZz7.38,4.08,3.52,2.49,1.02,0.67,0,…,0),D=diag(0.07,0.27,0.34,0.43,0.47,0.77,1,Zzz1.38,2.2,3,8.73,67.8,144.34),
from which we obtain the maximal value of Y^j to be 0.0534, with which, 19 pathways and one metabolite could be eliminated based on the proposed method while setting y=0.0534 in Tables 1 and 2 or the dashed part of Figure 1. The reduced dynamic model is given by
(22)x˙1=0.8282x1-0.038x6-0.0204x100.114x130.7x180.807x280.108×x310.0848x460.599x48-0.181-1.3423x1x300.915x330.0847,x˙2=1.3401x10.915x230.0848x300.915-17.166x20.706x30.0716x220.0848x240.915x47-0.0341,x˙3=0.3231x3-0.405x50.156x210.427x350.573x460.422x470.405x48-0.418-9.6952x20.376x30.489x240.554x410.446x47-0.00206,x˙4=x9x25-0.152x40.958x26x460.0348x48-0.862,x˙5=1.875x80.465x130.336x190.535x200.465-0.01923x3-0.717x50.413x60.306x8-0.29x10-0.0883×x210.756x290.244x460.748x470.718x48-0.741,x˙6=2.459x1-0.00921x6-0.0154x100.086x110.276×x280.813x320.276x390.6413-1.1528x50.0963x61.01x8-0.204×x10-0.062x140.0518x270.277x290.171x450.5x460.0222x48-0.0191,x˙8=0.5724x50.111x60.247x80.234x100.0713x290.197x380.803-1.9369x8x200.375x440.625,x˙9=16.242x20.679x30.0782x24x47-0.0372-x9x25,x˙10=0.156x40.724x50.106x60.259x8-0.223x10-0.0679x140.0568×x260.756x290.188x460.0506x48-0.672-0.8063x1-0.0101x6-0.0168×x100.99x11-0.879x150.911x280.0891x460.882x470.879x48-0.881,x˙11=2.0031x60.166x100.491x11-0.481x150.499×x270.166x360.335x460.483x470.481x48-0.483-2.4373x110.495x12-0.00542x160.574x320.166x420.261,x˙12=1.271x110.106x12-0.00836x160.885x370.115-9.1694x121.89x13-1.24x170.911x430.0893,x˙13=8.289x121.98x13-1.36x17-0.9387x1-0.0197×x130.775x180.618x190.382x460.458x48-0.139.Table 1
Result of model reduction byε in Case  1 (IMP as important part, UMP as unimportant part).
ε
y
IMP
UMP
Average residual error
0.1
0.0534
24
2
7.16
×
1
0
-
2
0.3
0.1602
20
6
337.0714
0.5
0.267
14
12
5.0268
×
1
0
17Table 2
Result of model reduction in Case  1 (ε=0.1).
Reduced pathways
Reduced concentration
υ
1
→
7
,
υ
2
→
7
,
υ
10
→
7
,
υ
13
→
7
,
υ
22
→
7
,
υ
33
→
7,
υ
34
→
7
,
υ
7
→
1
,
υ
7
→
2
,
υ
7
→
5
,
υ
7
→
6
,
υ
7
→
10,
x
7
υ
7
→
23
,
υ
7
→
28
,
υ
7
→
31
,
υ
7
→
40Graph of (a) reactive weight, (b) sensitivity analysis, (c) flux rates, and (d) compromise solution data in Case 1 (γ1=0.4,γ2=0.3,γ3=0.3).
(a)
(b)
(c)
(d)Note thatx7 has been removed from the original model. The effect of model reduction can be confirmed from observation of Figure 3. It is seen that the responses of concentrations for the original model and reduction model are quite consistent. With regard to the computational aspect, we conducted the simulation study using MATLAB code for Runge-Kutta 4th order on a PC featured with a 2.33 GHz microprocessor (Q8200). For this case, expenditure of the computation time for the original and reduction models is, respectively, 14.8 and 13.6 seconds. The latter saves the computation cost up to 9%.Dynamic responses of the noisy states in Case 1: (a) original model and (b) simplified model (xi as original concentration, xri as reduced concentration).
(a)
(b)


### 3.2. Case 2: Purine Metabolism

In Figure4, the concentration (dependent variable) possesses 16 species. We illustrate the procedure with the dynamics of PRPP, which is affected by 6 reactions. The synthesizing reaction has R5P as substrate, which is explicitly included as an independent variable, and is modulated by numerous inhibitors. One is PRPP itself, and four others of relevance are AMP, ADP, GDP, and GTP. One uses the aggregated pools x4 and x8 as substitutes and adjust the kinetic orders accordingly, because the latter are not represented individually in the model [3]:
(23)x˙1=0.898x1-0.03x4-0.45x8-0.04x170.65x180.7-31.335x11.27x2-0.165×x4-0.067x60.036x8-0.237x130.085x150.075x18-0.0092,x˙2=1.629x10.72x2-0.28x40.319x7-0.03x8-0.024x130.144x18-0.06-5.584x20.454x4-0.156x7-0.012x80.126x18-0.11,x˙3=3.5932x20.4x4-0.24x80.2x18-0.05-66544.7x30.99x4-0.95,x˙4=0.082082x10.00025x30.004x4-0.00418x50.0023x60.00037x110.988-614.27x40.054x5-0.0045x80.128x9-0.00003x100.00087x18-0.00028,x˙5=7.2067x40.2x5-0.6-9.006x50.368,x˙6=0.29x50.9-220.14x10.495x4-0.79x60.748,x˙7=1.2823x20.15x7-0.09x8-0.03-0.3738x40.12x70.16,x˙8=0.048842x10.0034x40.00014x70.00019x8-0.0034x110.996x150.001-408.11x2-0.000057x40.0498x7-0.00029x80.133x9-0.00009×x10-0.000029x18-0.001,x˙9=0.00229x40.002x9-0.0058x100.017x120.98-3.3x90.431x100.324,x˙10=0.00151x80.0058x9-0.017x10-0.0058x120.986-2.2655x90.413x100.34,x˙11=1024.195x40.05x80.13-0.115389x11,x˙12=5.5062x90.42x100.33-0.0032545x12,x˙13=0.1033x20.43x40.409x90.04x18-0.19-8.2827x10.817x2-0.66x130.528,x˙14=0.7234x130.343x150.236-0.9325x140.569,x˙15=0.2615x80.881x100.02x18-0.33-136.6x10.925x8-0.925x150.438,x˙16=0.949x140.55-0.00008744x162.21.Figure 4
Graph of purine metabolite in which curved arrows enter or leave the pathway corresponding to purine ring and ribose moieties.For the definitions of all state variables, see [3]. Simulation results of the network are illustrated in Figures 5 and 6. The entries of the matrices Λ and D are
(24)Λ=diag(296.47,264.8,94.1,85.69,69.59,16.66,9.95,ZZz7.38,4.08,3.52,2.49,1.02,0.67,0,…,0),D=diag(0.07,0.27,0.34,0.43,0.47,0.77,1,Zzz1.38,2.2,3,8.73,67.8,144.34),Graph of (a) reactive weight, (b) sensitivity analysis, (c) flux rates, and (d) compromise solution data in Case 2 (γ1=0.4,γ2=0.3,γ3=0.3).
(a)
(b)
(c)
(d)Figure 6
Dynamic responses of the noisy states in Case 2 wherexi denotes the original concentration and xri denotes the reduced concentration: (a) original model (in blue line) and (b) simplified model (in green line).As in the previous case, it can be observed from Tables3 and 4 that 10 pathways and 3 metabolites are ignorable when choosing y=0.0003. This makes it clear that three input fluxes and three output fluxes (six UMPs) are ignorable. Ignoring the six UMPs is equivalent to eliminating three state variables x6, x14, and x16. The six dashed lines in Figure 4 show the pathways to be ignorable after performing the analysis. The corresponding reduced model is given by
(25)x˙1=0.898x1-0.03x4-0.45x8-0.04x170.65x180.7-31.335x11.27x2-0.165x4-0.067x60.036x8-0.237x130.085×x150.075x18-0.0092,x˙2=1.629x10.72x2-0.28x40.319x7-0.03x8-0.024x130.144x18-0.06-5.584x20.454x4-0.156x7-0.012x80.126x18-0.11,x˙3=3.5932x20.4x4-0.24x80.2x18-0.05-66544.7x30.99x4-0.95,x˙4=0.082082x10.00025x30.004x4-0.00418x50.0023x60.00037x110.988-614.27x40.054x5-0.0045x80.128x9-0.00003x100.00087x18-0.00028,x˙5=7.2067x40.2x5-0.6-9.006x50.368,x˙7=1.2823x20.15x7-0.09x8-0.03-0.3738x40.12x70.16,x˙8=0.048842x10.0034x40.00014x70.00019x8-0.0034x110.996x150.001-408.11x2-0.000057x40.0498x7-0.00029x80.133x9-0.00009×x10-0.000029x18-0.001,x˙9=0.00229x40.002x9-0.0058x100.017x120.98-3.3x90.431x100.324,x˙10=0.00151x80.0058x9-0.017x10-0.0058x120.986-2.2655x90.413x100.34,x˙11=1024.195x40.05x80.13-0.115389x11,x˙12=5.5062x90.42x100.33-0.0032545x12,x˙13=0.1033x20.43x40.409x90.04x18-0.19-8.2827x10.817x2-0.66x130.528,x˙15=0.2615x80.881x100.02x18-0.33-136.6x10.925x8-0.925x150.438.Table 3
Result of model reduction byε in Case  2 (IMP as important part, UMP as unimportant part).
ε
y
IMP
UMP
Average residual error
0.0002
1.2002
×
1
0
-
4
30
2
0
0.0005
3.0005
×
1
0
-
4
26
6
1.041
×
1
0
-
5
0.0008
4.8008
×
1
0
-
4
24
8
0.0918Table 4
Result of model reduction in Case  2 (ε=0.0005).
Reduced pathways
Reduced concentration
υ
6
→
1
,
υ
5
→
6
,
υ
6
→
4
,
υ
13
→
14
,
x
6
,
x
14
,
x
16
υ
15
→
14
,
υ
1
→
6
,
υ
14
→
16The computational costs for the original model and reduction model are, respectively, 1 and 0.9 seconds. The reduced model saves the computation cost up to 10%. Note that this case possesses simpler interconnected structure than that of Case 1, therefore, the computation time is far less than that in Case 1; We summarize in Figure7 the residual errors between the original and reduced models for the two cases above. Again, the result confirms the proposed approach.Residual errors of model reduction for Cases 1 and 2.
(a)
Case 1
(b)
Case 2


### 3.3. Case 3: Saccharomyces cerevisiae Network

Organisms use glucose to produce ethanol, glycerol, glycogen, and trehalose in the yeast fermentation pathway. For the map illustrated in Figure8, it has an important point of regulation that external glucose is transported into the cell. Glucose-6-phosphate (G6P) inhibits the step. Internal glucose is phosphorylated to G6P by hexokinase and G6P flow into distinct pathways. Here, one only considers two terms. The first one is toward ethanol, and the next one is toward the polysaccharides glycogen and trehalose. This is insignificant for the flow into the oxidative and nonoxidative pentose pathways. The process is omitted under the experimental conditions of interest. The omission allows less kinetic information to establish the model. Thus, it is a major modeling step. However, it exhibits a drawback, because it has conditions to limit model predictions that the pentose branches can be ignored in comparison with the ethanol branch. In mathematical complexity, the gained reduction is weighed with the limitation in scope.Figure 8
Map of anaerobic fermentation pathway inSaccharomyces cerevisiae.G6P and fructose-6-phosphate can be converted into each other. Thus, G6P and fructose-6-phosphate are considered to be at equilibrium. In the next step, fructose-6-phosphate is phosphorylated by phosphofructokinase. If the enzyme fructose-1,6-diphosphate is active, it will lead to a futile cycle between fructose-6-phosphate and fructose-1, 6-diphosphate. However, this enzyme is typically inactive under the experimental conditions of interest. And the model cannot include the cycle. Glycerol or phosphoenolpyruvate is produced by fructose-1,6-diphosphate. Finally, the production of ethanol is catalyzed by pyruvate kinase. In addition to these major reaction pathways, there are a number of reactions between the ADP and ATP procedures. The model can be mathematically expressed as follows [3]:
(26)x˙1=0.9018x2-0.2344x6-3.1833x10.7464x50.0243x7,x˙2=3.1833x10.7464x50.0243x7-0.8187x20.664x5-0.354x80.899×x110.0008,x˙3=0.5234x20.7318x5-0.3941x8-0.0148x30.584x40.03x50.119x90.944×x120.056x14-0.575,x˙4=0.022x30.6159x50.1308x9x14-0.6088-0.0945x30.05x40.533x5-0.0822×x10,x˙5=0.088x30.3505x40.2665x50.0243x90.5x100.5x14-0.3044-3.327x10.2134x20.19x50.362x70.286x80.257x110.0002x130.457.In Figure9, some important messages from the graphs were discerned. The reactions of PEP were eliminated from the simulated result. However, the result obtained is not satisfactory as expected. To explore the cause of failure, we conduct a comparison among examples to identify the differences. One can use the value of average residual R- to determine suitability of the model simplification. For this case, R- is 2.8417. Obviously, it is too large to endorse suitability of the method for model simplification.Graph of (a) reactive weight, (b) sensitivity analysis, (c) flux rates, and (d) compromise solution data in Case 3 (γ1=0.4,γ2=0.3,γ3=0.3).
(a)
(b)
(c)
(d)


### 3.4. Case 4: Dictyostelium Network

Figure10 illustrates pulses of cAMP produced by ACA. ACA is activated by the binding of extracellular cAMP to the surface receptor CAR1. The protein kinase ERK2 is activated to transmit the signal to ACA via ligand-bound CAR1. The activation of ACA is independent of ERK2 in an alternate circuit. If cAMP is accumulated internally, then the protein kinase PKA is activated via binding to the regulatory subunit of PKA. PKA inhibits ERK2 and the cAMP phosphodiesterase REG is activated by phosphorylating it. REG A is activated by a protein phosphatase. The process makes internal cAMP hydrolyze via REG A. When REGA hydrolyzes the internal cAMP, its regulatory subunit inactivates PKA activity. And CAR1 is returned to its high-affinity state by protein phosphatase. The phosphodiesterase PDE degrades secreted cAMP, before it diffuses between cells. The double horizontal line expresses the membrane surface of the cell [16, 17].Figure 10
Dictyostelium network.The mathematical model is given by(27)x˙1=2x7-0.9x1x2,x˙2=2.5x5-1.5x2,x˙3=0.6x7-0.8x2x3,x˙4=1-1.3x3x4,x˙5=0.3x1-0.8x4x5,x˙6=0.7x1-4.9x6,x˙7=23x6-4.5x7.Figure11 shows the results for the respective analyses of this case in which the average residual error value is 0.9488. For Cases 1–4, the values of R- are 7.16×10-2, 1.041×10-5, 2.8417, and 0.9488, respectively. It is quite obvious that the value of R- in Cases 3 and 4 is large than that in other cases. Thus, the model cannot be reduced further. It is quite obvious that not every biochemical network can be simplified in its structural model. However, one can quantify the mathematical model of the metabolic network while comparing relative importance of each reaction pathway.Graph of (a) reactive weight, (b) sensitivity analysis, (c) flux rates, and (d) compromise solution data in Case 4 (γ1=0.4,γ2=0.3,γ3=0.3).
(a)
(b)
(c)
(d)


## 3.1. Case 1: Tricarboxylic Acid Cycle

In Figure1, the tricarboxylic acid cycle (TCA cycle), also called Krebs cycle and citric acid cycle, is a series of enzyme catalysed chemical reactions. TCA cycle is considered an important organism reaction, that pyruvate oxidation, carbohydrate, protein, and fatty acid are transformed into energy in the cell mitochondria. Acetyl-CoA requires Oxaloacetate (OAA) to continuously decompose and release CO2. If the carrier cannot be recycled, then the carrier will exhaust and the reaction will not be activated. Therefore, intermediate is constantly translated to carrier by the reaction cycle. The network possesses 13 metabolites (dependent variable), while the 35 enzymes are accounted as independent variables and coded as x14 through x48. The dynamic equations can be expressed as [15]
(20)x˙1=0.8282x1-0.038x6-0.0204x70.106x100.114x130.7x180.807×x280.108x310.0848x460.599x48-0.181-1.3423x1x300.915x330.0847,x˙2=1.3401x10.915x70.0848x230.0848x300.915-17.166x20.706x30.0716x220.0848x240.915x47-0.0341,x˙3=0.3231x3-0.405x50.156x210.427x350.573x460.422x470.405x48-0.418-9.6952x20.376x30.489x240.554x410.446x47-0.00206,x˙4=x9x25-0.152x40.958x26x460.0348x48-0.862,x˙5=1.875x70.0274x80.465x130.336x190.535x200.465-0.01923x3-0.717x50.413x60.306x8-0.29x10-0.0883×x210.756x290.244x460.748x470.718x48-0.741,x˙6=2.459x1-0.00921x6-0.0154x70.0162x100.086x110.276×x280.813x320.276x390.6413-1.1528x50.0963x61.01x8-0.204x10-0.062×x140.0518x270.277x290.171x450.5x460.0222x48-0.0191,x˙7=2.1167x10.129x20.129x220.129x330.129x340.741-3.4893x1-0.0187x6-0.0311x70.868x100.174x230.129×x280.165x310.129x400.577,x˙8=0.5724x50.111x60.247x80.234x100.0713x290.197x380.803-1.9369x8x200.375x440.625,x˙9=16.242x20.679x30.0782x24x47-0.0372-x9x25,x˙10=0.156x40.724x50.106x60.259x8-0.223x10-0.0679x140.0568×x260.756x290.188x460.0506x48-0.672-0.8063x1-0.0101x6-0.0168×x70.0177x100.99x11-0.879x150.911x280.0891x460.882x470.879x48-0.881,x˙11=2.0031x60.166x100.491x11-0.481x150.499x270.166×x360.335x460.483x470.481x48-0.483-2.4373x110.495x12-0.00542x160.574x320.166x420.261,x˙12=1.271x110.106x12-0.00836x160.885x370.115-9.1694x121.89x13-1.24x170.911x430.0893,x˙13=8.289x121.98x13-1.36x17-0.9387x1-0.0197x70.0196×x130.775x180.618x190.382x460.458x48-0.139.Figure 1
Tricarboxylic acid cycle network inDictyostelium discoideum.One is referred to [3] for the definitions of all state variables depicted above. Figure 2 illustrates the result of reactive weight, sensitivity analysis, flux rates, and compromise solution. The entries of the matrices Λ and D are given below:
(21)Λ=diag(296.47,264.8,94.1,85.69,69.59,16.66,9.95,ZZz7.38,4.08,3.52,2.49,1.02,0.67,0,…,0),D=diag(0.07,0.27,0.34,0.43,0.47,0.77,1,Zzz1.38,2.2,3,8.73,67.8,144.34),
from which we obtain the maximal value of Y^j to be 0.0534, with which, 19 pathways and one metabolite could be eliminated based on the proposed method while setting y=0.0534 in Tables 1 and 2 or the dashed part of Figure 1. The reduced dynamic model is given by
(22)x˙1=0.8282x1-0.038x6-0.0204x100.114x130.7x180.807x280.108×x310.0848x460.599x48-0.181-1.3423x1x300.915x330.0847,x˙2=1.3401x10.915x230.0848x300.915-17.166x20.706x30.0716x220.0848x240.915x47-0.0341,x˙3=0.3231x3-0.405x50.156x210.427x350.573x460.422x470.405x48-0.418-9.6952x20.376x30.489x240.554x410.446x47-0.00206,x˙4=x9x25-0.152x40.958x26x460.0348x48-0.862,x˙5=1.875x80.465x130.336x190.535x200.465-0.01923x3-0.717x50.413x60.306x8-0.29x10-0.0883×x210.756x290.244x460.748x470.718x48-0.741,x˙6=2.459x1-0.00921x6-0.0154x100.086x110.276×x280.813x320.276x390.6413-1.1528x50.0963x61.01x8-0.204×x10-0.062x140.0518x270.277x290.171x450.5x460.0222x48-0.0191,x˙8=0.5724x50.111x60.247x80.234x100.0713x290.197x380.803-1.9369x8x200.375x440.625,x˙9=16.242x20.679x30.0782x24x47-0.0372-x9x25,x˙10=0.156x40.724x50.106x60.259x8-0.223x10-0.0679x140.0568×x260.756x290.188x460.0506x48-0.672-0.8063x1-0.0101x6-0.0168×x100.99x11-0.879x150.911x280.0891x460.882x470.879x48-0.881,x˙11=2.0031x60.166x100.491x11-0.481x150.499×x270.166x360.335x460.483x470.481x48-0.483-2.4373x110.495x12-0.00542x160.574x320.166x420.261,x˙12=1.271x110.106x12-0.00836x160.885x370.115-9.1694x121.89x13-1.24x170.911x430.0893,x˙13=8.289x121.98x13-1.36x17-0.9387x1-0.0197×x130.775x180.618x190.382x460.458x48-0.139.Table 1
Result of model reduction byε in Case  1 (IMP as important part, UMP as unimportant part).
ε
y
IMP
UMP
Average residual error
0.1
0.0534
24
2
7.16
×
1
0
-
2
0.3
0.1602
20
6
337.0714
0.5
0.267
14
12
5.0268
×
1
0
17Table 2
Result of model reduction in Case  1 (ε=0.1).
Reduced pathways
Reduced concentration
υ
1
→
7
,
υ
2
→
7
,
υ
10
→
7
,
υ
13
→
7
,
υ
22
→
7
,
υ
33
→
7,
υ
34
→
7
,
υ
7
→
1
,
υ
7
→
2
,
υ
7
→
5
,
υ
7
→
6
,
υ
7
→
10,
x
7
υ
7
→
23
,
υ
7
→
28
,
υ
7
→
31
,
υ
7
→
40Graph of (a) reactive weight, (b) sensitivity analysis, (c) flux rates, and (d) compromise solution data in Case 1 (γ1=0.4,γ2=0.3,γ3=0.3).
(a)
(b)
(c)
(d)Note thatx7 has been removed from the original model. The effect of model reduction can be confirmed from observation of Figure 3. It is seen that the responses of concentrations for the original model and reduction model are quite consistent. With regard to the computational aspect, we conducted the simulation study using MATLAB code for Runge-Kutta 4th order on a PC featured with a 2.33 GHz microprocessor (Q8200). For this case, expenditure of the computation time for the original and reduction models is, respectively, 14.8 and 13.6 seconds. The latter saves the computation cost up to 9%.Dynamic responses of the noisy states in Case 1: (a) original model and (b) simplified model (xi as original concentration, xri as reduced concentration).
(a)
(b)


## 3.2. Case 2: Purine Metabolism

In Figure4, the concentration (dependent variable) possesses 16 species. We illustrate the procedure with the dynamics of PRPP, which is affected by 6 reactions. The synthesizing reaction has R5P as substrate, which is explicitly included as an independent variable, and is modulated by numerous inhibitors. One is PRPP itself, and four others of relevance are AMP, ADP, GDP, and GTP. One uses the aggregated pools x4 and x8 as substitutes and adjust the kinetic orders accordingly, because the latter are not represented individually in the model [3]:
(23)x˙1=0.898x1-0.03x4-0.45x8-0.04x170.65x180.7-31.335x11.27x2-0.165×x4-0.067x60.036x8-0.237x130.085x150.075x18-0.0092,x˙2=1.629x10.72x2-0.28x40.319x7-0.03x8-0.024x130.144x18-0.06-5.584x20.454x4-0.156x7-0.012x80.126x18-0.11,x˙3=3.5932x20.4x4-0.24x80.2x18-0.05-66544.7x30.99x4-0.95,x˙4=0.082082x10.00025x30.004x4-0.00418x50.0023x60.00037x110.988-614.27x40.054x5-0.0045x80.128x9-0.00003x100.00087x18-0.00028,x˙5=7.2067x40.2x5-0.6-9.006x50.368,x˙6=0.29x50.9-220.14x10.495x4-0.79x60.748,x˙7=1.2823x20.15x7-0.09x8-0.03-0.3738x40.12x70.16,x˙8=0.048842x10.0034x40.00014x70.00019x8-0.0034x110.996x150.001-408.11x2-0.000057x40.0498x7-0.00029x80.133x9-0.00009×x10-0.000029x18-0.001,x˙9=0.00229x40.002x9-0.0058x100.017x120.98-3.3x90.431x100.324,x˙10=0.00151x80.0058x9-0.017x10-0.0058x120.986-2.2655x90.413x100.34,x˙11=1024.195x40.05x80.13-0.115389x11,x˙12=5.5062x90.42x100.33-0.0032545x12,x˙13=0.1033x20.43x40.409x90.04x18-0.19-8.2827x10.817x2-0.66x130.528,x˙14=0.7234x130.343x150.236-0.9325x140.569,x˙15=0.2615x80.881x100.02x18-0.33-136.6x10.925x8-0.925x150.438,x˙16=0.949x140.55-0.00008744x162.21.Figure 4
Graph of purine metabolite in which curved arrows enter or leave the pathway corresponding to purine ring and ribose moieties.For the definitions of all state variables, see [3]. Simulation results of the network are illustrated in Figures 5 and 6. The entries of the matrices Λ and D are
(24)Λ=diag(296.47,264.8,94.1,85.69,69.59,16.66,9.95,ZZz7.38,4.08,3.52,2.49,1.02,0.67,0,…,0),D=diag(0.07,0.27,0.34,0.43,0.47,0.77,1,Zzz1.38,2.2,3,8.73,67.8,144.34),Graph of (a) reactive weight, (b) sensitivity analysis, (c) flux rates, and (d) compromise solution data in Case 2 (γ1=0.4,γ2=0.3,γ3=0.3).
(a)
(b)
(c)
(d)Figure 6
Dynamic responses of the noisy states in Case 2 wherexi denotes the original concentration and xri denotes the reduced concentration: (a) original model (in blue line) and (b) simplified model (in green line).As in the previous case, it can be observed from Tables3 and 4 that 10 pathways and 3 metabolites are ignorable when choosing y=0.0003. This makes it clear that three input fluxes and three output fluxes (six UMPs) are ignorable. Ignoring the six UMPs is equivalent to eliminating three state variables x6, x14, and x16. The six dashed lines in Figure 4 show the pathways to be ignorable after performing the analysis. The corresponding reduced model is given by
(25)x˙1=0.898x1-0.03x4-0.45x8-0.04x170.65x180.7-31.335x11.27x2-0.165x4-0.067x60.036x8-0.237x130.085×x150.075x18-0.0092,x˙2=1.629x10.72x2-0.28x40.319x7-0.03x8-0.024x130.144x18-0.06-5.584x20.454x4-0.156x7-0.012x80.126x18-0.11,x˙3=3.5932x20.4x4-0.24x80.2x18-0.05-66544.7x30.99x4-0.95,x˙4=0.082082x10.00025x30.004x4-0.00418x50.0023x60.00037x110.988-614.27x40.054x5-0.0045x80.128x9-0.00003x100.00087x18-0.00028,x˙5=7.2067x40.2x5-0.6-9.006x50.368,x˙7=1.2823x20.15x7-0.09x8-0.03-0.3738x40.12x70.16,x˙8=0.048842x10.0034x40.00014x70.00019x8-0.0034x110.996x150.001-408.11x2-0.000057x40.0498x7-0.00029x80.133x9-0.00009×x10-0.000029x18-0.001,x˙9=0.00229x40.002x9-0.0058x100.017x120.98-3.3x90.431x100.324,x˙10=0.00151x80.0058x9-0.017x10-0.0058x120.986-2.2655x90.413x100.34,x˙11=1024.195x40.05x80.13-0.115389x11,x˙12=5.5062x90.42x100.33-0.0032545x12,x˙13=0.1033x20.43x40.409x90.04x18-0.19-8.2827x10.817x2-0.66x130.528,x˙15=0.2615x80.881x100.02x18-0.33-136.6x10.925x8-0.925x150.438.Table 3
Result of model reduction byε in Case  2 (IMP as important part, UMP as unimportant part).
ε
y
IMP
UMP
Average residual error
0.0002
1.2002
×
1
0
-
4
30
2
0
0.0005
3.0005
×
1
0
-
4
26
6
1.041
×
1
0
-
5
0.0008
4.8008
×
1
0
-
4
24
8
0.0918Table 4
Result of model reduction in Case  2 (ε=0.0005).
Reduced pathways
Reduced concentration
υ
6
→
1
,
υ
5
→
6
,
υ
6
→
4
,
υ
13
→
14
,
x
6
,
x
14
,
x
16
υ
15
→
14
,
υ
1
→
6
,
υ
14
→
16The computational costs for the original model and reduction model are, respectively, 1 and 0.9 seconds. The reduced model saves the computation cost up to 10%. Note that this case possesses simpler interconnected structure than that of Case 1, therefore, the computation time is far less than that in Case 1; We summarize in Figure7 the residual errors between the original and reduced models for the two cases above. Again, the result confirms the proposed approach.Residual errors of model reduction for Cases 1 and 2.
(a)
Case 1
(b)
Case 2


## 3.3. Case 3: Saccharomyces cerevisiae Network

Organisms use glucose to produce ethanol, glycerol, glycogen, and trehalose in the yeast fermentation pathway. For the map illustrated in Figure8, it has an important point of regulation that external glucose is transported into the cell. Glucose-6-phosphate (G6P) inhibits the step. Internal glucose is phosphorylated to G6P by hexokinase and G6P flow into distinct pathways. Here, one only considers two terms. The first one is toward ethanol, and the next one is toward the polysaccharides glycogen and trehalose. This is insignificant for the flow into the oxidative and nonoxidative pentose pathways. The process is omitted under the experimental conditions of interest. The omission allows less kinetic information to establish the model. Thus, it is a major modeling step. However, it exhibits a drawback, because it has conditions to limit model predictions that the pentose branches can be ignored in comparison with the ethanol branch. In mathematical complexity, the gained reduction is weighed with the limitation in scope.Figure 8
Map of anaerobic fermentation pathway inSaccharomyces cerevisiae.G6P and fructose-6-phosphate can be converted into each other. Thus, G6P and fructose-6-phosphate are considered to be at equilibrium. In the next step, fructose-6-phosphate is phosphorylated by phosphofructokinase. If the enzyme fructose-1,6-diphosphate is active, it will lead to a futile cycle between fructose-6-phosphate and fructose-1, 6-diphosphate. However, this enzyme is typically inactive under the experimental conditions of interest. And the model cannot include the cycle. Glycerol or phosphoenolpyruvate is produced by fructose-1,6-diphosphate. Finally, the production of ethanol is catalyzed by pyruvate kinase. In addition to these major reaction pathways, there are a number of reactions between the ADP and ATP procedures. The model can be mathematically expressed as follows [3]:
(26)x˙1=0.9018x2-0.2344x6-3.1833x10.7464x50.0243x7,x˙2=3.1833x10.7464x50.0243x7-0.8187x20.664x5-0.354x80.899×x110.0008,x˙3=0.5234x20.7318x5-0.3941x8-0.0148x30.584x40.03x50.119x90.944×x120.056x14-0.575,x˙4=0.022x30.6159x50.1308x9x14-0.6088-0.0945x30.05x40.533x5-0.0822×x10,x˙5=0.088x30.3505x40.2665x50.0243x90.5x100.5x14-0.3044-3.327x10.2134x20.19x50.362x70.286x80.257x110.0002x130.457.In Figure9, some important messages from the graphs were discerned. The reactions of PEP were eliminated from the simulated result. However, the result obtained is not satisfactory as expected. To explore the cause of failure, we conduct a comparison among examples to identify the differences. One can use the value of average residual R- to determine suitability of the model simplification. For this case, R- is 2.8417. Obviously, it is too large to endorse suitability of the method for model simplification.Graph of (a) reactive weight, (b) sensitivity analysis, (c) flux rates, and (d) compromise solution data in Case 3 (γ1=0.4,γ2=0.3,γ3=0.3).
(a)
(b)
(c)
(d)


## 3.4. Case 4: Dictyostelium Network

Figure10 illustrates pulses of cAMP produced by ACA. ACA is activated by the binding of extracellular cAMP to the surface receptor CAR1. The protein kinase ERK2 is activated to transmit the signal to ACA via ligand-bound CAR1. The activation of ACA is independent of ERK2 in an alternate circuit. If cAMP is accumulated internally, then the protein kinase PKA is activated via binding to the regulatory subunit of PKA. PKA inhibits ERK2 and the cAMP phosphodiesterase REG is activated by phosphorylating it. REG A is activated by a protein phosphatase. The process makes internal cAMP hydrolyze via REG A. When REGA hydrolyzes the internal cAMP, its regulatory subunit inactivates PKA activity. And CAR1 is returned to its high-affinity state by protein phosphatase. The phosphodiesterase PDE degrades secreted cAMP, before it diffuses between cells. The double horizontal line expresses the membrane surface of the cell [16, 17].Figure 10
Dictyostelium network.The mathematical model is given by(27)x˙1=2x7-0.9x1x2,x˙2=2.5x5-1.5x2,x˙3=0.6x7-0.8x2x3,x˙4=1-1.3x3x4,x˙5=0.3x1-0.8x4x5,x˙6=0.7x1-4.9x6,x˙7=23x6-4.5x7.Figure11 shows the results for the respective analyses of this case in which the average residual error value is 0.9488. For Cases 1–4, the values of R- are 7.16×10-2, 1.041×10-5, 2.8417, and 0.9488, respectively. It is quite obvious that the value of R- in Cases 3 and 4 is large than that in other cases. Thus, the model cannot be reduced further. It is quite obvious that not every biochemical network can be simplified in its structural model. However, one can quantify the mathematical model of the metabolic network while comparing relative importance of each reaction pathway.Graph of (a) reactive weight, (b) sensitivity analysis, (c) flux rates, and (d) compromise solution data in Case 4 (γ1=0.4,γ2=0.3,γ3=0.3).
(a)
(b)
(c)
(d)


## 4. Conclusions

This paper proposes a multiobjective analysis approach based on the SVD to identify the major characteristics in biochemical networks and performs model reduction from its original model. We use theS-system model to capture the dynamic behaviours of the biochemical networks and integrate the criteria of reactive weight, sensitivity, and flux analyses into the multiobjective function to reduce the unimportant states and related parameters. It is demonstrated that the simplified model is achievable and explainable with information obtained via the systematic approach. Several representative examples have been tested and verified which endorse applicability of the presented approach. Despite preliminary success of the demonstrative cases, most biochemical reactions in real world may not be that simple as those presented here. However, the proposed approach does provide an option to simplify the task of model reduction for complex biochemical networks.



---

*Source: 263973-2013-05-08.xml*