This article was previously published in PLOS ONE 19(9): e0303674 (Kabus, Cloet, et al., 2024) and is a chapter of my dissertation.

Authors:

• Desmond Kabus^1,2 (ORCID: 0000-0002-6965-5211)

• Marie Cloet¹ (ORCID: 0000-0002-8974-6401)

• Christian Zemlin³ (ORCID: 0000-0001-5834-5544)

• Olivier Bernus⁴ (ORCID: 0000-0003-3917-5791)

• Hans Dierckx¹ (ORCID: 0000-0003-0899-8082)

Institutions:

1. Department of Mathematics, KU Leuven Campus Kortrijk (KULAK), Etienne Sabbelaan 53, 8500 Kortrijk, Belgium

2. Laboratory of Experimental Cardiology, Leiden University Medical Center (LUMC), Albinusdreef 2, 2333 ZA Leiden, The Netherlands

3. Division of Cardiothoracic Surgery, Department of Surgery, University of Washington School of Medicine, 660 South Euclid Avenue, 63110 St Louis, MO, United States of America

4. Univ. Bordeaux, Inserm, Centre de Recherche Cardio-Thoracique de Bordeaux U1045, IHU Liryc, Hôpital Xavier Arnozan, Avenue du Haut Lévêque, 33600 Pessac, France

Correspondence: [email protected]

DOI: 10.1371/journal.pone.0303674

Editor: Rafael Sachetto Oliveira, Universidade Federal de Sao Joao del-Rei, Brazil

Keywords:

• excitable media
• cardiac electrophysiology
• open-source
• finite differences
• anisotropic diffusion
• reaction-diffusion equation
• C++

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.

Popular summary:

We present Ithildin, an open-source library for reaction-diffusion systems such as the electrical waves in cardiac tissue controlling the heart beat. We demonstrate the versatility of Ithildin by example simulations in various cell models and geometries, from simple 2D simulations to detailed ones in ventricular geometry. Our simulations highlight key features of Ithildin, such as recording pseudo-electrograms or filament trajectories. We hope that our work will contribute to the growing understanding of cardiac electrophysiology and inform future experimental and theoretical studies.

With the Ithildin framework, we want to open up new gateways in the numerical simulation of reaction-diffusion systems, such as the electrical activation patterns in the heart.¹

Our motivation to write a reaction-diffusion solver comes from the numerical study of electrical patterns inside the heart (Clayton et al., 2011). These patterns, which are incompletely understood, are a main cause of death and even as a chronic disease, they complicate people’s lives. In the past decades, computer models of arrhythmia have allowed mechanistic insight in the origin and control of arrhythmias (Fenton & Karma, 1998). On the longer term, it is thought that digitized versions of patients' hearts could help offer better diagnostics and planning of procedures; such personalized heart models are called cardiac digital twins (Gillette et al., 2021; Koopsen et al., 2024; Niederer et al., 2019; Trayanova et al., 2020). In view of open science, we have decided to share the code that has been steadily developed in our group since 2007 with the scientific community.

A flowchart outlining the functionality of Ithildin can be found in Fig. 1. Ithildin is designed to comply with the 2011 version of the ISO-C++ standard (ISO, 2011), but it compiles with all newer versions, including the current 2023 ISO-C++ standard (ISO, 2014, 2017, 2020, 2023). The software facilitates forward Euler and Runge-Kutta finite-difference solutions for reaction-diffusion systems in $N$-dimensional space, such as the monodomain equation for cardiac electrophysiology, with specified boundary conditions (Clayton et al., 2011).

The framework offers quick computation through CPU parallelization using OpenMPI (Graham et al., 2006). It also boasts decent documentation of available features, made accessible through Doxygen (van Heesch, 2023). Ithildin writes easy-to-parse YAML log files to document simulation setups (döt Net et al., 2023). Additionally, it allows convenient output of frames of recorded variables at regular intervals in the form of NumPy NPY files (Harris et al., 2020). The software supports easy and powerful post-processing with the Python module for Ithildin (Kabus et al., 2022), including integration with Scientific Python (Virtanen et al., 2020), 2D visualization with Matplotlib (Hunter, 2007), and 3D visualization with ParaView (Ahrens et al., 2005).

Ithildin also allows the recording of pseudo-electrograms (EGMs) and state variables at full numerical time resolution, as well as the tracking of filaments, which represent the instantaneous rotation axes of rotors. The software features a flexible setup for in-silico experiments, also called simulations, through a simple class-based C++ interface.

Various types of geometries are implemented, ranging from a simple 1D cable and spirals in 2D tissue to whole-heart geometry and even 4D hyperspace. The space can be partitioned to use multiple cell models in the same experiment via Model_multi. Realistic stimulation protocols can be added as Stimulus objects and may be started by a Trigger. Ithildin also includes a logging system with minimal impact on computation speed and various levels of verbosity.

In this paper, we provide an overview of this framework, guiding the reader through its components. Results from several in-silico experiments are presented as the main components of Ithildin are introduced. The details of these so-called simulations are outlined towards the end of this paper in section 5, along with a tabular overview in Table 4.

The diffusion of electrical signals in cardiac tissue can be modelled as a reaction-diffusion system, where the diffusion tensor ${{\bm{{D}}}}$ represents the anisotropic properties of the medium. This tensor encapsulates the spatial orientation of fibers in the medium and the effects of inhomogeneities on signal propagation. With the local unit vectors along the fibers ${{\bm{{e}}}}_{\text{f}}$, normal to fibers in the sheet plane ${{\bm{{e}}}}_{\text{s}}$, and normal to both of these ${{\bm{{e}}}}_{\!\times\!} = {{\bm{{e}}}}_{\text{f}} \times {{\bm{{e}}}}_{\text{s}}$, forming a orthonormal basis, the fiber orientation is encoded using diffusivities $D_{{\text{f}}, {\text{s}}, {\!\times\!}}{{\left( {{\bm{{x}}}} \right)}}$ in each of these directions (Clayton et al., 2011): $$ {{\bm{{D}}}} = D_{\text{f}} {{\bm{{e}}}}_{\text{f}} {{{{\bm{{e}}}}_{\text{f}}}^\mathrm{T}} + D_{\text{s}} {{\bm{{e}}}}_{\text{s}} {{{{\bm{{e}}}}_{\text{s}}}^\mathrm{T}} + D_{\!\times\!} {{\bm{{e}}}}_{\!\times\!} {{{{\bm{{e}}}}_{\!\times\!}}^\mathrm{T}} \qquad{(1)}$$

The core equation governing the evolution of the state variable vector ${{\underline{{u}}}}$ is the reaction-diffusion equation: $$ \partial_t {{\underline{{u}}}} = {{\underline{{P}}}} \nabla \cdot {{\bm{{D}}}} \nabla {{\underline{{u}}}} + {{\underline{{r}}}}{{\left( {{\underline{{u}}}} \right)}} \qquad{(2)}$$ or in index notation: $$ \partial_t u_m{{\left( t,{{\bm{{x}}}} \right)}} = \textstyle\sum_{m'} P_{mm'} \; \textstyle\sum_{n} \partial_{n} \textstyle\sum_{n'} D_{nn'}{{\left( {{\bm{{x}}}} \right)}} \; \partial_{n'} u_{m'}{{\left( t, {{\bm{{x}}}} \right)}} + r_{m}{{\left( {{\underline{{u}}}};{{\bm{{x}}}} \right)}} \qquad{(3)}$$ for $m, m' \in {{\left\{1,...,M \right\}}}$ with the number of state variables $M$, and $n, n' \in {{\left\{1,...,N \right\}}}$ with the number of spatial dimensions $N$, using the notation $\partial_n = \partial_{x_n}$ for the spatial partial derivatives.

Here, ${{\underline{{u}}}}$ is the state variable vector and ${{\underline{{r}}}}({{\underline{{u}}}})$ accounts for the reaction term and is called the model. We refer to the first component of ${{\underline{{u}}}}$ as $u$, which for electrophysiological models is the transmembrane voltage $V_\text{m}$ or a rescaled version of it, see also section 4.2. For two-variable models, the second component of ${{\underline{{u}}}}$ is often referred to as the restitution or recovery variable $v$. In the term representing diffusion, ${{\bm{{D}}}}$ is determined by the geometry of the medium and the presence of inhomogeneities, see section 4.1. The projection matrix ${{\underline{{P}}}}$ is typically a diagonal matrix describing whether or not a variable is diffused. For instance, only the first variable of the AP96 model (Aliev & Panfilov, 1996) is diffused, such that ${{\underline{{P}}}} = {\operatorname{diag}{{\left( 1, 0 \right)}}}$.

Different notation is used to distinguish between vectors ${{\bm{{x}}}}$ and matrices ${{\bm{{D}}}}$ in physical space in bold font, and underlined vectors ${{\underline{{u}}}}$ and matrices ${{\underline{{P}}}}$ with respect to state variables. We use lowercase letters for vectors and uppercase for matrices. An overview of the most relevant quantities in Ithildin is given in Table 1.

Ithildin obtains approximate solutions of the reaction-diffusion equation (Eq. 2) via a finite-differences approach: Time $t$ and space ${{\bm{{x}}}}$ are discretized on a grid and values ${{\underline{{u}}}}{{\left( t, {{\bm{{x}}}} \right)}}$ are associated with the vertices of this grid. We choose a fixed temporal resolution, the time step $\Delta t$, and constant spatial grid spacing $\Delta {{\bm{{x}}}}$. The values ${{\underline{{u}}}}{{\left( t + \Delta t, {{\bm{{x}}}} \right)}}$ at a subsequent time-step are computed based on the previous ones, according to discretized versions of the governing equations, i.e., the reaction-diffusion equation (Eq. 2), together with boundary and initial conditions.

Starting from an initial state, the state variable vector ${{\underline{{u}}}}$ is integrated over time using a so-called time stepping scheme leading to an approximate solution of the reaction-diffusion system using finite differences. Ithildin implements two main stepping schemes to choose from: forward Euler and the classic Runge-Kutta method (RK4) (Euler, 1794; Press, 2007).

Defining ${{\underline{{f}}}}$ as the right hand side of the reaction-diffusion equation (Eq. 2), the forward Euler method takes the form (Press, 2007): $$ {{\underline{{u}}}} {{\left( t + \Delta t, {{\bm{{x}}}} \right)}} = {{\underline{{u}}}} {{\left( t, {{\bm{{x}}}} \right)}} + \Delta t \, {{\underline{{f}}}} {{\left( t, {{\bm{{x}}}}; {{\underline{{u}}}} \right)}} + O{{\left( \Delta t^2 \right)}} \qquad{(4)}$$ This method is the default time integration scheme in Ithildin. Despite its numerical error being of order $O{{\left( \Delta t^2 \right)}}$, with a sufficiently small time step, the accuracy of the Euler method is adequate for our use cases.

Due to the Courant-Friedrichs-Lewy condition (CFL), a stability criterion for the integration of the reaction-diffusion equation, $\Delta t$ needs to be chosen sufficiently small (Courant et al., 1928). Ithildin automatically chooses an appropriate time step based on the CFL condition for the different supported geometries, cf. section 4.1. For example, for the most simple implemented geometry contained in Ithildin, i.e., isotropic diffusion (see Geometry_Iso in section 4.1 and Table 2), the CFL condition is enforced by setting (Li et al., 1994): $$ \Delta t < {{\left[ 2\max{{{\underline{{P}}}}} \sum_{n=1}^{N} \frac{D_{nn}}{\Delta x_n^2} \right]}}^{-1} \qquad{(5)}$$ where $D_{nn}$ are the diagonal components of ${{\bm{{D}}}}$ and $\max{{{\underline{{P}}}}}$ is the maximum value of ${{\underline{{P}}}}$.

Higher accuracy at the cost of more computations per time step, can be achieved with the RK4 method (Press, 2007): $$\begin{aligned} {{\underline{{u}}}} {{\left( t + \Delta t, {{\bm{{x}}}} \right)}} &= {{\underline{{u}}}} {{\left( t, {{\bm{{x}}}} \right)}} + \frac{1}{6}\Delta{{\underline{{u}}}}_1 + \frac{1}{3}\Delta{{\underline{{u}}}}_2 + \frac{1}{3}\Delta{{\underline{{u}}}}_3 + \frac{1}{6}\Delta{{\underline{{u}}}}_4 + O{{\left( \Delta t^5 \right)}} \\ \Delta{{\underline{{u}}}}_1 &:= \Delta t\, {{\underline{{f}}}} {{\left( t, {{\bm{{x}}}}; {{\underline{{u}}}} \right)}} \\ \Delta{{\underline{{u}}}}_2 &:= \Delta t\, {{\underline{{f}}}} {{\left( t + \textstyle\frac{1}{2}\Delta t, {{\bm{{x}}}}; {{\underline{{u}}}} + \textstyle\frac{1}{2}\Delta{{\underline{{u}}}}_1 \right)}} \\ \Delta{{\underline{{u}}}}_3 &:= \Delta t\, {{\underline{{f}}}} {{\left( t + \textstyle\frac{1}{2}\Delta t, {{\bm{{x}}}}; {{\underline{{u}}}} + \textstyle\frac{1}{2}\Delta{{\underline{{u}}}}_2 \right)}} \\ \Delta{{\underline{{u}}}}_4 &:= \Delta t\, {{\underline{{f}}}} {{\left( t + \Delta t, {{\bm{{x}}}}; {{\underline{{u}}}} + \Delta{{\underline{{u}}}}_3 \right)}} \end{aligned}\qquad{(6)}$$ Note that ${{\underline{{f}}}}$ needs to be evaluated four times for the RK4 method and only once for the Euler method. While RK4 is still an explicit method subject to instability at too large $\Delta t$, a larger $\Delta t$ value than for the Euler method can typically be used.

The stepping scheme to be used can be chosen in Ithildin on a per-variable level via Model::steppings.

For the numerical solution of the reaction-diffusion equation (Eq. 2), the spatial derivative in its right hand side must be computed, i.e., the diffusion operator ${{\underline{{P}}}}\nabla\cdot{{\bm{{D}}}}\nabla{{\underline{{u}}}}$. This is implemented as weighted sums of the value of ${{\underline{{u}}}}$ at neighboring vertices on the grid of the discretized domain. The weights for the calculation of the stencil depend on the chosen type of diffusion (section 4.1). The two main types in this software are a first order stencil, including only the nearest neighbors, and a second order stencil, including also the next to nearest neighbors.

In the simplest case (see Geometry_Iso in section 4.1 and Table 2), we consider isotropic and homogeneous diffusivity. The diffusion operator can then be computed via the Laplacian operator $\nabla^2$. This is done with a 5-point stencil for the 2D case, a 7-point stencil for the 3D case, etc., see also Fig. 2, panel (a). Consequently, the weights are calculated as: $$ \begin{aligned} w_n &= \frac{1}{\Delta x_n} \quad \text{for }1\leq n\leq 2 N \\ w_0 &= -\sum_{n=1}^{2N}w_n \end{aligned} \qquad{(7)}$$ where $N$ is the number of dimensions and $\Delta x_n$ is the grid resolution in the direction of the neighbor corresponding to weight $w_i$, e.g., $\Delta x_5 = \Delta z$. Note that the indices correspond to those displayed in Fig. 2.

For the more general orthotropic diffusion, the stencil for numerical differentiation includes the nearest and diagonal neighbors of a grid point, see Fig. 2, panel (b). The weights for orthotropic diffusion are obtained by a combination of central differences approximations to derivatives and linear interpolation of values between two grid points. For more details, the reader is referred to the documentation of Geometry_OrtAniso.

The source code of Ithildin is written in the C++ programming language, following the 2011 version of the ISO-C++ standard (ISO, 2011). This was chosen to facilitate programming at a level relatively close to the hardware, but also using some of the useful data-structures that are contained in the standard template library (STL).

Ithildin is designed to run in parallel on multiple CPU cores using MPI, specifically OpenMPI (Graham et al., 2006). Using the mpirun command, multiple instances, also known as processes, of the same compiled executable are started that run across different processor cores. The processes are then coordinated such that there is one so-called manager process, that manages a bunch of so-called worker processes. The manager does an equal share of the computational work, just like a worker. The only difference is that the manager distributes and directs information to be exchanged from one process to another. In Ithildin, the memory is not shared across processes, instead each process works on its share of the computational domain. We split the domain in the $x$-direction, such that each process is responsible for computations on a roughly equal share of vertices inside the to-be-simulated medium. Even when some parts of the domain are classified as exterior points, i.e., points on which no calculations need to be performed, they are taken into account when splitting the domain between processes. This splitting is possible because the reaction-diffusion systems to be studied with Ithildin are local, meaning that the temporal evolution at each time $t$ and each point ${{\bm{{x}}}}$ in space depends only on the current state vector ${{\underline{{u}}}}{{\left( t, {{\bm{{x}}}} \right)}}$ at that position and its spatial derivatives (Eq. 2). Internally, a layer of so-called ghost points is added around each process’ part of the domain such that the spatial derivatives can still be calculated in the same way as for any other point in the domain (section 2.3). The values ${{\underline{{u}}}}$ on these ghost points are exchanged with the neighboring processes, as coordinated by the manager process. Additional ghost points are also used to enforce Neumann boundary conditions, which is done by setting the weights for the calculation of the numerical spatial derivatives accordingly.

To run a simulation in Ithildin, the user needs to define a main() function that is to be called by all subprocesses. This is usually done using a C++ file calling the required components of the Ithildin library defining the main() function, a so-called main file. Ithildin can be installed as a shared library on Unix-based systems that is required by executables obtained from compiling main files. Alternatively, it is also possible to statically compile the Ithildin library with a main file into a stand-alone executable.

Besides the necessary preparations for using MPI, running a simulation using Ithildin’s Sim class requires three main components, which are instances of three classes:

1. Model: the reaction term ${{\underline{{r}}}}{{\left( {{\underline{{u}}}} \right)}}$, typically a cell model,
2. Geometry: the discretized diffusion term ${{\underline{{P}}}} \nabla \cdot {{\bm{{D}}}} \nabla {{\underline{{u}}}}$ for a chosen geometry of the medium, and
3. Source: the stimulus protocol to use as well as inhomogeneities in the medium.

In the following section 4, an overview is given for each of these classes and their derived classes. Combining all of the components, an illustrative, minimal main file can be obtained, in which a planar wave crosses the medium in positive $x$-direction:

#include "ithildin.h"

int main(int argc, char** argv){
    // 0. initialize MPI
    Mpiclass mpi(argc, argv);

    // 1. select reaction term
    Model_SmooKa model;

    // 2. select diffusion term
    vector<int> size{30, 30, 1};
    vector<float> dx{1., 1., 1.};
    Geometry_Iso geom{&mpi, size, dx, &model};

    // (set up simulation)
    Sim sim{1., 10, &mpi, &model, &geom, "example"};

    // 3. define stimulus protocol
    Source sour{&model, &geom, &mpi, &sim};
    sour.stimulate({{0}, {1.}, Shape::Rect({0, 0}, {5, 0})});

    // (run simulation)
    return sim.run(&sour);
}

More detailed examples for main files can be found in the appendix, where we set up the numerical examples used throughout this paper.

We consider the geometry, the model and the source to be the inputs of Ithildin, see also Fig. 1. Upon running the simulation, Ithildin produces a variety of outputs, as files in easy-to-parse standardized data formats. The names of these files all begin with the so-called stem, consisting of a descriptive series name and a serial number, which defaults to a time-stamp. In the following, an overview of the usual output files of Ithildin by suffix appended to the stem is provided:

• _log.yaml: The log file contains metadata describing the setup of the simulation, as well as metadata about the conditions under which the simulation was run. This file is always output by Ithildin and is considered the central file of the results, as it points to the relevant other files that are only conditionally written during the simulation. While earlier versions of Ithildin used a non-standard format for the log files, in current versions, the YAML format is used, making it easy to parse by both: humans and machines (döt Net et al., 2023).
• _main.cpp: A copy of the C++ code in the main file may be included in the results for reproducibility.
• _git.diff: If run in a Git-repository with changes since the last commit, a patch file of these changes will be included in the results.
• _*.txyz.npy: For each of the state variables ${{\underline{{u}}}}$, a so-called var file in the NumPy NPY format will be written (Harris et al., 2020). These files contain the $N+1$-dimensional floating-point number array of the evolution of a state variable $u$ in the whole computational grid over time. Note that the order of indices is $(t, x, y, z)$, meaning that time is the slowest varying index and the $x$-axis the second slowest varying index. This order was chosen because the domain is split across processes along the $x$-axis, such that the processors can open and write to these files sequentially. For 2-dimensional simulations, the third spatial dimension, the $z$-axis, is one vertex thick, such that each var file still has four dimensions. For higher-dimensional simulations, more axes are added, leading to more dimensions in the var files.
• _inhom.txyz.npy: The inhom field, see details in section 4.3, is also stored in the NPY format, but with integer values, and only for the initial time-step. This file hence has the shape $(1, N_x, N_y, N_z)$ with $N_n$ denoting the number of vertices in each of the spatial dimensions.
• _hist*.csv: Comma-separated value (CSV) files describing the temporal evolution of ${{\underline{{u}}}}{{\left( t,{{\bm{{x}}}}_\text{s} \right)}}$ at a chosen sensor position ${{\bm{{x}}}}_\text{s}\in\Omega\subset\mathbb{R}^{N}$, see details in section 4.5.
• _egm*.csv: CSV files containing the recorded pseudo-EGM $\Phi{{\left( t,{{\bm{{x}}}}_\text{e} \right)}}$ at a chosen electrode location ${{\bm{{x}}}}_\text{e}\in\mathbb{R}^{N}$, see details in section 4.6.
• _tipdata.yaml: This YAML file is written if filament-tracking is turned on, see also section 4.7, and contains the detected phase singularities in regular time intervals.

While we have selected these file formats to be easy to read using a wide variety of software, we have also developed the Python module for Ithildin to facilitate interacting with the results of an Ithildin simulation and converting them to a variety of file formats, for instance writing files in the extensible data model and format (XDMF) that can be used to view simulation results in ParaView (Ahrens et al., 2005; Kabus et al., 2022; Kabus, De Coster, et al., 2024; Mark et al., 2007). The Python module also offers post-processing and analysis methods, for instance the computation of action potential duration (APD), conduction velocity (CV), various phases, phase defect detection, functions acting on filaments and filament trajectories, and several plotting functions (Kabus et al., 2022; Kabus, De Coster, et al., 2024).

The diffusion term in the reaction-diffusion equation (Eq. 2) is stated as ${{\underline{{P}}}} \nabla \cdot {{\bm{{D}}}} \nabla {{\underline{{u}}}}$. The conduction in cardiac tissue and hence the diffusion is stronger along the fiber direction than normal to the fibers (Clayton et al., 2011). This is encoded in the diffusion matrix ${{\bm{{D}}}}$.

As an example, in Fig. 3, the fiber direction is drawn on the surface of the ventricular geometry used in Sim. 4. The fibers are additionally colored by their fiber helix angle (Streeter et al., 1969). The voxel-based representation of this geometry was obtained by cutting a human heart into slices, digitizing and stacking them (Hren, 1996; ten Tusscher et al., 2007).

While the projection matrix ${{\underline{{P}}}}$ is managed by the Model class, see section 4.2, the diffusion matrix ${{\bm{{D}}}}$ and the handling of the spatial derivatives is implemented in the Geometry class. The most important things the Geometry class takes care of, are:

• Initialization of the computational grid, along with the strides and pointers, which are important for efficient computing;
• Computation of the entries of the diffusion tensor, based on the main directions of diffusion and the respective diffusion values;
• Computation and storage of the weights for the stencils of the numerical spatial derivatives (section 2.3) respecting the Neumann boundary conditions; and
• Handling the upper bound for the time step due to the CFL condition (section 2.2).

The Geometry class is a base class and should not be used directly to construct a computational domain. Instead, there are several subclasses, each representing a different type of diffusion or extrinsic shape. An overview of the Geometry subclasses is given in Table 2.

The Geometry_ND class is a subclass of Geometry_Iso and the Ellipsoid, Hyperboloid, and Paraboloid classes are subclasses of QuadricSurface.

The details on how the extrinsic and intrinsic curvature affect the diffusion tensor and hence the weights for the discretized differential operator, are discussed in the documentation of the code. There, the different ways to initialize the available domain types are given as well.

Note that some additional features are implemented in the Geometry class, that do not have a direct link with the diffusion term, such as inhomogeneities (section 4.3) and filament detection (section 4.7).

The Model class serves as the base class for all cardiac electrophysiology models in the Ithildin framework. It encapsulates common functions and variables used across different models. Key features of the Model class include:

• Implementation of the reaction term ${{\underline{{r}}}}{{\left( {{\underline{{u}}}} \right)}}$ in the reaction-diffusion equation.
• Storage of model-specific metadata such as relevant citations.
• Handling of variable-related information, including their names, indices, and resting values.
• Management of the values of the projection matrix ${{\underline{{P}}}}$.

Derived classes extend the Model class to implement specific cardiac electrophysiology models in the reactionterm function.

An overview of the cell models that are currently included in the source code of this project is provided in Table 3. More models may be added as additional classes.

The Ithildin framework also introduces several ModelWrapper classes that provide additional functionality and allow for the combination of models. These wrappers enable the recording of diffusion terms, reaction terms, local activation times (LAT), and local deactivation times (LDT) as additional state variables. This is done by adding code to the reaction term calculations of the underlying model. The wrappers inherit from the base ModelWrapper class, which is a wrapper that leaves the model unchanged. The primary ModelWrapper classes are:

• ModelWrapper_RecordDiffusion records diffusion terms of selected variables as additional variables.
• ModelWrapper_RecordReaction records reaction terms of selected variables as additional variables.
• ModelWrapper_RecordActivationTime records LAT for selected variables.
• ModelWrapper_RecordDeactivationTime records LDT for selected variables.
• ModelWrapper_RescaleTimeSpace linearly rescales the wrapped model in time and space.
• ModelWrapper_RescaleVars linearly rescales selected state variables of the wrapped model.

The Model_multi class enables the combination of multiple submodels into a single model. The behavior of the combined model may vary depending on the location ${{\bm{{x}}}}$ and is determined by one of its submodels. Which model is to be used depends on the integer value of the inhom field, which is used to describe spatial inhomogeneities, such as obstacles. Inhomogeneities will be further explained in the following, cf. section 4.3. This class is particularly useful for simulating scenarios where different regions of cardiac tissue exhibit distinct behaviors.

The Ithildin C++ framework provides a structured and modular approach to modeling cardiac electrophysiology. The Model class serves as the base for different models, while various ModelWrapper classes and the Model_multi class offer extended functionalities for recording and combining different model aspects.

In simulations of cardiac electrophysiology, accurately modeling the spatial properties of the cardiac tissue is essential. Inhomogeneities represent variations in the tissue’s characteristics, such as its electrical conductivity or cellular properties, that influence the propagation of electrical signals.

An inhomogeneity in the Ithildin framework is a distinct region within the simulation domain with different properties compared to its surroundings. In the context of cardiac electrophysiology, these properties could correspond to variations in the electrical conductivities of cells, cellular properties, or even the absence of excitable cells altogether. Inhomogeneities are defined by an integer field called inhom associated with each point in the domain. Grid points with a non-zero inhom value are considered interior points, indicating that the reaction-diffusion equation needs to be solved on these points. If the selected reaction term is a Model_multi (see also section 4.2), for an inhom value of $n$, the $n$th submodel will be used for this point. For example, in Fig. 4, the inhom field for Sim. 1 is visualized. Two different cell models are used for the values 1 and 2. The points where inhom has the value 0, are considered exterior points. The set of all interior points is the physical domain $\Omega$. At the boundary of the physical domain, Neumann boundary conditions are applied.

To incorporate inhomogeneities into the simulation, the framework provides methods to add them to the domain. These methods allow specifying the shape, location, and properties of each inhomogeneity. For example, a rectangular inhomogeneity could be added by specifying its width, height and location.

Ithildin provides a flexible way to define stimulation protocols using the Stimulus and Trigger classes, as well as the scheduling functionality of the Source class.

The Stimulus class enables the specification of temporal and spatial characteristics of voltage-based or current-based stimuli. It allows defining which state variables should be affected directly, the associated values, and whether the stimulus sets the variable directly or whether it is additive and hence behaving like a current source. Temporal modulation of stimulus strength is achieved through the amplitude function, while spatial constraints are managed by the shape function, an instance of the Shape class.

The Shape class describes a geometric shape via its characteristic function which can be used to define regions of interest in a simulation domain. The core principle is that the function evaluates a given position vector and returns a value: 1 if the position is inside the defined shape, and 0 if it is outside. However, values in the range ${{\left[ 0,1 \right]}}$ may also be used to create a smooth transition. A smoothly varying characteristic function may be useful to create more-realistic stimuli that deposit current in a smooth profile. This simple yet powerful concept forms the basis for constructing intricate spatial configurations.

The Shape class provides several pre-defined shape functions, though additional shapes can easily be added by defining a characteristic function:

• Ellipsoid: Defined by radii, a center and optionally the Euler angles, this shape represents a general three-dimensional ellipsoid with the specified orientation.
• Sphere: A special case of an ellipsoid where all radii are equal, forming a three-dimensional sphere.
• Ellipse in $xy$-plane: This two-dimensional shape resembles an ellipse lying on the $xy$-plane, defined by radii and a center.
• Cylinder along $z$-axis: Representing a three-dimensional cylinder centered along the $z$-axis, this shape is defined by a radius and a center. In 2D, it defines a disk.
• Rectangular cuboid: Defining a three-dimensional region, this shape is specified by two opposing corner points, creating a cuboid. In 2D, it defines a rectangle.
• Half plane: A plane defined by an origin and an outward normal vector splits the three-dimensional space into a half-space. In 2D, a straight line splits the plane in a similar way.

Characteristic functions can also be loaded from files in the NPY format (Harris et al., 2020). This feature facilitates the incorporation of custom shapes derived from external data sources.

The scheduling functionality via Source::schedule offers the ability to execute functions at specified points in time during the simulations. This feature greatly enhances experimental flexibility by allowing the execution of arbitrary code snippets at chosen moments during simulations. The scheduler is especially useful for introducing dynamic changes to the simulation environment, such as modifying stimuli or conditions mid-simulation.

The Trigger class provides a means to orchestrate actions based on specific conditions. Triggers encapsulate the decision-making process of when to execute a particular action, influenced by condition checks and coordination modes. Different coordination modes allow for the synchronization of trigger actions across multiple processes, facilitating complex simulations of activation waves. These modes are to trigger on each process individually once the condition is met, once the condition is met in a specific process, once the condition is met in any process, or once it is true in all processes.

Within the framework of cardiac electrophysiology, triggers are essential for defining stimulation protocols, e.g. for the S1S2 protocol, which is illustrated in Fig. 5: After a first excitation wave passes a sensor position, a second wave is triggered behind a part of the first waveback to stimulate spiral waves.

Besides the sensors for triggering stimuli, sensors in the context of this simulation framework are components that monitor and record the state variables of the simulated system at chosen positions during the simulation. We call this the history at a given sensor position.

These sensors are used to gather data about the behavior of the system at particular time intervals, regulated by the sensorlag parameter. This parameter controls the frequency at which sensor data is collected, allowing for flexibility in recording intervals. Notably, the recording frequency set by sensorlag need not align with the simulation’s time step or the duration between frames. The data will be recorded at the first time step after the specified sensorlag duration. This high-resolution temporal data can be used to study individual points in the medium in detail.

The recorded data are then written to designated comma separated value (CSV) output files associated with each sensor. These files are used to store the collected data over the course of the simulation.

The first four panels of Fig. 6 contain time traces of the transmembrane voltage $u$, the restitution variable $v$, the recorded value of the diffusion term, and the LAT at a given sensor position for Sim. 1. The times of the three stimuli are indicated by the black vertical lines. The other panels will be explained in the subsequent sections.

The EGM is a measurement of the potential generated by the charge distribution in cardiac tissue over time at a point in space, outside the tissue. In theory, this is a measurement of the extracellular potential ${\Phi_\text{e}}$. However, since Ithildin is a monodomain solver, the extracellular potential is not part of the model equations (Clayton et al., 2011). The code therefore calculates an approximation of the extracellular potential at points outside the mesh using the Egm class. To obtain this approximation, the pseudo-bidomain theory is used (Bishop & Plank, 2011). The approximation, referred to as a pseudo-EGM, uses the simplifications that the intracellular and extracellular conductivities are proportional, such that there is an explicit formula to calculate the extracellular potential, and that the bath conductivity is homogeneous.

The extracellular potential ${\Phi_\text{e}}$ at position ${{\bm{{x}}}}_\text{e}$ over time is then calculated as: $$ {\Phi_\text{e}}{{\left( {{\bm{{x}}}}_\text{e},t \right)}} = \frac{ \tau_\text{e} }{ 4\pi } \int_{\Omega}\mathrm{d}{{\bm{{x}}}}\; \frac{ \nabla \cdot {{\bm{{D}}}} \nabla u({{\bm{{x}}}},t) }{ {{\left\lVert {{\bm{{x}}}}_\text{e} - {{\bm{{x}}}} \right\rVert}} } \qquad{(8)}$$ where $\Omega$ denotes the computational domain, i.e., the simulated heart muscle tissue, $u$ is a state variable of the model representing the transmembrane potential, , which is usually encoded as the first variable in the state vector ${{\underline{{u}}}}$. Furthermore, ${{\bm{{D}}}}$ is the diffusion matrix from the reaction-diffusion equation (Eq. 2), and $\tau_\text{e}$ is a proportionality factor with units ${\mathrm{{m}{s}}}$.

Essentially, the used approximation for the EGM is a convolution of the diffusion term $\nabla \cdot {{\bm{{D}}}} \nabla u$ for the first variable with a kernel ${{\left\lVert {{\bm{{x}}}}_\text{e} - {{\bm{{x}}}} \right\rVert}}^{-1}$. Depending on the choice of the diffusion tensor, made by the user in the Geometry class (section 4.1), a different prefactor of the integral is required. Hence, the prefactor $\frac{\tau_\text{e}}{4\pi}$ is user-defined and can be changed using the function set_prefactor.

The integral in Eq. 8 is discretized and its calculation is implemented such that the additional amount of storage and number of calculations is limited as much as possible. For instance, as the required diffusion term $\nabla \cdot {{\bm{{D}}}} \nabla V_\text{m}$ is already computed during the forward stepping of the reaction-diffusion equation, it can be stored as an additional state variable using a ModelWrapper_RecordDiffusion and subsequently used in the pseudo-EGM calculation. This model wrapper is essential for the functioning of the Egm class and hence is a requirement when setting up a simulation with pseudo-EGM calculation.

The result is a CSV file with the pseudo-EGM data for each electrode that is defined by the user. The computed pseudo-EGM at a given position for Sim. 1 is displayed in the fifth panel of Fig. 6.

Formally, filaments can be understood as a line of wave break, i.e., a line where an activation and recovery surface come together (Clayton et al., 2006; Winfree, 1994). The activation surface can be seen as the wavefront, while the recovery surface can be seen as the waveback. When considering an excitable system in 2D, filaments become tips, being the point of intersection between the activation and recovery curve. Since a point cannot be excited and recovering at the same time, points on a filament are also called phase singularities.

Using this definition of filament points, detection algorithms have been designed. Our code relies on the one described by Fenton & Karma (1998). Additionally, this algorithm has been extended to grids in any dimension, where the generalization of filaments are called superfilaments (Cloet et al., 2023).

The filament point detection algorithm is included in the Geometry class. Its goal is to compute the points where the wavefront and waveback meet. While looping over all coordinate planes and grid points, it is checked whether there is an intersection of isolines in the adjacent voxel faces of a grid point. The location of the filament point is then estimated by bi-linear interpolation. An illustrative sketch of this method is given in Fig. 7.

The last panel in Fig. 6 displays the number of rotors over time for Sim. 1, which are found via filament detection. It can be seen that at the S2 stimulus, a single rotor is formed, and subsequently pairs of rotors as figure-of-eight spiral pairs.

More details of this process can be seen in Fig. 4 which shows the trajectories of these filaments in Sim. 1 tracked using the Python module for Ithildin. The trajectories are colored by time and the formation of a new tip is denoted by stars and their decay by crosses. Tips that still persist at the final frame of the simulation are indicated by points. The tip trajectory denoted with index 0 is formed by the S1S2 protocol and meanders around the medium. It persists until the end of the simulation. The figure-of-eight spiral wave pair denoted by indices 1 and 2 is formed by a conduction block breaking up. The spiral tip with index 2 runs into the boundary to disintegrate there.

While a phase singularity can be seen as points where all phases meet—in mathematics called a pole, a phase defect is a point at which there is a discrete jump from one phase value to another in an otherwise continuously varying phase. While phase defects are well known in physics, they were only recently identified in excitable media, within linear-core rotors and conduction block regions (Arno et al., 2021; Tomii et al., 2021). Phase defects are lines in 2D and surfaces in 3D.

In Ithildin, the phase defect detection is done with its Python module. For instance, during simulation, Ithildin can record the local activation times (LAT) which can then be used to compute the activation time phase in post-processing (Arno et al., 2021; Kabus et al., 2022). A variety of methods exist to then localize the phase defect (Kabus et al., 2022; Tomii et al., 2021).

Two examples for phase defect detection in Sim. 1 are given in Fig. 8 and Fig. 9. Both display four frames over time of the transmembrane voltage $V_\text{m}$, the activation time phase ${\varphi}$, and the phase defect $\rho$ computed via the cosine method (Kabus et al., 2022; Tomii et al., 2021). In Fig. 8, the phase defect of a single spiral wave is tracked. It can be seen that the phase defect extends due to conduction block such that the spiral wave moves across the domain. In Fig. 9, the break-up of a conduction block line into a figure-of-eight spiral wave pair can be seen. The conduction block line is a phase defect of zero topological charge (Arno et al., 2024a; Goryachev & Kapral, 1996; Mermin, 1979) which, in this case, reaches a critical length breaking apart into two oppositely charged spiral waves with much shorter phase defect lines.

In Fig. 10, we present the final frame of Sim. 4 in ventricular geometry, visualized with ParaView (Ahrens et al., 2005). It is colored by the normalized transmembrane voltage. Both, the classical tip and the phase defect surface are visualized. The spiral waves revolve around the phase defect surfaces.

More complete documentation of Ithildin can be generated with Doxygen (van Heesch, 2023) and can be also found online, see the data availability statement for details. The documentation contains detailed information on how to get started installing and working with ithildin.

Ithildin is used for numerical experiments in various use cases. For instance, it was used to study the structure of the core of rotors emerging in in-silico cardiac-electrophysiology cell models, leading to the description of phase defect lines (Arno et al., 2021, 2024a, 2024b; Kabus et al., 2022). Also, higher-dimensional rotors waves were simulated and the emerging super-filaments were detected using Ithildin (Cloet et al., 2023). In the creation of novel data-driven cell models using state space expansion, Ithildin was used to generate synthetic training data sets (Kabus, De Coster, et al., 2024).

Five simulations were conducted for this paper to illustrate the features of Ithildin. An overview of the simulations is given in Table 4, along with references to the figures that were generated with these data sets while details on their simulation setup can be found in their C++ code in the appendix.

Niederer et al. (2011) proposed a benchmark problem that is now used by the in-silico modelling community to validate and compare cardiac electrophysiology solvers (Campos et al., 2016; Niederer et al., 2011). Ithildin passes the benchmark as implemented in Sim. 5. In the benchmark problem an excitation wave travels through a cuboid-shaped medium following the cell model by ten Tusscher & Panfilov (2006). In Fig. 11, we present the results from the benchmark in a similar way as in the original publication introducing the benchmark (Niederer et al., 2011): Point $P_1={{{{\left[ 0,0,0 \right]}}}^\mathrm{T}}$ is the corner of the cuboid where the stimulus is applied and $P_8={{{{\left[ 20,7,3 \right]}}}^\mathrm{T}}{\mathrm{{m}{m}}}$ the furthest-away opposite corner. Consider the plane going through $P_1$ and $P_8$ as well as $P_{10}={{{{\left[ 0,7,1.5 \right]}}}^\mathrm{T}}{\mathrm{{m}{m}}}$. We call the distance along the long axis of the plane slicing through the cubioid $\xi_1$ and the short axis $\xi_2$. In the panels (a), the LAT on this plane is shown for the benchmark simulation at high and low resolution in space, $\Delta x = {0.1~\mathrm{{m}{m}}}$ and $\Delta x = {0.5~\mathrm{{m}{m}}}$ respectively, at the same temporal resolution $\Delta t = {0.005~\mathrm{{m}{s}}}$. Panel (b) also shows the LAT on the line from $P_1$ to $P_8$ for the different spatial resolutions $\Delta x \in {{\left\{{0.1~\mathrm{{m}{m}}}, {0.2~\mathrm{{m}{m}}}, {0.5~\mathrm{{m}{m}}} \right\}}}$ at the same temporal resolution. In panel (c), the LAT value at $P_8$ is compared across all the combinations of spatial and temporal resolution. Just like for most other solvers which are compared in the benchmark paper, it can be seen that at the coarsest spatial resolution, the excitation wave is slowed down significantly in the transversal direction (Niederer et al., 2011). Similarly to the other finite-difference solvers, the simulation fails at $\Delta t = {0.05~\mathrm{{m}{m}}}$ and $\Delta x = {0.1~\mathrm{{m}{m}}}$ (Niederer et al., 2011), due to numerical instability. When a check of the CFL condition is enabled (Eq. 5), Ithildin suggests to lower $\Delta t$ to given this spatial resolution, leading to a simulation at which no instability is observed.

Panel (d) of Fig. 11 shows the speed-up in computation time from parallelization by comparing the computation time for Sim. 5 at $\Delta x = {0.2~\mathrm{{m}{m}}}$ and $\Delta t = {0.01~\mathrm{{m}{s}}}$ on an 8-core Intel i7-10875H processor using 1, 2, 4, and 8 processes. In the double-logarithmic plot, it can be seen with linear regression that the computational speed increases almost linearly going from one to four processes, $t_\text{sim} \propto N_\text{proc}^{-p}$ with $p \approx 1$. Going to eight processes leads to diminishing returns, resulting in a smaller speed-up in computational time, as the overhead due to the boundary-exchange between processes grows. On this system using eight processes, at $\Delta x = {0.2~\mathrm{{m}{m}}}$ and $\Delta t = {0.01~\mathrm{{m}{s}}}$, solving the benchmark problem takes ${65.733~\mathrm{s}}$ of computation time in Ithildin, while it takes ${611.231~\mathrm{s}}$ in cbcbeat (Rognes et al., 2017). For both codes, we have turned off output to the disk and included the initial setup of the problem, such as compilation and memory allocation.

The Ithildin framework is a tool for the simulation of cardiac electrophysiology. The code offers a number of assets that allow for simulations targeting a variety of phenomena. For example, there are many instances available to set the local anisotropy of the myocardium via the Geometry class. Please refer to the example presented in section 4.1 for further details. The Geometry class also allows for the simulation and filament tracking in a domain with an arbitrary number of spatial dimensions (Arno et al., 2021, 2024a, 2024b; Cloet et al., 2023; Kabus et al., 2022).

These include the ability to define inhomogeneities, which can be used to model domains of any shape and with any kind of obstacles. This is used in Sim. 1 and Sim. 4. Furthermore, the user can define a variety of stimulation protocols, of which some are showcased in the example simulations. We also note the following options: EGM calculation, the possibility of recording the temporal evolution of variables at sensor positions, and the Python module for post-processing and analysis.

It is evident that Ithildin represents but one instance of software designed for the simulation of cardiac electrophysiology. Another finite-differences solver is BeatBox (Antonioletti et al., 2017), which also relies on domain splitting. A non-exhaustive list of examples of established simulation software based on the finite element method (FEM) are Chaste (Cooper et al., 2020), openCARP (Plank et al., 2021), lifex-ep (Africa et al., 2023), cbcbeat (Rognes et al., 2017), GEMS (Arens et al., 2018), CEPS (CARMEN, 2024), and simcardems (Finsberg et al., 2023). Besides the methods of finite differences, elements, and volumes, approaches from computational fluid dynamics such as the lattice Boltzmann method may be used (Campos et al., 2016). Several of these packages have more advanced features than our software. In addition, some have broader applications than cardiac electrophysiology. More software can be found in (Niederer et al., 2011).

Limitations of our software are the fact that, in the context of cardiac electrophysiology, Ithildin only supports the monodomain model as it is a reaction-diffusion problem conforming with Eq. 2, and the fact that the domain decomposition happens solely in one coordinate direction. Further limitations are that Ithildin is fundamentally voxel-based and hence does not support tetrahedral meshes, and that only Neumann boundary conditions are implemented. Additionally, it is only maintained by a small research team. This paper serves to enhance the visibility of the software and to invite fellow cardiac modelers to use and contribute to the project. We believe that the GitLab environment is an effective medium for facilitating interaction between users, identifying issues and suggesting improvements.

Our work with Ithildin has demonstrated its ability to study complex wave patterns in cardiac tissue (Arno et al., 2021, 2024a, 2024b; Cloet et al., 2023; Kabus et al., 2022; Kabus, De Coster, et al., 2024), but it is important to acknowledge the limitations of the fully explicit numerical scheme used: High resolution in time and space may be required for numerical stability. Nevertheless, there are use cases where Ithildin is practical in the context of clinical applications. For example, it can be used to study rotor waves and their implications for cardiac function, which have important implications for diagnosing and treating cardiac diseases. Additionally, Ithildin can be used to develop new models or simulations that are tailored to specific experimental conditions or clinical scenarios.

In this work, we introduced Ithildin, an open-source library that allows for numerical simulation and analysis of rotor waves. We demonstrated the versatility of Ithildin through a series of simulations, including spiral break-up in the Smooth-Karma model, the S1S2 protocol in the AP96 model, and 2D and 3D spiral waves in the BOCF model in ventricular geometry.

Our simulations highlighted several key features of Ithildin, such as the different implemented geometries and reaction terms, inhomogeneities, and stimuli, as well as recording data such as the pseudo-EGM or filament trajectories. These findings contribute to the growing understanding of rotor waves in cardiac electrophysiology and have the potential to inform future experimental and theoretical studies.

Overall, our work demonstrates the power of Ithildin as a tool for studying complex wave patterns in cardiac tissue. We hope that this library will be useful to researchers seeking to better understand the dynamics of rotor waves and their implications for cardiac function.

The source code of version 3.5.1 of the Ithildin software implemented for this paper is publicly available at https://gitlab.com/heartkor/ithildin and has been archived on Zenodo (DOI: 10.5281/zenodo.12799245). This archive also contains the data generated by the simulations used throughout this paper. Documentation of the code is publicly available at https://heartkor.gitlab.io/ithildin/.

This section contains C++ code to run the four simulations that are used throughout this paper. An overview of their parameters is given in Table 4.

This is the source code to the Ithildin simulation used as the main example throughout this paper.

#include "ithildin.h"

int main(int argc, char** argv){
    Mpiclass mpi(argc, argv);

    // define model
    Model_SmooKa model1;

    // define model for obstacles
    Model_SmooKa model2;
    model2.hidden_eps *= 1.5;

    // rescale time and space
    const float scale_time = 2.4; // ms/t.u.
    const float scale_space = 0.273; // mm/s.u.
    ModelWrapper_RescaleTimeSpace model1a{
        &model1, scale_time, scale_space, 2,
    };
    ModelWrapper_RescaleTimeSpace model2a{
        &model2, scale_time, scale_space, 2,
    };

    // rescale u to voltage, such that: u is in range [0, 1]
    const float ucrit = 0.5;
    const float ustim = 0.75;
    ModelWrapper_RescaleVars model1b{&model1a, "u", 0.25, 0.};
    ModelWrapper_RescaleVars model2b{&model2a, "u", 0.25, 0.};

    // record diffusion term for EGMs
    ModelWrapper_RecordDiffusion model1c{&model1b};
    model1c.record("u");
    ModelWrapper_RecordDiffusion model2c{&model2b};
    model2c.record("u");

    // record local activation time
    ModelWrapper_RecordActivationTime model1d{&model1c};
    model1d.record("u", ucrit);
    ModelWrapper_RecordActivationTime model2d{&model2c};
    model2d.record("u", ucrit);

    // combine models
    Model_multi model{{&model1d, &model2d}};

    // geometry parameters
    const size_t Nx = 100;
    const size_t Ny = 70;
    const float dx = 0.2;
    const float dy = dx;
    const float Lx = Nx*dx;
    const float Ly = Ny*dy;

    // set up geometry
    vector<int> size{Nx, Ny, 1};
    vector<float> deltas{dx, dy, 1.};
    Geometry_Iso geom{&mpi, size, deltas, &model};
    geom.add_inhom(
            Shape::Sphere(0.25*Ly, {0.7*Lx, 0.3*Ly}), 2
    );
    geom.add_inhom(
            Shape::Rect({0.2*Lx, 0.7*Ly}, {0.3*Lx, 0.8*Ly}), 2
    );
    geom.add_inhom(Shape([&](vector<float> p) {
        return p[1] < 0.1*Ly*(1. + 0.7*sin(14.*p[0]/Ly));
    }, "waves"), 0);

    // set up simulation and timing
    const int Nframe = 100; // number of frames
    const float framedur = 20.; // time between frames
    const float sampledur = 2.; // sampling time for time traces
    Sim sim{framedur, Nframe, &mpi, &model, &geom, "everything", 0};

    // filament tracking
    model.define_tip("u", ucrit, "v", 0.95);
    sim.enable_filament_recording(0, sampledur);

    // recording temporal data at a specific location
    sim.add_history_point({0.5*Lx, 0.5*Ly});
    sim.sensorlag = sampledur;

    // record EGMs
    sim.add_egm_electrode({0.5*Lx, 0.5*Ly, 1}, sampledur, "diffusu");
    const float lambda = 0.3; // proportionality factor between \
        // intra- and extracellular conduction matrix
    const float beta = 1e5; // 1/m; area-to-volume ratio
    const float capacitance = 1e-2; // F/m²; specific cell \
        // membrane capacitance
    const float conductivity = 0.5; // S/m; lumped conductivity
    for(Egm* egm : sim.egm_collection) {
        egm->construct_kernel(lambda, beta, capacitance, conductivity);
    }

    // set up stimulus protocol
    Source sour(&model, &geom, &mpi, &sim);
    sour.set_val("u", ustim);

    // S0: initial stimulus
    sour.stimulate(Shape::Rect({0, 0}, {0.1*Lx, Ly}));

    // wait until S1S2 protocol starts
    const float timeS1 = 600.; // ms
    sour.schedule.emplace(timeS1, [&]() {
        // S1: first stimulus of S1S2
        sour.stimulate(Shape::Rect({0, 0}, {0.1*Lx, Ly}));

        // S2: second stimulus of S1S2
        sour.setupS2(Point::Phys(&geom, &mpi, 0.5*Lx, 0.9*Ly), 0, ucrit,
            [&]() -> void {
                sour.stimulate(Shape::Rect({0, 0}, {Lx, 0.5*Ly}));
            }
        );
    });

    return sim.run(&sour);
}

This example illustrates the S1S2 protocol in an easy-to-interpret setting.

#include "ithildin.h"

int main(int argc, char** argv){
    Mpiclass mpi(argc, argv);

    Model_AP model0(0.15, 0.002, 8, 0.2, 0.3, 20);

    // rescale time and space [@aliev1996simple]
    const float scale_time = 12.9; // ms/t.u.
    const float scale_space = 1.; // mm/s.u.
    ModelWrapper_RescaleTimeSpace model{
        &model0, scale_time, scale_space, 2,
    };

    vector<int> size = {120, 120, 1};
    vector<float> dx = {1.0, 1.0, 1.0};
    float Lx = size[0]*dx[0];
    float Ly = size[1]*dx[1];
    Geometry_Iso geom{&mpi, size, dx, &model};

    float framedur = 1.*scale_time;
    int Nframe = 65;
    Sim sim{framedur, Nframe, &mpi, &model, &geom, "s1s2", 0};
    sim.enable_filament_recording(0, framedur);

    Source sour{&model, &geom, &mpi, &sim};
    sour.set_val("u", 0.75);

    // S1
    sour.stimulate(Shape::Rect({0, 0}, {0.1f*Lx, Ly}));

    // S2
    const Point sensor = Point::Phys(&geom, &mpi, 0.6f*Lx, 0.5f*Ly, 0);
    LOG("sensor", sensor.get_phys());
    sour.setupS2(sensor, 0, 0.5, [&]() -> void {
        sour.stimulate(Shape::Rect({0, 0}, {Lx, 0.2f*Ly}));
    });

    return sim.run(&sour);
}

This simulation is used to get an initial state for a 3D simulation in ventricular geometry in the following example.

#include "ithildin.h"

int main(int argc, char** argv){
    Mpiclass mpi(argc, argv);

    // define model
    Model_BO modelbo{1};
    const float ucrit = 0.5;
    const float ustim = 1.0;

    // record local activation time
    ModelWrapper_RecordActivationTime model{&modelbo};
    model.record("u", ucrit);

    // set up geometry
    vector<int> size{450, 450, 1};
    vector<float> dx{0.3, 0.3, 1.};
    const float Lx = size[0]*dx[0];
    const float Ly = size[1]*dx[1];
    Geometry_Iso geom = Geometry_Iso(&mpi, size, dx, &model);

    // set up simulation and timing
    const float framedur = 20.;
    const int Nframe = 34;
    Sim sim{framedur, Nframe, &mpi, &model, &geom, "bocf2d", 0};

    // set up stimulus protocol
    Source sour(&model, &geom, &mpi, &sim);
    sour.set_val("u", ustim);

    // S1: first stimulus of S1S2
    sour.stimulate(Shape::Rect({0, 0}, {0.1*Lx, Ly}));

    // S2: second stimulus of S1S2
    const Point sensor = Point::Phys(&geom, &mpi, 0.75*Lx, 0.5*Ly);
    LOG("sensor", sensor.get_phys());
    sour.setupS2(sensor, 0, ucrit, [&]() -> void {
        sour.stimulate(Shape::Rect({0, 0}, {0.75*Lx, 0.75*Ly}));
    });

    return sim.run(&sour);
}

This 3D spiral wave is stimulated by extending a 2D spiral wave and placing it on ventricular geometry.

#include "ithildin.h"

int main(int argc, char** argv){
    Mpiclass mpi(argc, argv);

    // choose model
    Model_BO modelbo{1};
    const float ucrit = 0.5;
    const float ustim = 1.0;

    // record local activation time
    ModelWrapper_RecordActivationTime model{&modelbo};
    model.record("u", ucrit);

    // set up geometry
    vector<float> Ds {1, 0.25, 0.25};
    vector<int> size {168, 208, 231};
    vector<float> deltas {0.43, 0.43, 0.5};
    Geometry_OrtAniso geom{
        size, deltas, "ventricles", &model, Ds, &mpi,
    };

    // set up simulation
    float framedur = 20.;
    int Nframe = 30;
    Sim sim{framedur, Nframe, &mpi, &model, &geom, "bocf3d", 0};

    // filament tracking
    sim.enable_filament_recording(0, framedur);

    // set up stimuli
    Source sour{&model, &geom, &mpi, &sim};

    // transform 2D simulation onto 3D ventricles, PD on RV
    NDArray<float> transform; transform.read("bocf2d3d.transform.npy");
    NDArray<float> vars = sour.read_from_stem("bocf2d_0", 34);
    const vector<size_t> shape{
        vars.get_shape(0), size[0], size[1], size[2],
    };
    vars = vars.interpolate(&mpi, shape,
        [&transform](std::vector<float> x){
            if(x.size() != 4) { ERROR("invalid shape"); }
            std::vector<float> z{1., x[1], x[2], x[3]};
            std::vector<float> y(4, 0);
            for(size_t i=0; i<4; i++) {
                for(size_t j=0; j<4; j++) {
                    y[i] += transform(i,j)*z[j];
                }
            }
            return std::vector<float>{x[0], y[1], y[2], y[3]};
        }
    );
    sour.set_frame(vars);

    return sim.run(&sour);
}

This is the benchmark problem as proposed by Niederer et al. (2011) using the cell model by ten Tusscher & Panfilov (2006).

#include "ithildin.h"

int main(int argc, char** argv){
    std::vector<float> arg = strings_to_floats(argc, argv);
    Mpiclass mpi(argc, argv);

    // arg 1: serial number
    const int serialnr = (arg.size() > 1) ? round(arg.at(1)) : 0;
    MASTER_INFO("serialnr", serialnr);

    // arg 2: discretization {0.1, 0.2, 0.5}
    const float dx = (arg.size() > 2) ? arg.at(2) : 0.2; // mm

    // arg 3: PDE time steps {0.005, 0.01, 0.05}
    const float dt = (arg.size() > 3) ? arg.at(3) : 0.05; // ms

    // define model
    Model_TP06 model_tp06;
    ModelWrapper_RecordActivationTime model{&model_tp06};
    model.record("V", 0);

    // constants
    const float C_m = 0.01; // µF/mm² (membrane capacitance)
    const float chi = 140; // 1/mm (surface-to-volume ratio)

    // diffusivities intra/extra longitudinal/transversal
    const float scale = 1 / (chi * C_m); // mm³/µF
    const float Dil = scale*0.17; // mm²/ms
    const float Dit = scale*0.019; // mm²/ms
    const float Del = scale*0.62; // mm²/ms
    const float Det = scale*0.24; // mm²/ms

    // diffusivities longitudinal/transversal
    const float Dl = Dil*Del / (Dil+Del); // mm²/ms
    const float Dt = Dit*Det / (Dit+Det); // mm²/ms

    // diffusivity tensor only scales, Pmat stores magnitude
    vector<float> D{1., Dt/Dl, Dt/Dl}; // 1
    model.set_Pmat(Dl); // mm²/ms

    // domain size
    const float Lx = 20; // mm
    const float Ly = 7; // mm
    const float Lz = 3; // mm

    // set up geometry
    vector<float> deltas{dx, dx, dx};
    vector<int> size{
        int(Lx/deltas[0])+1, int(Ly/deltas[1])+1, int(Lz/deltas[2])+1,
    };
    vector<float> tissue_angles{0., 0.};
    vector<float> diffpars{1, 1, 0};
    Geometry_OrtAniso geom{
        size, deltas, &model, tissue_angles, diffpars, D, &mpi,
    };

    // limit output
    for(size_t ivar=0; ivar<model.get_Nvar(); ivar++)
    { geom.set_bwritevar(ivar, 0); }
    geom.set_bwritevar(model.get_ivar("V"), 1);
    geom.set_bwritevar(model.get_ivar("latV"), 1);

    // set up simulation and timing
    const int Nframe = 15; // number of frames
    const float framedur = 10.; // time between frames in ms
    Sim sim{
        dt, framedur, Nframe, 0, &mpi, &model, &geom,
        "niederer2011benchmark", serialnr, false,
    };

    // set up stimulus
    Source sour(&model, &geom, &mpi, &sim);
    const float I_stim = 50; // mA/cm³ (stimulus current)
    const float i_stim = I_stim * scale; // A/F = V/s = mV/ms
    sour.set_val("V", i_stim); // mV/ms
    sour.stimulate({
            Shape::Rect({0, 0, 0}, {1.5, 1.5, 1.5}),
            1.0, // amplitude factor
            2.0, // ms; duration
            true, // flag: additive stimulus
    });

    return sim.run(&sour);
}

Acknowledgments:

We are grateful to Daniël A. Pijnappels, Antoine A.F. de Vries and our other collaborators at the LUMC, as well as Louise Arno, Lore Leenknegt, and Nathan Dermul at KU Leuven in Kortrijk for useful feedback and discussions.

While no generative AI has been used to write the Ithildin source code, generative AI has been used by the authors to aid in the writing of the text, specifically the Large Language Model (LLM) implementations GitHub Copilot, ChatGPT using GPT-3.5 and GPT-4o, as well as Mistral 0.2 via the Ollama software. The authors confirm that they have followed the current ethical publishing practices of PLOS One.

Funding:

DK is supported by KU Leuven grant GPUL/20/012. MC is supported by KU Leuven grant STG/19/007 and FWO-Flanders fellowship, grant 11PMS24N. HD is supported by KU Leuven grant STG/19/007. OB is supported by Agence Nationale de la Recherche grant ANR-IHUA-04. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests:

The authors have declared that no competing interests exist.

Copyright:

© 2024 Kabus et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Author contributions:

DK: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Data Curation, Writing – original draft, Writing – review & editing, Visualization, Project administration. MC: Methodology, Software, Validation, Formal analysis, Investigation, Data Curation, Writing – original draft, Writing – review & editing, Visualization. CZ: Software. OB: Software. HD: Conceptualization, Software, Resources, Writing – review & editing, Supervision, Project administration, Funding acquisition.