Complexity analysis and parameter estimation of dynamic metabolic systems.

A metabolic system consists of a number of reactions transforming molecules of one kind into another to provide the energy that living cells need. Based on the biochemical reaction principles, dynamic metabolic systems can be modeled by a group of coupled differential equations which consists of parameters, states (concentration of molecules involved), and reaction rates. Reaction rates are typically either polynomials or rational functions in states and constant parameters. As a result, dynamic metabolic systems are a group of differential equations nonlinear and coupled in both parameters and states. Therefore, it is challenging to estimate parameters in complex dynamic metabolic systems. In this paper, we propose a method to analyze the complexity of dynamic metabolic systems for parameter estimation. As a result, the estimation of parameters in dynamic metabolic systems is reduced to the estimation of parameters in a group of decoupled rational functions plus polynomials (which we call improper rational functions) or in polynomials. Furthermore, by taking its special structure of improper rational functions, we develop an efficient algorithm to estimate parameters in improper rational functions. The proposed method is applied to the estimation of parameters in a dynamic metabolic system. The simulation results show the superior performance of the proposed method.

1. Introduction

Living cells require energy and material for maintaining their essential biological processes through metabolism, which is a highly organized process. Metabolic systems are defined by the enzymes dynamically converting molecules of one type into molecules of another type in a reversible or irreversible manner. Modeling and parameter estimation in dynamic metabolic systems provide new approaches towards the analysis of experimental data and properties of the systems, ultimately leading to a great understanding of the language of living cells and organisms. Moreover, these approaches can also provide systematic strategies for key issues in medicine, pharmaceutical, and biotechnological industries [1]. The formulation and identification of metabolic systems generally includes the building of the mathematical model of biological process and the estimating of system parameters. Because the components of a pathway interact not only with each other in the same pathway but also with those in different pathways, most (if not all) of mathematical models of metabolic systems are highly complex and nonlinear. The widely used approaches for modeling inter- and intracellular dynamic processes are based on mass action law [1-4]. By mass action law, the reaction rates are generally polynomials in concentrations of metabolites with reaction constants or rational functions which are a fraction and whose denominator and numerators are polynomials in concentrations of metabolites with reaction constants [1-4]. As a result, the mathematical model is nonlinear not only in the states but also in the parameters. Estimation of these parameters is crucial to construct a whole metabolic system [5-7].

In general, all algorithms for nonlinear parameter estimation can be used to estimate parameters in metabolic systems, for example, Gauss-Newton iteration method, and its variants such as Box-Kanemasu interpolation method, Levenberg damped least squares methods and Marquardt's method [8, 9]. However, these iteration methods are initial-sensitive. Another main shortcoming is that these methods may converge to the local minimum of the least squares cost function and thus cannot find the real values of parameters. Furthermore, because of their highly complexity and nonlinearity, Gauss-Newton iteration method and its variants cannot efficiently and accurately estimate the parameters in metabolic systems [5–7, 10, 11].

In this paper, we propose a systematic method for estimating parameters in dynamic metabolic systems. Typically mathematical model of dynamic metabolic systems consists of a group of nonlinear differential equations, some of which contains several rational functions in which parameters are nonlinear. In Section 2, we propose a method for model complexity analysis via the stoichiometric matrix. As a result, we obtain a group of equations, each of which contains only one-rational function plus polynomial functions which we called an improper rational function. Then, based on the observation that in the improper rational functions both the denominator and numerator are linear in parameters while polynomials are also linear in parameters, we develop an iterative linear least squares method for estimating parameters in dynamic metabolic systems in Section 3. The basic idea is to transfer optimizing a nonlinear least squares objective function into iteratively solving a sequence of linear least squares problems. In Section 4, we apply our developed method to estimate parameters in a metabolism system. Finally we give conclusions and some directions of future work along with this study in Section 5.

2. Model Complexity Analysis for Parameter Estimation

A dynamic metabolic system consists of k substances (molecules), and m reactions can be described by a system of differential equations as follows: where x represents the concentrations of molecule i, r represents the reaction rate j, and c represents the stoichiometric coefficient of molecule i in reaction j. The mass action law in biochemical kinetics [2–4, 12] states that the reaction rate is proportional to the probability of a collision of the reactants. This probability is in turn proportional to the concentration of reactants. Therefore, reaction rate r is a function of the concentrations of molecules involved in reaction j and proportion constants.

The stoichiometric coefficient c assigned to molecule i and reaction j can be put into a so-called stoichiometric matrix C = [c ]. Let X = [x 1, x 2,…, x ] and r = [r 1, r 2,…, r ], and let β = [β 1,β 2,…,β ] represent the vector consisting of all independent proportion constants, and then (1) can be rewritten in the following vector-matrix format:

In principle, the stoichiometric coefficient c in matrix C is a constant integer and can be decided according to how molecule i is involved in reaction j. According to mass action law, the expression of reaction rates can be determined to be polynomials or rational functions with reaction constants [2–4, 12]. The challenge to build up the mathematic model of dynamic metabolic system (2) is to estimate the parameter vector β, especially when some reaction rates are in the form of rational functions in which parameters are nonlinear.

If each differential equation in (2) contains one-rational function without or with polynomial functions, the parameters in model (2) can be estimated by algorithms in [13, 14] or a new algorithm proposed in the next section of this paper. Unfortunately, each differential equation contains a linear combination of several rational functions, which makes the parameter estimation in those coupled differential equations more difficult. The stoichiometric matrix contains very important information about the structure of the metabolic systems and is widely used to analyze the steady state and flux balance of metabolic systems [2-4]. In this paper, via the stoichiometric matrix, we propose a systematic method to transfer a system of differential equations (2) into another system of differential equations, in which each differential equation contains at most one-rational function.

Running Example. To illustrate the proposed method, we use the upper part of glycolysis system as a running example, showing how the method is applied to this system step after step. The schematic representation of this system is shown in Figure 1. The model for this metabolic system is described by the system of differential equations (2) as follows:

Figure 1

Schematic representation of the upper part of glycolysis [4].

Based on the mass action law, the individual reaction rates can be expressed as

Model (3) has six ordinary differential equations (ODEs) and 15 parameters contained in eight reaction rates, three out of which are rational functions. Some ODEs contain more than one rational reaction rates, which makes the parameter more difficult.

Comparing (3) to (2) we have the state vector: X = [Gluc6P; Fruc6P; Fruc1,6P2; ATP, ADP, AMP] and stoichiometric matrix:

In the following, we describe our proposed method to analyze the complexity of model (2) through the running example.

Step 1

Collect the columns in the stoichiometric matrix corresponding to the rational reaction rates in model (2) to construct a submatrix C and collect other columns (corresponding to polynomial reaction rates) to construct a submatrix C . Therefore, we have

where r is the subvector of r and consists of all rational reaction rates while r is another subvector of r and consists of all polynomial reaction rates. In this step, we should make sure that the rank of matrix C equals the number of rational reaction rates. If the rank of matrix C does not equal the number of rational reaction rates, it means that some rational reaction rates are not independent. Then we combine dependent rational reaction rates together to create a new reaction rate such that all resulted rational reaction rates should be linearly independent [14]. As a result, the rank of matrix C will equal the number of rational reaction rates.

For the running example, we have

and r = [r 1, r 3, r 4] and r = [r 2, r 5, r 6, r 7, r 8]. The rank of matrix C equals 3, which is the number of rational reaction rates.

Step 2

Calculate the left inverse matrix of C . That is, calculate C − such that

As matrix C has the column full rank, matrix C − satisfying (8) exists although it is typically not unique. For a given matrix C , C − can be easily found by solving (8) which is a linear algebraic system. If it is not unique, any matrix satisfying (8) works for our proposed method.

For the running example, we can have

Step

Multiply (6) by matrix C − from the left to obtain or

From its expression, each differential equation in the system (11) contains only one-rational reaction rates plus a linear combination of polynomial reaction rates.

For the running example, we have

Step 4

Calculate matrix C ⊥ such that

where C ⊥ has the full row rank and rank⁡(C ⊥) + rank⁡(C −) = the  number  of  rows  in  C . Note that C ⊥ can be easily found by solving (13), which is a homogenous linear algebraic system. Again if it is not unique, any matrix satisfying (13) works for our proposed method.

Then multiply (6) by matrix C ⊥ from the left to obtain

or

For the running example, we can have

Step 5

Let D = C ⊥ C . If rank⁡(D) ≥ the  number  of  columns, then solving (15) yields

If rank⁡(D) < the  number  of  columns, it means that some polynomial reaction rates in (15) are linearly dependent. Then combine the linearly dependent rates and construct a new reaction rate vector and full column rank matrix such that

and then solving (18) yields

For the running example, we have rank⁡(D) < the  number  of  columns. As the first four columns are linearly dependent, we can have a new reaction rates −2r 2 − 2r 5 + r 6 − r 7. Therefore, we have and furthermore, noting that (d/dt)(ATP + ADP + AMP) = 0, from (19) we have

After these five steps, dynamic metabolic system (2) is transferred into a system of differential equations, in which each differential equation contains one-rational function plus polynomial functions ((11) or (12)) or only polynomial function ((19) or (21)). Parameters in (19) can be analytically estimated by well-known least squares methods. In the next section, we describe an algorithm to estimate parameters in (11).

3. Parameter Estimation Algorithm

After its complexity analysis, estimating parameters in dynamic metabolic system is reduced to mainly estimating parameters in a rational function plus polynomial, which we call the improper rational function. These functions are nonlinear in both parameters and state variables. Therefore, estimation of parameters in these models is a nonlinear estimation problem. In general, all algorithms for nonlinear parameter estimation can be used to estimate parameters in the improper rational functions, for example, Gauss-Newton iteration method and its variants such as Box-Kanemasu interpolation method, Levenberg damped least squares methods, Marquardt's method [9–12, 15], and more sophisticated methods [16]. However, these iteration methods are initial sensitive. Another main shortcoming is that most of these methods may converge to the local minimum of the least squares cost function and thus cannot find the real values of parameters. In the following, we describe an iterative linear least squares method to estimate parameters in the improper rational functions. The basic idea is to transfer optimizing a nonlinear least squares objective function into iteratively solving a sequence of linear least squares problems.

Consider the general form of the following improper rational functions: where the vector X consists of the state variables and the p-dimensional vector β consists of all parameters in the improper rational function (22), which can naturally be divided into three groups: those in the numerator of the rational functions β   (i = 1,…, p ), those in the denominator of the rational function β   (j = 1,…, p ), and those in the polynomial β   (k = 1,…, p ), where we have that p + p + p = p. N (X)  (i = 0,1,…, p ), D (X)  (j = 0,1,…, p ), and P (X)  (k = 1,…, p ) are the known functions nonlinear in the state variable X and do not contain any unknown parameters. Either N 0(X) or D 0(X) must be nonzero, and otherwise from sensitivity analysis [9, 16] the parameters in model (22) cannot be uniquely identified.

If there is no polynomial part, model (22) is reduced to a rational function. Recently, several methods have been proposed for estimating parameters in rational functions [5, 6, 13, 14]. The authors in [5, 6] have employed general nonlinear parameter estimation methods to estimate parameters in rational functions. As shown in their results, the estimation error is fairly large. We have observed that in rational functions both the denominator and numerator are linear in the parameters. Based on this observation, we have developed iterative linear least squares methods in [13, 14] for estimating parameters in rational functions. Mathematically, improper rational function (22) can be rewritten as the following rational function:

However, in the numerator of the model above, there are p p + p + p coefficients while there are p + p + p unknown parameters. When p = 1, the number of parameters is equal to the numbers of coefficients, and the methods developed in [13, 14] can be applied. However, when p > 1, those methods are not applicable as the number of parameters is less than the number of coefficients in the numerator.

In order to describe an algorithm to estimate parameters in the improper rational function (22) for n given groups of observation data y and X   (t = 1,2,…, n), we introduce the following notation:

To estimate parameters in the improper rational function (22), as in [11], we form a sum of the weighted squared errors (the cost function) with the notions above as follows:

Minimizing J(β) with respect to β = [β ,β ,β ] can give the nonlinear least squares estimation of parameters β , β , and β . We rewrite the objective function (22) as follows:

In the objective function (26), for a given parameters in the first term, we have where

Then for given parameters , we can estimate the parameters β = [β ,β ,β ] by linear least squares method as follows:

Based on the above discussion, we propose the following iterative linear least squares method.

Choose the initial guess for β 0.

Iteratively construct matrix A(β ) and vector b by (28) and (29), respectively, and then solve the linear least squares problem:

which gives the solution until the stopping criterion is met, where β = [β ,β ,β ] is the estimation of parameters β at step s.

From (31), if the estimation sequence β 1, β 2,… is converged to β*, the objective function (26) reaches its minimum value at β*. That is, β*is the estimation of parameters in model (22).

There are several ways to set up a stopping criterion. In this paper the stopping criteria are chosen as

where ||·|| is the Euclidean norm of the vector and ε is a preset small positive number, for example, 10−5.

4. Application

To investigate the method developed in previous sections, this study generates artificial data from the dynamic metabolic system in the running example with the biochemically plausible parameter values [4] listed in column 2 of Table 1 and initial values: Gluc6P(0) = 1 mM, Fruc6P(0) = 0 mM, Fruc1,6P2(0) = 0 mM, ATP(0) = 2.1 mM, ADP(0) = 1.4 mM, and AMP (0) = 0.1 mM. The trajectory of this system is depicted in Figure 2. From Figure 2, the concentrations of all molecules except for Frucose-1,6-biphosphate reach their its steady states after about 0.1 minutes while Frucose-1,6-biphosphate after 0.5 minutes. Therefore, we do not use the data simulated after 0.5 minutes.

Table 1

The true value (from [4]), estimated value, and relative estimation errors.

Parameter name	True value	Estimated value	REE (%)
V max⁡,2 (mM·min−1)	50.2747	50.2447	0.0001
K ATP,1 (mM)	0.10	0.10000	0.0399
k 2 (mM−1·min−1)	2.26	2.2599	0.0049
V max⁡,3 f (mM·min−1)	140.282	139.4917	0.5633
V max⁡,3 r (mM·min−1)	140.282	141.3623	0.7701
K Gluc6P,3 (mM)	0.80	0.7999	1.3884
K Fruc6P,3 (mM)	0.15	0.1499	0.0930
V max⁡,4 (mM·min−1)	44.7287	44.6664	0.1372
K Fruc6P,4 (mM2)	0.021	0.0206	1.8457
k	0.15	0.1526	1.7447
k 5 (min−1)	6.04662	6.0466	0.0007
k 6 (min−1)	68.48	68.4837	0.0054
k 7 (min−1)	3.21	3.20797	0.0078
k 8f (min−1)	432.9	432.8408	0.0137
k 8r (min−1)	133.33	133.314	0.0120

Figure 2

Trajectory of system (3).

Although no noise is added to the artificial data in the simulation, noises are introduced in numerically calculating the derivatives by finite difference formulas. In general, the higher the sampling frequency and more data points are used, the more accurate the numerical derivatives are. On the other hand, we may not obtain data with the high frequency because of experimental limitations in practice. In this study, the sampling frequency is 100 data points per minute. In numerically calculating the concentration change rate at each time point from concentration x, we adopt the five-point central finite difference formula as follows:

The estimation accuracy of the proposed method is investigated in terms of relative estimation error which is defined as

As all parameters to be estimated are nonnegative, initial values are chosen as 0 or 1 in this study. The experimental results are listed in columns 3 and 4 in Table 1. From column 3 in Table 1, the estimated parameter values are very close to the corresponding true values. Actually the relative estimation errors calculated from (29) for all estimated parameters except for two are less than 1%. This indicates that the proposed method can accurately estimate the parameters in this system.

5. Conclusions and Future Work

In this study, we have first described a method to analyze the complexity of metabolic systems for parameter estimation, based on the stoichiometric matrix of the metabolic systems. As a result, the estimation of parameters in the metabolic systems has been reduced to the estimation of parameters in the improper rational functions or polynomial functions. Then we have developed an iterative linear least squares method for estimating parameters in the improper rational models. The results from its application to a metabolism system have shown that the proposed method can accurately estimate the parameters in metabolic systems.

We do not consider the noises in the data except those introduced by numerical derivatives in this study. One direction of future work is to investigate the influence of noises in the data to the estimation accuracy. In addition, low sampling frequency is expected, particularly for molecular biological systems as in practice measurements from them may be very expensive or it is impossible to sample measurements with high frequencies. Another direction of future work is to improve the estimation accuracy of the proposed method with low sampling frequencies.