Identification of Dynamic Metabolic Flux Balance Models Based on

Article pubs.acs.org/IECR

Identification of Dynamic Metabolic Flux Balance Models Based on Parametric Sensitivity Analysis Ricardo Martinez Villegas, Hector Budman,* and Ali Elkamel Department of Chemical Engineering, University of Waterloo, Waterloo, Ontario, Canada N2L 3G1

ABSTRACT: This paper deals with robust identification of dynamic metabolic flux model parameters based on parametric sensitivity analysis. First, the parameters in the model are ranked based on a global parametric sensitivity analysis to assess whether a subset of the parameters can be eliminated from further analysis. Then, the remaining significant parameters are identified based on the maximization of an overall parametric sensitivity measure subject to set based constraints that are derived from available data. The proposed method is illustrated on a case study of a model describing the diauxic growth of Escherichia coli in a batch culture.

1. INTRODUCTION This study discusses a new mathematical approach to perform identification of a dynamic metabolic flux model based on a global sensitivity analysis. Dynamic mathematical models involve sets of physicochemical parameters, which are mathematical variables that have to be calibrated to provide a good match between model predictions and available data. The current work focuses on parameter estimation for models based on dynamic flux balance analysis (DFBA), that is a particular method used to study microbial metabolic networks.1 This type of mathematical modeling approach can serve to simulate the metabolism of an individual cell by describing the flux distribution inside a cellular network. The method is based on maximizing a biological based objective such as growth rate subject to constraints on the rate of change of certain metabolites.1−4 Although other approaches have been developed to describe the dynamics of cell cultures, DFBA is generally advantageous in terms of the relatively smaller number of parameters that need to be calibrated to fit the data. However, the models resulting from DFBA may still need a large number of kinetic parameters to properly describe intricate kinetic regulations of the metabolic reactions of an organism. Generally, it is desirable to reduce the set of parameters that need to be calibrated for obtaining a desirable fit between model predictions and data to avoid overfitting of noisy data5 that can result in poor model prediction ability. Parametric sensitivity analysis can be used to reduce the number of model parameters that need to be calibrated. © XXXX American Chemical Society

Parametric sensitivity is often used to study the effect of changes in model parameters on model outputs of interest.6 When a system is operated within a region of high sensitivity, small changes in the parameter values will drastically affect the output. For this reason, it is crucial to quantify parametric sensitivity when designing, operating, or optimizing a system based on a model.7,8 For example, for biological systems, parametric sensitivity has also been used to guide the optimization of genetic circuits9 and for optimal design of biological experiments.10 Local parametric sensitivity analysis examines the individual effects of each parameter on the model outputs independently from each other at one specific moment. However, such analysis is often inaccurate, since it ignores correlations between the parameters. Consequently, sensitivity measures that explicitly account for these correlations11 are preferred. Parameter estimation must always be carried out with noisy data. In addition, unmeasured process disturbances are continuously occurring, thus creating further variability in the measurements.12 Hence, the model parameter values estimated from noisy data and in the presence of disturbances can only be identified within bounds with varying levels of confidence. Received: Revised: Accepted: Published: A

August 30, 2016 February 3, 2017 February 3, 2017 February 3, 2017 DOI: 10.1021/acs.iecr.6b03331 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research This work will focus on two main objectives related to the calibration of DFBA models: (i) the identification of a minimal number of DFBA model parameters that need to be calibrated for obtaining a desired quality of fit between model predictions and data and (ii) the calibration of the resulting subset of model parameters in the presence of noise and disturbances. Both objectives will be approached using a global parametric sensitivity analysis proposed in the literature.13 To tackle these two objectives, a novel approach is proposed that involves the maximization of a parametric sensitivity based metrics subject to set based constraints in the data. The solution of this optimization problem serves to identify the subset of model parameters with respect to which the sensitivity of model outputs is the highest. In addition, it is shown that the resulting parameters’ values have the largest associated probability under certain statistical assumptions regarding the outputs’ measurements. The paper is organized as follows: Section 2 briefly gives background information regarding dynamic metabolic flux models and the parametric sensitivity method used in the study. Section 3 shows the parameter estimation methodology based on a parametric sensitivity metric. Section 4 presents the results of the proposed method and comparison to parameter estimation based on least-squares minimization for a case study involving the diauxic growth of E. coli in a batch culture, followed by the conclusions in Section 5.

dgOxy(t ) dt

d g X (t ) dt

g 0 = g(t = 0) = [10.8, 0.4, 0.21, 0.001] V *Glucose dGlucose ≤ m dt K m + Glucose dOxygen ≤ Vo dt ⎡ v1 Acetate ⎢ −39.43 Glucose ⎢ ⎢ St = 0 Oxygen ⎢ ⎢−35 Biomass ⎢ ⎣1

v4

v2 9.46Glc + 12.92O2 → X v3 9.84Glc + 12.73O2 → 1.24Ac + X v4 −19.23Glc → 12.12Ac + X

(1)

where St represents the stoichiometric matrix m × n, corresponding to the simplified reactions listed above, with m being the metabolite number and n the number of fluxes, v’s represent the metabolic flux vectors in the system, gi and Sti are respectively the concentration of metabolite i and the row of the stoichiometric matrix corresponding to metabolite i (i = Oxy-oxygen, Glu-glucose, Ace-acetate, X-biomass), km is the substrate saturation constant ([mM]), kLa the oxygen transfer coefficient ([h−1]), Vo is the maximal oxygen uptake ([mmol/ gdw h]), and Vm is the maximal glucose uptake ([mmol/gdw h]). The nominal values of the parameters are presented in Table 1. Table 1. Nominal Parameter Values for Diauxic Growth Model Nominal Parameter Values

Km

kLa

Vo

Vm

0.015

7.5

15

10

It should be noticed that this is a highly simplified model that only describes the main nutrients (glucose, acetate, and oxygen) and the main byproduct. If the modeling of other metabolites is required, it is possible to couple the dynamic model in eq 1 to a steady state flux balance of other metabolites by using the known stoichiometric relations for E. coli. This has not been done in this study, since the steady state balance will not contain additional kinetic parameters and thus it is not directly relevant for the current study that focuses on parameter estimation. 2.2. Sensitivity Analysis Based Parameter Estimation Method. The sensitivity analysis used for the present research

subject to

dt

v3

v1 39.43Ac + 35O2 → X

vi

dgGlu(t )

⎤ ⎥ 0 1.24 12.12 ⎥ −9.46 − 9.84 19.23 ⎥ ⎥ −12.92 − 12.73 0 ⎥ ⎥⎦ 1 1 1 v2

Corresponding to reactions:

max (v1 + v2 + v3 + v4)

dt

= StXvX

g (t ) ≥ 0

2. BACKGROUND This section briefly reviews the parametric sensitivity method applied in the current study and the formulation of dynamic metabolic flux models. 2.1. Dynamic Flux Balance Algorithm. The model is based on using the metabolic reactions of the organism that define the metabolic network. The objective of the mathematical model is to describe the consumption and production of metabolites of a normal bioprocess. At each time interval, the model first solves for the fluxes (mol/h/unit biomass) associated with each one of the metabolic reactions composing the network. These fluxes are calculated from the maximization or minimization of a meaningful biological objective subject to constraints on the rate of change of the metabolites. Then, the calculated fluxes are used to predict the evolution of the metabolite concentrations at the next time interval by numerical integration. The specific case study to be analyzed in the present article is the diauxic growth of E. coli. For this case, as generally assumed for bacteria, the cell allocates resources so as to maximize the growth rate subject to the mass balance constraints and the positivity of the fluxes. The direction of the fluxes is assumed to be known a priori based on information from available databases such as KEGG (The Kyoto Encyclopedia of Genes and Genomes). A dynamic metabolic flux model for the E. coli diauxic growth model has been previously reported:1

dgAce(t )

= StOxy vX + kLa(0.21 − Oxygen)

= StAcevX = StGluvX B

DOI: 10.1021/acs.iecr.6b03331 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

The Uim describes the correlations among the parameters,11 since they involve contributions from each of the original parameters ηj. Then, the magnitudes of the singular values multiplying the corresponding Uim can be effectively used to rank the relative contributions of the principal components Uim to the outputs. Transformed parameters labeled as λi are defined that are related to the normalized model parameters ηj according to an orthogonal linear transformation as follows:

is based on the method proposed by Rand that explicitly accounts for correlations among parameters (2008).13 The first step consists of constructing a matrix made of the partial derivatives of each metabolite with respect to each of the parameters of interest. Normalization of the parameters is done to correctly account for differences of orders of magnitude of the parameters. The matrix is denoted by M and is defined as follows: ⎡ ∂g ∂g1 ⎢ 1 (t1) (t1) ∂η2 ⎢ ∂η1 ⎢ ∂g1 ⎢ ∂g1 ⎢ ∂η (t 2) ∂η (t 2) 2 ⎢ 1 ⎢ ⋮ ⎢ ⎢ ∂g1 ⋮ (t T) M=⎢ ⎢ ∂η1 ⎢ ⎢ ∂g2 (t ) ⋮ ⎢ ∂η 1 1 ⎢ ⎢ ⋮ ⎢ ∂gmp ⎢ ∂gmp (t T ) (t T) ⎢ ∂η2 ⎣ ∂η1

⎤ (t1) ⎥ ∂ηn ⎥ ⎥ ∂g1 ⎥ ... (t 2 ) ⎥ ∂ηn ⎥ ⎥ ⎥ ⎥ ⋱ ⋮ ⎥ ⎥ ⎥ ⎥ ⋮ ⎥ ⎥ ⎥ ⎥ ∂gmp ⎥ ... (t T)⎥ ∂ηn ⎦ ...

∂g1

λi =

δkj kj

= δ log kj

where the index i corresponds to the number of the transformed parameter λ and the index j refers to the original parameters, i.e. j = 1, ..., n. The Wij are elements of the matrix W defined before. Based on the singular value decomposition above, the deviations in the output variables can be expressed as a function of the singular values as follows: δgm(t ) =

⎛ Δt ⎞ ⎜ ⎟ M; 0 ≤ Δt ≤ T ⎝T ⎠

(7)

Combining and reordering eqs 6 and eq 7, the parametric sensitivity of each output i with respect to each parameter j can be calculated as follows: Sij = σiWij

(2)

(8)

and

∂gm ∂ηj

(3)

S

=

∑ Sij Uim (9)

i=1

A plot referred to as the sensitivity heat plot can be constructed using the information obtained from eqs 8 and 9. This plot is used to facilitate the visualization of the sensitivity, a time-dependent unit vector ( f i,m(t)) given as follows: fi ,m (t ) = σi(max|Wij|)|Ui ,m(t )| j

(10)

This unit vector is used to quantify the contribution of each of the transformed parameters λi (eq 6) toward the outputs m. A threshold is defined, e.g. 5% of the maximum value of f i,m(t), to distinguish between the significant versus the insignificant contributions of the parameters to the outputs of the problem. Finally, for visualization, the sensitivity coefficients Sij, given by eq 8, are used to construct an additional graph to be referred to as the Sensitivity Spectrum. This spectrum consists of the plot of the 3D bar graph log10|Sij| function of the transformed parameters λi and the original normalized parameters ηj. From this plot it is possible to simultaneously identify the original parameters ηj, which have the largest effects on the outputs. The effect of a parameter ηj on a particular output can be neglected if it results in low sensitivity coefficients for the significant λi, where the latter are determined according to the magnitude of the corresponding singular values. In summary, the two plots, i.e. the sensitivity heat plot and the sensitivity spectrum, are used together to assess the effect of the parameters on the outputs of interest.

(4)

For the special case of a linear dependency of the outputs with respect to the parameters where the parametric sensitivities remain constant over time, it is possible to show, using the definition of singular values, that the scaling in eq 4 leaves the singular values invariant with respect to the choice of Δt. Then a singular value decomposition of the resulting normalized matrix in eq 4 is derived:

M1 = Uσ VT

∑ λiσi Uim(t ) + O(|| δη ||2 ) i

where kj is the parameter to be examined for its effect, ηj is the normalized term of the parameter, n is the total number of parameters, gm is the model output at a specific time interval, and mp is the total number of outputs. The rows correspond to the different output variables from the model at each time interval, and the columns correspond to each of the parameters of the model. To mitigate the effect of the choice of the sampling interval used for discretization of the singular value decomposition results that follow, the M matrix is normalized as follows: M1 =

(6)

j

Using the following normalization of parameter deviations:

δηj =

∑ Wijδηj

(5)

where U has the same size as M1 (m × n) and these two matrices are orthogonal to each other, σ is a diagonal matrix (n × n) whose elements are the singular values. VT is an orthogonal matrix (n × n) that for clarity is further denoted as W(VT = W). The decomposition is needed to obtain the sensitivity principal components Uim(t), which are unit vectors and are columns of U, for which the subscript i refers to the number of principal components and m refers to the output m.

3. PARAMETER ESTIMATION APPROACH BASED ON PARAMETRIC SENSITIVITY ANALYSIS The key idea is to use the parametric sensitivity measures presented in the previous section to accomplish two main tasks C


Article

Industrial & Engineering Chemistry Research related to the identification of dynamic metabolic flux models’ parameters as follows: (i) To identify the parameters that do not significantly affect the outcomes of the model. (ii) To alleviate the computational effort of the parameter estimation problem based on the parametric sensitivity analysis. To accomplish these tasks, we propose to maximize a lumped measure of the individual parametric sensitivity coefficients defined in the previous section subject to set based constraints identified from data as follows:

representations have recently gained much attention in describing biological data for the purpose of identification of biological models.14 In biological systems it can be typically assumed that, because of measurement noise and/or unmeasured disturbances, the metabolites’ concentrations are always measured between upper and lower limits at each time interval for which data is collected, as shown in Figure 1. Since experimental runs are generally scarce due to experimental costs, it is difficult to assign a particular probability to trajectories within these sets. Instead, a uniform probability is assumed for each trajectory within the sets. The rationale for using the parametric sensitivity as the objective function in eq 11 is directly related to the use of set based constraints as follows. From eq 9 is possible to write that the deviations in the output gm can be approximated as follows:

max ∑ ∑ |Sij Ui,m| η̂

i

j

s. t. min g(t ) ≥ g(t ) ≥ max g(t ) ηU ≥ η̂ ≥ ηL

δgm(t ) ≈ Sg δη

(12)

m

(11)

where

where η̂ is the vector of estimates of the parameters η obtained from the solution of eq 11. The key objective of this optimization is to maximize the parametric sensitivity for all the data available for a particular experiment with respect to the values of the parameters. The combinations of two sources of variability were assumed for the data: unmeasured disturbances and measurement noise. Variations in model parameters from one batch to the next were used to simulate the effect of unmeasured disturbances on the process. For example, minor changes in growth media composition that may cause significant changes in the evolution of a cell culture may be suitably represented by equivalent variations in the maximal growth rate, e.g. parameter Vm given in model 1 given before. It should be noticed that these variations were used to produce in silico (simulated) data from the process while the model assumes fixed parameters’ values that must be identified. The heat sensitivity plot explained in section 2 can be used to find a time interval at which maximum parametric sensitivity is obtained. Then, based on the resulting sensitivity spectrum at the time interval of maximum sensitivity, it is possible to identify parameters that are below a certain threshold of significance and therefore do not need to be accurately calibrated. The constraints on the concentrations of metabolites gm at each sampling interval in eq 11 are referred to as set based constraints (Figure 1). These are used to describe the data by convex sets that can be easily handled by the optimization problem, while at the same time they are a suitable representation of the uncertainty in the data. Set based

S

S

S

Sg = [∑ Si1Ui m ∑ Si2 Ui m ... ...∑ Si n Ui m] m

i=1

i=1

i=1

(13)

and

δη T = [δη1 δη2...δηn]

(14)

Then, from eqs 12−14, a bound on the sum of absolute values of derivatives c

∑ i

δgm δηi

≤ || Sg (η)|| m

1

(15)

Considering the upper bound of eq 15 and integrating both sides gives: ηU

∫η ∑ L

i

δgm δηi

=

∫η

ηU L

|| Sg (η)|| |δη| m

1

(16)

where the LHS of eq 16 is bounded, since both gm and the η values are bounded according to eq 11. Then, if the integration on the LHS is also nonzero, ∥Sgm(η)∥1 can be properly normalized and used as a probability density function for η for any distribution of g that satisfies eq 16. Then, it is possible to define a joint probability for the occurrence of any combination of values of g and η as P(g ∩ η) = P(g |η)P(η)

(17)

Since, as stated above, it is assumed that for set based constraints the values of g are drawn from a uniform distribution, P(g|η) is consequently constant for any value of g within the sets, and thus: P(g ∩ η) = constant *P(η)

(18)

and since the probability P(η) can be calculated from the pdf(η) = constant × Sgm(η)1 as shown above, then we conclude from eq 18 that the joint probability of P(g∩η)can also be calculated from the pdf(η) = constant × Sgm(η)1. Accordingly, the parameter value of η with the maximal probability to occur is the one corresponding to the maximum of the pdf, i.e. the sensitivity function ||Sgm(η)1||, where the latter is identical to the objective function in eq 11 according to the definition in eq 13.

Figure 1. Set based constraint for glucose. D


Article

Industrial & Engineering Chemistry Research The arguments in eqs 12−18 are illustrated first for a very simple example using one value of g and one value of η. Let’s assume that δg as a function of different values of η̂ is given in Figure 2, which exhibits regions of variable sensitivity (δg/δη)

Figure 4. Probability distribution of η.

parameter value that has the largest probability to occur. The corresponding distribution of η will satisfy the set based constraint assumption of uniformly distributed changes of δgm within the sets. A resulting benefit from using parameters’ values obtained in regions of high parametric sensitivity is that they have smaller confidence regions based on Fisher Information Matrix arguments.15 This point is illustrated in the case study in the Results section. Finally, in order to speed-up the parameter estimation procedure, a two-step strategy is proposed as follows: (i) A search of nonsignificant parameters was first identified using a genetic algorithm to solve problem 11 with a coarse grid of the parameter space. (ii) The parameters that were found to be significant were then further calibrated using a genetic algorithm to solve eq 11 with a finer grid than the one used in step 1.

Figure 2. δgm(η) function of different solution η values.

where the latter is determined by the tangent to the graph in Figure 2. The specific function used in the example is: ⎡ (η − 0.5)2 ⎤ δg(η) = exp −⎢ ⎥ ⎣ 0.01 ⎦

(19)

It can be easily proved by derivation that the maximum sensitivity, i.e. the largest slope, occurs at η = 0.43. Then, the question is which particular distribution of η will result in a uniform distribution in the values of δgm as implied by the assumption of set based constraints on the data. If, for example, η would be uniformly distributed, it will result in a relatively low probability of occurrences of high values of δgm, whereas the probability of occurrences of small values of δgm will be very large, thus resulting in a nonuniform distribution of δgm. As shown by eqs 12−18, to obtain a uniform probability of δg, presented in Figure 3, the pdf of η should resemble the

4. RESULTS 4.1. Description of the Case Study. The parametric sensitivity based approach described above was used to estimate parameters in the model in the regions of high parametric sensitivity. The problem was solved using set based constraints that were obtained with simulated data. Accordingly, the upper and lower bounds of the sets were defined by positive or negative 0−10% parameter values’ changes with respect to the nominal values used in the model of Mahadevan et al. (2002)1 given in Table 1. These 0−10% fluctuations in model parameters were used to represent unmeasured disturbances in the process. As stated above, these disturbances may arise in a bioreactor due to changes in growth media, variability in the inoculums used to start each batch, or errors in initial conditions.16 In practice, only few trajectories may be available for formulating the set based constraints, and the upper and lower bounds of the available trajectories will be used to formulate the sets. In addition, Gaussian noise was added to all the outputs to simulate sensor noise. From the resulting simulations with combinations of parameter changes and noise, a family of trajectories that are presented in Figure 5A−5D for biomass, glucose, acetate, and oxygen, respectively, is obtained. Upper and lower bounds at regular time intervals were used as set based constraints in the optimization problem. The optimization problem was solved using a genetic algorithm from the Matlab optimization toolbox (The MathWorks Inc., Natwick, MA), since gradient descent based methods were found inefficient due to the presence of local maxima. The problem was solved for four parameters and with the following lower bounds, [0.0135, 6.75, 13.5, 9], and upper

Figure 3. Probability distribution of δgm(η).

function ∥Sgm(η)∥1, where η = 0.43 (the largest slope to the graph in Figure 2) will have the largest probability to occur while other values of η will have less probability. The probability distribution of η that results in a uniform distribution of δgm for this example is calculated numerically from inversion of the function δgm(η) above: η = 0.5 −

− 0.01*ln(δg + noise)

(20)

by assuming 1000 uniformly distributed samples of δgm (Figure 3) corrupted by 10% of Gaussian noise, and the resulting distribution is shown in Figure 4. Thus, the maximization of the parametric sensitivity function in eq 11 will result in η = 0.43, which corresponds to the E


Article


parameter, but the search region was approximately 10% of the coarser grid. The numbers of iterations for convergence were 340 and 95 for the coarse and fine grids, respectively. The solution of problem 11 was conducted in two sequential steps as follows: (i) a coarse optimization for eq 11 was done with the four parameters of the problem and (ii) a finer optimization of eq 11was conducted with a finer grid using only the parameters that were found to be relevant according to the parametric sensitivity analysis, fixing the nonsensitive parameters with a constant value from within the lower and upper limits previously defined. 4.2. Solution of Problem 11 with Coarse Grid To Find Significant Parameters. The heat plot and the sensitivity spectrum for the parameters’ values resulting from the maximization in eq 11 are presented in Figures 6A and 6B,

Figure 6. Global sensitivity analysis: (A) sensitivity heat map [acetate, glucose, oxygen, biomass]; (B) parametric sensitivity spectrum. Figure 5. Set based constraints for: (A) acetate, (B) glucose, (C) oxygen, and (D) biomass.

respectively. The heat plot shows the effect of the transformed parameters λi on the outputs of the problem, where the λi values capture the correlations between the parameters of the problem. The 3D sensitivity spectrum in Figure 6B shows the sensitivity coefficients given by eq 6 as a function of the parameters of the problem and the transformed parameters λi. From Figure 6A it is possible to observe that the transformed parameters (λi) with the highest influence on all outputs are λ1 followed by λ2 and λ3. The effect of these parameters is varying with time, with a maximum occurring at the time of glucose

bounds, [0.0165, 8.25, 16.5, 11]. These upper and lower limits were defined from the 0−10% fluctuation previously mentioned. The coarse optimization used a population size given by the combinations of 210 values for each parameter within the range between the upper and lower bounds defined above. The fine optimization used a population size corresponding to the combinations of 25 values for each F


Article


squares fitting is significantly smaller, thus resulting in larger confidence intervals, as shown below. Approximate confidence intervals were calculated for both estimation methods based on the FIM (Fisher Information Matrix) as follows:

depletion. Knowledge about the effect of the transformed parameters λion the outputs can be used in conjunction with Figure 6B to rank the significance of the contributions of the actual parameters ηj on the model. η1 (km) has a negligible contribution for each of the λi values, meaning that the substrate saturation constant has a negligible overall effect on the model outputs. Also, we observe from Figure 6B that the 4 (Vm) is the one with the highest contribution in the model since the variable λ1 that involves a large contribution from Vm is present in almost all the metabolite outputs (see Figure 6A). The second parameter in importance is the maximal uptake rate, η3 (Vo), followed by the oxygen transfer coefficient, η2 (kLa). It is clear from Figure 6B that for λ2 the most influential parameters are η2 and η3 (kLa and Vo). However, these two parameters have only a significant effect on the acetate concentration and a minor effect on glucose. This is expected due to the fact that in this model oxygen concentration controls the acetate uptake. In summary, the coarse level optimization is done with all four parameters but km can be kept constant at a constant value in the finer level optimization presented in the next section, thus saving computational time. If one is interested only in the effect of the parameters on the biomass at the end of the batch, it is clear from Figure 6A that only λ1 is important, in which case all parameters except for km are important, since they significantly contribute to λ1, according to Figure 6B. 4.3. Solution of Problem 11 with Fine Grid of Significant Parameters and Comparison to LeastSquares Regression. The quality of fitting was compared for the two methods in terms of the sum of squared errors and the parametric sensitivity coefficients for two cases: (i) parameters obtained from the maximization of sensitivity (problem in eq 11) solved with a coarse grid followed by solution with a finer grid using only three parameters as explained in the previous section; and (ii) parameters obtained from a standard least-squares nonlinear regression. The results for the sum of squared errors and the sensitivity norm ∥Sgm∥1 defined by eq 15 for these two sets of parameters are presented in Table 2, and the parameters sets corresponding to each method are presented in Table 3.

F = ST ∑−1 S

where S is a matrix whose elements are the partial differences of the outputs at each time interval with respect to each of the parameters, and it can be summed for all of the times analyzed, and ∑ is the covariance matrix of the measured noise. The FIM inverse matrix provides a lower bound for the covariance matrix of the errors in the parameter set,16where the elements of the inverse FIM are defined as: dii = F −1(i , i)

Sum of Squared Error

Sensitivity Coefficient

PS Method Least-Squares

1133.4 856.4

14.3402 12.1615

[θî ± tn ‐ p, α /2 dii ]

Km

kLa

Vo

Vm

Parameter Sets

7.6958 7.0853

13.4186 15.6051

9.0719 9.1247

PS Method Least-Squares

(23)

where t n − p, α /2 represents the t-distribution with n−p degrees of freedom and a confidence interval of 100 (1 − α)%.17 The resulting confidence intervals for both regression methods are presented in Table 4. Clearly, there is a substantial reduction in the magnitude of the confidence intervals: 8% for the oxygen transfer coefficient (kLa), 20% for the maximal oxygen uptake (Vo), and 65% for the maximal glucose uptake (Vm). Smaller confidence intervals are not only important for assessing the reliability of the parameter estimate but they may have great relevance for the purpose of robust optimization. For example, an optimization function relevant to the current case study may be the maximization of biomass at the end of the batch. In that case, the robust optimization function will involve the weighted sum of the expectation and the variance of the final biomass value. To illustrate the effect of the confidence intervals on a robust optimization formulation, we calculated the probability density function of the biomass at the end of the batch resulting from the model parameters’ uncertainty, where the latter is quantified by the confidence intervals. The parameters were assumed to be normally distributed with a variance equal to the confidence intervals obtained from the FIM based calculations in eq 18 and with mean values obtained either from the least-squares regression or from the maximization of the parametric sensitivity and shown in Table 3 for both methods. The mean and standard deviations in biomass for both regression methods are presented in Table 5, and the probability density functions of biomass for both methods are shown in Figure 7. As shown in Figure 7, the normal distribution histogram for the biomass based on the maximization of sensitivity shows a significantly narrower distribution as compared to the one obtained using least-squares. This considerable reduction in the variance in biomass shows the potential for reducing the conservatism of a robust optimization problem that optimizes a cost that involves this variance. In order to test the ability of the method to deal with a larger number of parameters, we studied a case with 5 parameters where the additional parameter was added as a constant to the numerator in the Monod expression for the glucose consumption constraint used in model 1. An additional

Table 3. Sum of Squared Error and Sensitivity Coefficient for Both Estimation Methods 0.0146 0.0165

(22)

The two sided confidence interval for the i-th estimated parameter is calculated as follows:

Table 2. Sum of Squared Error and Sensitivity Coefficient for Both Estimation Methods Curve Fitting Method

(21)

As expected, the sum of squared errors is smaller with the parameters obtained by least-squares regression as compared to the parametric sensitivity based regression, since the leastsquares regression uses the sum as the cost to be minimized. However, the parametric sensitivity coefficient for the leastG


Article


Table 4. Mean, Variance, and Confidence Intervals of the Significant Parameters for Both Parameter Estimation Methods PS Method

Least-Squares

Parameters

Mean

Variance

Confidence Intervals with tα,v = 1.746

Mean

Variance

Confidence Intervals with tα,v = 1.746

kLa Vo Vm

7.6958 13.4186 9.0719

0.0013 0.0041 0.0003

±0.0636 ±0.1114 ±0.0300

7.0853 15.6051 9.1247

0.0016 0.0065 0.0024

±0.0695 ±0.1406 ±0.0859

⎞ ⎛ Least‐Squares min⎜ ⎟ η ⎝ Parametric Sensitivity ⎠

Table 5. Mean and Standard Deviation of Biomass Concentration at Time 8 h for Both Parameter Estimation Methods Mean Biomass

Standard Deviation

Mean Biomass

Standard Deviation

S. t. Set based constraints from the metabolic flux model

0.8045

4.901e−05

0.8056

6.230e−05

ηU ≥ η ≥ ηL

Sensitivity Analysis

Least-Squares

(24)

The parameter estimates obtained from this optimization problem are shown in Table 6. It is evident that some of the Table 6. Parameter Sets of the Estimation Methods Parameter Set

Figure 7. Normal distribution of the biomass concentrations at time 8 h with initial concentrations of glucose 10.8 (mM) for both methods.

Case

km

kLa

Vo

Vm

Nominal Least-Squares Max PS Min (LS/PS)

0.015 0.016 0.014 0.016

7.5 7.085 7.695 7.176

15 15.61 13.42 14.35

10 9.124 9.072 9.180

parameters’ estimates are closer to the values calculated by least-squares regression as compared to the parameters obtained with the maximization of sensitivity. Also the sensitivities of the estimates were increased, thus resulting in smaller confidence intervals as compared to the estimates obtained by least-squares regression as shown in Table 7.

objective of this study was to test the ability of the algorithm to detect the lack of sensitivity to the artificially introduced fifth parameter. The 5-parameter case resulted as expected in correct recognition during the coarse optimization step of the two parameters that have no significant effect on outputs: the substrate saturation constant Km (as in the 4-parameter case solve above) and the artificially added parameter to the numerator of the Monod expression. In addition, it was found that the computational time in the 5-parameter case increased by only 4% as compared to the 4-parameter case. This relatively small increase in computations is due to the use of sensitive parameters in the optimization search, thus avoiding unnecessary iterations related to parameters with an insignificant effect on outputs. Also, it should be noticed that in terms of computational time for increasing number of parameters the proposed genetic algorithm based approach can be easily parallelized by partitioning the grid of parameters’ values among different processors. 4.4. Hybrid Estimation Approach. The least-squares methodology has been widely applied to estimate the parameter sets of many mathematical models used in engineering and for identifying biological models in particular.18 However, this method is only accurate if certain statistical assumptions are made.19,20 On the other hand, the minimal sum of square errors is often considered a good indicator of model accuracy. For this reason, a modification of the approach represented by Problem 11 can be used to achieve a trade-off between a low sum of square errors and a high parametric sensitivity, where the latter provides, as shown above, higher confidence in the parameter estimates. To achieve such a trade-off, the objective function in eq 11 can be modified as follows:

Table 7. Mean and Confidence Intervals of the Significant Parameters for the Minimization of the LS over the PS Method and Least-Squares Fitting Min (LS/PS)

Least-Squares

Param.

Mean

C. I.

Mean

C. I.

kLa Vo Vm

7.1764 14.348 9.1807

±0.0664 ±0.1014 ±0.0953

7.0853 15.605 9.1247

±0.0695 ±0.1406 ±0.0859

5. CONCLUSIONS This study investigates the application of parametric sensitivity analysis as a tool for identifying the parameters of a dynamic metabolic flux model. A novel approach to parameter estimation is proposed whereby the parametric sensitivity is maximized subject to set based constraints that are identified from data. It is shown that the proposed method served to find which parameters are not significant and to calculate parameters’ estimates with a large level of confidence. It was also shown that the estimates that correspond to maximal parametric sensitivity are the most probable ones following the assumption of uniformly distributed output changes within the set based constraints. Although the current paper illustrated the parametric sensitivity approach on E. coli growth, the proposed procedure is general. Dynamic metabolic flux models usually involve a H


Article


(18) Feng, X.; Tang, Y. J.; Dolan, K. D. Construction of a parsimonious kinetic model to capture microbial dynamics via parameter estimation. Inverse Probl. Sci. Eng. 2014, 22 (2), 309−324. (19) Englezos, P.; Kalogerakis, N. Applied parameter estimation for chemical engineers; CRC Press: 2000. (20) Van den Bos, A. Parameter estimation for scientists and engineers; John Wiley & Sons: 2007.

large number of metabolic reactions, thus making the identification of these models a very challenging task. The parametric sensitivity based approach is effective for such problems, since recognizing a priori the parameters with negligible effect on the outputs can serve to reduce the computational effort.

■

AUTHOR INFORMATION

Corresponding Author

*Tel.: +1 519 888 4567 ext 36980. Fax: +1 519 746 4979. Email: [email protected]. ORCID

Hector Budman: 0000-0002-0773-7457 Ali Elkamel: 0000-0002-6220-6288 Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS Financial support from CONACYT is acknowledged. REFERENCES

(1) Mahadevan, R.; Edwards, J. S.; Doyle, F. J. Dynamic flux balance analysis of diauxic growth in escherichia coli. Biophys. J. 2002, 83 (3), 1331−1340. (2) Budman, H.; Patel, N.; Tamer, M.; Al-Gherwi, W. A dynamic metabolic flux balance based model of fed-batch fermentation of bordetella pertussis. Biotechnol. Prog. 2013, 29 (2), 520−531. (3) Hjersted, J. L.; Henson, M. A. Optimization of fed-batch saccharomyces cerevisiae fermentation using dynamic flux balance models. Biotechnol. Prog. 2006, 22 (5), 1239−1248. (4) Orth, J. D.; Thiele, I.; Palsson, B. Ø. What is flux balance analysis? Nat. Biotechnol. 2010, 28 (3), 245−248. (5) Schaffer, C. Overfitting avoidance as bias. Machine learning 1993, 10 (2), 153−178. (6) Varma, A.; Morbidelli, M.; Wu, H. Parametric sensitivity in chemical systems; Cambridge University Press: 2005. (7) Iman, R. L.; Helton, J. C. An investigation of uncertainty and sensitivity analysis techniques for computer models. Risk Anal. 1988, 8 (1), 71−90. (8) Hamby, D. A comparison of sensitivity analysis techniques. Health Phys. 1995, 68 (2), 195−204. (9) Feng, X.; Hooshangi, S.; Chen, D.; Li, G. Optimizing Genetic Circuits by Global Sensitivity Analysis. Biophys. J. 2004, 87 (4), 2195. (10) Gadkar, K.; Gunawan, R.; Doyle, F. J., III. Modeling and analysis challenges in systems biology. In Proc. 9th Computer Applications in Biotech, Nancy, France, 2004. (11) Ades, A. E.; Lu, G. Correlations between parameters in risk models: estimation and propagation of uncertainty by Markov Chain Monte Carlo. Risk Anal. 2003, 23 (6), 1165−1172. (12) Michiels, W.; Engelborghs, K.; Roose, D.; Dochain, D. Sensitivity to infinitesimal delays in neutral equations. SIAM Journal on Control and Optimization 2002, 40 (4), 1134−1158. (13) Rand, D. A. Mapping global sensitivity of cellular network dynamics: Sensitivity heat maps and a global summation law. J. R. Soc., Interface 2008, 5, S59−S69. (14) Findeisen, R.; Imsland, L.; Allgower, F.; Foss, B. A. State and output feedback nonlinear model predictive control: An overview. European Journal of Control 2003, 9 (2−3), 190−206. (15) Walter, E.; Prozanto, L. Qualitative and quantitative experiment design for phenomenological modelsa survey. Automatica 1990, 26 (2), 195−213. (16) Dewasme, L.; Richelle, A.; Dehottay, P.; Georges, P.; Remy, M.; Bogaerts, P.; et al. Linear robust control of S. cerevisiae fed-batch cultures at different scales. Biochem. Eng. J. 2010, 53 (1), 26−37. (17) Draper, N. R.; Smith, H. Applied regression analysis; John Wiley & Sons: 2014. I


Identification of Dynamic Metabolic Flux Balance Models Based on

Recommend Documents