Automatic time-step adaptation of the forward Euler method in simulation of models of ion channels and excitable cells and tissue

https://doi.org/10.1016/j.simpat.2008.03.007Get rights and content

Abstract

The forward Euler method is commonly used to simulate models of ion channels and excitable cells and tissue. We show that the computational efficiency of these simulations could be significantly improved with automatic time-step adaptation. For this purpose a new easy-to-implement time-step adaptation method was developed.

Introduction

Solving a system of ordinary or combined ordinary and partial differential equations (ODE and PDE) is a common task when simulating models of excitable cells and tissues. ODEs usually describe the function of the ion channels in the cell membrane [1], [2] and the PDEs are used to describe spatial phenomena such as spreading of the action potential in a tissue [1] or diffusion of ions inside cells [3]. The forward Euler method with a fixed time-step is often used to simulate these models [1], [2], [4]. The advantage of the Euler method is its simple implementation especially in large scale problems which need parallel computing to be solved [1]. The weakness of the Euler method is the demand for a small integration time-step causing slow computation of the solution. In the models of ion channels and excitable cells and tissue, this problem is even enhanced; these models are often relatively stiff as they consist of rapidly changing components, such as action potential during upstroke [1] and gating and state-variables of the ion channel models, and thus they demand the use of an extremely small time-step [1], [2], [4]. In this paper we show that the simulation times of complex ion channel, cell and tissue models could be significantly improved without loss of the accuracy of the solution by adapting the time-step in the Euler method during simulation. For this purpose we developed a new easy to implement method for automatic step size detection in the Euler method which is suitable for stiff problems.

Section snippets

Methods

When applying the Euler method to solve a system of ODEs and PDEs, the PDEs are approximated to systems of ODEs by approximating the partial derivates [5]. Thus, the solving of an ODE or combined ODE and PDE problem is in both cases solving an ODE system. The time-step of the ODE system integration has to be chosen according to the ODE requiring the smallest time-step. This small time-step is used in the whole time-scale of the solution even though at some parts of the solution a longer

Results

We first investigated if we could reduce the CPU times by using the limited method with a certain value for B in (6) instead of using the non-limited method to estimate the proper time-step. We ran the test simulation at the error tolerance limit of 10−5 (0.001%) with the non-limited method and with the limited method with B values ranging from 2.5 to 1015. We found that between B values of 25–109 the CPU time was smaller with the limited method compared to the non-limited method (Fig. 1b). The

Discussion and conclusions

In this paper we present a method referred to as the limited method which automatically adapts the used time-step for Euler integrations of stiff differential equation systems. The step size determination based on the local relative error of the solution (3) (the non-limited method) can be derived from the estimate of the local absolute error [5]. The novel finding of this study was that the step size determination could be made more suitable to stiff problems if we prevent the underestimation

Acknowledgements

We thank SL. Hänninen for valuable comments on the manuscript. This study was supported by the Finnish Heart Research Foundation and Sigrid Juselius Foundation.

References (8)

  • J. Sneyd et al.

    A model of calcium waves in pancreatic and parotid acinar cells

    Biophys. J.

    (2003)
  • K.H.W.J. ten Tusscher et al.

    A model for human ventricular tissue

    Am. J. Physiol. Heart Circ. Physiol.

    (2004)
  • J. Sneyd et al.

    A dynamic model of the type-2 inositol trisphosphate receptor

    Proc. Natl. Acad. Sci. USA

    (2002)
  • P.S. Haddock et al.

    Subcellular [Ca2+]i gradients during excitation-contraction coupling in newborn rabbit ventricular myocytes

    Circ. Res.

    (1999)
There are more references available in the full text version of this article.

Cited by (5)

  • Discretization for a spray deposition model: Criteria for temporal and spatial differencing

    2013, Computers and Electronics in Agriculture
    Citation Excerpt :

    Fung (2003) also studied the effect of using predetermined coefficients with the weighted residual method in developing Δt integration algorithms for second order differential equations. Krohonen and Tavi (2008) developed a Δt adaptation method, the ‘Limited’ method, for performing Euler integrations of stiff differential equation systems at faster rates. Spatial discretization, achieved by finite differencing, is also essential in numerically approximating the solution of an otherwise continuous space-dependent model.

  • The most precise computations using Euler's method in standard floating-point arithmetic applied to modelling of biological systems

    2013, Computer Methods and Programs in Biomedicine
    Citation Excerpt :

    Table 1 summarizes a comparative study of different numerical methods for numerical integration over the interval [0, 20, 000] for the initial-value problem (9). ML: different Matlab's ODE solvers; KT: two methods proposed in [2] (limited and unlimited); EM: method suggested in this paper; E: Euler's method with constant step size; R: Rosenbrock classical method. The programs for Euler's method and method derived in this paper are written on C++ and worked in a single precision float.

  • Numerical analysis of heat transfer characteristics in microwave heating of magnetic dielectrics

    2012, Metallurgical and Materials Transactions A: Physical Metallurgy and Materials Science
View full text