Numerical solvers have various theoretical formulations such as single- vs. multi-step, low- vs. high-order, explicit vs. implicit solvers [6–8]. To systematically compare the numerical solvers, we have examined qualitative performance measures such as convergence and accuracy, and quantitative performance measures such as computational cost in terms of function evaluations, partial derivatives, LU decompositions, solutions of linear systems, 'take-off' points, period and amplitude of oscillation.

Adaptive schemes allow savings in computational cost without compromising accuracy, i.e., the solvers take smaller step-sizes when the results change rapidly, or larger step-sizes when the results move slowly. With a specified tolerance interval for the magnitude of step errors, the step-size may be adjusted so that smaller steps are taken where the step error is large and vice versa. This adaptive step-size feature is studied by performing simulations under 2 conditions: (a) relaxed relative error tolerance (RET), i.e. in MATLAB, 'RelTol' = 10^{-3}, and (b) strict relative error tolerance (SET), 'RelTol' = 10^{-6}.

Single-step numerical solvers depend only on one preceding time-point i.e. only y(t_{n-1}) is required to obtain y(t_{n}). Hence, trends from additional preceding time-points do not influence the solution of the immediate step. On the other hand, multi-step solvers require multiple preceding time-points to determine y(t_{n}). As such, multi-step methods are less sensitive to initial conditions compared to single-step methods in general. In the case of single-step Runge-Kutta (RK) methods – that form the basis for ode23, ode45 and ode23tb – the interval between t_{n-1} and t_{n} are further divided into subintervals based on the order of the RK method. This subdivision also helps to control error propagation at successive steps, leading to increased stability in the solution. Thus, with the exception of ode23 where the solution quickly diverges, RK methods in general are a popular choice for numerical solvers.

The explicit solvers compute y(t_{n}) using other known values, while the implicit solvers require an iterative process to solve y(t_{n}) instead. Convergence is not guaranteed in this iterative process and depends heavily on the termination criteria. On the other hand, because this iterative process refines the solution at each step, implicit solvers are generally more suited to solve stiff problems. An unexpected result occurs with ode23s, where the solution becomes unstable under RET.

The tradeoff for obtaining a high order of accuracy is increased computational cost. For example, the variable-order solver ode113 is capable of accuracy up to the 13th order. However, to achieve this accuracy, it performs extensive function evaluations depending on the highest order of the derivatives required. On the other hand, low-order solvers such as ode23 perform simpler function evaluations that are limited to computing the 3rd-order derivative.

From this comparison, our results suggest ode45 and ode15s are better numerical solvers when dealing with stiff systems of nonlinear ODEs. ode45 is an explicit RK method of order 4 while ode15s is variable-order capable of 1st–5th order accuracy. Furthermore, low- to medium-order methods appear to provide more accurate solutions compared to high-order methods. In addition, results suggest that under relaxed relative error tolerance, implicit solvers may handle stiff problems relatively better than explicit solvers; this is consistent with established guidelines in handling stiff problems. However, when the relative error tolerance is strict, all numerical solutions studied here converge to limit cycle behavior. Under this condition, explicit solvers have an advantage over implicit solvers in terms of computational cost because of the lesser number of operations required. In addition, LC behavior is dependent on the parameter set {s, q, w, f} and initial conditions for the species involved [9–12]. The choice of the parameter set and initial conditions is based on a previous solution [13] that has been verified to exhibit LC behavior. This implies that the parameter set used for this study may not be the best choice. Thus, the characterization of the parameter space for the Oregonator model remains a worthwhile problem.

The key advantage in selecting MATLAB numerical solvers for comparison is that they are widely used in engineering and science and these solvers' reliability is of critical interest to the community. For those experienced users who may implement these solvers individually, it is also critical for them to understand the intrinsic differences in all types of solves, so as to be able to tailor solver algorithms to specific problems. For instance, Cocherová modified the ode23 solver to solve the Hodge-Huxley model of electrical conduction through a nerve fiber [14]. Results from her work demonstrate improved accuracy and convergence of solutions even under RET.

The next step is to determine whether there exists an analytical system-model dependency, that is, to extract and map features of biological systems to concrete analytical models. This extraction process for biomolecular systems emulates the modeling process of conventional systems such as mechanical, electrical and fluid flow systems. The numerical solver selection results from this work should help in validating the analytical model and also in translating systems biology to clinical treatrment and therapeutic engineering applications.