Quantitative analysis of numerical solvers for oscillatory biomolecular system models

Chang F Quo,May D Wang

doi:10.1186/1471-2105-9-s6-s17

Abstract

BackgroundThis article provides guidelines for selecting optimal numerical solvers for biomolecular system models. Because various parameters of the same system could have drastically different ranges from 10-15 to 1010, the ODEs can be stiff and ill-conditioned, resulting in non-unique, non-existing, or non-reproducible modeling solutions. Previous studies have not examined in depth how to best select numerical solvers for biomolecular system models, which makes it difficult to experimentally validate the modeling results. To address this problem, we have chosen one of the well-known stiff initial value problems with limit cycle behavior as a test-bed system model. Solving this model, we have illustrated that different answers may result from different numerical solvers. We use MATLAB numerical solvers because they are optimized and widely used by the modeling community. We have also conducted a systematic study of numerical solver performances by using qualitative and quantitative measures such as convergence, accuracy, and computational cost (i.e. in terms of function evaluation, partial derivative, LU decomposition, and "take-off" points). The results show that the modeling solutions can be drastically different using different numerical solvers. Thus, it is important to intelligently select numerical solvers when solving biomolecular system models.ResultsThe classic Belousov-Zhabotinskii (BZ) reaction is described by the Oregonator model and is used as a case study. We report two guidelines in selecting optimal numerical solver(s) for stiff, complex oscillatory systems: (i) for problems with unknown parameters, ode45 is the optimal choice regardless of the relative error tolerance; (ii) for known stiff problems, both ode113 and ode15s are good choices under strict relative tolerance conditions.ConclusionsFor any given biomolecular model, by building a library of numerical solvers with quantitative performance assessment metric, we show that it is possible to improve reliability of the analytical modeling, which in turn can improve the efficiency and effectiveness of experimental validations of these models. Also, our study can be extended to study a variety of molecular-level system models for human disease diagnosis and therapeutic treatment.

Highlights

This article provides guidelines for selecting optimal numerical solvers for biomolecular system models
One major hurdle is that when describing a non-observable biomolecular system, the set of ordinary differential equations (ODEs) often do not have closedform analytical solutions
Two sets of simulations were performed under different conditions of error tolerance – (a) relaxed relative error tolerance (RET) i.e. in MATLAB, 'RelTol' = 10-3 and (b) strict relative error tolerance (SET) i.e. 'RelTol' = 10-6 ; other invariant conditions are a parameter set {s = 100, w = 3.835, q = 10-6, f = 1.1} and initial conditions {α = 20, η = 1.1, ρ = 20}

Summary

Introduction

This article provides guidelines for selecting optimal numerical solvers for biomolecular system models. Previous studies have not examined in depth how to best select numerical solvers for biomolecular system models, which makes it difficult to experimentally validate the modeling results. To address this problem, we have chosen one of the well-known stiff initial value problems with limit cycle behavior as a test-bed system model. We have chosen one of the well-known stiff initial value problems with limit cycle behavior as a test-bed system model Solving this model, we have illustrated that different answers may result from different numerical solvers. There is an urgent need to study numerical solver behavior in biomolecular systems modeling using both qualitative and quantitative performance measures.

Objectives

Methods

Results

Discussion

Conclusion