An Approach for Model Reduction of Biochemical Networks

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.

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. The S-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 typical S-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 𝜐⁺_𝑖∈ℜ and efflux 𝜐⁻_𝑖∈ℜ represented, respectively, by 𝑥̇_𝑖=𝜐⁺_𝑖−𝜐⁻_𝑖=𝛼_𝑖^𝑟+𝑚_𝑗=1𝑥^𝑔𝑖𝑗_𝑗−𝛽_𝑖^𝑟+𝑚_𝑗=1𝑥^ℎ𝑖𝑗_𝑗,𝑖=1,…,𝑟, (1) where 𝑥_𝑗,𝑗=1,…,𝑟 are dependent variables such as intermediate metabolites and products, 𝑥_𝑗,𝑗=𝑟+1,…,𝑟+𝑚 are independent variables such as substrates and enzymes, 𝑥_𝑖 is metabolic concentration, the positive numbers 𝛼_𝑖 and 𝛽_𝑖 refer to rate constants, and 𝑔_𝑖𝑗 and ℎ_𝑖𝑗 are kinetic order variables. Positive kinetic orders indicate activating influences while negative values mean inhibition. If the kinetic order 𝑔_𝑖𝑗  (ℎ_𝑖𝑗) equals zero, then the network will not contain the pathway from 𝑥_𝑗 to 𝑥_𝑖. This model can then be casted into a canonical matrix form regardless of the degree of its nonlinearity: 𝑑𝑥(𝑡)(𝑥(𝑡,𝑑𝑡=𝑆𝑓)) (2) where 𝑥(𝑡)∈ℜ^𝑟 is the concentration of reactants, 𝑆∈ℜ^𝑟×2𝑟 denotes the stoichiometric matrix, and 𝑓(⋅)∈ℜ^2𝑟 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.

If rank(𝑆)=𝑟, then 𝑆 can be decomposed by SVD as 𝑆=𝑈Σ𝑉^𝑇, (3) where 𝑈∈ℜ^𝑟×𝑟 and 𝑉∈ℜ^2𝑟×2𝑟 are unimodular matrices, Σ=[diag(𝜎₁,𝜎₂,…,𝜎_𝑟)0_𝑟×𝑟]∈ℜ^𝑟×2𝑟 is the matrix consisting of singular values (𝜎_𝑖) and null entries, the elements of 𝑈 are stoichiometric coefficients, and the elements of 𝑉 are reaction participation numbers in model.

To further investigate the effect of the flux rate on the biochemical network, the matrix 𝑆^𝑇𝑆 is diagonalized as 𝑆^𝑇𝑆=𝑉Σ^𝑇Σ𝑉^𝑇=𝑉Λ𝑉^𝑇, (4) where 𝑉=[𝑣_𝑗,𝑝]∈ℜ^2𝑟×2𝑟 and Λ=diag(𝜆₁,𝜆₂,…,𝜆_𝑟,0,…,0)∈ℜ^2𝑟×2𝑟 with 𝜆_𝑖=𝜎²_𝑖 are the eigenvector matrix and eigenvalue matrix, respectively. From the representation of (4), it is clear that 𝑉 involves reaction participation numbers and Λ consists of the square of singular values. Relative importance of the 𝑗th reaction can thus be measured via the weighted parameter 𝜂_𝑗 defined as 𝜂_𝑗=∑^𝑟_𝑝=1𝜆_𝑝𝑣_𝑗,𝑝∑^𝑟_𝑝=1𝜆_𝑝,𝑗=1,2,…,𝑟, (5) where 𝜆_𝑝 and 𝑣_𝑗,𝑝 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 𝛼_𝑖^𝑟+𝑚_𝑗=1𝑥^𝑔𝑖𝑗_𝑗=𝛽_𝑖^𝑟+𝑚_𝑗=1𝑥^ℎ𝑖𝑗_𝑗,𝑖=1,2,…,𝑟 (6) or ^𝑟+𝑚_𝑗=1𝑔_𝑖𝑗−ℎ_𝑖𝑗ln𝑥_𝑗=ln𝛽_𝑖−ln𝛼_𝑖,𝑖=1,2,…,𝑟. (7)

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 of S-system and can be employed in the evaluation of typical biochemical experiments. Consider the sensitivity of rate constants as 𝑧^𝛼_𝑖𝑗=𝜕ln𝑥_𝑖𝜕ln𝛼_𝑗,𝑧^𝛽_𝑖𝑗=𝜕ln𝑥_𝑖𝜕ln𝛽_𝑗, (8)(9) where 𝑧^𝛼_𝑖𝑗=−𝑧^𝛽_𝑖𝑗 in the steady state; 𝑖,𝑗=1,2,…,𝑟 by [3]. Define 𝑍=𝑧^𝛼_{𝑖𝑗𝑟×𝑟}=𝜕ln𝑥_𝑖𝜕ln𝛼_{𝑗𝑟×𝑟} (10) which can be decomposed as 𝑍^𝑇𝑍=𝐿𝐷𝐿^𝑇, (11) where 𝐿=[𝑙_𝑞,𝑙]∈ℜ^𝑟×𝑟 is the similarity transformation matrix and 𝐷=diag(𝑑₁,𝑑₂,…,𝑑_𝑟).

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 𝑗th reaction as 𝑤^𝛼_𝑞=𝑤^𝛽_𝑞=∑^𝑟_𝑙=1𝑑_𝑙𝑙_𝑞,𝑙∑^𝑟_𝑙=1𝑑_𝑙,𝑞=1,2,…,𝑟 (12) and let 𝑤_𝑗=𝑤^𝛼_𝑞,𝑞=𝑗+12𝑤,𝑗=odd_𝛽𝑞𝑗,𝑞=2,𝑗=even, (13) where 𝑤^𝛽_𝑞∈ℜ equals to 𝑤^𝛼_𝑞∈ℜ by 𝑍^𝑇𝑍 in (12). The factor 𝑤_𝑗 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 𝜐_𝑗. 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 𝑓_𝑗∈ℜ, as𝑓_𝑗=1𝑡_𝑓−𝑡₀^𝑡𝑓_𝑡0𝜐_𝑗𝑑𝑡,𝑗=1,2,…,2𝑟, (14) where 𝑡_𝑓 and 𝑡₀ are, respectively, the terminal and initial times for simulation; 𝜐_𝑗 represents the input flux of the ((𝑗+1)/2)th state when 𝑗 is odd; otherwise, it represents the output flux of the (𝑗/2)th state when 𝑗 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 𝜂_𝑗, consider 𝜂_min≤𝜂_𝑗≤𝜂_max with 0≤(𝜂_𝑗−𝜂_min)/(𝜂_max−𝜂_min)≤1, where 𝜂_min and 𝜂_max are, respectively, the lower and upper limits of 𝜂_𝑗. Similar treatments are applied for 𝑤_𝑗 and 𝑓_𝑗. Let 𝛾𝑌=¹𝜂_𝑗−𝜂_min𝜂_max−𝜂_min𝛾₂𝑤_𝑗−𝑤_min𝑤_max−𝑤_min𝛾₃𝑓_𝑗−𝑓_min𝑓_max−𝑓_min3×1, (15) where 𝛾_𝑖 is the weight of the 𝑖th subobjective; 𝑤_𝑗 and 𝑓_𝑗 are defined, respectively, in (13) and (14). Define 𝑌_𝑗=𝛾₁𝜂_𝑗−𝜂_min𝜂_max−𝜂_min^𝑠+𝛾₂𝑤_𝑗−𝑤_min𝑤_max−𝑤_min^𝑠+𝛾₃𝑓_𝑗−𝑓_min𝑓_max−𝑓_min^𝑠,^1/𝑠 (16) where 𝑗=1,2,…,2𝑟. The exponent 𝑠∈[1,∞) is used to reflect the general vector length. If 𝑠=1, the weighting represents the relative importance among each subobjective. Next, consider combination of the three objective functions by finding the cost value 𝑌_𝑗∈ℜ. It can help to distinguish important and unimportant parts: 𝑦=𝜀max_𝑗𝑌_𝑗, (17) where 𝑦 serves as a threshold and 𝜀 is a normalized factor to normalize the value of 𝑦. If 𝑌_𝑗<𝑦, 𝑌_𝑗 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 𝑔_𝑖𝑗 and ℎ_𝑖𝑗 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 error 𝑅_𝑖∈ℜ via 𝑅_𝑖=∫^𝑡𝑓_𝑡0𝑥_𝑖(𝑡)−𝑥^𝑟_𝑖(𝑡)²𝑑𝑡∫^𝑡𝑓_𝑡0𝑥²_𝑖(𝑡),𝑑𝑡 (18) where 𝑥^𝑟_𝑖(𝑡) denotes concentration of the reduced network at the sampling time instant 𝑡. If the error concentration is larger than the original system concentration, then the reduced system is unfeasible and could not be reduced. Furthermore, define 1𝑅=𝑟^𝑟_𝑖=1𝑅_𝑖. (19) When 𝑅→0⁺, it refers to the ideal situation and there is no residual error. On the other hand, when 𝑅→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 Figure 1, 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 CO₂. 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 𝑥₁₄ through 𝑥₄₈. The dynamic equations can be expressed as [15] 𝑥̇₁=0.8282𝑥₁^−0.038𝑥₆^−0.0204𝑥₇^0.106𝑥^0.114₁₀𝑥^0.7₁₃𝑥^0.807₁₈×𝑥^0.108₂₈𝑥^0.0848₃₁𝑥^0.599₄₆𝑥^−0.181₄₈−1.3423𝑥₁𝑥^0.915₃₀𝑥^0.0847₃₃,𝑥̇₂=1.3401𝑥₁^0.915𝑥₇^0.0848𝑥^0.0848₂₃𝑥^0.915₃₀−17.166𝑥₂^0.706𝑥₃^0.0716𝑥^0.0848₂₂𝑥^0.915₂₄𝑥^−0.0341₄₇,𝑥̇₃=0.3231𝑥₃^−0.405𝑥₅^0.156𝑥^0.427₂₁𝑥^0.573₃₅𝑥^0.422₄₆𝑥^0.405₄₇𝑥^−0.418₄₈−9.6952𝑥₂^0.376𝑥₃^0.489𝑥^0.554₂₄𝑥^0.446₄₁𝑥^−0.00206₄₇,𝑥̇₄=𝑥₉𝑥₂₅−0.152𝑥₄^0.958𝑥₂₆𝑥^0.0348₄₆𝑥^−0.862₄₈,𝑥̇₅=1.875𝑥₇^0.0274𝑥₈^0.465𝑥^0.336₁₃𝑥^0.535₁₉𝑥^0.465₂₀−0.01923𝑥₃^−0.717𝑥₅^0.413𝑥₆^0.306𝑥₈^−0.29𝑥^−0.0883₁₀×𝑥^0.756₂₁𝑥^0.244₂₉𝑥^0.748₄₆𝑥^0.718₄₇𝑥^−0.741₄₈,𝑥̇₆=2.459𝑥₁^−0.00921𝑥₆^−0.0154𝑥₇^0.0162𝑥^0.086₁₀𝑥^0.276₁₁×𝑥^0.813₂₈𝑥^0.276₃₂𝑥^0.6413₃₉−1.1528𝑥₅^0.0963𝑥₆^1.01𝑥₈^−0.204𝑥^−0.062₁₀×𝑥^0.0518₁₄𝑥^0.277₂₇𝑥^0.171₂₉𝑥^0.5₄₅𝑥^0.0222₄₆𝑥^−0.0191₄₈,𝑥̇₇=2.1167𝑥₁^0.129𝑥₂^0.129𝑥^0.129₂₂𝑥^0.129₃₃𝑥^0.741₃₄−3.4893𝑥₁^−0.0187𝑥₆^−0.0311𝑥₇^0.868𝑥^0.174₁₀𝑥^0.129₂₃×𝑥^0.165₂₈𝑥^0.129₃₁𝑥^0.577₄₀,𝑥̇₈=0.5724𝑥₅^0.111𝑥₆^0.247𝑥₈^0.234𝑥^0.0713₁₀𝑥^0.197₂₉𝑥^0.803₃₈−1.9369𝑥₈𝑥^0.375₂₀𝑥^0.625₄₄,𝑥̇₉=16.242𝑥₂^0.679𝑥₃^0.0782𝑥₂₄𝑥^−0.0372₄₇−𝑥₉𝑥₂₅,𝑥̇₁₀=0.156𝑥₄^0.724𝑥₅^0.106𝑥₆^0.259𝑥₈^−0.223𝑥^−0.0679₁₀𝑥^0.0568₁₄×𝑥^0.756₂₆𝑥^0.188₂₉𝑥^0.0506₄₆𝑥^−0.672₄₈−0.8063𝑥₁^−0.0101𝑥₆^−0.0168×𝑥₇^0.0177𝑥^0.99₁₀𝑥^−0.879₁₁𝑥^0.911₁₅𝑥^0.0891₂₈𝑥^0.882₄₆𝑥^0.879₄₇𝑥^−0.881₄₈,𝑥̇₁₁=2.0031𝑥₆^0.166𝑥^0.491₁₀𝑥^−0.481₁₁𝑥^0.499₁₅𝑥^0.166₂₇×𝑥^0.335₃₆𝑥^0.483₄₆𝑥^0.481₄₇𝑥^−0.483₄₈−2.4373𝑥^0.495₁₁𝑥^−0.00542₁₂𝑥^0.574₁₆𝑥^0.166₃₂𝑥^0.261₄₂,𝑥̇₁₂=1.271𝑥^0.106₁₁𝑥^−0.00836₁₂𝑥^0.885₁₆𝑥^0.115₃₇−9.1694𝑥^1.89₁₂𝑥^−1.24₁₃𝑥^0.911₁₇𝑥^0.0893₄₃,𝑥̇₁₃=8.289𝑥^1.98₁₂𝑥^−1.36₁₃𝑥₁₇−0.9387𝑥₁^−0.0197𝑥₇^0.0196×𝑥^0.775₁₃𝑥^0.618₁₈𝑥^0.382₁₉𝑥^0.458₄₆𝑥^−0.139₄₈. (20)

Figure 1 

Tricarboxylic acid cycle network in Dictyostelium 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 𝐷 are given below: (),(),Λ=diag296.47,264.8,94.1,85.69,69.59,16.66,9.95,7.38,4.08,3.52,2.49,1.02,0.67,0,…,0𝐷=diag0.07,0.27,0.34,0.43,0.47,0.77,1,1.38,2.2,3,8.73,67.8,144.34 (21) from which we obtain the maximal value of 𝑌_𝑗 to be 0.0534, with which, 19 pathways and one metabolite could be eliminated based on the proposed method while setting 𝑦=0.0534 in Tables 1 and 2 or the dashed part of Figure 1. The reduced dynamic model is given by 𝑥̇₁=0.8282𝑥₁^−0.038𝑥₆^−0.0204𝑥^0.114₁₀𝑥^0.7₁₃𝑥^0.807₁₈𝑥^0.108₂₈×𝑥^0.0848₃₁𝑥^0.599₄₆𝑥^−0.181₄₈−1.3423𝑥₁𝑥^0.915₃₀𝑥^0.0847₃₃,𝑥̇₂=1.3401𝑥₁^0.915𝑥^0.0848₂₃𝑥^0.915₃₀−17.166𝑥₂^0.706𝑥₃^0.0716𝑥^0.0848₂₂𝑥^0.915₂₄𝑥^−0.0341₄₇,𝑥̇₃=0.3231𝑥₃^−0.405𝑥₅^0.156𝑥^0.427₂₁𝑥^0.573₃₅𝑥^0.422₄₆𝑥^0.405₄₇𝑥^−0.418₄₈−9.6952𝑥₂^0.376𝑥₃^0.489𝑥^0.554₂₄𝑥^0.446₄₁𝑥^−0.00206₄₇,𝑥̇₄=𝑥₉𝑥₂₅−0.152𝑥₄^0.958𝑥₂₆𝑥^0.0348₄₆𝑥^−0.862₄₈,𝑥̇₅=1.875𝑥₈^0.465𝑥^0.336₁₃𝑥^0.535₁₉𝑥^0.465₂₀−0.01923𝑥₃^−0.717𝑥₅^0.413𝑥₆^0.306𝑥₈^−0.29𝑥^−0.0883₁₀×𝑥^0.756₂₁𝑥^0.244₂₉𝑥^0.748₄₆𝑥^0.718₄₇𝑥^−0.741₄₈,𝑥̇₆=2.459𝑥₁^−0.00921𝑥₆^−0.0154𝑥^0.086₁₀𝑥^0.276₁₁×𝑥^0.813₂₈𝑥^0.276₃₂𝑥^0.6413₃₉−1.1528𝑥₅^0.0963𝑥₆^1.01𝑥₈^−0.204×𝑥^−0.062₁₀𝑥^0.0518₁₄𝑥^0.277₂₇𝑥^0.171₂₉𝑥^0.5₄₅𝑥^0.0222₄₆𝑥^−0.0191₄₈,𝑥̇₈=0.5724𝑥₅^0.111𝑥₆^0.247𝑥₈^0.234𝑥^0.0713₁₀𝑥^0.197₂₉𝑥^0.803₃₈−1.9369𝑥₈𝑥^0.375₂₀𝑥^0.625₄₄,𝑥̇₉=16.242𝑥₂^0.679𝑥₃^0.0782𝑥₂₄𝑥^−0.0372₄₇−𝑥₉𝑥₂₅,𝑥̇₁₀=0.156𝑥₄^0.724𝑥₅^0.106𝑥₆^0.259𝑥₈^−0.223𝑥^−0.0679₁₀𝑥^0.0568₁₄×𝑥^0.756₂₆𝑥^0.188₂₉𝑥^0.0506₄₆𝑥^−0.672₄₈−0.8063𝑥₁^−0.0101𝑥₆^−0.0168×𝑥^0.99₁₀𝑥^−0.879₁₁𝑥^0.911₁₅𝑥^0.0891₂₈𝑥^0.882₄₆𝑥^0.879₄₇𝑥^−0.881₄₈,𝑥̇₁₁=2.0031𝑥₆^0.166𝑥^0.491₁₀𝑥^−0.481₁₁𝑥^0.499₁₅×𝑥^0.166₂₇𝑥^0.335₃₆𝑥^0.483₄₆𝑥^0.481₄₇𝑥^−0.483₄₈−2.4373𝑥^0.495₁₁𝑥^−0.00542₁₂𝑥^0.574₁₆𝑥^0.166₃₂𝑥^0.261₄₂,𝑥̇₁₂=1.271𝑥^0.106₁₁𝑥^−0.00836₁₂𝑥^0.885₁₆𝑥^0.115₃₇−9.1694𝑥^1.89₁₂𝑥^−1.24₁₃𝑥^0.911₁₇𝑥^0.0893₄₃,𝑥̇₁₃=8.289𝑥^1.98₁₂𝑥^−1.36₁₃𝑥₁₇−0.9387𝑥₁^−0.0197×𝑥^0.775₁₃𝑥^0.618₁₈𝑥^0.382₁₉𝑥^0.458₄₆𝑥^−0.139₄₈. (22)

Table 1 

Result of model reduction by 𝜀 in Case  1 (IMP as important part, UMP as unimportant part).

Table 2 

Result of model reduction in Case  1 (𝜀=0.1).

(a)

(b)

(c)

(d)

Figure 2 

Graph 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).

Note that 𝑥₇ 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%.

(a)

(b)

Figure 3 

Dynamic responses of the noisy states in Case 1: (a) original model and (b) simplified model (𝑥𝑖 as original concentration, 𝑥𝑟𝑖 as reduced concentration).

3.2. Case 2: Purine Metabolism

In Figure 4, 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 𝑥₄ and 𝑥₈ as substitutes and adjust the kinetic orders accordingly, because the latter are not represented individually in the model [3]: 𝑥̇₁=0.898𝑥₁^−0.03𝑥₄^−0.45𝑥₈^−0.04𝑥^0.65₁₇𝑥^0.7₁₈−31.335𝑥₁^1.27𝑥₂^−0.165×𝑥₄^−0.067𝑥₆^0.036𝑥₈^−0.237𝑥^0.085₁₃𝑥^0.075₁₅𝑥^−0.0092₁₈,𝑥̇₂=1.629𝑥₁^0.72𝑥₂^−0.28𝑥₄^0.319𝑥₇^−0.03𝑥₈^−0.024𝑥^0.144₁₃𝑥^−0.06₁₈−5.584𝑥₂^0.454𝑥₄^−0.156𝑥₇^−0.012𝑥₈^0.126𝑥^−0.11₁₈,𝑥̇₃=3.5932𝑥₂^0.4𝑥₄^−0.24𝑥₈^0.2𝑥^−0.05₁₈−66544.7𝑥₃^0.99𝑥₄^−0.95,𝑥̇₄=0.082082𝑥₁^0.00025𝑥₃^0.004𝑥₄^−0.00418𝑥₅^0.0023𝑥₆^0.00037𝑥^0.988₁₁−614.27𝑥₄^0.054𝑥₅^−0.0045𝑥₈^0.128𝑥₉^−0.00003𝑥^0.00087₁₀𝑥^−0.00028₁₈,𝑥̇₅=7.2067𝑥₄^0.2𝑥₅^−0.6−9.006𝑥₅^0.368,𝑥̇₆=0.29𝑥₅^0.9−220.14𝑥₁^0.495𝑥₄^−0.79𝑥₆^0.748,𝑥̇₇=1.2823𝑥₂^0.15𝑥₇^−0.09𝑥₈^−0.03−0.3738𝑥₄^0.12𝑥₇^0.16,𝑥̇₈=0.048842𝑥₁^0.0034𝑥₄^0.00014𝑥₇^0.00019𝑥₈^−0.0034𝑥^0.996₁₁𝑥^0.001₁₅−408.11𝑥₂^−0.000057𝑥₄^0.0498𝑥₇^−0.00029𝑥₈^0.133𝑥₉^−0.00009×𝑥^−0.000029₁₀𝑥^−0.001₁₈,𝑥̇₉=0.00229𝑥₄^0.002𝑥₉^−0.0058𝑥^0.017₁₀𝑥^0.98₁₂−3.3𝑥₉^0.431𝑥^0.324₁₀,𝑥̇₁₀=0.00151𝑥₈^0.0058𝑥₉^−0.017𝑥^−0.0058₁₀𝑥^0.986₁₂−2.2655𝑥₉^0.413𝑥^0.34₁₀,𝑥̇₁₁=1024.195𝑥₄^0.05𝑥₈^0.13−0.115389𝑥₁₁,𝑥̇₁₂=5.5062𝑥₉^0.42𝑥^0.33₁₀−0.0032545𝑥₁₂,𝑥̇₁₃=0.1033𝑥₂^0.43𝑥₄^0.409𝑥₉^0.04𝑥^−0.19₁₈−8.2827𝑥₁^0.817𝑥₂^−0.66𝑥^0.528₁₃,𝑥̇₁₄=0.7234𝑥^0.343₁₃𝑥^0.236₁₅−0.9325𝑥^0.569₁₄,𝑥̇₁₅=0.2615𝑥₈^0.881𝑥^0.02₁₀𝑥^−0.33₁₈−136.6𝑥₁^0.925𝑥₈^−0.925𝑥^0.438₁₅,𝑥̇₁₆=0.949𝑥^0.55₁₄−0.00008744𝑥^2.21₁₆. (23)

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 𝐷 are (),(),Λ=diag296.47,264.8,94.1,85.69,69.59,16.66,9.95,7.38,4.08,3.52,2.49,1.02,0.67,0,…,0𝐷=diag0.07,0.27,0.34,0.43,0.47,0.77,1,1.38,2.2,3,8.73,67.8,144.34 (24)

(a)

(b)

(c)

(d)

Figure 5 

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).

Figure 6 

Dynamic responses of the noisy states in Case 2 where 𝑥𝑖 denotes the original concentration and 𝑥𝑟𝑖 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 Tables 3 and 4 that 10 pathways and 3 metabolites are ignorable when choosing 𝑦=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 𝑥₆, 𝑥₁₄, and 𝑥₁₆. The six dashed lines in Figure 4 show the pathways to be ignorable after performing the analysis. The corresponding reduced model is given by 𝑥̇₁=0.898𝑥₁^−0.03𝑥₄^−0.45𝑥₈^−0.04𝑥^0.65₁₇𝑥^0.7₁₈−31.335𝑥₁^1.27𝑥₂^−0.165𝑥₄^−0.067𝑥₆^0.036𝑥₈^−0.237𝑥^0.085₁₃×𝑥^0.075₁₅𝑥^−0.0092₁₈,𝑥̇₂=1.629𝑥₁^0.72𝑥₂^−0.28𝑥₄^0.319𝑥₇^−0.03𝑥₈^−0.024𝑥^0.144₁₃𝑥^−0.06₁₈−5.584𝑥₂^0.454𝑥₄^−0.156𝑥₇^−0.012𝑥₈^0.126𝑥^−0.11₁₈,𝑥̇₃=3.5932𝑥₂^0.4𝑥₄^−0.24𝑥₈^0.2𝑥^−0.05₁₈−66544.7𝑥₃^0.99𝑥₄^−0.95,𝑥̇₄=0.082082𝑥₁^0.00025𝑥₃^0.004𝑥₄^−0.00418𝑥₅^0.0023𝑥₆^0.00037𝑥^0.988₁₁−614.27𝑥₄^0.054𝑥₅^−0.0045𝑥₈^0.128𝑥₉^−0.00003𝑥^0.00087₁₀𝑥^−0.00028₁₈,𝑥̇₅=7.2067𝑥₄^0.2𝑥₅^−0.6−9.006𝑥₅^0.368,𝑥̇₇=1.2823𝑥₂^0.15𝑥₇^−0.09𝑥₈^−0.03−0.3738𝑥₄^0.12𝑥₇^0.16,𝑥̇₈=0.048842𝑥₁^0.0034𝑥₄^0.00014𝑥₇^0.00019𝑥₈^−0.0034𝑥^0.996₁₁𝑥^0.001₁₅−408.11𝑥₂^−0.000057𝑥₄^0.0498𝑥₇^−0.00029𝑥₈^0.133𝑥₉^−0.00009×𝑥^−0.000029₁₀𝑥^−0.001₁₈,𝑥̇₉=0.00229𝑥₄^0.002𝑥₉^−0.0058𝑥^0.017₁₀𝑥^0.98₁₂−3.3𝑥₉^0.431𝑥^0.324₁₀,𝑥̇₁₀=0.00151𝑥₈^0.0058𝑥₉^−0.017𝑥^−0.0058₁₀𝑥^0.986₁₂−2.2655𝑥₉^0.413𝑥^0.34₁₀,𝑥̇₁₁=1024.195𝑥₄^0.05𝑥₈^0.13−0.115389𝑥₁₁,𝑥̇₁₂=5.5062𝑥₉^0.42𝑥^0.33₁₀−0.0032545𝑥₁₂,𝑥̇₁₃=0.1033𝑥₂^0.43𝑥₄^0.409𝑥₉^0.04𝑥^−0.19₁₈−8.2827𝑥₁^0.817𝑥₂^−0.66𝑥^0.528₁₃,𝑥̇₁₅=0.2615𝑥₈^0.881𝑥^0.02₁₀𝑥^−0.33₁₈−136.6𝑥₁^0.925𝑥₈^−0.925𝑥^0.438₁₅. (25)

Table 3 

Result of model reduction by 𝜀 in Case  2 (IMP as important part, UMP as unimportant part).

Table 4 

Result of model reduction in Case  2 (𝜀=0.0005).

The 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 Figure 7 the residual errors between the original and reduced models for the two cases above. Again, the result confirms the proposed approach.

(a) Case 1

(b) Case 2

Figure 7 

Residual errors of model reduction for Cases 1 and 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 Figure 8, 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 in Saccharomyces 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]: 𝑥̇₁=0.9018𝑥₂^−0.2344𝑥₆−3.1833𝑥₁^0.7464𝑥₅^0.0243𝑥₇,𝑥̇₂=3.1833𝑥₁^0.7464𝑥₅^0.0243𝑥₇−0.8187𝑥₂^0.664𝑥₅^−0.354𝑥₈^0.899×𝑥^0.0008₁₁,𝑥̇₃=0.5234𝑥₂^0.7318𝑥₅^−0.3941𝑥₈−0.0148𝑥₃^0.584𝑥₄^0.03𝑥₅^0.119𝑥₉^0.944×𝑥^0.056₁₂𝑥^−0.575₁₄,𝑥̇₄=0.022𝑥₃^0.6159𝑥₅^0.1308𝑥₉𝑥^−0.6088₁₄−0.0945𝑥₃^0.05𝑥₄^0.533𝑥₅^−0.0822×𝑥₁₀,𝑥̇₅=0.088𝑥₃^0.3505𝑥₄^0.2665𝑥₅^0.0243𝑥₉^0.5𝑥^0.5₁₀𝑥^−0.3044₁₄−3.327𝑥₁^0.2134𝑥₂^0.19𝑥₅^0.362𝑥₇^0.286𝑥₈^0.257𝑥^0.0002₁₁𝑥^0.457₁₃. (26)

In Figure 9, 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 𝑅 to determine suitability of the model simplification. For this case, 𝑅 is 2.8417. Obviously, it is too large to endorse suitability of the method for model simplification.

(a)

(b)

(c)

(d)

Figure 9 

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).

3.4. Case 4: Dictyostelium Network

Figure 10 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 𝑥̇₁=2𝑥₇−0.9𝑥₁𝑥₂,𝑥̇₂=2.5𝑥₅−1.5𝑥₂,𝑥̇₃=0.6𝑥₇−0.8𝑥₂𝑥₃,𝑥̇₄=1−1.3𝑥₃𝑥₄,𝑥̇₅=0.3𝑥₁−0.8𝑥₄𝑥₅,𝑥̇₆=0.7𝑥₁−4.9𝑥₆,𝑥̇₇=23𝑥₆−4.5𝑥₇. (27)

Figure 11 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 𝑅 are 7.16×10⁻², 1.041×10⁻⁵, 2.8417, and 0.9488, respectively. It is quite obvious that the value of 𝑅 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.

(a)

(b)

(c)

(d)

Figure 11 

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).

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 the S-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.

Acknowledgment

This research was sponsored by the National Science Council, Taiwan, under Grant NSC-98-2221-E-005-087-MY3.