Analysis of time-stepping methods for the monodomain model

TL;DR

This paper compares various time-stepping methods for solving the monodomain model with different ionic models, analyzing stability, convergence, and efficiency to guide method selection based on model stiffness and complexity.

Contribution

It provides a comprehensive analysis of explicit, semi-implicit, exponential, and deferred correction methods for the monodomain model, including stability criteria and performance evaluation.

Findings

Explicit and semi-implicit methods have specific stability limits depending on ionic model complexity.

All methods achieve optimal convergence order despite ionic model non-differentiability.

Efficiency varies with method and model stiffness, guiding optimal method choice.

Abstract

To a large extent, the stiffness of the bidomain and monodomain models depends on the choice of the ionic model, which varies in terms of complexity and realism. In this paper, we compare and analyze a variety of time-stepping methods: explicit or semi-implicit, operator splitting, exponential, and deferred correction methods. We compare these methods for solving the bidomain model coupled with three ionic models of varying complexity and stiffness: the phenomenological Mitchell-Schaeffer model, the more realistic Beeler-Reuter model, and the stiff and very complex ten Tuscher-Noble-Noble-Panfilov (TNNP) model. For each method, we derive absolute stability criteria of the spatially discretized monodomain model and verify that the theoretical critical time-steps obtained closely match the ones in numerical experiments. We also verify that the numerical methods achieve an optimal order of convergence on the model variables and derived quantities (such as speed of the wave, depolarization time), and this in spite of the local non-differentiability of some of the ionic models. The efficiency of the different methods is also considered by comparing computational times for similar accuracy. Conclusions are drawn on the methods to be used to solve the monodomain model based on the model stiffness and complexity, measured respectively by the eigenvalues of the model's Jacobian and the number of variables, and based on strict stability and accuracy criteria.