Living biological systems are dynamic. By dynamic, we mean that the biological variables change with time. Quantitative description of such systems are often expressed in terms of Ordinary Differential Equations (ODEs) that relate the dynamic variables to their rates of change. These relations are mathematical models for the system. To obtain explicit time-dependent behavior of the dynamic variables predicted by a model, one needs to solve the ODEs.

XSIM allows a model to be written as a set of ODEs that are solved (i.e., integrated) by the XSIM software. See section 3.3 of the Model Writer Reference Manual for information about writing ODE models.

Two solvers for ODEs have been implemented under XSIM:

DOPRI5: An explicit fourth order Runge-Kutta method with fifth order correction (Dormand & Price) with step size control. This solver solves non-stiff and mildly-stiff differential equation.

RADAU: An implicit Runge-Kutta method (RADAU IIA) of variable order that switches automatically between orders 5, 9, and 13. This solver will solve ODEs for stiff systems. It treats the Jacobian matrix, df/ dy, as either a full or a banded matrix and as either user-supplied or internally approximated by difference quotients.

The code for these solvers was obtained from Drs. E. Hairer and G. Wanner. We deeply appreciate Dr. Hairer's permission to distribute this code as part of XSIM. For additional information about the solvers, refer to the second edition of Solving Ordinary Differential Equations (I and II) by Ernst Hairer and Gerhard Wanner published by Springer-Verlag (ISBN 0-387-56670-8 and 0-387-56675-9). The code for the solvers can be obtained from Dr. Hairer's web page at:

http://www.unige.ch/math/folks/harier/software.html

Two levels of user control over the solvers are offered. The basic user level uses default values for most of the ODE control parameters. Users with a working knowledge for numerical methods for solving ODEs can have greater control over the performance of the ODE solvers by selecting the advanced user level and setting values for the advanced parameters for the selected ODE solver.

The window for basic ODE control and status is shown below.

User level

The choices of user level are Basic and Advanced. If the basic level is chosen, only the parameters shown in the upper box of the window need to be set. As explained in section 3.3, if the advanced level is chosen the user should access the advanced parameters for the selected solver using the appropriate button at the bottom of the window.

ODE solver selection

The choices of solver are DOPRI5 and RADAU. The stiffness of the model equations is the major consideration when choosing the solver to use. A set of equations is said to be stiff if the ratio of the rates of its fastest process (maximum eigenvalue) to its slowest process (minimum eigenvalue) is greater than 10. If the ratio is between 2 and 10, the equations are mildly-stiff.

If the derivatives can be computed in a few arithmetic calculations and the equations are, at most, mildly stiff, then DOPRI5 is an appropriate solver. If the equations are more complex and known to be stiff, RADAU is an appropriate solver. If the user does not know the stiffness of the equations, both DOPRI5 and RADAU solvers should be used and the solutions should be carefully compared.

Since it takes many more calculations to compute a model solution using RADAU than with DOPRI5, the user should, especially for complex models, balance the time required to calculate solutions against the need for the increased solution accuracy obtained with RADAU.

Error tolerances

Both relative and absolute error tolerances for the error tests can be specified. The overall tolerance is relative_tolerance * |y_i| + absolute_tolerance. The default values of 0.00001 and 0.000001, respectively, ensure that the global solution has an error of less than 0.001%.

Note that when any dynamic variable has an initial condition of zero or passes through zero at any later value of the independent variable, it is important to set the absolute tolerance to a non-zero value. Otherwise, the step will have zero overall tolerance, and the solver will remain at that value of the independent variable indefinitely. Thus, the absolute error tolerance should be set to a non-zero value when the absolute noise of the solution or when any of the solution values are may have a zero value.

Number of function evaluations: The total number of times the ODE function was evaluated during the current simulation time step.

Number of computed steps: The number of steps the ODE solver took the integrate the model over the current simulation time step.

Number of accepted steps: For the current simulation time step, the number of steps that were accepted by the ODE solver based on the error tolerances set for the simulation run.

Number of rejected steps: For the current simulation time step, the number of steps that were rejected by the ODE solver based on the error tolerances set for the simulation run. (Note: A rejection during the first step is not counted.)

Users with a working knowledge for numerical methods for solving ODEs can have greater control over the performance of the ODE solvers. Those not familiar with terms like ``Jacobian'' and ``step size stabilization'' should consult a reference on ODE methods before using the advanced control parameters.

Short descriptions of the advanced control parameters for DOPRI5 and RADAU are given in the following sections. For more detailed information, see Drs. Hairer and Wanner's book and the code for the solvers available on Dr. Hairer's web page. More information about these resources is given in section 2.

3.3.2. Advanced DOPRI5 control

The window for advanced control of the DOPRI solver is shown below.

Rounding unit: The rounding unit used in testing calculated values. Default value is 2.3 x 10^-16.

Step size safety factor: The safety factor used in step size prediction. Default is 0.9.

Step size selection parameters: Limits for step size selection. The new step is chosen subject to the restriction that . Defaults are 0.2 and 10.0 respectively.

Beta for step size control: Larger values of b, greater than 0.1, make the step size control more stable. Negative values provoke b = 0. Default value is 0.04.

Maximum step size: The maximum allowable step size. A value of 0 invokes a maximum step size equal to the simulation step size for the independent variable. Default value is 0.

Initial step size: The initial step size used for solving the ODEs. A value of 0 causes the initial step size to be estimated by the solver. Default value is 0.

Maximum # of allowed steps: The maximum number of steps to be used by the solver for this step of the independent variable. A value of 0 invokes a maximum number of steps of 10⁵. Default value is 10⁵.

Stiffness testing step: The test for stiffness is invoked after the number of steps specified. If a negative value is specified, the test for stiffness is never invoked. Default value is 1000.

3.3.3. Advanced RADAU control

The window for advanced control of the RADAU solver is shown below.

External Jacobian parameters: The external Jacobian switch allows the user to invoke an external routine to calculate the Jacobian. Note that this is only useful if such a routine has been provided as part of the model. If no such routine has been provided, the switch has no effect on the model output. If a routine has been provided and is invoked, the upper and lower bandwidths may be set.

Rounding unit: The rounding unit used in testing calculated values. Default value is 1 x 10^-16.

Step size safety factor: The safety factor used in step size prediction. Default value is 0.9.

Jacobian recomputation: Controls when the Jacobian is recomputed. When the Jacobian evaluation is costly, a value of 0.1 is recommended. For small systems, the value should be smaller (i.e. around 1 x 10^-3). Default value is 1 x 10^-3.

Limits for no step change: In the range , the step size is not changed. Together with a large value of the Jacobian recomputation parameter, this saves LU-decompositions and computing time for large systems. For small systems, values of 1 and 1.2 are recommended. For large systems, 0.99 and 2 should be tried. Default values are 1 and 1.2.

Maximum step size: The maximum allowable step size. A value of 0 invokes a maximum step size equal to the simulation step size for the independent variable. Default value is 0.

Step size selection parameters: Limits for step size selection. The new step is chosen subject to the restriction that . Defaults are 0.2 and 8.0 respectively.

Contractivity bounds: The order is increased if the computed contractivity factor is smaller than bound 1 (default value: 0.002). The order is decreased if the computed factor is larger than bound 2 (default value: 0.8). Additionally, the order is decreased only if the stepsize ratio is in the range . Default values are 0.8 and 1.4

Jacobian to Hesselberg switch: If the ``Yes'' option is chosen, the Jacobian matrix is transformed to Hessenberg form. This is particularly advantageous for large systems with a full Jacobian. This does not work for banded systems of implicit systems. Default value is ``No.''

Maximum number of allowed steps: This is the maximum number of allowed integration steps for each simulation step. A value of 0 invokes the default value. Default value is 100000.

Maximum number of Newton iterations: This parameter affects the number of iterations of Newton's method of interpolation used when changing the integration stepsize. The actual number of iterations used is Nmax + 2.5*(Ns - 3), where Nmax is the maximum number of iterations and Ns is the number of stages (discussed below).

Extrapolated collocation switch: If the ``Yes'' option is chosen, the extrapolated collocation solution is taken as the starting value of Newton's method. If the ``No'' option is chosen, starting values of 0 are used. The latter is recommended if Newton's method has difficulties with convergence which is indicated the number of steps taken is larger than the number of steps accepted plus the number of steps rejected. The default value is ``Yes.''

Step size strategy switch: Two options, Gustafsson and Classical, are available. The former seems to produce safer results. For simple problems, the latter often produces slightly faster runs.

M values: If the differential system has the special structure: Yi' = Yi+m2 for i = 1, ... , m1, where m1 is a multiple of m2, a substantial gain in computation efficiency can be achieved by setting values for m1 and m2. For example, for second order systems P' = V, V' = G(P,V), where P and V are vectors of dimension n/ 2, one should set m1 = m2 = n/2. Default values are m1 = 0 and m2 = m1.

Setting m1 > 0 impacts the structure of a user-supplied Jacobian changes and the meaning of the upper and lower Jacobian bandwidth parameters. See Dr. Hairer's books or website for details.

Ns stages parameters: Stages are sub-steps within the integration step. Possible values for the minimum and maximum number of steps are 1, 3, 5, and 7. Default values are 3 and 7 respectively. The number of initial stages is also limited to 1, 3, 5 and 7. By default it is evaluated from and set equal to the minimum number of stages.

Number of Jacobian evaluations: The number of evaluations of the Jacobian, either analytically or numerically.

Number of LU decompositions: The number of LU-decompositions of both matrices.

Number of FB substitutions: The number of forward-backward substitutions of both matrices. The number of substitutions used for stepsize selection is not included.

Copyright 1995-1998, University of Washington

Last modified: 01:22pm PDT October 17, 1998