Dynamic sensitivity analysis of biological systems

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][2][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 secondorder 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][8][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][12][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 cal-culation 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][15][16][17][18][19][20][21][22][23][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][29][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 t i to t i + η 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 t i + η and the same step size are used to propagate the sensitivity equations from t i to t i + η. 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: where [x] is the concentration of species x and k i is the rate constant. The initial concentration of C 2 H 6 is 5.951 × 10 -6 mol/cm 3 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  [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%.

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

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: where x is the concentration of cell mass, s is the concentration of glucose, p 1 is the concentration of ethanol, p 2 is the concentration of glycerol, V is the working volume of the fermenter, t f is the final fermentation time,

Reaction
Rate constants The initial concentration of C 2 H 6 is 5.951 × 10 -6 mol/cm 3 and all other initial concentrations are zeros. Units for rate constants are mol, cm, s and the temperature is 923 K.  normalized fermentation time, s F is the feed concentration of glucose, F is the feed rate, is the ethanol yield factor, and is the glycerol yield factor. The unstructured kinetic models for the specific cell growth and product formation are respectively expressed as follows: 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 = p 1 V/t f under some physical constraints, e.g., the residual glucose restriction s(t f ) ≤ s r for reducing the separation cost, s r 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: Figure 1 shows the computational time profile of the ethanol fed-batch fermentation model with the optimal feed rate.
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 s F , F and t f ) 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 , , , and s F . 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)   . 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.
The relative sensitivity with respect to t f is shown in Figure  3(a). As expected, an increase in t f 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).

J, it is clear that an increase in
or s F will have more impact than an equal relative increase in or μ m . The negative value of relative sensitivity for J with respect to t f 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 s F is higher than that with respect to at the final fermentation time, by increasing s F 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 conclu-sion that to increase J by increasing will be a better choice than increasing s F or , and decreasing the fermentation time.
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 F i , i = 1,..., 10, are shown  Figure 5. The effects on J are in decreasing order from F 1 to F 5 . 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.

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.  to t f , x(0), s(0), and V(0). (a) relative sensitivities with respect to t f ; (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/t f ).
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: 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 var-iables, N ∈ ‫ޒ‬ (n+m) × q is the stoichiometric matrix describing the interconnecting fluxes, x 0 and y 0 are initial concentrations of x and y respectively.
,..., , 1 where x j ∈ x for j ≤ n, x j ∈ y for j > n, γ i is the rate constant, and g ij is the kinetic order for each x j . Equation (2) can be expressed concisely as: 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: we let x n+1 = t and dx n+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 = {t i |i = 1,..., k}. An ODE solver is to find the value of x(t i ), t i ∈ 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 nonoverlapping intervals [t i-1 , t i ], i = 1,..., k. The unified formulas of the modified collocation method for the subinterval [t j-1 , t j ], t i-1 ≤ t j-1 <t j ≤ t i , can be expressed as where η j is the step size in time t j , Î is an identity-like matrix, and the coefficient matrices D and 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 piece-wise linear polynomials transformed from equation (5) for each subinterval [t j-1 , t j ], t i-1 ≤ t j-1 <t j ≤ t i have the same formulas as the modified Euler method: 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 [t i-1 , t i ], 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 (x n ) 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(x m , x n ) 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. [32] Let (X, d) be a non-empty complete metric space. A mapping ψ:

Theorem 4. Banach fixed point theorem
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). D 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 where c = x(t j-1 )+0.5η j f(x(t j-1 ), y(t j-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(t j ) is obtained using the iterative rule. If x i (t j ) tends to a limit x*(t j ) 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 x i (t j ), i = 1,..., ∞ and d(a, b) = ||a -b|| p where a and b are arbitrary x i (t j ) 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.
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 Comparing equations (11) with (12), we obtain 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 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 t j . 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: where N is the stoichiometric matrix, v i is the i th element of v ∈ ‫ޒ‬ q and g ij is the kinetic order for each x j ∈ x in v i . 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 3. An order set T = {t 1 ,..., t k } of sampling points, t i , 1 ≤ i ≤ k is the sampling time of the solution of each ODE, k is the number of sampling points.

A tolerance ε
Output: The set of solutions of dependent variables at each sampling time.
• For each sampling time t i in T.
2. Repeat the following steps until d t = t i -t i-1 .  (e) If μ * η j > 1, then update η j with 0.9/μ (f) Call iteration algorithm to compute the value of x c stepped forward η j from x p .
(g) If the iteration algorithm succeeds in computing x c , then d t ← d t + η j and η j ← t i -t i-1 -d t , otherwise exist this algorithm. 3.

Algorithm Iteration
Input: 1. A set of n ordinary differential equations = f(x, y).
4. Repeat the following steps until the iteration limitation is reached or the value of x(t + η t ) converges.
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(x i , θ j ) of dependent variable x i ∈ x with respect to parameter θ j ∈ θ is defined as where x i (t; θ j + Δθ j ) is the i th component of the solution of equation (5) with an increment Δθ j on the j th parameter.
The function x i (t; θ j + Δθ j ) can be expanded into a Taylor series as follows: where 0 <ξ < 1. If Δθ j is sufficiently small, the last term of equation (16) can be truncated, leading to a linear approximation of x i (t; θ j + Δθ j ). Substituting the linear approximation of x i (t; θ j + Δθ j ) into equation (15) leads to This is defined as the first-order local sensitivity of x i with respect to θ j [33]. The relative parameter sensitivity S(x i , θ j ) of x i with respect to θ j is defined as Similar to the parameter sensitivity, the absolute log gain l(x i , y j ) and log gain L(x i , y j ) of x i ∈ x with respect to y j ∈ y are expressed respectively as follows: 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 x i with respect to θ j is given as where f i is the i th element of f [34]. The absolute dynamic log gain of x i with respect to y j is similar to equation (17) by replacing θ j , s(x i , θ j ) with y j , l(x i , y j ) respectively when y j (t) is a constant. In the case which y j (t) is a continuous time-dependent function, the whole time domain of y j (t) is partitioned into N u time intervals (t k-1 , t k ), k = 1,..., N u . Function y j (t) is parameterized by the piecewise constant functions ω k (t), k = 1,..., N u , as follows: where u k , k = 1,..., N u , are constant input control parameters and The continuous time-dependent function y j (t) is approximated by equation (18). The dynamic log gain l(x i , y j ) of infinite dimension is transferred into N u dynamic log gains of dimension one with respect to the input control parameters, and expressed as.
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 N r 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 + N r dimensional vector which contains model parameters θ, initial conditions of dependent variables x 0 and constant input in r. Φ indicates a matrix of size n × (n + p + N r ) 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: The model equations are rewritten as The sensitivity equations in matrix form can be derived by applying the chain rule to the derivative of Φ: where η j is the step size in time t j . This equation can be rewritten as: where the constant vector The value of w(t j ) can be evaluated if the values of x(t j ) have been computed when solving equation (8) by the ODE solver. The value of the M(t j ) 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 t j is known, the value of the M(t j ) 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 According to the definition of the matrix norm, it is easy to verify that ||M|| p = ||∂f/∂x|| p and we have 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.