Automatic time-step adaptation of the forward Euler method in simulation of models of ion channels and excitable cells and tissue
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)
- et al.
A model of calcium waves in pancreatic and parotid acinar cells
Biophys. J.
(2003) - et al.
A model for human ventricular tissue
Am. J. Physiol. Heart Circ. Physiol.
(2004) - et al.
A dynamic model of the type-2 inositol trisphosphate receptor
Proc. Natl. Acad. Sci. USA
(2002) - et al.
Subcellular [Ca2+]i gradients during excitation-contraction coupling in newborn rabbit ventricular myocytes
Circ. Res.
(1999)
Cited by (5)
Discretization for a spray deposition model: Criteria for temporal and spatial differencing
2013, Computers and Electronics in AgricultureCitation 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 BiomedicineCitation 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 ScienceMathematical model of mouse embryonic cardiomyocyte excitation-contraction coupling
2008, Journal of General Physiology