A specialized ODE integrator for the efficient computation of parameter sensitivities
© Gonnet et al.; licensee BioMed Central Ltd. 2012
Received: 27 April 2011
Accepted: 22 March 2012
Published: 20 May 2012
Skip to main content
© Gonnet et al.; licensee BioMed Central Ltd. 2012
Received: 27 April 2011
Accepted: 22 March 2012
Published: 20 May 2012
Dynamic mathematical models in the form of systems of ordinary differential equations (ODEs) play an important role in systems biology. For any sufficiently complex model, the speed and accuracy of solving the ODEs by numerical integration is critical. This applies especially to systems identification problems where the parameter sensitivities must be integrated alongside the system variables. Although several very good general purpose ODE solvers exist, few of them compute the parameter sensitivities automatically.
We present a novel integration algorithm that is based on second derivatives and contains other unique features such as improved error estimates. These features allow the integrator to take larger time steps than other methods. In practical applications, i.e. systems biology models of different sizes and behaviors, the method competes well with established integrators in solving the system equations, and it outperforms them significantly when local parameter sensitivities are evaluated. For ease-of-use, the solver is embedded in a framework that automatically generates the integrator input from an SBML description of the system of interest.
For future applications, comparatively ‘cheap’ parameter sensitivities will enable advances in solving large, otherwise computationally expensive parameter estimation and optimization problems. More generally, we argue that substantially better computational performance can be achieved by exploiting characteristics specific to the problem domain; elements of our methods such as the error estimation could find broader use in other, more general numerical algorithms.
In systems biology, mathematical models often take the form of system of ordinary differential equations (ODEs). These are approximations of the underlying mechanisms such as enzyme-catalyzed biochemical reactions that are applicable when molecule numbers are sufficiently high, and when the spatial distributions of components in a cell can be neglected. More specifically, ODE models consider the rate of change in a set of states (e.g. species concentrations) as a function of the system’s current state, its inputs, and its inherent kinetic parameters that capture, for instance, affinities of molecular interactions .
In contrast to systems modeling in domains such as physics, however, model parameters and initial conditions for systems biology models are often not known, or they can only be roughly approximated. As few kinetic parameters can be measured directly, parametric uncertainty often prevails . How the system variables depend on these system parameters can therefore be of interest, e.g. to help find parameters such that the simulated system matches some observed or desired behavior. Dependencies between system variables and parameters are captured by the local parameter sensitivities that describe to what extent the state of the system changes when parameter values are perturbed from a reference value. Formally, local parameter sensitivities comprise the set of derivatives of all system variables with respect to the system parameters. As with the state dynamics, the parameter sensitivities’ time evolution follows a system of ODEs .
From a computational point of view it is important to note that in all but the simplest cases, starting from a set of initial conditions, there is no direct way to compute the solutions of a system of ODEs (states or parameter sensitivities in our case) for an arbitrary time. The variables are therefore integrated numerically in small steps over time, until the desired end time is reached. Consequently, efficient and accurate numerical integration methods are critical for many applications.
The computational effort for numerical integration is linked to the system size, and over time mathematical models have become increasingly detailed to achieve better predictions. Nevertheless, even models of moderate complexity result in numerical challenges when parameter sensitivities are needed. For instance, the parameter sensitivities can be integrated naively alongside the system variables, but this implies integrating a system of size n x ×(1 + n p ), where n x and n p are the number of system variables and system parameters respectively .
Additionally, the solution of a system of ODEs is often used in system identification processes where global optimization or probabilistic inference are required [4, 5]. In such cases, thousands, if not millions of trajectories need to be computed. Assessing the quality of the identified model, for instance with respect to the uncertainty in parameter values, again requires computing the local parameter sensitivities . Although local sensitivity information can often help improve the overall estimation process, sensitivity computations are rarely included for performance reasons. Specific efficient methods exist for cases in which only scalar valued functionals are optimized  or oscillatory systems are considered . Yet in many other cases, such as optimal control [9, 10], the identification of relevant parameters , model reduction and simplification  or parameter training , the full parameter sensitivities need to be computed. Consequently, improvements with respect to the speed with which the original ODE systems and their parameter sensitivities can be reliably integrated may affect the entire process significantly.
These issues are not new and they concern many application domains. Considerable efforts have been invested in establishing reliable and efficient general-purpose ODE solvers for dynamic systems and—to a lesser extent—for the associated parameter sensitivities. Here, however, we are concerned with solving systems of ODEs as they typically occur in the simulation of biochemical reaction networks in systems biology . We show that rather general characteristics of such systems allow for the development of application domain-oriented ODE solvers with novel numerical features (with potential broader applicability), and with superior performance compared to state-of-the-art, widely employed general-purpose solvers. To provide some context for this claim, we first briefly review key characteristics of systems biology models in the form of ODEs, and general methods for the numerical integration of ODEs.
where f(x(t),p) is a system of functions f i (x(t),p) modelling the conversion rate of each respective variable x i (t,p) at time t, and p is a vector of n p system parameters.
where J f (x(t)) is the Jacobian matrix of f(x(t)) with respect to x(t); note that we drop explicit dependencies on p to simplify notation. Initial conditions for Eq. (3) are set according to whether the initial conditions for the states in Eq. (1) depend on the parameters or not .
Such problems are often well solved by general purpose ODE solvers, but (bio)chemical reaction networks offer a number of features that may be exploited by more specialized solvers, resulting in faster and/or more precise simulations. For instance, in enzyme kinetics, reversible association and dissociation processes are usually much faster than product formation. The resulting stiffness severely limits the types of numerical methods that can be used for ODE integration.
with closed and open circles indicating non-zero and zero elements, respectively. Even in this dense sub-network, the number of non-zeros implies that we do not need to compute a substantial number of terms to determine the Jacobian.
Many large-scale biological networks have a scale-free structure, that is, most of their nodes have few interactions, but a small number of hubs with many interactions exist . This prevents an easy decomposition of a large network into subsystems that can be handled (and integrated) independently. Therefore, despite the sparsity of the Jacobian, model size remains a major issue for numerical performance.
Two more general aspects also need to be considered. Firstly, due to the growing use of abstract modeling software, the reactions and the underlying reaction equations are usually available to us as abstract models, such as smbl, which we can analyze and manipulate analytically. Secondly, since parametric uncertainty is abundant in biology, sensitivity analysis, i.e. the integration of the parameter sensitivities s k (t), requires particular attention. Hence, an ideal ODE solver for our application domain would efficiently and reliably handle large, stiff dynamic systems including their parameter sensitivities, and optimally exploit the systems’ non-trivial sparsity and analytic access.
If the factors ∂ k x(t) / ∂ t k / k! decrease sufficiently quickly and the higher-order terms become insignificant as of some degree n, then we can reliably approximate the new solution by a polynomial of degree n in h. In explicit integrators, previously computed values of x(t) and the derivatives ∂ x(t) / ∂t are used to construct a polynomial g n (t) of degree n and to extrapolate the value of x(t + h) ≈ g n (t + h). In implicit integrators, a solution x(t + h) is sought such that it matches that of a polynomial g n (t) of degree n interpolated through previous values of x(t) and/or their derivatives, and the derivative at the solution x(t + h) itself. In general, implicit integrators are more accurate for stiff ODEs, where the derivatives in Eq. (5) do not decay sufficiently quickly.
Common ODE integration schemes and the values that are used to approximate the polynomial in Eq. (5)
Backward Differentiation Formula (BDF)
, x(t−h k )
x(t), , x(t−h k )
Second-derivative rule (this work)
x(t), , , ,
Despite the commensurate degree of freedom in designing ODE integrators, and the number of algorithms for the numerical integration of ODEs that have been published over the past 40 years, only very few of them have found wide-spread application. Practical considerations—any method should be easily accessible to its end users, who are usually not interested in manipulating or even formulating the underlying equations themselves—are certainly major causes for this convergence . However, a closer analysis of the most popular solvers for stiff ODE systems reveals another cause, namely incremental evolution.
In this area, the first major piece of software was the GEAR package , which by 1996 evolved into cvode, a part of the Sundials suite of nonlinear and differential/algebraic equation solvers . The default integrators in Matlab (The MathWorks, Natick, MA) such as ode15s employ similar integration rules and error estimates. Both the Sundials suite and Matlab are used increasingly in systems biology , but it is not evident that they are optimal for this application domain.
(for notational simplicity, we will write J f (t) and f(t) instead of J f (x(t)) and f(x(t)), respectively). Note that this second derivative with respect to the time t should not be confused with the second-order sensitivities described and used in [9, 10, 23], which are the second derivatives of the system variables with respect to the system parameters.
The use of second derivatives was first suggested in  and later studied in detail in ,  and , and the resulting formulas were shown to have good stability properties. A more recent study  reinforces the stability and potential efficiency gains for stiff systems through second-derivative methods. However, despite several published implementations [27, 28], these methods have not yet found wide acceptance because, despite being A-stable, they are only stable at infinity if only the second derivative at t + h is used  (see Section S3 in Additional file 1 and Additional file 2 for details).
In most cases, the evaluation of the second derivatives is not much more expensive than the evaluation of f(t).
The Jacobians and are evaluated at the initial guess .
in the case of the BDF and the Adams-Moulton formula of degree four, respectively. These truncation errors are 72 times and 19 times larger than the error of our second-derivative formula (assuming the fifth derivative x (5)(t) is approximately constant in [t n−3,t n + 1], see Section S2 and Figure S1A in Additional file 1 and Additional file 2 for details). The large difference stems from the dependence of the interpolation error on the width of the interpolation interval, e.g. for the BDF and the Adams-Moulton formula, this interval is four times larger.
In any ODE integration scheme, the local error estimate and the step-size adjustment are crucial to both its accuracy and its efficiency. The step-size adjustment uses the error estimate of a previous integration step to predict the largest possible next step h satisfying the required tolerance. With imprecise error estimates, the step-size adjustment has to be conservative to preserve accuracy, or it risks producing an imprecise result.
In most implicit ODE solvers, the local error is either estimated from the difference between the initial estimate , usually computed with an explicit rule, and the final converged step x(t + h), or as the difference between two rules of different degree over the previous x(t) and the converged step x(t + h). These approaches mainly consider computational efficiency because, ideally, to estimate the error of a formula of degree d 1, we need to compute a better approximation of degree d 2 > d 1. The difference between both converged solutions x 1(t + h) and x 2(t + h) can then be used to approximate the difference between the lower-degree estimate and the exact solution x ⋆(t + h). However, this requires two Newton iterations to compute both solutions and if both rules have different weights for the values of and , we would need to invert or decompose two different matrices to compute a Newton iteration (Eq. (8)).
We propose a different approach that may better reconcile accuracy with computational cost. We first compute the converged lower-degree solution x 1(t + h) and use it as an initial estimate for the Newton iteration of the higher-degree solution. Since we are not actually interested in the exact solution x 2(t + h), but only in an approximation of the difference between the two solutions, it suffices to compute just one Newton step to get a first-oder approximation of that difference. Note that, in principle, this still requires the inversion or decomposition of a different matrix for the Newton iteration.
Note that if the assumptions on x (5)(t) do not hold, the error estimate in the next time step will fail, causing the step size to be reduced automatically. Furthermore, if we adjust h to fulfill the requested tolerance τ exactly, the error estimate will be larger than τ approximately half of the time. Therefore, in practice, we choose σ such that the next error will be τ/2. This gives us a recipe to adjust the step size from one integration step to the next and, hence, the last key ingredient of a functional second-derivative ODE solver (see Section S2 and Figure S1B of the Additional file 1 and Additional file 2 for details).
For n x variables and n p parameters, the naive approach to sensitivity calculation implies integrating a system of n x ×(1 + n p ) variables and, by consequence, inverting or decomposing matrices of that size within the Newton iteration. However, the system variables x(t) do not depend on the parameter sensitivities, yet the sensitivities depend on x(t). Hence, we can, in each step, first compute the values x(t + h) and, once they have converged, compute the s k (t + h) in a separate step using the same integration rule. This staggered approach was first introduced in Caracotsis & Stewart , then extended by Maly & Petzold , and finally implemented in the Sundials cvodes ODE solver, a modified version of cvode capable of sensitivity analysis .
We then have two alternatives to compute s k (t + h): either iteratively using Newton’s method to solve Eq. (14), or directly by inverting or decomposing the matrix on the left-hand side of Eq. (15). The iterative approach via Eq. (14) is equivalent to the one suggested in , and we can re-use the inverted or decomposed iteration matrix used to compute the variables x(t + h) in Eq. (8). However, one has to re-evaluate the Jacobians at each iteration to compute and because and are evaluated at the initial estimate , and not at the converged solution x(t + h). To compute s k (t + h) directly, which corresponds to the original approach in , we need to re-compute the Jacobians and the inverse or decomposition of the left-hand side of Eq. (15). For small n x , however, this extra matrix computation may be advantageous over running the Newton iteration for each parameter p k .
Furthermore, if the Jacobians do not vary significantly over time, they can be re-used as the matrices and for the Newton iteration of the next step. Such an approach offers an advantage if the cost of running an additional n p Newton iterations to compute the parameter sensitivities iteratively outweighs the cost incurred by the slower convergence due to using older Jacobians in the next step. In our second-derivative integrator, we therefore compute the parameter sensitivities directly as per Eq. (15).
In order to generate the matrices J f (t) and J Jf (t), as well as the second derivative automatically, we established a framework that automatically translates arbitrary models from the standard SBML format  to Matlab functions or C-language code. The framework also generates routines to compute the parameter derivatives ∂ f(t)/∂ p and necessary for the parameter sensitivity computations. This conversion, which needs to be done only once per model, exploits the sparsity of the corresponding matrices by generating compact expressions for their non-zero entries only, making them efficient to evaluate. It uses the Matlab Symbolic Toolbox to manipulate, differentiate and simplify the resulting expressions automatically (see Sections S1.3 and S1.4 in the Additional file 1 and Additional file 2 for details).
We implemented the second-derivative ODE integrator as odeSD in Matlab and as odeSD mex in the C programming language, using the Matlab mex interface with calls to the LAPACK and BLAS libraries for the linear algebra operations. Both solvers provide an interface similar to that of the Matlab default integrators. Additionally, a native C-language version, odeSD c , was implemented for use outside of the Matlab programming environment. All implementations could operate on any type of ODE-based model, but the overall implementation is targeted to systems biology models in standard SBML format, for which we developed an automatic model conversion framework (see Methods). The implementation details are described in Sections S1.1 and S1.2 of the Additional file 1 and Additional file 2.
Matlab’s default integrator for stiff systems, ode15s, which uses a 5-point Numerical Differentiation Formula (NDF), a more stable variant of the BDF integration rules; it is used as the default integrator in SBToolbox2 ,
The Matlab interface supplied by the sundialsTB toolbox  served to run the Sundials integrators which are implemented in C.
Systems biology test models and their key characteristics, namely number of states (n x ), number of parameters (n p ), number of non-zeros of the Jacobians nnz(J f ), integration time interval (t), and biological system described by the model
nnz(J f )
nnz(J Jf )
Hornberg et al.
Steady state, ERK1 phosphorylation and kinase/phosphatase control.
Kholodenko et al.
Steady state, short term signaling by the EGF receptor.
Singh et al.
Steady state, IL-6 signal transduction in hepatocytes.
Borisov et al.
Steady state, insulin-EGF network interactions in mitogenic signaling.
Ung et al.
Steady state, regulation of EGFR endocytosis and EGFR-ERK signaling by crosstalk.
Elowitz & Leibler 
Oscillatory, synthetic network of transcriptional regulators.
Leloup & Goldbeter 
Oscillatory, circadian oscillations of PER and TIM proteins in Drosophila.
Wolf et al.
Oscillatory, autonomous metabolic oscillations in continuous culture of S. cerevisiae.
Goldbeter & Pourquié 
Oscillatory, segmentation clock by crosstalk of Notch, Wnt, and FGF pathways.
Xie & Kulasiri 
Oscillatory, circadian rhythms in Drosophila with interlocked feedback loops.
The performance comparison with sensitivity calculations requires two additional considerations: Since ode15s and radau5 do not provide any special functionality for computing parameter sensitivities, we used an augmented system of size n x ×(n p + 1) including an analytic sparse Jacobian for each model, and the sensitivities s k (t), k=1…m were integrated alongside the system variables. cvodes, the sensitivity analysis-enabled version of cvode from the sundials package, uses a simultaneous integrator based on the method of Maly & Petzold . Optionally, the staggered integrator of Feehery et al. can be selected, but this did not produce better results.
In all cases, parameter sensitivities were integrated to the same precision as the system variables. As with the integration without sensitivities, precision/work diagrams were computed for all models with sensitivities, omitting the cases in which ode15s failed completely.
To explain the performance, consider that although odeSD and odeSD mex also evaluate the Jacobians in each step, both the much smaller number of steps required (Figure 3D), which is of no particular advantage when the parameter sensitivities are not computed, and the re-use of the Jacobians for sensitivity computations lead to significantly shorter execution times since substantially fewer evaluations of f(·) and of the Jacobians J (·)(·) are needed (Figure 3D). Since the latter dominates the integration cost, all three versions of odeSD outperform cvodes in all but the smallest systems. These results are not a consequence of the sparsity of the systems per se, but of the more precise integration rule which can be computed efficiently thanks to sparsity. The better error estimate and the computation of local parametric sensitivities are therefore particular strengths of our ODE solver based on second derivatives.
We have presented an integrator for ODE systems resulting from the modeling of chemical and biological reaction networks, which are often stiff and sparse. For the realistic systems biology models tested, the new integrator outperforms commonly used state of the art integrators when parameter sensitivities are required. It is competitive in integrating the system equations alone, despite limitations for specific models near the steady state. The improvements with respect to sensitivity calculations are critical for many applications to drive highly compute-intensive (global) optimization and estimation processes.
The improvements themselves are due to a combination of several factors: The more accurate second-derivative rule allows us, in combination with a better error estimate, to take larger steps, which in turn allows us to reduce the number of otherwise expensive sensitivity calculations. The re-use of the Jacobians from the sensitivity calculations further reduces the total computational costs. Although each integration step is more expensive than in first-derivative methods, due to the additional second-derivative information that needs to be computed, far less steps are required in total, resulting in a more efficient method.
To be of practical relevance for applications in systems biology, odeSD and odeSD mex are accessible via Matlab interfacesa, and we plan to make them more easily available through integrated modeling environments such as the SBToolbox2  and COPASI , e.g. via the native C-language interface. To accelerate larger optimization and estimation processes, further efficiency improvements through the use of sparse matrix routines and by more elaborate step size control schemes are possible.
In terms of numerical algorithms, to our knowledge, this is the first practical application of a second-derivative integration method with good performance. Key, novel features of our integrator, such as a more precise error estimator and direct computation of parameter sensitivities with re-use of the Jacobians, may well be suited for other problems or types of integrators. Importantly, while our integrator has been developed for (bio)chemical and reaction networks, it is still quite general in targeting stiff and sparse ODE systems. Overall, we feel that a lot is to be gained by adapting general algorithms to specific problem domains, and that results from the work on specific problems can spill over to the broader field.
The integrator (Matlab and C-language versions) and the model conversion framework are available via http://www.csb.ethz.ch/tools.
aWe note that it is not possible to execute odeSD mex in the 64-bit Windows version of Matlab R2010a as well as in 64-bit Linux versions prior to R2010a in our testing environments, for reasons of memory allocation problems in Matlab.
We thank Ernst Hairer for comments and discussions. This work was supported by the EU FP6 project BaSysBio (LSHG-CT-2006-037469), the EU FP7 project UniCellSys (HEALTH-F4-2008-201142), the Swiss Initiative for Systems Biology SystemsX.ch (RTD project YeastX), by an ETH Excellence Scholarship (LW), and by Swiss National Science Foundation’s Individual Support Fellowships Nr. PBEZP2-127959 and Nr. PA00P2-134146 (PG).
This article is published under license to 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.