  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.
