Quantitative comparison of numerical solvers for models of oscillatory biochemical systems

The goal of this work is to find optimal numerical solvers to study models of biological systems of interest. We report 2 simple rules to select the optimal numerical solver(s) for solving stiff, complex oscillatory systems. Numerical solvers currently used to solve mathematical models in systems biology can be ill-conditioned due to stiffness. As a case study, we choose the classic Belousov-Zhabotinskii (BZ) reaction, described by the Oregonator model. We determine the optimal numerical solver(s) to handle stiff initial-value problems that lead to limit cycle behavior, by systematically comparing qualitative and quantitative performance measures i.e. convergence, accuracy and computational cost of numerical solvers that are widely used by the engineering and modeling community. This is a cornerstone for our long-term research objective to systematically study a variety of molecular-level models for biomedical systems for human disease diagnosis and therapeutic treatment in order to understand and predict disease mechanism and progression.