Computer Programs in Biomedicine 10 (1979) 1-15 © Elsevier/North-Holland Biomedical Press
Section I: Programs, program packages and systems THE EFFICIENT NUMERICAL SOLUTION OF BIOLOGICAL SIMULATION PROBLEMS Richard E. PLANT Department of Mathematics, University of California, Davis, CA 95616, USA
A subroutine called DESOL for the numerical solution of ordinary differential equations of the type arising in biological simulation problems is described. DESOL is about as efficient as current high quality integrators, but because of its compactness it can be easily used on small computers. The subroutine has excellent stability properties and functions with very little required input from the user. In addition, it has features which aid in debugging associated programs. Several test and example problems are given, as is the derivation of the major formulae used in the package. Numerical integration
Biological models
Subroutines
1. Introduction Experience with the numerical solution of biological simulation problems (e.g., [1,2]) has lead to the conclusion that many such problems are not well suited for the more classical numerical methods. On the other hand, the particular properties of these problems are such that many of the features provided by recently developed, high quality integration subroutines are not useful in this application. Since each of these features entails a certain amount of computational overhead, the subroutines are rendered less efficient in computing the solution. In addition, the use of these features may result in a sacrifice of other desirable properties such as compactness, ease of pro-
gramming and numerical stability. Accordingly, a subroutine called DESOL which incorporates the useful features of the integrators such as GEAR [3], EPISODE [4] and ODE [5], but which is better suited to biological simulation modelling, has been developed. The subroutine package consisting of DESOL and two auxilliary subroutines is described in this report. As an example of the capabilities of the package, fig. 1 represents the numerical solution of the voltage component of the Hodgkin-Huxley equations. In order to obtain a smooth, accurate curve, this solution was plotted every 0.02 ms, so, for example, a fixed step size fourth order Runge-Kutta method would require 1750 steps and 7000 function evaluations to compute the curve. The DESOL package obtained the soluUon in 261 steps and 553 function evaluations with no loss in accuracy.
2. General description :~_ i
,
T (MSEC)
Fig. 1. Numerical solution of the Hodgkin-Huxley equations. The figure actually contains two curves: a solid curve representing the solution obtained by DESOL and a dashed curve representing the 'exact' solution, obtained as described in the text. A slight separation of the curves may be noted in the second peak.
DESOL is designed to solve simulation problems which may be written in the form of a vector differential equation as: dy =f(v, t),Y(to) =Yo dt
(1)
where y = (yl, y : ..... yn) and f(y, t) = (fl, f2 ..... fn), the f i each being functions o f y and t. Such simulation problems are often dismissed as requiring only a
2
R.E. Plant, Efficient numerical solution
small amount of computer time for their solution. In many cases, however, the problems must be integrated repeatedly with many small changes in the parameters, for example, in parameter estimation procedures. The total time required by the problem, and the savings provided by an efficient solution package, may therefore be substantial. The DESOL package actually consists of three subroutines: DESOL itself, a driver subroutine which organizes and supervises the integration, and two auxilliary subroutines which are called by DESOL: ADAMS, the basic integrator which performs one step of numerical integration, and INTERP, which interpolates the solution to the outer point. The package is designed to be most efficient for problems which have the following general characteristics: (i) The accuracy requirements are not too stringent (e.g., a global error on the order of a few percent is acceptable). (ii) The problem may be mildly 'stiff' (i.e., have natural time constants which may differ by one or two orders of magnitude). In addition, the solutions may oscillate or move from one value to another rapidly, or the function f o f eq. (1) may be discontinuous or badly behaved. (iii) Values of the solution are required at many points to facilitate graphical output. (iv) The functions f/(y, t) may be relatively expensive to evaluate. (v) The program may be run on a minicomputer, and thus should be very compact. Properties (i) through (if) are similar to those discussed by Shampine [7], while property (iv) is the exact opposite. This is because many biological problems involve experimentally derived functions which are fairly complex. For example, on the Burroughs B6700 computer each function evaluation in the Hodgkin-Huxley equations of fig. 1 requires roughly 125 multiplications and divisions and about the same number of additions. This is considerably above the dividing line for operations between 'expensive' and 'cheap' functions (about 25 operations) given by Hull et al. [8]. The properties (i) through (v) are the main consideration in deciding which numerical method to use in the program (the various methods available are discussed, e.g., in [9]). Properties (iii) and (iv) imply that the method of choice is the Adams-Bashforth-Moulton predictor-
corrector method, and this is the method used in DESOL. In the form used (cf. the Appendix) the method requires only two function evaluations per accepted step and one per rejected step, so that the steps are relatively inexpensive. With regard to (iii), DESOL makes use of the following property of the Adams methods (see [5], ch. 4 for discussion and proofs). We may write the Adams-Moulton method as; tn+ l Yn + 1 = Yn + f
P(t) dt
(2)
tn where fin and Yn + 1 are the approximations to Y(tn) and Y(tn+l)and P(t) is a vector of interpolating polynomials. If tout is any output point such that tn ~ tout ~< tn +1, then we may evaluate Yout according to the formula: tout
Your =Yn+l + J
P(t) dt
(3)
tn+ l
and, roughly speaking,Yout will have the same accuracy aSYn andyn+l. This is very convenient because, for example, if we wish to print out y every 0.01 units of t and the most efficient step size is 0.25 units of t, when we may compute 25 values based on a single step. Properties (i) and (v) imply that the method should be a low order one, and that the order should be fixed. The DESOL package uses a second order predictor and a third order corrector on most steps (in starting and when it encounters a presumed discontinuity it uses a first order predictor and second order corrector). This order was chosen for the following reasons: (a) The absolute stability properties of this combination are as good as those of lower order methods and better than those of higher order methods ([5], p. 135); (b) The cost of computing the integration formulas is small; (c) Variable order methods which were tested on biological simulation problems almost always operated at that order. Property (ii) implies that the method should be a variable step size one. This again favors the Adams methods which permit a very simple estimation of the
R.E. Plant, Efficient numerical solution
local error using a formula known as Milne's device [9]. In addition, property (ii) implies that standard functional iteration rather than, say, Newton iteration [3] is generally acceptable, and this is the type used. The method used in DESOL is a truly variable-step Adams method, as opposed to the more common technique which is to use a fixed step formula and at each change in step size to either restart the integration using a Runge-Kutta scheme or to interpolate to the new starting values. The variable step formulas make use of the divided difference form of the interpolating polynomial in eq. (2) and in the analogous predictor equation (Appendix). This method was originally discussed by Krogh, for example, in [10], and has been implemented by Shampine and Gordon [5]. The DESOL package does not make use of the algorithm derived by Krogh to advance a step since at the second and third order it is easier to compute the required coefficients directly. The DESOL package is different from many others in that it requires very little from the user. Discussions between the author and potential users have revealed that many of these potential users are not very familiar with recent developments in numerical analysis and would be best served by an integration package which requires an absolute minimum of information from the user, and which contains many debugging aids. On the other hand, many potential users are quite sophisticated in numerical methods and may even wish ~o make minor modifications in the program. The package has been written in an attempt to accomodate both types of users. Detailed instructions for use of the package are given in the comment cards at the beginning of subroutine DESOL (section 4); the following is an outline of what the user must provide. First, the user must provide a main program which calls subroutine DESOL. The program must contain arrays Y, SAVE, and ISAVE, dimensioned at least NEQN (the number of equations in the system), 12 + 4 X NEQN, and 4, respectively. Before the first call, the user must set NEQN, TO, TOUT, and IFLAG. TO is the initial value of t (i.e., to in eq. (1)), TOUT is the value of t at which the solution is desired, and IFLAG must be set at zero as an indicator that this is the initial call to DESOL. The user must also set each value in the Y array equal to the corresponding initial value
3
y~, and should set each value in the arrays SAVE and ISAVE equal to zero. If all goes well, DESOL will return to the main program with the values if(tout), computed with an accuracy of about one percent, in the array Y and with IFLAG equal to 1. (If IFLAG > 1, a failure has occurred as described below.) To con. tinue the integration, the user simply sets TO to the current T O U T , provides a new TOUT, and calls DESOL again. No modification in any of the other parameters should be made, with the following single exception. Some problems involve functions fO', t) which are discontinuous in t. An example would be:
dy
--~ = a ( t ) - y ,
y(O) = 1
where
a(t)
= [ O: t ~< 1
t 1:
t> 1
The DESOL package attempts to detect such discontinuities and to deal with them, but the process is made more efficient if the main program informs DESOL of such a discontinuity. This is done by setting IFLAG = 0 on each call for which TO has a value at which the function f i s discontinuous. An example is given in section 5. If DESOL is called with values TOUT ~< TO, then a single integration step is performed and the value of TOUT is set equal to the value of t at the end of this step, with the contents of the Y array containing the values o f y at this value of t. Thus, if the user does not wish to make use of the interpolation feature but would rather have the values of the solution at the natural steps, he or she may accomplish this by simply setting TOUT = TO after each call. This is illustrated in section 5. This feature has been found to be especially useful as a debugging tool since it allows the user to see the output at the actual computation points. In addition to a main program, the user must provide a subroutine FTN (NEQN, T, Y, F). It is in this subroutine that the values of the functions fi(.v, t), i = 1 ..... NEQN, are computed. The parameters NEQN and T are the number of equations and the current value of t, respectively, the array Y contains the values i f , i = 1 ..... NEQN, and the computed values f O ' , t) must be placed in the array F. Both Y and
4
R.E. Plant, Efficient numerical solution
F must be dimensioned NEQN. Only values in array F should be altered in subroutine FTN, neither NEQN, T, nor the array Y should be changed. Certain of the parameters have values which have been 'tuned' on the basis of experience, and which the user may wish to change. These parameters are all assigned values in subroutine ADAMS. The most important are as follows: (i) ERRMAX, the maximum accepted magnitude of the relative local truncation error of a component of the solution ERRMAX is set at 5 × 10 -3 in the program. Experience indicates that this is sufficient to guarantee a global error of ~
Section I: Programs, program packages and systems THE EFFICIENT NUMERICAL SOLUTION OF BIOLOGICAL SIMULATION PROBLEMS Richard E. PLANT Department of Mathematics, University of California, Davis, CA 95616, USA
A subroutine called DESOL for the numerical solution of ordinary differential equations of the type arising in biological simulation problems is described. DESOL is about as efficient as current high quality integrators, but because of its compactness it can be easily used on small computers. The subroutine has excellent stability properties and functions with very little required input from the user. In addition, it has features which aid in debugging associated programs. Several test and example problems are given, as is the derivation of the major formulae used in the package. Numerical integration
Biological models
Subroutines
1. Introduction Experience with the numerical solution of biological simulation problems (e.g., [1,2]) has lead to the conclusion that many such problems are not well suited for the more classical numerical methods. On the other hand, the particular properties of these problems are such that many of the features provided by recently developed, high quality integration subroutines are not useful in this application. Since each of these features entails a certain amount of computational overhead, the subroutines are rendered less efficient in computing the solution. In addition, the use of these features may result in a sacrifice of other desirable properties such as compactness, ease of pro-
gramming and numerical stability. Accordingly, a subroutine called DESOL which incorporates the useful features of the integrators such as GEAR [3], EPISODE [4] and ODE [5], but which is better suited to biological simulation modelling, has been developed. The subroutine package consisting of DESOL and two auxilliary subroutines is described in this report. As an example of the capabilities of the package, fig. 1 represents the numerical solution of the voltage component of the Hodgkin-Huxley equations. In order to obtain a smooth, accurate curve, this solution was plotted every 0.02 ms, so, for example, a fixed step size fourth order Runge-Kutta method would require 1750 steps and 7000 function evaluations to compute the curve. The DESOL package obtained the soluUon in 261 steps and 553 function evaluations with no loss in accuracy.
2. General description :~_ i
,
T (MSEC)
Fig. 1. Numerical solution of the Hodgkin-Huxley equations. The figure actually contains two curves: a solid curve representing the solution obtained by DESOL and a dashed curve representing the 'exact' solution, obtained as described in the text. A slight separation of the curves may be noted in the second peak.
DESOL is designed to solve simulation problems which may be written in the form of a vector differential equation as: dy =f(v, t),Y(to) =Yo dt
(1)
where y = (yl, y : ..... yn) and f(y, t) = (fl, f2 ..... fn), the f i each being functions o f y and t. Such simulation problems are often dismissed as requiring only a
2
R.E. Plant, Efficient numerical solution
small amount of computer time for their solution. In many cases, however, the problems must be integrated repeatedly with many small changes in the parameters, for example, in parameter estimation procedures. The total time required by the problem, and the savings provided by an efficient solution package, may therefore be substantial. The DESOL package actually consists of three subroutines: DESOL itself, a driver subroutine which organizes and supervises the integration, and two auxilliary subroutines which are called by DESOL: ADAMS, the basic integrator which performs one step of numerical integration, and INTERP, which interpolates the solution to the outer point. The package is designed to be most efficient for problems which have the following general characteristics: (i) The accuracy requirements are not too stringent (e.g., a global error on the order of a few percent is acceptable). (ii) The problem may be mildly 'stiff' (i.e., have natural time constants which may differ by one or two orders of magnitude). In addition, the solutions may oscillate or move from one value to another rapidly, or the function f o f eq. (1) may be discontinuous or badly behaved. (iii) Values of the solution are required at many points to facilitate graphical output. (iv) The functions f/(y, t) may be relatively expensive to evaluate. (v) The program may be run on a minicomputer, and thus should be very compact. Properties (i) through (if) are similar to those discussed by Shampine [7], while property (iv) is the exact opposite. This is because many biological problems involve experimentally derived functions which are fairly complex. For example, on the Burroughs B6700 computer each function evaluation in the Hodgkin-Huxley equations of fig. 1 requires roughly 125 multiplications and divisions and about the same number of additions. This is considerably above the dividing line for operations between 'expensive' and 'cheap' functions (about 25 operations) given by Hull et al. [8]. The properties (i) through (v) are the main consideration in deciding which numerical method to use in the program (the various methods available are discussed, e.g., in [9]). Properties (iii) and (iv) imply that the method of choice is the Adams-Bashforth-Moulton predictor-
corrector method, and this is the method used in DESOL. In the form used (cf. the Appendix) the method requires only two function evaluations per accepted step and one per rejected step, so that the steps are relatively inexpensive. With regard to (iii), DESOL makes use of the following property of the Adams methods (see [5], ch. 4 for discussion and proofs). We may write the Adams-Moulton method as; tn+ l Yn + 1 = Yn + f
P(t) dt
(2)
tn where fin and Yn + 1 are the approximations to Y(tn) and Y(tn+l)and P(t) is a vector of interpolating polynomials. If tout is any output point such that tn ~ tout ~< tn +1, then we may evaluate Yout according to the formula: tout
Your =Yn+l + J
P(t) dt
(3)
tn+ l
and, roughly speaking,Yout will have the same accuracy aSYn andyn+l. This is very convenient because, for example, if we wish to print out y every 0.01 units of t and the most efficient step size is 0.25 units of t, when we may compute 25 values based on a single step. Properties (i) and (v) imply that the method should be a low order one, and that the order should be fixed. The DESOL package uses a second order predictor and a third order corrector on most steps (in starting and when it encounters a presumed discontinuity it uses a first order predictor and second order corrector). This order was chosen for the following reasons: (a) The absolute stability properties of this combination are as good as those of lower order methods and better than those of higher order methods ([5], p. 135); (b) The cost of computing the integration formulas is small; (c) Variable order methods which were tested on biological simulation problems almost always operated at that order. Property (ii) implies that the method should be a variable step size one. This again favors the Adams methods which permit a very simple estimation of the
R.E. Plant, Efficient numerical solution
local error using a formula known as Milne's device [9]. In addition, property (ii) implies that standard functional iteration rather than, say, Newton iteration [3] is generally acceptable, and this is the type used. The method used in DESOL is a truly variable-step Adams method, as opposed to the more common technique which is to use a fixed step formula and at each change in step size to either restart the integration using a Runge-Kutta scheme or to interpolate to the new starting values. The variable step formulas make use of the divided difference form of the interpolating polynomial in eq. (2) and in the analogous predictor equation (Appendix). This method was originally discussed by Krogh, for example, in [10], and has been implemented by Shampine and Gordon [5]. The DESOL package does not make use of the algorithm derived by Krogh to advance a step since at the second and third order it is easier to compute the required coefficients directly. The DESOL package is different from many others in that it requires very little from the user. Discussions between the author and potential users have revealed that many of these potential users are not very familiar with recent developments in numerical analysis and would be best served by an integration package which requires an absolute minimum of information from the user, and which contains many debugging aids. On the other hand, many potential users are quite sophisticated in numerical methods and may even wish ~o make minor modifications in the program. The package has been written in an attempt to accomodate both types of users. Detailed instructions for use of the package are given in the comment cards at the beginning of subroutine DESOL (section 4); the following is an outline of what the user must provide. First, the user must provide a main program which calls subroutine DESOL. The program must contain arrays Y, SAVE, and ISAVE, dimensioned at least NEQN (the number of equations in the system), 12 + 4 X NEQN, and 4, respectively. Before the first call, the user must set NEQN, TO, TOUT, and IFLAG. TO is the initial value of t (i.e., to in eq. (1)), TOUT is the value of t at which the solution is desired, and IFLAG must be set at zero as an indicator that this is the initial call to DESOL. The user must also set each value in the Y array equal to the corresponding initial value
3
y~, and should set each value in the arrays SAVE and ISAVE equal to zero. If all goes well, DESOL will return to the main program with the values if(tout), computed with an accuracy of about one percent, in the array Y and with IFLAG equal to 1. (If IFLAG > 1, a failure has occurred as described below.) To con. tinue the integration, the user simply sets TO to the current T O U T , provides a new TOUT, and calls DESOL again. No modification in any of the other parameters should be made, with the following single exception. Some problems involve functions fO', t) which are discontinuous in t. An example would be:
dy
--~ = a ( t ) - y ,
y(O) = 1
where
a(t)
= [ O: t ~< 1
t 1:
t> 1
The DESOL package attempts to detect such discontinuities and to deal with them, but the process is made more efficient if the main program informs DESOL of such a discontinuity. This is done by setting IFLAG = 0 on each call for which TO has a value at which the function f i s discontinuous. An example is given in section 5. If DESOL is called with values TOUT ~< TO, then a single integration step is performed and the value of TOUT is set equal to the value of t at the end of this step, with the contents of the Y array containing the values o f y at this value of t. Thus, if the user does not wish to make use of the interpolation feature but would rather have the values of the solution at the natural steps, he or she may accomplish this by simply setting TOUT = TO after each call. This is illustrated in section 5. This feature has been found to be especially useful as a debugging tool since it allows the user to see the output at the actual computation points. In addition to a main program, the user must provide a subroutine FTN (NEQN, T, Y, F). It is in this subroutine that the values of the functions fi(.v, t), i = 1 ..... NEQN, are computed. The parameters NEQN and T are the number of equations and the current value of t, respectively, the array Y contains the values i f , i = 1 ..... NEQN, and the computed values f O ' , t) must be placed in the array F. Both Y and
4
R.E. Plant, Efficient numerical solution
F must be dimensioned NEQN. Only values in array F should be altered in subroutine FTN, neither NEQN, T, nor the array Y should be changed. Certain of the parameters have values which have been 'tuned' on the basis of experience, and which the user may wish to change. These parameters are all assigned values in subroutine ADAMS. The most important are as follows: (i) ERRMAX, the maximum accepted magnitude of the relative local truncation error of a component of the solution ERRMAX is set at 5 × 10 -3 in the program. Experience indicates that this is sufficient to guarantee a global error of ~
