The Ithildin library for efficient numerical solution of anisotropic reaction-diffusion problems in excitable media

Abstract

Ithildin is an open-source library and framework for efficient parallelized simulations of excitable media, written in the C++ programming language. It uses parallelization on multiple CPU processors via the message passing interface (MPI). We demonstrate the library's versatility through a series of simulations in the context of the monodomain description of cardiac electrophysiology, including the S1S2 protocol, spiral break-up, and spiral waves in ventricular geometry. Our work demonstrates the power of Ithildin as a tool for studying complex wave patterns in cardiac tissue and its potential to inform future experimental and theoretical studies. We publish our full code with this paper in the name of open science.

Fig 1. Ithildin can be used to solve reaction-diffusion problems in excitable media.

Required inputs for the software are: the diffusion tensor and geometry of a medium—such as the heart muscle, a reaction term—the so-called model, and source terms—typically a stimulation protocol. Ithildin can then calculate the evolution of the model variables in the medium over time. During calculation, Ithildin records relevant spatio-temporal data and metadata, as well as detecting rotor cores as filaments. The data visualized here are taken from several simulations which will be discussed below.

Fig 2. Stencils for isotropic diffusion (a) and orthotropic diffusion (b) with numbering of the involved grid points.

Fig 3. Ventricle geometry with fiber direction, colored by the fiber helix angle.

Fig 4. Inhomogeneities in Sim 1.

The field inhom describes where unexcitable obstacles or exterior points are located with the value inhom = 0 and which cell model is to be used inside if inhom > 0. This figure also visualizes the filament trajectories colored by time, cf. section 4.7. Birth of a tip is denoted by stars, their deaths by crosses, and tips that are still around at the end of the simulation by points. The three trajectories with the longest lifetime are numbered by indices 0 through 2.

Fig 5. The S1S2 protocol illustrated for Sim 2.

The first and third row display the transmembrane voltage u at selected frames in time, and the second and fourth row the state space phase [17, 47]. The first stimulus is applied in the first frame in the hatched region at the left edge. The sensor location is marked with a cross. Right after the third frame, the sensor triggers the second stimulus in the hatched region at the bottom edge. A spiral wave forms. The phase singularity at the center of the spiral is tracked as the white curve, which is dotted for the entire trajectory of the phase singularity, and solid for its trajectory since the previous frame.

Fig 6. Extracted time traces of Sim 1.

The first two panels (a, b) contain the model variables u and v, followed by data computed by model wrappers, namely the diffusion term ∇ ⋅ D∇u (c) and the LAT (d), at an interior point, cf. section 4.5. Panel (e) contains the pseudo-EGM at an exterior point, next to the 2D domain, cf. section 4.6. The last panel (f) contains the number of rotors detected via phase singularities, cf. section 4.7. The black vertical lines indicate times at which a stimulus was applied.

Fig 7. Illustration of the filament point tracking algorithm.

Each coordinate plane adjacent to a grid point (i, j, k) is checked for an intersection of two surfaces (purple and orange lines) which are usually isosurfaces of state variables in $\underset{\_}{u}$. In this example, a filament point (white dot) will be found in the yz-face.

Fig 8. Phase defect detection for a single meandering spiral in Sim 1.

By recording the activation times of the transmembrane voltage V_m (top row), the activation time phase φ [17, 51] can be computed (middle row). Where this phase is discontinuous, a phase defect is localized. This is visualized as the phase defect density $\varrho$ (bottom row). In these four frames, a single rotor is tracked shortly after its formation. Due to conduction block, which can be seen as an extended phase defect line, the spiral moves through the medium. Afterwards, the spiral remains mostly stationary, leading to a shorter phase defect line.

Fig 9. Phase defect detection during figure-of-eight spiral wave pair creation in Sim 1.

Visualization in the same style as Fig 8. A long conduction block line breaks apart into two rotors.

Fig 10. Final frame of Sim 4 of the BOCF model in ventricle geometry.

The classical filament is plotted as the purple lines. The phase defect surface is contained in the gray contour.

Fig 11. Results of the cardiac electrophysiology benchmark.

The benchmark was proposed by Niederer et al. (2011) [56] and is implemented in Sim 5. To ease comparisons with the paper introducing the benchmark, we present our results in a similar style. In panels (a)-(c), we present the LAT in the medium at different spatial and temporal resolutions. Coloring according to LAT is consistent across these panels. In panel (d), we measure the computational speed of Ithildin at different number of used processes. A linear fit is shown in the double-logarithmic plot which excludes the data point at N_proc = 8.