A specialized ODE integrator for the efficient computation of parameter sensitivities

Background 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. Results 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. Conclusions 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.

To achieve reliable timing data, all models were integrated 10 times. Evaluation numbers for derivatives and Jacobians were obtained by designing appropriate function wrappers. Numerical precision was evaluated as the maximal relative deviation between a reference solution (obtained with radau5 at relative numerical tolerance of 10 −15 ) and the approximation with variable tolerance at the end time of integration over all system states. States with values below machine precision were excluded from this analysis. To determine the average performance as a function of numerical precision, precisions were binned (with centers as shown in the corresponding figures), and performance indicators were averaged over these intervals and over all models. For the automatic function generation from SBML models, we used Symbolic Toolbox V5.2 and the current version of SBtoolbox2 that supports all features of SBML Level 1 and 2.

S1.1 Details of the Matlab Integrator
To retain a certain level of compatibility, odeSD uses the same calling sequence as the default integrators in Matlab: [ T, X, TE, YE, IE, S, SE ] = odeSD ( f, t, x0, opts, params, varargin{:} ) where the return value T is a vector containing the time steps and X is a matrix whose columns are the variables at each time step. The odeSD mex integrator uses the same interface, yet is called using odeSD wrapper instead of odeSD. The return value S contains a three-dimensional array containing the values of the parameter sensitivities at each time step. The return values TE, YE, IE and SE contain, analogously to the Matlab integrators, the data at points where user-specified events (see the Events option below) were triggered. The parameter sensitivities are only computed if the return value S is requested. The parameter t is a vector of length at least two containing the time interval over which to integrate. If t is of length larger than two, the values of x(t) are only stored and returned for those times. Otherwise, all computed values of x(t) are returned. The optional parameter opts contains an ODE-options structure as is available for the other Matlab integrators in which the options RelTol, AbsTol, NormControl, NonNegative, OutputFcn, OutputSel, Refine, Stats, InitialStep, MaxStep, Events and Jacobian can be set. The optional parameter params, which is required when the output S is requested, contains a vector of parameter values to be passed to the right-hand side f for sensitivity analysis and the variable parameters varargin{:} are passed to the right-hand side f.
The right-hand side f is expected to be either a function or a handle on a function of the form where dxdt and d2xdt2 are the first and second derivativesẋ(t) andẍ(t) respectively. The output arguments J f and J Jf are the Jacobian matrices of the first and second derivatives, respectively. The outputs dfdp and d2fdpdt, which are only required if the parameter sensitivities need to be computed, are matrices containing the derivatives of the first two outputs with respect to the parameters where the kth column of dfdp contains the derivatives of the variables with respect to the kth parameter. Similarly, the function specified with the Jacobian option in opts should have the form where the first derivative dxdt is supplied as it may simplify the computation of the Jacobians. If no initial step size is specified in opts, a default value of 1/100th of the integration interval is assumed, which is repeatedly scaled by a factor of 0.7 until the extrapolated initial guessx(t 0 + h) = x 0 +hf (x 0 ) +h 2 /2J(x 0 )f (x 0 ) satisfies all non-negativity constraints specified by the option NonNegative.
In every step, the algorithm computes an initial guessx(t + h) by constructing a BDF over the last three points and extrapolating to the new time t + h. In case no previous steps are available, a linear extrapolation using the zeroth and first derivative at t is used, and if only one step is available, the zeroth and first derivatives of the last two steps are employed.
If the step size has changed from the previous step, the iteration matrix in (8) is reassembled using the most recent copies of the JacobiansJ f (t + h,x(t + h)) andJ Jf (t + h,x(t + h)) and is decomposed using an LU factorization. If after two iterations of Newton's method the solution diverges or is not expected to converge within five steps (see [2] for details on how convergence is estimated), the Jacobians are re-evaluated, the iteration matrix is reassembled and decomposed, and the iteration is restarted. If the Jacobians are up to date, the iteration is abandoned and the step size h is reduced by a factor of 0.7.
Once the Newton iteration has converged, the error is approximated as in (11). The maximum relative component-wise error is used to compute the scaling factor σ for the next step as in (12) with a tolerance τ of half of the requested tolerance. Note that since the scaling σ is computed to achieve half of the required tolerance, a step only fails if σ < (1/2) 1/5 ≈ 0.87. If the error estimate is below the requested tolerance, the algorithm proceeds to compute, if requested, the parameter sensitivities, as per (15). For the parameter sensitivities, the same error estimator and step size scaling are applied as for the system variables.
The fifth-degree rule used in the error estimate of both the system variables and the parameter sensitivities is where h −1 is the size of the previous step.

S1.2 Details of the C-language integrator
The interface of the native c-language odeSD c is similar to that of the Matlab version. The calling sequence, as defined in odeSD.h is int odeSD ( f_t *f , dfdx_t *dfdx , int nr_tspan , double *tspan , int N , const double *x0 , struct odeSD_opts *opts , int nr_params , const double *params , void *varargin , double **t_out , double **x_out , double **s_out ); where f t is a function of the type and dfdx t is a function of the type int dfdx ( double t , const double *x , const double *f , const double *p , void *varargin , double *J , double *dJ , double *dfdp , double *d2fdtdp ); The functions f and dfdx compute the first derivatives and Jacobians, respectively. The resulting vectors and matrices are stored in column-major order in the output variables f (first derivative), dfdt (second derivative), J (Jacobian of first derivative), dJ (Jacobian of second derivative), dfdp (derivative of f with respect to the parameters) and dfdp (derivative of dfdt with respect to the parameters). These variables point to memory allocated by odeSD and are NULL if the value is not required. The values x and p contain the N system variables and nr params parameters respectively. The variable varargin passed to odeSD is passed on to f and dfdx.
The output variables t out, x out and s out are analogous to the return values of the Matlab integrator and will contain pointers to these arrays in column-major ordering, allocated by odeSD using malloc.
The argument opts is a pointer to a structure of the type odeSD opts defined in odeSD.h which contains options analogous to those of odeset in Matlab. The global variable odeSD opts default contains the default settings. The arguments x0 and params contain the initial system variables and the parameters respectively and tspan is a pointer to a vector containing nr tspan time steps at which to evaluate the system variables and parameter sensitivities. If nr tpsan is two, the outputs are stored for each computed step.
The function odeSD returns the number of output values stored or any value < 0 on error. A stack of any errors can be displayed using errs dump( FILE *out ), defined in errs.h.

S1.3 Automatic Generation of Functions
An important part of this work was the extraction of the necessary mathematical information used by the solvers from each model in SBML representation [3]. To this end, we developed a Matlab (MathWorks, Nantucket / MA) interface which, given an SBML model, fully automatically generates a series of suitable Matlab files that contain both the analytical representations of the system's derivatives f (t) andẍ(t) as well as the Jacobians J f (t), J Jf (t), ∂f (t)/∂p, and ∂ẍ(t)/∂p.
The entries of the derivativeẍ(t) and of the matrix ∂ẍ(t)/∂p are generated explicitly, i.e. without evaluating J f (t)f (t). The matrix J Jf (t) is constructed as per (9), where the entries of (∂J f (t)/∂x)f (t) are computed explicitly and added to the square of the previously computed J f (t).
The framework makes use of SBToolbox2 [4] for reading the initial SBML model and of Matlab's Symbolic Toolbox for performing various differentiations on the system's equations. The latter was preferred over parsing and differentiating the expressions in Matlab, as the Symbolic Toolbox provides the ability to simplify the sometimes clumsy or redundant expressions resulting from a straight-forward differentiation. The generated function files are written in both standard Matlab and c-language formats and can be used, in principle, with appropriate wrappers, in any solver. For example, given an SBML model in Matlab's current directory, we can write: modelInfo = xml2odefun(modelname,{'c files'}).
All the necessary function files in c-language format that describe the system will be generated in the current directory, and a structure providing useful information of the system, modelInfo, is returned.
Producing the system description files for high-dimensional (both in terms of state and parametric space) systems is a non-trivial task because the computational cost grows quadratically with the system dimension. In a fully interconnected system we will never be able to avoid such a computational cost. Remember, however, that biological systems are generally poorly interconnected, resulting in sparse Jacobians. Parsing the system's equations prior to differentiation allows us to pinpoint the elements of the various Jacobians that are non-zero and perform all computations only on these elements. In all the systems under study, the computational cost grows only linearly with the system size.

S1.4 Framework Usage and Example
The following short Matlab code shows how to construct the integrator input files using our framework: We first pass the SBML model through the converter function xml2odefun (line 3) such that the necessary Matlab files, named MyModel.m and jac MyModel.m, analogously to the name of the SBML file, are generated. In line 12, we extract the initial conditions for both the variables x0 and the parameters p0 by calling the right hand side without any arguments. The Jacobian generated by the converter is then added to an options structure (line 18) so that it can be passed to the integrator (line 21) odeSD, which then computes the integration and sensitivity analysis. In lines 24 and 27 we plot the system variables and the parameter sensitivities of the third variable respectively.
Analogously, the c-language files generated by the converter when called with the 'c files' option, can be passed to odeSD c as follows 1 /* S t a n d a r d headers .

S2 Accuracy of the second derivative rule and of the error estimate
To illustrate the effect of the interval width, we determined the truncation error in each step for our second-derivative rule as well as a BDF and Adams-Bashforth rule of the same degree when integrating the simple oscillator We employed a fixed step size of h = 0.1 and used the exact values of x(t) = (sin t, cos t) T for the previous steps in each case. As predicted, the resulting error of the second-derivative method for this test case is at least one order of magnitude lower compared to the other methods (Supporting Fig. S1A). Hence, we expect a substantial improvement of accuracy or efficiency compared to standard first-derivative methods.
Using the same test function (1) as before, we compared our new error estimate with the-computationally more efficient (O(n) vs. O(n 2 ))-estimate based on the difference of g 4 (t + h) and g 5 (t + h) over the converged solution x 1 (t + h). The ratio of predicted to actual errors shown in Supporting Fig. S1B indicates a superior accuracy of our estimate. Furthermore, the almost uniform spacing of the errors along the x-axis in Figure S4 is a good indication of the robustness of the error estimate.

S3 Results without parameter sensitivities
The results show that our integrator requires less steps than the first-derivative methods. The smaller number of steps, however, does not always translate into a speed advantage, especially for the larger models with steady-state behavior (see Fig. 1 and Supporting Fig. S2 for details). On average, the second-derivative integrator requires an approximately equal number of function evaluations (Fig. 1C), but a much larger number of evaluations of the Jacobians (Fig. 1D). This effect is a result of the instability of the second-degree rules at infinity (see Section "Methods: A second-derivative integrator") and was already observed in [5] and in [6] when studying the "rober" model therein. It can be understood by considering the explicit formulation of the Newton iteration matrix M(t + h) (8) in terms of the Jacobian J Jf (9). The matrixJ f (t + h) is not the exact Jacobian and whatever perturbation it contains will be amplified in (J f (t + h)) 2 , causing the Newton iteration to converge less often than in the first-derivative rule, i.e. the step size is limited by variation in the second-derivative Jacobian J Jf .
This, however, is only a serious disadvantage for the models with a large number of variables and it is no disadvantage at all if the Jacobians need to be evaluated at every step, e.g. for very stiff systems or when computing parameter sensitivities. Integration of such models is not infeasible, but inefficient when approaching the steady state.   Figure S2: Detailed performance comparison without parameter sensitivities. All models were integrated for the time spans shown in Table 2. (A-H) Computation times for the individual models listed in Table 2 with varying relative tolerances (see X-axis labels) using odeSD (white bars), odeSD mex (black), ode15s (red), radau5 (green), and cvodes (blue).  Figure S3: Detailed performance comparison with parameter sensitivities. All models were integrated for the time spans shown in Table 2. (A-H) Computation times for the individual models listed in Table 2 with varying relative tolerances (see X-axis labels) using odeSD (white bars), odeSD mex (black), ode15s (red), radau5 (green), and cvodes (blue).  Table 2.