An Explicit Local Space-Time Adaptive Framework for Monodomain Models

TL;DR

This paper tackles the bottleneck of slow monodomain cardiac electrophysiology simulations by introducing an explicit local space-time adaptive framework that combines discontinuous Galerkin spatial discretization with tree-based adaptive mesh refinement and synchronous local time stepping. The method rests on a primal symmetric interior penalty DG formulation, yielding per-element decoupled ODEs and enabling localized refinement in space and time driven by local error indicators, including a Kelly-type spatial indicator and an RV-T temporal indicator. Key contributions include the DG-based formulation, an efficient local time stepping strategy with barrier time stepping and CFL-based substep selection, and a comprehensive evaluation across conduction velocity, spiral wave, and idealized left ventricle benchmarks, reporting wall-clock speedups from 2× to 20× in serial while maintaining accuracy. The results demonstrate substantial practical potential for accelerating cardiac simulations, particularly in scenarios with localized wavefronts, and suggest pathways for further optimization and parallelization, as well as applicability to electromechanical models with moving domains.

Abstract

We present a new explicit local space-time adaptive framework to decrease the time required for monodomain simulations for cardiac electrophysiology. Based on the localized structure of the steep activation wavefront in solutions to monodomain problems, the proposed framework adopts small time steps and a tree-based adaptive mesh refinement scheme only in the regions necessary to resolve these localized structures. The time step and mesh adaptation selection process is fully controlled by a combination of local error indicators. The main contributions of this work consist in the introduction of a primal symmetric interior penalty formulation of the monodomain model and an efficient algorithmic strategy to manage local time stepping for its temporal discretization. In a first serial implementation of this framework, we report decreases in wall-clock time between 2 and 20 times with respect to an optimized implementation of a commonly used numerical scheme, showing that this framework is a promising candidate to accelerate monodomain simulations of cardiac electrophysiology.

Figures (5)

• Figure 1: Left: Example of tree-based adaptive mesh refinement on a single quadrilateral element. In this case, refinement operations add four children (i.e., reading the figure left to right) while coarsening operations collapse them into their parent elements (i.e., reading the figure right to left). Right: Adaptive mesh refinement with hexahedral elements during the conduction velocity benchmark NieKerBenBerBerBraCheClaFenGarHeiLanMalPatPlaRodRoySacSeeSkaSmi:2011:vct (rendered with GLVis). It can be clearly seen that the mesh is properly refined around the wavefront and coarsened everywhere else.
• Figure 2: Schematic for one barrier time step on a mesh with three root elements. Here the center element is refined at the beginning of the time step. Lines highlighted in yellow illustrate the active elements where the solution has been computed in the last time substep. Solid red arrows indicate the advance in time whereas dashed red arrows correspond to computing the interpolated values (red circular markers).
• Figure 3: Left: Schematic illustration of the benchmark setup. The stimulus is applied in the region marked with 'S' and the red diagonal represents the line along which local activation times are measured. Center: The number of time steps for all elements in the mesh over time for the newly proposed S-LTS scheme. The root mesh contains 420 elements and the element number peaks at 25159 elements at 27ms. This corresponds to the point in time when the wavefront size peaks. Hence it can be observed that the AMR follows the wavefront closely since the number of element evaluations directly correlates with the number of elements through the enforcement of the CFL condition. Right: Local activation times along the measurement line obtained with the LTG (reference) and S-LTS solutions. Note that the reference solution is in good quantitative agreement with the solution reported in the original conduction velocity benchmark study of NieKerBenBerBerBraCheClaFenGarHeiLanMalPatPlaRodRoySacSeeSkaSmi:2011:vct. We computed the solution for 50ms. The proposed S-LTS method required 81s while the reference simulation, using an optimized implementation of the classical operator splitting scheme from KriSarKlu:2013:nqoQuGar:1999:aas (see text for more details), required 1237s ($\approx 21min$), resulting in a speed up of $\approx 15$. IO operations as well as the time for the LAT computation were excluded from these time measurements.
• Figure 4: Spiral wave benchmark on a 16cm by 16cm domain. A) Initial condition for the transmembrane potential. B) Initial condition for the h-gate. This setup is known to induce spiral waves across a large number of ionic models, (see, e.g., PatGalCorKabFenGra:2020:duq). C) Transmembrane potential at $t=1s$ for the Lie-Trotter-Godunov operator splitting of the monodomain model with time step length $\Delta t=0.01ms$, grid size $0.156mm$ and order $1$ Lagrange ansatz space. The legend to interpret the magnitude of the transmembrane potential is the same as the one reported in panel A. The black line represents the transmembrane isopotential line at $-30mV$ for the newly proposed scheme. For the reference simulation and our scheme, we can observe that the wavefront shape is similar and the location of the wavefronts is close. Furthermore, no lattice pinning is observed.
• Figure 5: Benchmark on idealized left ventricle modeled as a truncated ellipsoid. Left: Reference activation time on a fine grid with element size between $0.15mm$ and $0.3mm$, and a small time step length of $\Delta t=0.001ms$ together with the isocontours at $t=1, 6, 12, 18, 24 ms$. Right: Difference in activation times between the reference solution shown on the left and the newly proposed S-LTS scheme. The small artifacts results from a combination of a low sampling rate for the S-LTS scheme (the S-LTS solution is sampled only at the barrier time step of 0.1ms), as well as a slightly varying wave speed in our scheme. For example, we notice a slight increase of the wavespeed for the S-LTS scheme toward the apical region. Nevertheless, the difference in activation times remains below $0.5ms$.