Email updates

Keep up to date with the latest news and content from BMC Bioinformatics and BioMed Central.

This article is part of the supplement: Seventh International Conference on Bioinformatics (InCoB2008)

Open Access Research

Dynamic sensitivity analysis of biological systems

Wu Hsiung Wu1, Feng Sheng Wang2* and Maw Shang Chang1

Author Affiliations

1 Department of Computer Science and Information Engineering, National Chung Cheng University, Chiayi 62102, Taiwan

2 Department of Chemical Engineering, National Chung Cheng University, Chiayi 62102, Taiwan

For all author emails, please log on.

BMC Bioinformatics 2008, 9(Suppl 12):S17  doi:10.1186/1471-2105-9-S12-S17


The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1471-2105/9/S12/S17


Published:12 December 2008

© 2008 Wu et al; licensee BioMed Central Ltd.

This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

A mathematical model to understand, predict, control, or even design a real biological system is a central theme in systems biology. A dynamic biological system is always modeled as a nonlinear ordinary differential equation (ODE) system. How to simulate the dynamic behavior and dynamic parameter sensitivities of systems described by ODEs efficiently and accurately is a critical job. In many practical applications, e.g., the fed-batch fermentation systems, the system admissible input (corresponding to independent variables of the system) can be time-dependent. The main difficulty for investigating the dynamic log gains of these systems is the infinite dimension due to the time-dependent input. The classical dynamic sensitivity analysis does not take into account this case for the dynamic log gains.

Results

We present an algorithm with an adaptive step size control that can be used for computing the solution and dynamic sensitivities of an autonomous ODE system simultaneously. Although our algorithm is one of the decouple direct methods in computing dynamic sensitivities of an ODE system, the step size determined by model equations can be used on the computations of the time profile and dynamic sensitivities with moderate accuracy even when sensitivity equations are more stiff than model equations. To show this algorithm can perform the dynamic sensitivity analysis on very stiff ODE systems with moderate accuracy, it is implemented and applied to two sets of chemical reactions: pyrolysis of ethane and oxidation of formaldehyde. The accuracy of this algorithm is demonstrated by comparing the dynamic parameter sensitivities obtained from this new algorithm and from the direct method with Rosenbrock stiff integrator based on the indirect method. The same dynamic sensitivity analysis was performed on an ethanol fed-batch fermentation system with a time-varying feed rate to evaluate the applicability of the algorithm to realistic models with time-dependent admissible input.

Conclusion

By combining the accuracy we show with the efficiency of being a decouple direct method, our algorithm is an excellent method for computing dynamic parameter sensitivities in stiff problems. We extend the scope of classical dynamic sensitivity analysis to the investigation of dynamic log gains of models with time-dependent admissible input.

Background

A mathematical model to understand, predict, control, or even design a real biological system is a central theme in systems biology. The most often used mathematical models for dynamic biological systems are formulated as nonlinear ordinary differential equations (ODEs). The critical challenges to get an ODE model are structure identification and parameter estimation of the model. To identify the structure and parameters of a dynamic model, the most important and essential job is to find the solution of ODEs efficiently and accurately. This job can be treated by analytical and numerical methods. Analytical methods are limited to certain special forms of ODEs, but numerical methods have no such limitations. There are several numerical methods can be used to solve ODEs [1-3], e.g., Taylor-series methods, modified Euler methods, and Runge-Kutta methods with variable step size control. The Taylor-series method takes more computation time to solve ODEs if the various derivatives are complicated, and the error is difficult to determine for arbitrary functions. The modified Euler method is a special case of a second-order Runge-Kutta method, and is more efficient compared to the Taylor-series method [4]. Fourth-order Runge-Kutta method is the most widely used ODE solver to meet requirements on both efficiency and accuracy. Collocation methods [5] are another common algorithms for solving ODEs and have been used for more than forty years. Wang [6] proposed a modified collocation method to transform ODEs into algebraic equations, and solved them by the Newton-Raphson method or an iteration method with step length restriction. The restricted step size is fixed and computed by trial and error. To overcome this drawback, we propose an adaptive step size control approach based on the Banach fixed point theorem for the modified collocation method in this paper.

There are different types of gains and dynamic sensitivities defined for sensitivity analysis [7-9]. The relative change of a dependent variable in response to a relative change in an independent variable is called a logarithmic gain, or a log gain. Log gains describe the change of dependent variables due to an environment change and are very useful for the assessment of the robustness and parameter estimation of a model. The change of a dependent variable in response to a change in a parameter is called a parameter sensitivity. In contrast to log gains, parameter sensitivities are the change of dependent variables correspond to a structure change in the model. The Biochemical Systems Theory (BST) [10] and Metabolic Control Analysis (MCA) [11-13] have achieved a great success in addressing the sensitivities at a steady state. However, the transient or periodic behavior is the primary interest in many systems (e.g., oscillation systems and fermentation systems that do not have a steady state). In these systems, the parameter sensitivities and log gains change with time, therefore the calculation methods for the steady state responses can not function. Dynamic sensitivity analysis is used in studying time-varying sensitivities in dynamic biological systems. A dynamic biological system can be characterized using logarithmic gains, sensitivities with respect to parameters and initial conditions. Several methods have been published to evaluate dynamic sensitivities [14-24]. They can be divided into the indirect methods (IDMs) and the direct methods (DMs). In the IDMs, the value of one dedicated parameter is varied infinitesimally while the values of other parameters are fixed. The model equations are solved anew for these sets of values of the parameters that differ in the value of the dedicated parameter only. The sensitivity of each variable with respect to this dedicated parameter is computed using the difference between the solutions of that variable for the two sets of values of the parameters divided by the infinitesimal difference of the dedicated parameter. In the DMs, using an ODE solver to solve the model equations and sensitivity equations simultaneously is the most used method for computing dynamic sensitivities. Shiraishi et al. [25] published a variable-order, variable-step Taylor-series method that can be used as an ODE solver providing a highly accurate calculation to compute dynamic sensitivities. This method is limited to general mass action (GMA) models described by power-law differential equations. Runge-Kutta methods with variable step size control can be used to compute dynamic sensitivities for most of the nonlinear differential equations, but is inefficient to determine the step size in a large dimensional system including the model differential equations and sensitivity differential equations. Due to the efficiency, Dunker [15] proposed the decoupled direct method (DDM), in which the sensitivity equations are solved separately from the model equations. He said: "the decoupled method has advantages in simplicity, stability, accuracy, efficiency, storage requirements, and program size over other methods". Although the DDM approach has so many advantages, the step size for the time profile determined by the error control of model equations is unable to be used for the sensitivities when the sensitivity equations are more stiff than the model equations and will generate inaccurate results.

Dynamic sensitivity analysis evaluates the influences on dependent variables due to variations of parameters, initial conditions and independent variables. In many practical applications, e.g., the fed-batch fermentation systems, the system admissible input (corresponding to independent variables of the system) can be time-dependent. The main difficulty for investigating the dynamic log gains of these systems is the infinite dimension due to the time-dependent input. Shiraishi et al. [26] proposed an efficient method, based on a combination of the recasting technique and the Taylor-series method, for calculating the time courses of log gains to investigate the dynamic behavior of log gains for oscillation models with a limit cycle. The method is limited to the computations of dynamic log gains with respect to constant independent variables. We extend the computations of dynamic log gains to a model with continuous time-varying admissible input based on the finite parameterization method (PM). The classical PM was created for numerical solutions of optimal control problems [27]. The central idea of the method relies on a simple approximation mechanism. The whole time domain of a continuous admissible input is partitioned into several subintervals, and the input for each subinterval is approximated by a piecewise constant function. The dynamic log gains with respect to the continuous admissible input can be computed based on the partial derivations of dependent variables with respect to the piecewise constant input [28-30].

In this paper, we present an algorithm with an adaptive step size control that can be used for computing the solution and dynamic sensitivities of an autonomous ODE system simultaneously. This algorithm is the modified collocation method, proposed by Wang [6], with an adaptive step size control approach. Although our algorithm is one of the decouple direct methods in computing dynamic sensitivities of an ODE system, the step size determined by model equations can be used on the computations of the time profile and dynamic sensitivities with moderate accuracy even when sensitivity equations are more stiff than model equations. In the algorithm, the modified collocation method is used to transform model and sensitivity equations into algebraic equations, and the approximated solution is solved by an iteration method. This algorithm can be extended easily to solve problems of mixed differential and algebraic equations (DAEs) by combing algebraic equations with that transferred from differential equations. In our algorithm for computing dynamic sensitivities of an ODE system, the model equations and sensitivity equations are solved alternatively in two stages. First, the model equations are advanced from ti to ti + η using the iteration method, where η is the step size decided by model equations based on the fixed-point theorem. Second, the solution of model equations at ti + η and the same step size are used to propagate the sensitivity equations from ti to ti + η. For dynamic systems with continuous time-dependent admissible input, the dynamic log gains are computed based on the parameterization techniques. The PM is used to approximate the original infinite-dimensional problem by a finite dimensional one with piecewise constant input. The dynamic log gain for this approximation problem is defined as the percentage change of a dependent variable in response to an infinitesimal percentage change for each piecewise constant input.

To show this algorithm can perform dynamic sensitivity analysis on very stiff ODE systems with moderate accuracy, it is implemented and applied to two sets of chemical reactions: pyrolysis of ethane and oxidation of formaldehyde. The accuracy of this algorithm is demonstrated by comparing the dynamic parameter sensitivities from this new method and that from the direct method with Rosenbrock stiff integrator based on the indirect method. The same dynamic sensitivity analysis is performed on an ethanol fed-batch fermentation system with a time-varying feed rate to evaluate the applicability of the algorithm to realistic models with time-dependent admissible input.

Results and discussion

To illustrate the accuracy of our algorithm, it is implemented and applied to stiff chemical mechanisms for the pyrolysis of ethane as well as the oxidation of formaldehyde. These systems have been shown to be unstable using both the DM and the Green's function method [15]. The same dynamic sensitivity analysis is performed on an ethanol fed-batch fermentation system with a time-varying feed rate to evaluate the applicability of the algorithm to realistic models with time-dependent admissible input.

Pyrolysis of ethane

The chemical mechanism for the pyrolysis of ethane is a very stiff system and consists of seven species in five reactions. The chemical reactions and rate constants are shown in Table 1 and are described by GMA model equations as follows:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M1">View MathML</a>

where [x] is the concentration of species x and ki is the rate constant. The initial concentration of C2H6 is 5.951 × 10-6 mol/cm3 and all other initial concentrations are zeros. All sets of sensitivity coefficients with respect to all rate constants and initial conditions are computed simultaneously without any difficulty using our algorithm with a tolerance of 10-7. The normalized sensitivity coefficients for the pyrolysis of ethane at 1 s and 20 s calculated by our algorithm are shown in Table 2. The results obtained by the indirect method (IDM) according to the finite difference approximation and the direct method with Rosenbrock stiff integrator (R/DM) [24] are also given in Table 2 for comparison. The results of our algorithm are of equal accuracy to R/DM in comparison to the IDM and the maximum relative error is 0.58%.

Table 1. The mechanism for ethane pyrolysis.

Table 2. Sensitivity coefficients for ethane pyrolysis.

Oxidation of formaldehyde

The formaldehyde oxidation mechanism is a larger system, involves 15 species in 25 reactions. The chemical reactions and rate constants are shown in Table 3 and are described by GMA model equations as follows:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M2">View MathML</a>

where [x] is the concentration of species x and ki is the rate constant. The initial concentrations in mol/cm3 are [CH2O] = 1.124 × 10-7, [O2] = 2.109 × 10-6, [CO] = 4.699 × 10-6, [M ] = 1.1772 × 10-5, and all other initial concentrations are zeros. Sensitivity coefficients with respect to all rate constants and initial conditions are computed successfully using our algorithm with a tolerance of 10-9. The normalized sensitivity coefficients for O and H2O at 0.005 s calculated by our algorithm are presented in Table 4. The results obtained by IDM and the direct method with Rosenbrock stiff integrator (R/DM) are also given in Table 4 for comparison. Our results are in good agreement with the R/DM in comparison to the IDM, and the maximum relative error is 0.25%. The discrepancies between the results of our algorithm and the R/DM method are sufficiently small to prove that this new method is capable of performing dynamic sensitivity analysis for stiff differential equations as accurate as direct methods.

Table 3. The mechanism for formaldehyde oxidation.

Table 4. Sensitivity coefficients for formaldehyde oxidation.

Ethanol fed-batch fermentation

The dynamic sensitivity analysis of an ethanol fed-batch fermentation process, a real dynamic biological system never reaching a steady state, is used to elucidate the applicability of our algorithm. Wang et al. [31] built a mathematical kinetic model of fermentation for ethanol and glycerol production using Saccharomyces diastaticus LORRE 316, which is a high ethanol tolerance yeast. The mathematical kinetic model for the fed-batch process consists of the dynamic behavior of biomass, glucose, ethanol and glycerol, and its dynamic mass balance equations are expressed as follows:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M3">View MathML</a>

where x is the concentration of cell mass, s is the concentration of glucose, p1 is the concentration of ethanol, p2 is the concentration of glycerol, V is the working volume of the fermenter, tf is the final fermentation time, τ = t/tf is the normalized fermentation time, sF is the feed concentration of glucose, F is the feed rate, <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M4">View MathML</a> is the ethanol yield factor, and <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M5">View MathML</a> is the glycerol yield factor. The unstructured kinetic models for the specific cell growth and product formation are respectively expressed as follows:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M6">View MathML</a>

Using a batch fermentation model, Wang et al. obtained the optimal values of 19 parameters [31]. The initial and feed concentrations of glucose are set to 10 and 200 g/L, the initial concentration of biomass is set to 2 g/L, and the starting working volume is set to 1.5 L in the computations of optimal feed rate and optimal fermentation time to maximize the ethanol production rate J = p1V/tf under some physical constraints, e.g., the residual glucose restriction s(tf) ≤ sr for reducing the separation cost, sr is the concentration of the desired residual glucose. The optimal final fermentation time is 12.836 hours and the optimal feed rate F* for the fed-batch fermentation model [31] is as follows:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M7">View MathML</a>

(1)

Figure 1 shows the computational time profile of the ethanol fed-batch fermentation model with the optimal feed rate.

thumbnailFigure 1. Dynamic behavior of the ethanol fed-batch fermentation model. Time profiles of cell mass (x), glucose (s), ethanol (p1), glycerol (p2), and the working volume (V) for the ethanol fed-batch fermentation computed using the optimal feed rate F* and the feed glucose of 200 g/L. The horizontal scale is in normalized fermentation time (t/tf).

Our algorithm is applied to the ethanol fed-batch fermentation model using the initial conditions as described above. All dynamic sensitivities with respect to 22 parameters (including sF, F and tf) and initial conditions, and the dynamic log gains with respect to time-varying feed rate are computed simultaneously without any difficulty. Figure 2 shows the dynamic relative sensitivities with respect to μm, <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M4">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M5">View MathML</a>, and sF. When the maximum specific growth rate μm is increasing, the rate of consuming glucose is increasing such that the concentration of residue glucose is decreasing. This situation is compatible with the trend in Figure 2(a). The increases in the ethanol and glycerol yield factor cause the increases in the production of ethanol and glycerol, and more glucose remains at the final time. As Figures 2(b) and 2(c) show, to increase the production of ethanol and glycerol by improving the ethanol yield factor is better than by increasing the glyverol yield factor. Figure 2(d) shows that if the feed concentration of glucose is increasing, the cell growth and the production of ethanol and glycerol are increasing. Under this condition, S. diastaticus LORRE 316 is unable to completely consume glucose to produce ethanol during the fermentation time and more glucose remains at the final time.

thumbnailFigure 2. Dynamic relative sensitivities with respect to μm, <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M4">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M5">View MathML</a>, and sF. (a) relative sensitivities with respect to μm; (b) relative sensitivities with respect to <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M4">View MathML</a>; (c) relative sensitivities with respect to <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M5">View MathML</a>; (d) relative sensitivities with respect to sF. The horiziontal scale is in normalized fermentation time (t/tf).

The relative sensitivity with respect to tf is shown in Figure 3(a). As expected, an increase in tf causes a low relative increase in the concentration of cell mass and a high relative decrease in the concentration of residue glucose. Figures 3(b), 3(c) and 3(d) show the dynamic relative sensitivities with respect to the initial conditions x, s, and V. When the initial concentration of cell mass increases, the residue glucose decreases, and the production of ethanol will increase a little, but the production of glycerol will decrease a little at the final fermentation time. Starting the fermentation process with more glucose will cause more glucose to remain and the production of ethanol and glycerol to increase a little at the final fermentation time as shown in Figure 3(c). Figure 3(d) shows that if the initial working volume is increasing, all concentrations of cell mass, glucose, ethanol, and glycerol are decreasing at the final fermentation time.

thumbnailFigure 3. Dynamic relative sensitivities with respect to tf, x(0), s(0), and V(0). (a) relative sensitivities with respect to tf; (b) relative sensitivities with respect to the initial value of x; (c) relative sensitivities with respect to the initial value of s; (d) relative sensitivities with respect to the initial value of V. The horizontal scale is in normalized fermentation time (t/tf).

We are interested in the ethanol production rate J in the fermentation process. The effects on J with respect to μm, <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M4">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M5">View MathML</a>, sF, and tf are shown in Figure 4. To increase J, it is clear that an increase in <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M4">View MathML</a> or sF will have more impact than an equal relative increase in <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M5">View MathML</a> or μm. The negative value of relative sensitivity for J with respect to tf means a decrease in the fermentation time will get a higher J at the expense of more residual glucose. Though the relative sensitivity of J with respect to sF is higher than that with respect to <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M4">View MathML</a> at the final fermentation time, by increasing sF to increase J will cause more glucose left at the final time and increase the cost to separate the residue glucose and the ethanol product. We can make a conclusion that to increase J by increasing <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M4">View MathML</a> will be a better choice than increasing sF or <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M5">View MathML</a>, and decreasing the fermentation time.

thumbnailFigure 4. Dynamic relative sensitivities of J. The relative sensitivities of ethanol production rate with respect to μm, <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M4">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M5">View MathML</a>, sF, and tf. The horizontal scale is in normalized fermentation time (t/tf).

The feed rate F(t) of the fed-batch fermentation model is a time-dependent input control variable, so that the computation of the effect on J with respect to F(t) is an infinite dimensional problem. The fermentation time is divided into ten equal time partitions, and the optimal feed rate F* for the fed-batch fermentation model [31] is approximated by ten piecewise constant functions. The ten input control parameters, denoted by Fi, i = 1,..., 10, are shown in equation (1). The proposed algorithm computes the dynamic log gains based on the parameterization method. All dynamic log gains with respect to Fi, i = 1,..., 10, are computed (data not shown here) with the parameter sensitivities simultaneously. The dynamic log gains of J with respect to Fi, i = 1,..., 10, are computed by

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M8">View MathML</a>

Due to the optimal values of Fi, i = 6,..., 10 are equal to or very close to 0, The dynamic log gains of J with respect to them are small and can be ignored. The dynamic log gains of J with respect to Fi, i = 1,..., 5, are shown in Figure 5. The effects on J are in decreasing order from F1 to F5. Increasing the feed rate at an early stage will get a higher J at the final fermentation time than that at a later stage without considering the residual glucose.

thumbnailFigure 5. Dynamic log gains of J w.r.t. F(t). The dynamic log gains of ethanol production rate with respect to Fi, i = 1,..., 5. The horizontal scale is in normalized fermentation time (t/tf).

Conclusion

To deeply study the dynamic behavior of a biological system, one of the methods is to model it as a mathematical model. The most used mathematical model for simulating biological systems is the ODE model. The essential task for modeling and simulating a biological system is to find the solution of an ODE model efficiently and accurately. We present an algorithm with an adaptive step size control that can be used for computing the solution and dynamic sensitivities of an ODE system simultaneously. Instead of using error control to decide the step size in solving the model equations, our algorithm computes the step size based on the fixed-point theorem and the same step size can be used in solving the sensitivity equations.

Dynamic sensitivity analysis is a useful tool to investigate the behavior of dynamic systems. In the direct methods for solving the dynamic sensitivities, sensitivity equations and model equations are coupled and solved together at the expense of more computation time. In contrast, sensitivity equations and model equations are solved separately in the decouple direct methods. The DDMs are more efficient than the DMs due to the dimension of ODEs. The chief disadvantage of DDMs is the requirement of error control on both model equations and sensitivity equations. Our algorithm with an efficient step control approach based on the fixed-point theorem is used to address the disadvantage of DDMs. Analogous to the DMs, the same step size obtained by model equations is used on both model and sensitivity equations. It has been implemented and applied to well-known stiff problems with the same accuracy compared to the direct method with Rosenbrock stiff integrator (R/DM). As our algorithm is one of the DDMs, it has the efficiency of the DDMs and the same accuracy of the DMs as presented in the section describing the results. By combining the efficiency and accuracy, our algorithm is an excellent method for computing dynamic parameter sensitivities in stiff problems.

We extend the scope of classical dynamic sensitivity analysis to the investigation of dynamic log gains of models with time-dependent admissible input. The parameterization method is used to approximate the infinite-dimensional computation problem for dynamic log gains in models with time-dependent admissible input by a classical finite-dimensional computation problem of dynamic log gains. Then, all dynamic log gains and parameter sensitivities can be obtained simultaneously from our algorithm. Appropriate parameterization allows one to obtain a more efficient way to compute the dynamic log gains with respect to a continuous time-dependent input than that by finite difference approximation. Finally, the new proposed algorithm is applied to the ethanol fed-batch fermentation system, a real dynamic biological system which never reaches a steady state, with a time-varying feed rate for elucidating the applicability to realistic models with time-dependent admissible input. Through the dynamic sensitivity analysis of the ethanol fed-batch fermentation model, we conclude that to get a higher ethanol production rate by increasing the ethanol yield factor is a good choice.

Methods

A dynamic biological system is always modeled as a nonlinear ODE system:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M9">View MathML</a>

(2)

where x(t) ∈ ℝn is a vector of dependent variables, y(t) ∈ ℝm is a vector of independent variables, θ ∈ ℝp is a vector of parameters, v ∈ ℝq is a vector of fluxes between the variables, N ∈ ℝ(n+m) × q is the stoichiometric matrix describing the interconnecting fluxes, x0 and y0 are initial concentrations of x and y respectively. Using a different kinetic description for v results in a different mathematical model. In the Michaelis-Menten model, each element vi of v is of the form

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M10">View MathML</a>

(3)

where xj is the substrate, <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M11">View MathML</a> is the maximum flux for vi, and Ki is the half-saturated flux for vi. For the GMA systems, the kinetic equation for vi is expressed as a power-law function

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M12">View MathML</a>

(4)

where xj x for j n, xj y for j >n, γi is the rate constant, and gij is the kinetic order for each xj. Equation (2) can be expressed concisely as:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M13">View MathML</a>

(5)

where the function f is assumed to be continuous and differentiable in all its arguments x, y,and θ. This assumption on f is satisfied for both of equations (3) and (4).

For a model described by a nonautonomous dynamic system as follows:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M14">View MathML</a>

(6)

we let xn+1 = t and dxn+1/dt = 1, equation (6) can be rewritten as equation (5) with x(t) ∈ ℝn+1. This is an (n + 1)-dimensional autonomous dynamic system. Similarly, an n-dimensional time dependent equation is a special case of an (n+ 1)-dimensional autonomous dynamic system. Using this trick, we can always remove any time dependence by adding an extra dimension to the system. Thus, without loss of generality, we will consider the autonomous dynamic systems expressed by equation (5) unless stated otherwise.

ODE solver

Given a set of ODEs expressed as equation (5) and a set of time points T = {ti|i = 1,..., k}. An ODE solver is to find the value of x(ti), ti T for a given θ. Many ODE solvers with variable step size control can be used to solve equation (5). Wang [6] proposed a modified collocation method with Lagrange polynomials as shape functions to transform ODEs into algebraic equations. The whole time domain of the problem is divided into a number of non-overlapping intervals [ti-1, ti], i = 1,..., k. The unified formulas of the modified collocation method for the subinterval [tj-1, tj], ti-1 tj-1 <tj ti, can be expressed as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M15">View MathML</a>

(7)

where ηj is the step size in time tj, Î is an identity-like matrix, and the coefficient matrices D and <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M16">View MathML</a> depend on the shape functions. The accuracy and efficiency of collocation methods depend largely on the degree of shape functions. The modified collocated equations with piecewise linear polynomials transformed from equation (5) for each subinterval [tj-1, tj], ti-1 tj-1 <tj ti have the same formulas as the modified Euler method:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M17">View MathML</a>

(8)

Instead of solving the ODEs in equation (5) directly, we find the solution of algebraic equations in equation (8) step-by-step for each time interval [ti-1, ti], i = 1,..., k. The solution obtained from equation (8) is a good approximation solution of ODEs in equation (5) when the step size ηj is small enough. How to decide the step size is an important problem for all ODE solvers. A larger step size can cause the solution to be inaccurate and divergent, and a smaller step size is inefficient for the computations. Wang [6] uses the Newton-Raphson method with step length restriction to solve equation (8). The restricted step size is fixed and computed by trial and error. To overcome this drawback, we propose an adaptive step size control approach based on the Banach fixed point theorem for the modified collocation method in this paper. We will show that this approach can determine the step size automatically and efficiently when computing the solutions and dynamic sensitivities of equation (5) simultaneously.

The Banach fixed point theorem and some terminologies for describing the theorem are defined below.

Definition 1. Metric space [32]

A metric space (X, d) is a set X where a notion of distance d (called a metric) between elements of the set is defined.

Definition 2. Cauchy sequence [32]

A sequence (xn) in a metric space (X, d) is said to be Cauchy if for every ε > 0 there is a positive integer N such that for all natural numbers m, n >N, the distance d(xm, xn) is less than ε.

Definition 3. Complete metric space [32]

A metric space (X, d) in which every Cauchy sequence has a limit in X is called complete.

Theorem 4. Banach fixed point theorem [32]

Let (X, d) be a non-empty complete metric space. A mapping ψ: X X is called a contraction on X if there is a nonnegative real number q < 1 such that for all a, b in X

d(ψ(a), ψ(b)) ≤ q·d(a, b).(9)

The contraction ψ on X admits one and only one fixed point x* ∈ X such that ω (x*) = x*.

We now show how to decide the step size ηj in equation (8) based on the Banach fixed point theorem. Equation (8) is an implicit expression of x, and it can be rewritten as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M18">View MathML</a>

(10)

where c = x(tj-1)+0.5ηjf(x(tj-1), y(tj-1); θ) is a constant vector. When the values of independent variable y and θ are given, the solution of x = g(x) can be found by an iteration process. A sequence of values of x(tj) is obtained using the iterative rule. If xi(tj) tends to a limit x*(tj) when i → ∞, it is the answer of x = g(x) and called a fixed point of the function g(x). Let X be the set of xi(tj), i = 1,..., ∞ and d(a, b) = ||a - b||p where a and b are arbitrary xi(tj) and ||·||p is the p-norm. Then (X, d) forms a non-empty complete metric space. If g is a contraction on X, the Banach fixed point theorem guarantees the existence of a fixed point and the convergence of the iteration process to that fixed point. By the equation (9) and the definition of distance function d, for a, b X we obtain

||g(a) - g(b)||p q||a - b||p, q < 1.(11)

We suppose that g(x) is a continuous and differentiable function on X. By the generalized mean value theorem and the definition of matrix norm, we have

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M19">View MathML</a>

(12)

Comparing equations (11) with (12), we obtain

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M20">View MathML</a>

(13)

where ||·||p is the p-norm of a matrix. By substitution of the term on the right of the equal sign in equation (10) for g, equation (13) can be rewritten as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M21">View MathML</a>

(14)

This equation is used to compute the maximum ηj with p = 2 when the process of finding the solution of equation (8) is in progress. The Jacobian matrix ∂f/∂x in equation (14) must be evaluated at each time tj. The computation of the Jacobian matrix can be done by evaluating the analytic formula of the partial derivative of f with respect to x which is user-provided, or by the finite difference approximation. For the GMA systems, the model equations are expressed in power-law and the value of the Jacobian matrix can be straightforwardly obtained using the analytic formula as follows:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M22">View MathML</a>

where N is the stoichiometric matrix, vi is the ith element of v ∈ ℝq and gij is the kinetic order for each xj x in vi. For efficiency, we approximate the value of 2-norm of the Jacobian matrix ∂f/∂x by ||∂f/∂x||Δ, where n is the dimension of x and ||∂f/∂x||Δ is the maximum absolute value of the element of the Jacobian matrix. The proposed algorithm AMCM, Adaptive Modified Collocation Method, is shown as follows

Algorithm AMCM

Input:

1. A set of n ordinary differential equations <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M23">View MathML</a> = f(x, y) with n dependent variables xi, i = 1,..., n and m independent variables yi, i = 1,..., m.

2. Two order sets x0 = {xi(t0)|i = 1,..., n} of initial values of x and y0 = {yi(t0)|i = 1,..., m} of initial values of y.

3. An order set T = {t1,..., tk} of sampling points, ti, 1 ≤ i k is the sampling time of the solution of each ODE, k is the number of sampling points.

4. A tolerance ε

Output: The set of solutions of dependent variables at each sampling time.

• For each sampling time ti in T.

1. ηj ti - ti-1, dt ← 0, xc x(ti-1), yc y(ti-1)

2. Repeat the following steps until dt = ti - ti-1.

(a) xp xc, yp yc

(b) Evaluate the Jacobin matrix <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M24">View MathML</a>

(c) Compute the upper bound μ of the value of ||A||2 by <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M25">View MathML</a>

(d) If μ * ε ≥ 1, it means the ODEs are stiff, then exit this algorithm.

(e) If μ * ηj > 1, then update ηj with 0.9/μ

(f) Call iteration algorithm to compute the value of xc stepped forward ηj from xp.

(g) If the iteration algorithm succeeds in computing xc, then dt dt + ηj and ηj ti - ti-1 - dt, otherwise exist this algorithm.

3. x(ti) ← xc, y(ti) ← yc.

• return x(ti), i = 1,..., k.

End of Algorithm AMCM

Algorithm Iteration

Input:

1. A set of n ordinary differential equations <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M23">View MathML</a> = f(x, y).

2. x(t), y(t), ηj and the iteration limitation.

Output: x(t + ηj).

1. Evaluate the value of f(x(t), y(t)).

2. x(t + ηj) ← x(t) + f(x(t), y(t)) * ηj.

3. y(t + ηj) ← y(t) + f(x(t), y(t)) * ηj.

4. Repeat the following steps until the iteration limitation is reached or the value of x(t + ηt) converges.

(a) Evaluate the value of f(x(t + ηj), y(t + ηj)).

(b) x(t + ηj) ← x(t) + 0.5 * ηj * (f(x(t), y(t)) + f(x(t + ηj), y(t + ηj)))

5. If the iteration limitation is reached, then exit this algorithm; otherwise, return x(t + ηj).

End of Algorithm Iteration

Dynamic sensitivity Solver

For a model described by equation (5), the absolute parameter sensitivity s(xi, θj) of dependent variable xi x with respect to parameter θj θ is defined as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M26">View MathML</a>

(15)

where xi(t; θj + Δθj) is the ith component of the solution of equation (5) with an increment Δθj on the jth parameter. The function xi(t; θj + Δθj) can be expanded into a Taylor series as follows:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M27">View MathML</a>

(16)

where 0 <ξ < 1. If Δθj is sufficiently small, the last term of equation (16) can be truncated, leading to a linear approximation of xi(t; θj + Δθj). Substituting the linear approximation of xi(t; θj + Δθj) into equation (15) leads to

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M28">View MathML</a>

This is defined as the first-order local sensitivity of xi with respect to θj [33]. The relative parameter sensitivity S(xi, θj) of xi with respect to θj is defined as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M29">View MathML</a>

Similar to the parameter sensitivity, the absolute log gain l(xi, yj) and log gain L(xi, yj) of xi x with respect to yj y are expressed respectively as follows:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M30">View MathML</a>

Once the local sensitivity is known, the calculation of the relative sensitivity is straightforward. So, for briefing, we limit our explanation on the absolute sensitivity only below.

The absolute dynamic sensitivity of xi with respect to θj is given as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M31">View MathML</a>

(17)

where fi is the ith element of f [34]. The absolute dynamic log gain of xi with respect to yj is similar to equation (17) by replacing θj, s(xi, θj) with yj, l(xi, yj) respectively when yj(t) is a constant. In the case which yj(t) is a continuous time-dependent function, the whole time domain of yj(t) is partitioned into Nu time intervals (tk-1, tk), k = 1,..., Nu. Function yj(t) is parameterized by the piecewise constant functions ωk(t), k = 1,..., Nu, as follows:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M32','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M32">View MathML</a>

(18)

where uk, k = 1,..., Nu, are constant input control parameters and

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M33','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M33">View MathML</a>

The continuous time-dependent function yj(t) is approximated by equation (18). The dynamic log gain l(xi, yj) of infinite dimension is transferred into Nu dynamic log gains of dimension one with respect to the input control parameters, and expressed as.

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M34">View MathML</a>

where j = 1,..., Nu.

We extended the proposed algorithm AMCM to compute the dynamic sensitivities. The dynamic sensitivities with respect to parameters (include the initial conditions) and absolute dynamic log gains can be computed simultaneously. Let u be the vector of input control parameters and r be an Nr dimensional vector of constants which contains constant independent variables in y and the input control parameters in u. When all components of y are constant, the vector r is equal to y. Let z be an n + p + Nr dimensional vector which contains model parameters θ, initial conditions of dependent variables x0 and constant input in r. Φ indicates a matrix of size n × (n + p + Nr) which contains the absolute sensitivities with respect to model parameters and initial conditions, and the absolute log gains with respect to constant independent variables and input control parameters. z and Φ have the following form:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M35','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M35">View MathML</a>

The model equations are rewritten as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M36','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M36">View MathML</a>

The sensitivity equations in matrix form can be derived by applying the chain rule to the derivative of Φ:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M37','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M37">View MathML</a>

(19)

Let ϕi be the ith column vector of Φ and zi be the ith element of z. The matrix sensitivity equations in equation (19) are rearranged into a vector of linear ODEs:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M38','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M38">View MathML</a>

(20)

Where

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M39','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M39">View MathML</a>

In order to find the solution of equation (20), the whole time domain is divided into a number of non-overlapping time intervals [ti-1, ti], i = 1,..., k. The sensitivity equations in equation (20) are transformed to algebraic equations using modified collocation method with piecewise linear polynomials as shape functions for each subinterval [tj-1, tj], ti-1 tj-1 <tj ti:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M40','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M40">View MathML</a>

(21)

where ηj is the step size in time tj. This equation can be rewritten as:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M41','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M41">View MathML</a>

(22)

where the constant vector

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M42','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M42">View MathML</a>

The value of w(tj) can be evaluated if the values of x(tj) have been computed when solving equation (8) by the ODE solver. The value of the M(tj) in equation (22) must be calculated before the solution of equation (21) can be obtained by an iteration method. When the value of the Jacobian matrix ∂f/∂x in time tj is known, the value of the M(tj) can be obtained straightforwardly. There is no requirement for computing the value of the Jacobian matrix here due to it has been computed for step size control when solving equation (8) by the ODE solver.

Equation (22) is similar to equation (10), we can apply the fixed point theorem to ϕ = h(ϕ) to get the maximum ηj satisfying the following condition

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M43','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M43">View MathML</a>

According to the definition of the matrix norm, it is easy to verify that ||M||p = ||∂f/∂x||p and we have

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/S12/S17/mathml/M44">View MathML</a>

(23)

Equation (23) is the same as equation (14), so the same criterion is used to determine the step size ηj for computing the time course and sensitivity profile with the guarantee of convergence and existence of a fixed point. The dynamic sensitivities can be obtained directly by solving equation (21) after the time profile of x has been computed and the step size ηj has been decided by the ODE solver.

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

WHW developed and implemented the algorithm and drafted the manuscript. FSW conceived of the study, participated in its design and coordination, and helped to draft the manuscript. MSC assisted in developing the algorithm and finalizing the manuscript. All authors read and approved the final manuscript.

Acknowledgements

The financial support from the National Science Council, Taiwan, ROC (Grant NSC97-2221-E-194-010-MY3 and NSC97-2627-B-194-001), is highly appreciated.

This article has been published as part of BMC Bioinformatics Volume 9 Supplement 12, 2008: Asia Pacific Bioinformatics Network (APBioNet) Seventh International Conference on Bioinformatics (InCoB2008). The full contents of the supplement are available online at http://www.biomedcentral.com/1471-2105/9?issue=S12.

References

  1. Albrecht P: A New Theoretical Approach to Runge-Kutta Methods.

    SIAM Journal on Numerical Analysis 1987, 24:391-406. Publisher Full Text OpenURL

  2. Albrecht P: The Runge-Kutta Theory in a Nutshell.

    SIAM Journal on Numerical Analysis 1996, 33:1712-1735. Publisher Full Text OpenURL

  3. Inc M, Bildik N, Bulut H: A comparison of numerical ODE solvers based on Euler methods.

    Math Comput Appl 1998, 3:153-159. OpenURL

  4. Gerald CF, Wheatley PO: Applied numerical analysis. Addison-Wesley; 1994. OpenURL

  5. Villadsen JV, Stewart WE: Solution of boundary-value problems by orthogonal collocation.

    Chem Eng Sci 1967, 22:1483-1501. Publisher Full Text OpenURL

  6. Wang FS: A modified collocation method for solving differential-algebraic equations.

    Applied Mathematics and Computation 2000, 116:257-278. Publisher Full Text OpenURL

  7. Savageau MA: Parameter sensitivity as a criterion for evaluating and comparing the performance of biochemical systems.

    Nature 1971, 229:542-544. PubMed Abstract | Publisher Full Text OpenURL

  8. Savageau MA: The behavior of intact biochemical control systems.

    Current Topics in Cellular Regulation 1972, 6:63-129. OpenURL

  9. Voit EO: Computational analysis of biochemical systems: a practical guide for biochemists and molecular biologists. Cambridge, UK: Cambridge University Press; 2000. OpenURL

  10. Savageau M: Biochemical Systems Analysis. A Study of Function and Design in Molecular Biology. Reading, MA: Addison-Wesley; 1976. OpenURL

  11. Kacser H, Burns JA: The control of flux.

    Symp Soc Exp Biol 1973, 27:65-104. PubMed Abstract OpenURL

  12. Fell DA: Metabolic control analysis: a survey of its theoretical and experimental development.

    Biochem J 1992, 286:313-330. PubMed Abstract | PubMed Central Full Text OpenURL

  13. Heinrich R, Schuster S: The Regulation of Cellular Systems. New York: Chapman & Hall; 1996. OpenURL

  14. Caracotsios M, Stewart WE: Sensitivity Analysis of Initial Value Problems with Mixed Ode's and Algebraic Equations.

    Comput Chem Eng 1985, 9:359-365. Publisher Full Text OpenURL

  15. Dunker AM: The decoupled direct method for calculating sensitivities coefficients in chemical kinetics.

    J Chem Phys 1984, 81:2385-2393. Publisher Full Text OpenURL

  16. Dougherty EP, Hwang JT, Rabitz H: Further developments and applications of the Green's function method of sensitivity analysis in chemical kinetics.

    J Chem Phys 1979, 71:1794-1808. Publisher Full Text OpenURL

  17. Dickinson RP, Gelinas RJ: Sensitivity analysis of ordinary differential equation systems. A direct method.

    J Comput Phys 1976, 21:123-143. Publisher Full Text OpenURL

  18. Leis JR, Kramer MA: The simultaneous solution and sensitivity analysis of systems described by ordinary differential equations.

    ACM Trans Math Softw 1988, 14:45-60. Publisher Full Text OpenURL

  19. Mauch K, Arnold S, Reuss M: Dynamic sensitivity analysis for metabolic systems.

    Chemical Engineering Science 1997, 52:2589-2598. Publisher Full Text OpenURL

  20. Ghosh A, Miller D, Zou R, Pais H, Sokhansanj B, Kriete A: Integrated Spatio-temporal Model of Cell Signaling.

    FOSBE 2005. OpenURL

  21. Hwang JT, Dougherty EP, Rabitz S, Rabitz H: The Green's function method of sensitivity analysis in chemical kinetics.

    J Chem Phys 1978, 69:5180-5191. Publisher Full Text OpenURL

  22. Ingalls BP, Sauro HM: Sensitivity analysis of stoichiometric networks: an extension of metabolic control analysis to non-steady state trajectories.

    Journal of Theoretical Biology 2003, 222:23-36. PubMed Abstract | Publisher Full Text OpenURL

  23. Shen J: A direct method of calculating sensitivity coefficients of chemical kinetics.

    J Chem Phys 1999, 111:7209-7214. Publisher Full Text OpenURL

  24. Zou R, Ghosh A: Automated sensitivity analysis of stiff biochemical systems using a fourth-order adaptive step size Rosenbrok integration method.

    AIEE Proc Sys Biol 2006, 153:79-90. Publisher Full Text OpenURL

  25. Shiraishi F, Takeuchi H, Hasegawa T, Nagasue H: A Taylor-series solution in Cartesian space to GMA-system equations and its application to initial-value problems.

    Applied Mathematics and Computation 2002, 127:103-123. Publisher Full Text OpenURL

  26. Shiraishi F, Hatoh Y, Irie T: An efficient method for calculation of dynamic logrithmic gains in biochemical systems theory.

    Journal of Theoretical Biology 2005, 234:79-85. PubMed Abstract | Publisher Full Text OpenURL

  27. Gorbunov VK: The parameterization method for optimal control problems.

    Comput Math Math Phys 1979, 19:212-224. OpenURL

  28. Gorbunov VK, Lutoshkin IV: The parameterization method in optimal control problems and differential-algebraic equations.

    Journal of Computational and Applied Mathematics 2006, 185:377-390. Publisher Full Text OpenURL

  29. Li R, Teo KL, Wong KH, Duan GR: Control parameterization enhancing transform for optimal control of switched systems.

    Mathematical and Computer Modelling 2006, 43:1393-1403. Publisher Full Text OpenURL

  30. Teo KL, Jennings LS, Lee HW, Rehbock V: The control parameterization enhancing transform for constrained optimal control problems.

    J Austral Math Soc Ser B 1999, 40:314-335. OpenURL

  31. Wang FS, Su TL, Jang HJ: Hybrid differential evolution for problems of kinetic parameter estimation and dynamic optimization of an ethanol fermentation process.

    Ind Eng Chem Res 2001, 40:2876-2885. Publisher Full Text OpenURL

  32. Kreyszig E: Introductory Functional Analysis with Applications. Wiley; 1989. OpenURL

  33. Varma A, Morbidelli M, Wu H: Parameter sensitivity in chemical systems. Cambridge: Cambridge University Press; 1999. OpenURL

  34. Tomovic R, Vukobratovic M: General Sensitivity Theory. New York: American Elsevier; 1972. OpenURL