Skip to main content

lifex-ep: a robust and efficient software for cardiac electrophysiology simulations

Abstract

Background

Simulating the cardiac function requires the numerical solution of multi-physics and multi-scale mathematical models. This underscores the need for streamlined, accurate, and high-performance computational tools. Despite the dedicated endeavors of various research teams, comprehensive and user-friendly software programs for cardiac simulations, capable of accurately replicating both normal and pathological conditions, are still in the process of achieving full maturity within the scientific community.

Results

This work introduces \(\texttt {life}^{\text{x}}\)-ep, a publicly available software for numerical simulations of the electrophysiology activity of the cardiac muscle, under both normal and pathological conditions. \(\texttt {life}^{\text{x}}\)-ep employs the monodomain equation to model the heart’s electrical activity. It incorporates both phenomenological and second-generation ionic models. These models are discretized using the Finite Element method on tetrahedral or hexahedral meshes. Additionally, \(\texttt {life}^{\text{x}}\)-ep integrates the generation of myocardial fibers based on Laplace–Dirichlet Rule-Based Methods, previously released in Africa et al., 2023, within \(\texttt {life}^{\text{x}}\)-fiber. As an alternative, users can also choose to import myofibers from a file. This paper provides a concise overview of the mathematical models and numerical methods underlying \(\texttt {life}^{\text{x}}\)-ep, along with comprehensive implementation details and instructions for users. \(\texttt {life}^{\text{x}}\)-ep features exceptional parallel speedup, scaling efficiently when using up to thousands of cores, and its implementation has been verified against an established benchmark problem for computational electrophysiology. We showcase the key features of \(\texttt {life}^{\text{x}}\)-ep through various idealized and realistic simulations conducted in both normal and pathological scenarios. Furthermore, the software offers a user-friendly and flexible interface, simplifying the setup of simulations using self-documenting parameter files.

Conclusions

\(\texttt {life}^{\text{x}}\)-ep provides easy access to cardiac electrophysiology simulations for a wide user community. It offers a computational tool that integrates models and accurate methods for simulating cardiac electrophysiology within a high-performance framework, while maintaining a user-friendly interface. \(\texttt {life}^{\text{x}}\)-ep represents a valuable tool for conducting in silico patient-specific simulations.

Peer Review reports

Background

Cardiac electrophysiology focuses on the heart conduction system from both the normal and pathological perspectives: e.g., it involves the study, diagnosis and treatment planning of cardiac arrhythmias [1, 2].

Nowadays, several clinical tools are widely employed to address these rhythm disorders. The Electrocardiogram (ECG) provides a recording of the electrical activity of the heart. Various deviations from sinus rhythm, such as Atrial Fibrillation (AF), Ventricular Tachycardia (VT), Left Bundle Branch Block (LBBB), can be monitored and identified based on the distinctive shape and morphology of the ECG on an individualized patient-specific level [3]. The ECG may be combined with other imaging data, such as Magnetic Resonance Imaging (MRI) and Computed Tomography (CT), or electroanatomical maps that are directly recorded on the internal and external surfaces of the heart. Thanks to the aforementioned tools, clinicians can successfully reconstruct both the heart anatomy of a patient and multiple electrophysiology properties of the cardiac tissue. These information lay the foundations for traditional decision-making in cardiology. Indeed, medications and surgeries, such as the implantation of pacemakers or cardioverter defibrillators, are planned accordingly [4].

In recent years, the advent of computer models and in-silico simulations in cardiology has enabled the integration of novel scientific tools in standard clinical practice [5, 6]. Physics-based mathematical models and data-driven methods are combined to generate digital replicas of human hearts containing a detailed electrophysiology description, both at the cellular level and at the organ scale [7]. Patient-specific clinical data are embedded inside numerical simulations of the cardiac function for precision medicine [8]. For instance, personalized electrophysiology simulations are employed to assess risk stratification of arrhythmias, to define optimal catheter-based ablation targets, and to perform Cardiac Resynchronization Therapy (CRT) [9,10,11,12].

In this work, we present \(\texttt {life}^{\text{x}}\)-ep, a publicly available software specifically designed for conducting numerical simulations of cardiac electrophysiology, encompassing both normal and pathological conditions. \(\texttt {life}^{\text{x}}\)-ep is built upon the foundation of \(\texttt {life}^{\text{x}}\) [13], an open-source, high-performance C++ Finite Element (FE) numerical solver capable of tackling multi-physics, multi-scale, and multi-domain differential problems. Leveraging the deal.II [14] FE core, \(\texttt {life}^{\text{x}}\) was conceived as part of the iHEART project (refer to Section “Funding”) with the primary aim of providing the scientific community with a cutting-edge FE solver for cardiac modeling.

As part of the \(\texttt {life}^{\text{x}}\) ecosystem, the software released in \(\texttt {life}^{\text{x}}\)-ep has already been widely employed in combination with other modules for cardiac simulations (some of which publicly available: \(\texttt {life}^{\text{x}}\) [13], \(\texttt {life}^{\text{x}}\)-fiber [15] and \(\texttt {life}^{\text{x}}\)-cfd [16], see also Fig. 1), including electrophysiology [17,18,19,20], mechanics [21,22,23], electromechanics [24,25,26,27], fluid dynamics [28,29,30,31,32,33,34,35], fluid–structure interaction [36], electro-mechano-fluid interaction [37, 38] and myocardial perfusion [39, 40]. This wide range of applications stands as a proof of the flexibility and usability of \(\texttt {life}^{\text{x}}\)-ep.

Overall, developing a personalized computer model of the human heart electrophysiology that accurately represents all the underlying biological aspects has garnered significant attention. Crucial to achieving this goal is the development of efficient parallel numerical algorithms that can handle computationally intensive tasks with high accuracy. In this context, \(\texttt {life}^{\text{x}}\)-ep proposes a software that offers a user-friendly yet very detailed and highly customizable interface, a variety of modeling and numerical options to choose from, while also ensuring accuracy and computational efficiency, representing a valuable tool for conducting in silico patient-specific simulations.

The present release focuses on the modeling of cardiac electrophysiology. In the following paragraphs, we provide a concise overview of cardiac electrophysiology and the prevailing mathematical methods employed to model it.

Cardiac electrophysiology: physiology and modeling

The heart wall consists of three distinct layers: the internal thin endocardium, the external thin epicardium and the thick muscular cardiac tissue known as the myocardium. The latter is predominantly composed of cardiomyocytes, which are specialized, striated excitable muscle cells responsible for the essential cardiac function. When these cardiomyocytes are stimulated by an electrical impulse, a change in the electro-chemical balance of the cell membrane results in a series of biochemical reactions that determine a large variation of the transmembrane potential, namely the voltage differential between the intra and extracellular spaces of the cell. The rapid depolarization and subsequent slow repolarization mechanism is known as the action potential. This process is triggered and controlled by the opening and closing of voltage-gated ion channels that make the cell membrane permeable to specific ionic species, like sodium, potassium, and calcium. The transmembrane potential changes as a result of the ionic fluxes, which, in turn, are driven by the voltage difference itself.

The cells of the heart tissue are connected to each other through gap junctions, i.e. intercellular low-resistance ionic channels, which allow the electric signal to travel from cell to cell across the whole cardiac muscle (known as myocardium). The propagation of the electric signal within the myocardium is highly anisotropic, with preferential directions for conduction determined by the presence of muscle fibers [17]. These fibers are organized in sheets, that determine a second preferential direction of propagation. For further details on the biophysical mechanisms of cardiac electrophysiology, we refer the interested reader to [1, 2, 41]. A visual representation of the different spatial scales involved in this complex chain of physical processes is given in Fig. 2.

The mathematical description of these processes is made of two building blocks: a ionic model, describing the chemical processes taking place at the cellular scale, and an action potential propagation model, representing the spatial propagation of the action potential wavefront at the tissue scale.

Ionic models

Several ionic models with different levels of biophysical detail and suitable for specific types of cells have been proposed in the literature [41]. Most of them are written in the form of the following Ordinary Differential Equations (ODEs)

$$\begin{aligned} \left\{ \begin{aligned}&\frac{du(t)}{dt}+{{\mathcal {I}}_{\text{ion}}^{}}(u(t),{\varvec{w}}^{}(t)) = {{\mathcal {I}}_{\text{app}}}(t),&\qquad t \in (0, T],\\&\frac{d{\varvec{w}}^{}(t)}{dt}={\varvec{H}}^{}(u(t),{\varvec{w}}^{}(t)),&\qquad t \in (0, T],\\&u(0)=u_0^{}, \, {\varvec{w}}^{}(0)={\varvec{w}}_0^{},&\end{aligned} \right. \end{aligned}$$
(1)

where the unknowns are \(u=u(t)\), the transmembrane potential, and the vector \({\varvec{w}}^{}={\varvec{w}}^{}(t)=(w_1,\ldots ,w_{M_{}})\), collecting \(M_{}\) ionic variables. Ionic variables may include concentrations of different ionic species and gating variables (describing the opening probability of ionic channels). The term \({{\mathcal {I}}_{\text{ion}}^{}}\) describes the electric current generated by the flux of ionic species across the cell membrane, while \({{\mathcal {I}}_{\text{app}}}\) represents an externally applied current. The dynamics of ionic variables is modeled by the function \({\varvec{H}}^{}\). Different ionic models are characterized by the number of ionic variables \(M_{}\) and the definition of the two functions \({{\mathcal {I}}_{\text{ion}}^{}}\) and \({\varvec{H}}^{}\).

The \(\texttt {life}^{\text{x}}\)-ep release includes two phenomenological ionic models (Aliev-Panfilov (APf) [42] and Bueno-Orovio (BO) [43]) and two physiological ones, one for ventricular cells (ten Tusscher-Panfilov 2006 (TTP06) [44]) and one for atrial cells (Courtemanche-Ramirez-Nattel (CRN) [45]).

Action potential propagation models

One of the most popular action potential propagation models is the monodomain equation, a Partial Differential Equation (PDE) in which the spatial propagation of the electric signal through gap junctions is accounted for by a diffusion term [41].

Referring to Fig. 3, we consider a computational domain \(\Omega _{} \subset {\mathbb {R}}^3\), representing the region of the myocardium of interest (e.g., a slab of cardiac tissue, atria, ventricles or whole heart geometries), see Fig. 3(a). In the whole domain, we consider a local orthonormal triplet of vectors \({{\textbf{f}}_0}\), \({{\textbf{s}}_0}\) and \({{\textbf{n}}_0}\), defining the fiber, the sheetlet and the sheet-normal directions see Fig. 3(b), respectively [17].

A notable feature of \(\texttt {life}^{\text{x}}\)-ep is its capability to divide the domain \(\Omega _{}\) into a generic number \(N\) of subdomains, thereby allowing distinct electrophysiology properties to be assigned to each individual region (e.g. different ionic models or different electrical conductivities). In practical implementation, this goal is realized by providing the simulation with a computational mesh, wherein each subdomain is distinctly labeled. More formally, we introduce a partition of \(\Omega _{}\) into \(N\) disjoint subdomains, namely \(\Omega _{1}, \Omega _{2}, \dots , \Omega _{N}\) (more precisely, we assume \({\overline{\Omega }}_{} = \cup _{i=1}^{N} {\overline{\Omega }}_{i}\) and \(\Omega _{i} \cap \Omega _{j} = \emptyset\) for \(i \ne j\)). We denote by \(\Gamma _{i} = \partial \Omega _{} \cap \partial \Omega _{i}\) the external boundary of the i-th subdomain, and by \(\Gamma _{ij} = \partial \Omega _{i} \cap \partial \Omega _{j}\) the interface between the i-th and the j-th subdomains, see Fig. 3(a).

For each subdomain \(\Omega _{i}\), we consider a ionic model, featuring \(M_{i}\) ionic variables, and characterized by the functions \({\varvec{H}}^{i}\) and \({{\mathcal {I}}_{\text{ion}}^{i}}\) and initial conditions \(u_0^{i}\) and \({\varvec{w}}_0^{i}\). Moreover, we introduce the unknown \({\varvec{w}}^{i} :\Omega _{i}\times [0,T] \rightarrow {\mathbb {R}}^{M_{i}}\). The transmembrane potential, namely \(u:\Omega _{}\times [0,T] \rightarrow {\mathbb {R}}\), is instead defined in the whole computational domain.

The monodomain model reads as follows [5, 41]

$$\begin{aligned} \left\{ \begin{aligned}&\frac{\partial u}{\partial t}+{{\mathcal {I}}_{\text{ion}}^{}}(u,{\varvec{w}}^{i}) -\nabla \cdot ({\varvec{D}}^{i}\nabla u)={{\mathcal {I}}_{\text{app}}}({{\textbf {x}}},t),{} & {} \, \text{ in } \Omega _{i}\times (0,T], \\&\frac{d{\varvec{w}}^{i}}{dt}={\varvec{H}}^{i}(u,{\varvec{w}}^{i}),{} & {} \, \text{ in } \Omega _{i}\times (0,T], \\&({\varvec{D}}^{i}\nabla u) \cdot {\textbf{n}}_{i} = 0,{} & {} \, \text{ on } \Gamma _{i}\times (0,T], \\&({\varvec{D}}^{i}\nabla u) \cdot {\textbf{n}}_{i} + ({\varvec{D}}^{j}\nabla u) \cdot {\textbf{n}}_{j} = 0,{} & {} \, \text{ on } \Gamma _{ij}\times (0,T], \, \\&u=u_0^{i},{} & {} \, \text{ in } \Omega _{i}\times \{0\}, \\&{\varvec{w}}^{i}={\varvec{w}}_0^{i},{} & {} \, \text{ in } \Omega _{i}\times \{0\}, \end{aligned} \right. \end{aligned}$$
(2)

where \(i,j = 1,\dots ,N\), and the diffusion tensor \({\varvec{D}}^{i}\) is defined as

$$\begin{aligned} {\varvec{D}}^{i} = \sigma _{\text {l}}^{i} {{\textbf{f}}_0}\otimes {{\textbf{f}}_0}+ \sigma _{\text {t}}^{i} {{\textbf{s}}_0}\otimes {{\textbf{s}}_0}+ \sigma _{\text {n}}^{i} {{\textbf{n}}_0}\otimes {{\textbf{n}}_0}, \end{aligned}$$
(3)

where \(\sigma _{\text {l}}^{i}, \sigma _{\text {t}}^{i}, \sigma _{\text {n}}^{i}\) being positive coefficients denoting the longitudinal, transversal and normal conductivities, respectively [17]. The operator \(\otimes\) denotes the outer product acting on the vectors \({{\textbf{f}}_0}\), \({{\textbf{s}}_0}\) and \({{\textbf{n}}_0}\). Given two vectors \({\textbf{v}}\) and \({\textbf{w}}\), their outer product is a tensor, component-wise defined as \(({\textbf{v}} \otimes {\textbf{w}})_{ij} = {\textbf{v}}_i {\textbf{w}}_j\), for \(i, j = 1, 2, 3\). Homogeneous Neumann boundary conditions are prescribed on the whole boundary \(\partial \Omega\) to impose the condition of electrically isolated domain, with \({\textbf{n}}_{i}\) denoting the outward normal unit vector to the boundary. On the interface between subdomains, we impose continuity of flux conditions.

The action potential is triggered by an external applied current \({{\mathcal {I}}_{\text{app}}}({{\textbf {x}}},t)\), that mimics the presence of a natural or artificial pacemaker, or the role of specialized conduction systems, such as the Purkinje network [1, 2] or other relevant bundles. Since these anatomical entities are not explicitly modeled in this release, \(\texttt {life}^{\text{x}}\)-ep provides great flexibility in the definition of \({{\mathcal {I}}_{\text{app}}}({\textbf {x}},t)\), allowing for the definition of multiple stimuli of various shape and with customized duration and delay. The applied current is assumed to be in the following form:

$$\begin{aligned} {{\mathcal {I}}_{\text{app}}}({\textbf {x}}, t) = \sum _{i = 1}^{N_\text {stim}}g^t_i(t) g^{{\textbf {x}}}_i({\textbf {x}}), \end{aligned}$$

wherein \(N_\text {stim} \in {\mathbb {N}}\) is the number of applied stimuli, \(g^t_i:(0, T) \rightarrow \{0, 1\}\) is an indicator function that only activates the i-th stimulus during a prescribed time interval, and \(g^{{\textbf {x}}}_i: \Omega _{} \rightarrow [0, \infty )\) is a space-dependent function that describes the shape of the stimulus, defined either as the indicator function of a given subset of \(\Omega _{}\) (for cubic, spherical or planar impulses) or as a Gaussian function. Location, shape and timing of each stimulus can be customized by the user.

The possibility of prescribing different properties to different regions of the domain enables several use cases. For instance, it allows multi-chamber simulations, in which different ionic models are assigned to atria and ventricles [24]. Moreover, it also gives the possibility to simulate pathological scenarios, such as the presence of scars or fibrosis [18, 25, 46, 47]. Non-conductive zones can be introduced in the model, by labelling a given subregion of the domain as non-conductive. In this case, the corresponding portion of the domain will be excluded from the simulation, thus creating a block in the conduction [24]. Furthermore, grey zones can be modeled as well, by suitably adjusting the tissue conductivities and the ionic model properties in grey zone and fibrotic subregions [18, 25, 48, 49].

We remark that the monodomain equation is often expressed as

$$\begin{aligned} \chi _\text{m}\left( C_\text{m}\frac{\partial u}{\partial t}+{\widehat{{\mathcal {I}}}_{\text{ion}}^{}}(u,{\varvec{w}}^{})\right) -\nabla \cdot (\widehat{{\varvec{D}}}^{}\nabla u)=\chi _\text{m}{\widehat{{\mathcal {I}}}_{\text{app}}}({{\textbf {x}}},t), \end{aligned}$$

where \(C_\text{m}\) denotes the membrane cell capacitance and \(\chi _\text{m}\) is the membrane surface-to-volume ratio. This formulation can clearly be recast to the one of (2) by the rescaling \({{\mathcal {I}}_{\text{ion}}^{}} = {\widehat{{\mathcal {I}}}_{\text{ion}}^{}} / C_\text{m}\), \({{\mathcal {I}}_{\text{app}}}{} = {\widehat{{\mathcal {I}}}_{\text{app}}}{} / C_\text{m}\) and \({\varvec{D}}^{} = \widehat{{\varvec{D}}}^{} / (\chi _\text{m}C_\text{m})\). These equivalences are helpful when comparing different literature sources.

Cardiac electrophysiology: numerical discretization

We partition the temporal domain (0, T] into \(N_T\) subintervals with a time step \(\Delta t=t_{n+1} - t_{n}\), where \(t_n = n\Delta t\) for \(n = 0, 1, \dots , N_T\). We denote with a subscript n the approximation of a variable at time \(t_n\) (e.g., \(u_n \approx u(t_n)\)). Time derivatives are approximated by means of a Backward Differentiation Formula (BDF) scheme of order \(\sigma\), with \(\sigma \in \{1, 2, 3\}\). Given a generic time-dependent function f, its time derivative is thus approximated by

$$\begin{aligned} \frac{\partial f}{\partial t} \approx \frac{\alpha _\text {BDF}f_{n+1} - f_{\text {BDF},n}}{\Delta t}, \end{aligned}$$

where \(f_{\text {BDF},n}\) is a linear combination of \(f_n\), \(f_{n-1}\), \(\dots\), and \(\alpha _\text {BDF}\) is a coefficient depending on the order of the scheme [50]. We also define explicit extrapolations of order \(\sigma\), denoted by the subscript EXT (e.g., \(f_{\text {EXT},n+1} \approx f(t_{n+1})\)).

The numerical discretization of the system (2) is obtained by means of the implicit-explicit time advancing scheme described in [26]. Given the solution up to time step \(t_n\), to compute the solution at time \(t_{n+1}\):

  1. 1

    solve the time-discrete ionic model equations: for \(i = 1, \dots , N\),

    $$\begin{aligned} \frac{\alpha _\text {BDF}{\varvec{w}}^{i}_{n+1} - {\varvec{w}}^{i}_{\text {BDF},n}}{\Delta t} = {\varvec{H}}^{i}(u_{\text {EXT},n+1},{\varvec{w}}^{i}_{\text {EXT},n+1},{\varvec{w}}^{i}_{n+1}), \quad \text{ in } \Omega _{i}. \end{aligned}$$
    (4)

    We treat some of the ionic variables appearing in \({\varvec{H}}^{i}\) with an implicit formulation, and others with an explicit formulation, so that the resulting problem can be solved by directly inverting the equations [26].

  2. 2

    solve the time-discrete monodomain equation: for \(i,j = 1, \dots , N\),

    $${ \left\{ \begin{aligned} & \begin{aligned} \frac{{\alpha _{{\text{BDF}}} u_{{n + 1}} - u_{{\text{BDF},n}} }}{{\Delta t}} &+ {\mathcal{I}}_{{\text{ion}}}^{{}} (u_{{\text{EXT},n + 1}}, \boldsymbol{w}_{{n + 1}}^{i} ) \\&- \nabla \cdot (\boldsymbol{D}^{i} \nabla u_{{n + 1}} ) = {\mathcal{I}}_{{\text{app}}} ({\bf x},t), \end{aligned} &&\text{in } \Omega _{i} ,\\ &(\boldsymbol{D}^{i} \nabla u_{{n + 1}} ) \cdot \mathbf{n}_{i} = 0, &&\text{on } \Gamma _{i}, \\ &(\boldsymbol{D}^{i} \nabla u_{{n + 1}}) \cdot \mathbf{n}_{i} + (\boldsymbol{D}^{j} \nabla u_{{n + 1}} ) \cdot \mathbf{n}_{j} = 0,&&\text{on }\Gamma _{{ij}} .\end{aligned}\right.}$$
    (5)

    Notice that the ionic current term is treated explicitly by using the extrapolated potential \(u_{\text {EXT},n+1}\), so that the resulting problem is linear in \(u_{n+1}\).

For the spatial discretization, we introduce a tetrahedral or hexahedral mesh over \(\Omega\), and use the FE method [51] to approximate the solution variables \(u\) and \({\varvec{w}}^{i}\) as piecewise polynomials of order \(p\). \(\texttt {life}^{\text{x}}\)-ep supports polynomials of orders 1 and 2 on tetrahedral meshes, and polynomials of arbitrary degree on hexahedral meshes.

The ionic model (4) is solved at each support point of the degrees of freedom of the FE space. The ionic current \({{\mathcal {I}}_{\text{ion}}^{}}\) is also evaluated at every support point, and then interpolated onto quadrature nodes on the interior of the mesh elements, in the approach known as Ionic Current Interpolation (ICI) [52, 53]. The time-discrete monodomain equation (5) is discretized in space using the FE method, leading to the following algebraic linear system of equations

$$\begin{aligned} \left( \frac{\alpha _\text {BDF}}{\Delta t}{\textsf{M}}+ {\textsf{K}}\right) {\textsf{u}}{}_{n+1} = \frac{1}{\Delta t}{\textsf{M}}{} {\textsf{u}}{}_{\text {EXT},n+1} - {\textsf{s}}_{n+1} + {\textsf{f}}_{n+1}, \end{aligned}$$
(6)

where, denoting by \(\varphi _j\) the basis functions of the FE space, \({\textsf{M}}\) is the mass matrix of entries \({\textsf{M}}_{jk} = \int _\Omega \varphi _k \varphi _j d{\textbf{x}}\), \({\textsf{K}}\) is the stiffness matrix of entries \({\textsf{K}}_{jk} = \sum _{i = 1}^{N} \int _{\Omega _{i}} {\varvec{D}}^{i}\nabla \varphi _k \cdot \nabla \varphi _j d{\textbf{x}}\) and \({\textsf{f}}_{n+1}\) is the vector arising from the applied current term \({{\mathcal {I}}_{\text{app}}}\), of entries \(({\textsf{f}}_{n+1})_i = \int _{\Omega _{}} {\mathcal {I}}_{\text {app},n+1}\,\varphi _i d{\textbf{x}}\). The vector \({\textsf{s}}_{n+1}\) arises from the ionic current terms \({{\mathcal {I}}_{\text{ion}}^{i}}\). Exploiting the ICI formulation, it can be computed as

$$\begin{aligned} {\textsf{s}}_{n+1} = \sum _{i = 1}^{N} {\textsf{M}}^i {\textsf{I}}_{\text {ion},n+1}^i, \end{aligned}$$

where \({\textsf{M}}^i\) is the mass matrix for subdomain \(\Omega _{i}\), whose entries are \({\textsf{M}}_{jk}^i = \int _{\Omega _{i}} \varphi _k \varphi _j d{\textbf{x}}\), and \({\textsf{I}}_{\text {ion},n+1}^i\) is the vector of the evaluations of \({{\mathcal {I}}_{\text{ion}}^{i}}\) at the support points of the FE space. The integrals that arise from the FE discretization are numerically approximated using the Gauss–Legendre quadrature rule, with the minimum number of points to ensure exact integration of the mass matrix.

Overview of existing software

One of the first landmarks in computational cardiology has been defined by the pioneering work of Hodgkin and Huxley in the mid-1950s [54]. The first cardiac action potential model has been developed by D. Noble in 1962 [55]. Since then, a plethora of cardiac electrophysiology models has been proposed in the literature [6, 41, 56, 57]. These mathematical models feature a different degree of biophysical complexity and act either at the microscopic level, by describing the behavior of single cardiomyocytes, or at the organ scale, where an ensemble of many myocardial cells is considered.

The development of software to perform electrophysiology simulations is still mainly steered by academia and public institutions, despite this variety of mathematical models, the presence of increasingly elaborated numerical methods, and the evolution of computer hardware. Indeed, several research tools for multi-physics and multi-scale cardiac simulations have been proposed in the last two decades. Among them, an important role is played by openCARP [58], an open-source C++ simulation environment integrated with cellML, a public repository that encompasses many cell-based mathematical models [59], Chaste [60], an open-source C++ library mainly developed at the University of Oxford, Cardioid [61], a highly efficient and scalable tool for high resolution electrophysiology simulations mainly developed at the Lawrence Livermore National Laboratory, and Alya [62], a high-performance code from the Barcelona Supercomputing Center. Prominent examples of industrial projects that demonstrate translational research efforts for the heart are given by the Living Heart Project by Dassault Systémes [63, 64] and the services and software provided by NumeriCor GmbHFootnote 1 [65].

The aforementioned tools allow to perform single-chamber, bi-atrial, bi-ventricular or four-chamber heart electrophysiology simulations by means of accurate, yet computationally expensive physiologically-based models, such as the bidomain or monodomain equation coupled with the TTP06 and CRN ionic models [41, 44, 66], or the more efficient reaction-eikonal equation [67]. Model parameters can be calibrated on a patient-specific basis to match ECG or Body Surface Potential Mapping (BSPM) [68]. Different pathological scenarios involving AF [69, 70], VT [9, 12] and LBBB  [11, 71] have been addressed by employing these software tools. Moreover, heterogeneity in the tissue and cellular properties can be prescribed to incorporate the presence of scar, grey zones and fibrosis importing measurements from clinical data, such as Late Gadolinium Enhancement-Magnetic Resonance Imaging (LGE-MRI) [10], contrast enhanced CT [72], or the Imaging Itensity Ratio (IIR) [73]. This allows to locally vary Conduction Velocity (CV) and Action Potential (AP) morphology. Another application is related to the in-silico assessment of drugs efficacy by means of numerical simulations [74, 75], where model parameters can be tuned to replicate the effects of pharmacological therapies. Electrophysiology simulations are also used to evaluate gender differences in healthy and pathological conditions involving arrhythmias [75, 76].

Compared to the aforementioned software and libraries, the \(\texttt {life}^{\text{x}}\)-ep solver stands out with several distinctive features and advantages, mostly inherited from the \(\texttt {life}^{\text{x}}\) core structure [13]. It is designed to be user-friendly and easy to use, even for biomedical researchers without extensive experience in numerical methods. The solver is implemented in C++ using advanced programming paradigms and leverages Message Passing Interface (MPI) for distributed memory parallelism. Moreover, it supports the possibility to import arbitrary meshes with either hexahedral or tetrahedral elements, and incorporates advanced numerical solvers based on the Trilinos linear algebra backend, thus ensuring precise control over the numerical setting and accuracy. The solver also exhibits ideal scalability up to thousands of cores, as demonstrated in Section Strong scalability test, allowing to efficiently simulate large-scale scenarios. In addition to the numerical and programming features stemming from its foundation on \(\texttt {life}^{\text{x}}\), \(\texttt {life}^{\text{x}}\)-ep offers two options for prescribing myocardial fibers, which can be either imported from a file or generated online taking advantage of the previous release \(\texttt {life}^{\text{x}}\)-fiber [15], based on the LDRBMs presented in [17]. Moreover, it supports spatial heterogeneity in the choice of both models and physical coefficients, easily configurable through a convenient parameter file, without the need to access and modify the source code.

In general, \(\texttt {life}^{\text{x}}\)-ep stands out in its ease of use, performance, and compatibility with common I/O (input/output) file formats, as well as its comprehensive and self-contained infrastructure, achieved by combining sophisticated mathematical models with accurate numerical schemes. To the best of our knowledge, none of the packages mentioned above exhibits similar features altogether.

Implementation

In this section, we present the technical specifications of \(\texttt {life}^{\text{x}}\)-ep and provide a comprehensive documentation of its user interface. The aim is to guide users through the entire process, from downloading the software to successfully running a full cardiac electrophysiology simulation.

\(\texttt {life}^{\text{x}}\)-ep offers a numerical solver tailored for cardiac electrophysiology, leveraging the mathematical models and numerical algorithms discussed in the previous section. The linear algebra backend is provided by Trilinos [77], integrated into deal.II. This incorporates the implementation of various linear solvers (CG and GMRES) and flexible black-box preconditioners (AMG, additive Schwarz, block Jacobi) supported by \(\texttt {life}^{\text{x}}\). The code is inherently parallel and designed to run efficiently on a diverse range of architectures, spanning from personal laptops to High-Performance Computing (HPC) facilities and cloud platforms. To ensure reliability and performance, \(\texttt {life}^{\text{x}}\)-ep has been thoroughly tested on multiple systems, including a cluster node equipped with 192 cores based on Intel Xeon Gold 6238R (2.20 GHz) at MOX, Dipartimento di Matematica, Politecnico di Milano, as well as the GALILEO100 supercomputer available at CINECA (Intel CascadeLake 8260, 2.40GHz, technical specifications available at https://wiki.u-gov.it/confluence/display/SCAIUS/UG3.3

For further insights into the core functionalities of \(\texttt {life}^{\text{x}}\), we recommend referring to [13]. Additionally, in the following sections, we provide a concise guide on how to quickly get started and run simulations using \(\texttt {life}^{\text{x}}\)-ep.

Running simulations in \(\texttt {life}^{\text{x}}\) -ep

The distribution and installation process of \(\texttt {life}^{\text{x}}\)-ep is designed to be user-friendly and platform-independent. The software is conveniently provided as a binary AppImageFootnote 2 executable, which can be obtained from https://doi.org/10.5281/zenodo.8085266. Along with the executable, all the necessary input files required to reproduce the numerical results presented in the subsequent sections are included.

The adoption of the AppImage format ensures a universal package that is compatible with x86-64 Linux operating systems, eliminating the need for multiple distribution-specific versions. From the user’s perspective, this translates to a seamless and straightforward download-then-run experience, without the hassle of manually managing system dependencies.

As an AppImage, \(\texttt {life}^{\text{x}}\)-ep has been built on Debian Buster,Footnote 3 which corresponds to the current oldoldstable version. This follows the principle of “Build on old systems, run on newer systems”.Footnote 4 Consequently, the software is expected to function on virtually any recent x86-64 Linux distribution, provided that glibcFootnote 5 version 2.28 or higher is installed.

Once downloaded and extracted the \(\texttt {life}^{\text{x}}\)-ep archive, the AppImage file needs to be made executable by typing the following command in a terminal:

figure a

Finally, lifex_electrophysiology-1.5.0-x86_64.AppImage can be executed with:

figure b

No root permissions are required for the commands mentioned above to run successfully. However, it is important to note that the AppImage relies on the userspace filesystem framework called FUSE.Footnote 6 Please ensure that FUSE is installed on your system. If you encounter any errors, the following commands may be helpful in resolving the issue:

figure c

Additionally, we recommend referring to the AppImage troubleshooting guide.Footnote 7

The following command will provide an inline help that includes detailed information about all the available command line options and their purpose:

figure d

The executable allows to run test cases with an arbitrary number of disjoint subdomains \(\Omega _i,\,i=1,\dots ,N\), which are also referred to as volumes. The configuration of the simulation is supplied through a parameter file. The user can generate a template parameter file using the following command:

figure e

To match different user needs, the level of detail in the parameter files can be adjusted using the optional minimal or full option after the -g flag. The minimal option reduces the level of detail, making it suitable for initial usage of \(\texttt {life}^{\text{x}}\)-ep. On the other hand, the full option increases the level of detail, exposing advanced options, such as parameter choices and detailed options on linear algebra and preconditioning, among others. If the user does not specify either the minimal or full option, an intermediate verbosity level is selected by default. In all cases, parameters that are not present in the file will retain their default values. This flexibility allows users to customize the level of detail in the parameter files according to their specific requirements and familiarity with the software.

The parameters are written in a plain text file, organized as a list of key-value pairs grouped in subsections, which describe the configuration for the simulation to be run. Each parameter is accompanied by a brief documentation within the parameter file itself, explaining its meaning.

In the provided command, the optional argument -vol \(\texttt {<}\)volume labels\(\texttt {>}\) allows to specify a list of user-defined labels. These labels are used to differentiate each subdomain, enabling the selection of heterogeneous model options such as ionic model type, coefficients, and electrical conductivities for each subdomain. If not specified, a single subdomain characterized by its global Volumetric parameters is assumed to exist.

The following example illustrates a parameter file specifying three subdomains: Healthy, Fibrosis, and Scar:

figure f

The input mesh is expected to include at least one volumetric tag corresponding to each of the subdomains that need to be differentiated. These subdomain tags are specified in the Material IDs list, located under each respective subdomain section.

The parameter file includes a dedicated section called Fiber generation, which serves the purpose of enabling the importing from a file of the myocardial fibers or the online generation of them on various geometry types, such as slabs, ventricles, and atria. To accomplish this, \(\texttt {life}^{\text{x}}\)-ep incorporates the functionalities of its predecessor, \(\texttt {life}^{\text{x}}\)-fiber [15], which utilizes the Laplace-Dirichlet Rule-Based Methods presented in [17]. This unique feature of \(\texttt {life}^{\text{x}}\)-ep sets it apart from other existing software alternatives.

Once the user has edited the parameter file, the simulation can be started using the following command:

figure g

When executing the command, the volume labels provided must match the ones used for generating the parameter file. It is essential to use consistent volume labels throughout the process to ensure proper identification and configuration of the subdomains within the simulation.

A parallel simulation is started prepending the command with the mpirun or mpiexec wrapper (which may vary depending on the MPI implementation available), e.g.:

figure h

where N represents the desired number of parallel processes. The binary package supports parallel execution using MPICH (https://www.mpich.org/) version 4.0 or higher.

The parameter file also includes options that enable the serialization of the solution, allowing the simulation to be paused or stopped at any point and then resumed at a later time using the serialized data. This feature is particularly useful when dealing with long-running simulations or when unexpected interruptions occur.

License and third-party software

This work is copyrighted by the \(\texttt {life}^{\text{x}}\)-ep authors and licensed under the Creative Commons Attribution Non-Commercial No-Derivatives 4.0 International License.Footnote 8

It should be noted that \(\texttt {life}^{\text{x}}\)-ep incorporates several third-party libraries, which are separately copyrighted by their respective authors and whose use is covered by various permissive licenses.

Third-party software bundled with (in binary form), required by, copied, modified, or explicitly used in \(\texttt {life}^{\text{x}}\)-ep include:

\(\texttt {life}^{\text{x}}\)Footnote

https://lifex.gitlab.io/

:

[13]: the open-source, high-performance software providing the core functionalities for the numerical solution of the Finite Element problems described in the previous section;

deal.IIFootnote

https://www.dealii.org/

:

[14]: it provides support to mesh handling, assembling and solving Finite Element problems (compiled with enabled support to TrilinosFootnote 11 for linear algebra data structures and solvers) and to input/output functionalities;

BoostFootnote

https://www.boost.org/

:

[78]: its modules Filesystem and Math are used for manipulating files/directories and for advanced mathematical functions and interpolators, respectively;

VTKFootnote

https://vtk.org/

:

[79]: it is used for importing external surface or volume input data and coefficients appearing in the mathematical formulation.

Some of the packages listed above, as stated by their respective authors, rely on additional third-party dependencies that may also be bundled (in binary form) with \(\texttt {life}^{\text{x}}\)-ep, although not used directly. These dependencies include: ADOL-C,Footnote 14ARPACK-NG,Footnote 15BLACS,Footnote 16Eigen,Footnote 17FFTW,Footnote 18GLPK,Footnote 19HDF5,Footnote 20HYPRE,Footnote 21METIS,Footnote 22MPICH,Footnote 23MUMPS,Footnote 24NetCDF,Footnote 25OpenBLAS,Footnote 26PETSc,Footnote 27ParMETIS,Footnote 28ScaLAPACK,Footnote 29Scotch,Footnote 30SuiteSparse,Footnote 31SuperLU,Footnote 32oneTBB,Footnote 33p4est.Footnote 34

These libraries are free software and have relatively few restrictions on their use. However, please note that different terms may apply. For detailed information on the licenses and copyright statements for these packages, please refer to the content of the folder share/doc/licenses/ in the \(\texttt {life}^{\text{x}}\)-ep archive.

Results and discussion

To highlight the versatility of \(\texttt {life}^{\text{x}}\)-ep, this section presents the results obtained in various numerical examples conducted on a range of idealized and realistic geometries, encompassing both normal and pathological scenarios. To facilitate the reproduction of these test cases, this section provides a getting started guide, detailing the following steps:

  • Generating and importing the input data (e.g., computational meshes and fibers);

  • Configuring the parameter files specific to each test case and executing the corresponding simulation;

  • Performing post-processing of the results and visualizing the output.

Finally, a strong scalability test is presented in Section Strong scalability test.

Input data

All parameter files and meshes related to the numerical simulations described below can be downloaded from the \(\texttt {life}^{\text{x}}\)-ep release archive https://doi.org/10.5281/zenodo.8085266.

This guide provides various pre-configured hexahedral and tetrahedral meshes, including:

  • four idealized geometries: cardiac slab tissue, left atrium, left ventricle, and ventricular slab, see Fig. 4(a–d, g);

  • two realistic geometries: left atrium and left ventricle, see Fig. 4(e–f, h–i).

We emphasize that the provided example meshes are solely used to illustrate the \(\texttt {life}^{\text{x}}\)-ep features, as users have the flexibility to input any (idealized or patient-specific) meshes into \(\texttt {life}^{\text{x}}\)-ep.

The cardiac tissue slab, idealized left atrial, and left ventricular meshes are characterized by a single volumetric tag representing the entire myocardium, Fig. 4(a–c). On the other hand, the ventricular slab, realistic left atrial, and left ventricular meshes have multiple distinct volumetric tags, Fig. 4(d–i). The ventricular slab is divided into sub-endocardial, myocardial, and sub-epicardial layers, Fig. 4(d, g). The realistic left atrium and left ventricle include regions with scars, grey zones, and fibrosis, Fig. 4(e, f, h, i). The former meshes are used for single volume simulations, while the latter for multi-volume simulations. In all cases, volumetric tags must be defined during the mesh generation process [80].

The geometrical models for the idealized cardiac meshes were created using the built-in CAD engine of gmsh,Footnote 35 an open-source 3D FE mesh generator. The gmsh scripts used to generate these meshes are included in the \(\texttt {life}^{\text{x}}\)-ep release archive (https://doi.org/10.5281/zenodo.8085266). For detailed information on the syntax of the scripts, we refer to the online documentation of gmsh. Tetrahedral mesh generation for the slab tissue, idealized left atrium, and ventricular slab are also performed using gmsh. On the other hand, the idealized and realistic left ventricular hexahedral meshes were instead generated using the cubitFootnote 36 mesh generation software toolkit. Finally, the realistic left atrial mesh was perfomed using the semi-automatic meshing tools developed in [80], based on the Vascular Modelling Toolkit (vmtk) software [81].

The realistic left atrium and left ventricle, also containing the scar and fibrotic regions, have been produced starting from the openly available meshes used in [82] (for the left atriumFootnote 37) and in [83] (for the left ventricleFootnote 38). For the latter, we used the (vmtk) software [81], in conjunction with the cubit mesh generator.

Mesh files can be specified in the \(\texttt {life}^{\text{x}}\)-ep parameter file for electrophysiology simulation within the section named Mesh and space discretization, by setting the proper element type (tetrahedra or hexahedra), the polynomial order to be used in the FE discretization and an appropriate mesh rescaling factor (e.g., if the coordinates in the input mesh file are given in millimeters, Scaling factor = 1e-3, since \(\texttt {life}^{\text{x}}\)-ep internally handles all physical quantities in the International System of Units.

figure i

Regarding the prescription of the myocardial fiber architecture, an essential building block for cardiac electrophysiology simulations, \(\texttt {life}^{\text{x}}\)-ep provides two options. Users can either import the myofibers from a file or generate them online using LDRBMs [17] by incorporating the \(\texttt {life}^{\text{x}}\)-fiber release package, which was recently published in [15].

To generate the myocardial fibers using LDRBMs, users can select the appropriate Geometry type within the parameter file under the section named Fiber generation. Each Geometry type corresponds to the specific LDRBM applicable for different geometries, such as (ventricular and spherical) slabs, (cut at base and complete) left ventricular and left atrial geometries. The parameters related to the fiber generation for a particular Geometry Type are located within a subsection with the same name. We refer to [15] for a more comprehensive guide on configuring LDRBMs, and to [17] for an in-depth mathematical description.

figure j

As an alternative, users can choose to import myofibers from a file by setting Geometry type = Import from file. The myofibers are imported from a VTU file format with unstructured grid data. In this file, the three fiber directions (\({{\textbf{f}}_0}\), \({{\textbf{s}}_0}\), and \({{\textbf{n}}_0}\) representing the fiber, sheet, and sheet-normal directions, respectively) must be normalized and embedded as point-data arrays. Additionally, users need to set an appropriate geometry rescaling factor (e.g., if the coordinates in the input fiber file are given in millimeters, Scaling factor = 1e-3). This ensures that the imported myofibers align properly with the geometry of the cardiac electrophysiology model.

figure k

The myocardial fibers, generated using LDRBM, for the electrophysiology simulations presented hereafter are illustrated in Fig. 5.

Simulations of normal cardiac electrophysiology

We present normal electrophysiology simulations applied to the set of idealized geometries, namely a rectangular slab of cardiac tissue, an idealized left atrium, an idealized left ventricle and a layered ventricular slab.

Slab benchmark

To perform software verification, we consider the N-version slab benchmark proposed in [84]. This benchmark involves a rectangular slab domain of size \((3 \times 7 \times 20) \times 10^{-3}\) m, as depicted in Fig. 4(a). The fiber directions within the domain are oriented along the long axis (0.02 m), as shown in Fig. 5(a). For comprehensive modeling and geometrical information regarding the benchmark definition, we refer to the original paper [84].

The simulation can be run by first generating the parameter file using

figure l

then configuring the simulation by editing the lifex_electrophysiology.prm default parameter file, and finally running the simulation by typing

figure m

The same can be obtained by directly using the already prepared parameter file params_slab.prm uploaded in the release archive:

figure n

To match with [84], precise specifications are provided for the ionic model employed (in this test case, the TTP06 model [44]), including initial conditions and conductivity values along the myofiber directions (longitudinal, transversal, and normal).

figure o

The applied stimulus current, delivered to a volume of \((1.5 \times 1.5 \times 1.5) \times 10^{-3}\) m, situated at one corner of the slab, is prescribed in the subsection Applied current. The stimulus has a duration of \(2 \times 10^{-3}\) s and an amplitude of 35.714 V/s.

figure p

In Fig. 6(a) and Fig. 7(a), we display a snapshot of the transmembrane potential and the total activation time computed as output of the numerical simulation, respectively. The computation of the activation time is evaluated in a given point in the cardiac muscle as the time when the transmembrane potential derivative \(\frac{\partial u}{\partial t}\) reaches its maximum value. This can be enabled in the parameter file under the subsection named Activation time.

The problem was solved using both tetrahedral and hexahedral conforming meshes and the BDF2 scheme, using eight combinations of spatial resolutions (\(dx = [0.5, 0.2, 0.1, 0.05] \times 10^{-3}\) m) and time steps (\(\Delta t = [0.05, 0.01, 0.005, 0.001] \times 10^{-3}\) s), as reported in Fig. 8. The results showcase the activation time at points along the diagonal line of the slab, as also shown in Fig. 9(a). For both the tetrahedral and hexahedral meshes, the numerical solutions converge towards the finer space-time discretization (\(\Delta t = 0.001 \times 10^{-3}\) s, \(dx = 0.05 \times 10^{-3}\) m), yielding a latest activation time of 41.8 ms and 42.0 ms, respectively. These findings align with the values reported in the original N-version benchmark paper [84].

To further evaluate the \(\texttt {life}^{\text{x}}\)-ep results in comparison to the other eleven codes participating in the N-version benchmark [84], we report in Fig. 9 the activation time (for \(\Delta t = 0.005 \times 10^{-3}\) s, \(dx = 0.1 \times 10^{-3}\) m) at specific points along the slab diagonal (namely P1-P8-P9, see Fig. 9(a)) for all the codes, including \(\texttt {life}^{\text{x}}\)-ep. The activation times at points P1-P8-P9 for all the other eleven codes are available in the electronic supplementary material of [84]. As shown in Fig. 9(b), the \(\texttt {life}^{\text{x}}\)-ep results fall within the range of activation time values computed by the other codes. Moreover, at such refinement level (\(dt = 0.005 \times 10^{-3}\) s, \(dx = 0.1 \times 10^{-3}\) m), the \(\texttt {life}^{\text{x}}\)-ep numerical solutions (in both tetrahedral and hexahedral meshes) are in excellent agreement with the reference converged activation time values of 42-43 ms computed at point P8, where the accumulation of errors tends to occur as the wave propagates across the cuboid.

Idealized left atrium

We simulate the propagation of the electrophysiology wavefront in an idealized left atrial geometry (see Fig. 4(b)) using the APf ionic model [42]. To run the simulation, the user can modify the generic parameter file lifex_electrophysiology.prm or use the predefined file params_ideal_la.prm available in the release archive.

figure q

We employ a second-order BDF temporal scheme with a time step \(\Delta t=10^{-4}\) s and a final time \(T=0.15\) s. Furthermore, we initialize the ionic model by running a 1000-cycle long single-cell simulation applying a stimulus with period 0.8 s.

figure r

The simulation is initiated using a single spherical impulse with a radius of \(r=3~\times ~10^{-3}\) m. The parameters related to the impulse site, amplitude, duration, and initial time are specified in the Applied current/Spherical subsection of the parameter file. Notice that, when utilizing the APf model, it is essential to properly rescale the impulse amplitude and also conductivity values embedded within the diffusion tensor (3) in accordance with the model’s formulation [42]. Finally, the myofiber architecture is prescribed using the atrial LDRBM described in [17], as depicted in Fig. 5(b), by specifying Geometry type = Left atrium in the subsection Fiber generation. It is important to remark that in the APf model formulation, the transmembrane potential u is a dimensionless variable ranging from 0 to 1. However, for visualization purposes, the actual transmembrane potential \(u_{\text {mV}}\) in millivolts is obtained by postprocessing the numerical results using the formula \(u_{\text {mV}} = (100u - 80)~\)[mV] [42]. This conversion is applied to the data shown in Fig. 6(b).

Idealized left ventricle

We simulate the electrophysiology wavefront propagation in an idealized left ventricular hexahedral mesh, shown in Fig. 4(c), using the Bueno-Orovio (BO) ionic model [43]. The simulation can be run either by modifying the generic parameter file lifex_electrophysiology.prm or by utilizing the dedicated parameter file params_ideal_lv.prm.

figure s

This test case can be executed at various levels of hierarchical grid refinements by modifying the parameter value Number of refinements in the Mesh and space discretization subsection. We remark that this is possible for all \(\texttt {life}^{\text{x}}\) simulations conducted with hexahedral meshes, but is not available for tetrahedral meshes, due to lack of support in deal.II.

figure t

For this simulation, we use a second-order BDF temporal scheme with a time step of \(\Delta t=5 \times 10^{-5}\) s and a final time of \(T=0.15\) s. The ionic model is initialized using 1000-cycle long single-cell simulations with a cardiac period of 0.8 s. We employ a pacing protocol where three ventricular endocardial areas are activated using Gaussian impulses [85]. The characteristics of these impulses, such as amplitude, duration, and initial time, can be set in the file parameter subsection Applied current/Gaussian.

figure u

The myofiber architecture is prescribed using the Rossi-Lassila (RL) ventricular LDRBM [17], see Fig. 5(c), by specifying Geometry type = Left ventricle and setting Algorithm type = RL in the subsection Fiber generation.

The simulation results are reported in Fig. 6(c) and Fig. 7(c), where a snapshot of the transmembrane potential and the total activation time are displayed, respectively. Note that the transmembrane potential \(u_{\text {mV}}\) shown in Fig. 6(c) is obtained by postprocessing the numerical results using the formula \(u_{\text {mV}}=(85.7u-84)~[mV]\) [43].

Ventricular slab

This test case serves as an explanatory example for the multi-volume simulation framework of \(\texttt {life}^{\text{x}}\)-ep. We simulate a portion of an idealized left ventricular geometry, also referred to as ventricular slab [15], where the computational domain is divided into three volumetric regions: Sub-endocardial, Myocardial, and Sub-epicardial layers, as shown in Fig. 4(d).

To perform multi-volume simulations, the parameter template is created by specifying the volume labels on the command line as follows:

figure v

Doing so, the file lifex_electrophysiology.prm will contain three subsections named Sub Endocardium, Myocardium and Sub Epiucardium, where the volumetric tags (Material IDs) and all volume-specific parameters are prescribed. We use the TTP06 ionic model with a different Cell type for each layer (Endocardium, Myocardium and Epicardium) [44].

figure w

The simulation can then be run by providing on the command line the same volume labels used when generating the parameter file:

figure x

A second-order BDF temporal scheme is employed in this simulation, with a time step of \(\Delta t=5 \times 10^{-5}\) s and a final time of \(T=0.12\) s. The ionic model is initialized in each volume by conducting 1000-cycle long single-cell simulations with a cardiac period of 0.8 s. One single-cell simulation is run for every volume, so that the initialization is consistent with the spatially-varying parameters. The simulation onset is obtained with a single spherical impulse, of radius \(r=2.5 \times 10^{-3}\) m. The myofiber architecture is prescribed in the Fiber generation subsection by utilizing the Slab LDRBM, as shown in Fig. 5(d). The fiber orientations are set to exhibit a linear transmural variation from \(60^{\circ }\) to \(-60^{\circ }\), passing from the endocardial to the epicardial surface.

The simulation results are shown in Fig. 6(d) and Fig. 7(d), which display a snapshot of the transmembrane potential and the total activation time, respectively.

Simulations of pathological cardiac electrophysiology

In the following section, we present the pathological electrophysiology simulations applied to realistic left atrial and left ventricular geometries.

Realistic pathological left atrium

We perform a simulation of reentrant drivers typical of AF in a realistic left atrial tetrahedral mesh [45, 82], as shown in Fig. 4(e, h).

In this pathological scenario, four volumes have been introduced, named Healthy, Fibrosis Mild, Fibrosis and Scar. We use the CRN ionic model adopting different ionic conductances and conductivity values to model the pathophysiological behaviour in the fibrotic regions, that are labelled according to their IIR [82]. Specifically, in the Fibrosis Mild volumetric region, we changed, with respect to the default CRN values, the parameters related to the transient outward current conductance gto, the L-type intracellular \(\text {Ca}^{2+}\) current conductance gCaL, and the delayed rectifier current represented by gKur_fix and gKur_var, to model the effects of chronic AF. In the Fibrosis region, we set specific current conductance values for gK1, gCaL and gNa, which stand for the inward rectifier, the L-type calcium, potassium and sodium current conductances, respectively. Conductances and conductivities were adjusted according to [48]. An additional region, labeled as Healthy, models healthy tissue, without modification to ionic parameters. Finally, we consider a region labelled as Scar, which we model as non-conductive, thus excluding that region from the domain of the monodomain and ionic models and creating a block in the conduction (see Section Action potential propagation models).

figure y

For the spatial discretization we employ second-order FEs by setting FE space degree = 2. We use a first-order BDF temporal scheme with a time step \(\Delta t=5~\times ~10^{-5}\) s and a final time \(T=1.5\) s. We initialize the ionic model by running a 24-cycle long single-cell simulation using an impulse period of 0.5 s [48]. We use a pacing protocol with a sequence of four multiple spherical impulses delivered every \(220 \times 10^{-3}\) s. This timing is set under the parameter set Impulse initial times in the Applied current/Spherical subsection of the parameter file. Finally, the myofiber architecture is prescribed using the atrial LDRBM [17], as depicted in Fig. 5(e), by specifying Geometry type = Left atrium in the subsection Fiber generation.

The simulation can be run using

figure z

The simulation results are presented in Fig. 6(e) and Fig. 7(e), which display a snapshot of the transmembrane potential and the total activation time, respectively. In this videoFootnote 39 (online version) we show the evolution of the transmembrane potential.

Realistic pathological left ventricle

We simulate a Ventricular Tachycardia (VT) macro-reentrant circuit in a realistic left ventricle [83], as shown in Fig. 4(f, i).

To account for the presence of grey zone and scar, we consider three volumes: Healthy, Border Zone and Scar. Different ionic conductances in the TTP06 ionic model and conductivity values in the monodomain equation are used to characterize the electrophysiology properties within each region. In particular, in the Border Zone, we reduce the conductances of the peak sodium, L-type calcium, and potassium currents GNa, GCaL, Gkr and Gks_myo, according to [49]. The Healthy region has a normal configuration, while the Scar region is modeled as non-conductive (setting Disable conduction = true).

figure aa

We use second-order FEs for the spatial discretization, with one level of hierarchical grid refinement (Number of refinements = 1), and a second-order BDF temporal scheme with a time step of \(\Delta t=5 \times 10^{-5}\) s and a final time of \(T=2.5\) s. The ionic model is initialized using 1000-cycle long single-cell simulations with a cardiac period of 0.8 s. We employ a pacing protocol with a sequence of multiple spherical impulses delivered in a specific ventricular endocardial area. The myofiber architecture is prescribed using the RL ventricular algorithm [17], see Fig. 5(f), by specifying Geometry type = Left ventricle and setting Algorithm type = RL in the subsection Fiber generation.

The simulation can be run using

figure ab

The simulation results are presented in Fig. 6(f) and Fig. 7(f), which display a snapshot of the transmembrane potential and the total activation time, respectively. In this videoFootnote 40 (online version) we show the evolution of the transmembrane potential.

Strong scalability test

We consider the slab benchmark of [84], and discretize the domain with a structured hexahedral mesh of 47185920 elements and 47744577 nodes, corresponding to a mesh size of approximately \(dx = {3.6 \times 10^{-5}}{\hbox {m}}\). We set \(\Delta t = {1 \times 10^{-4}}{\hbox {s}}\) and \(T = {5\times 10^{-2}}{\hbox {s}}\), and we use linear FEs and the second-order BDF scheme for time discretization.

Fig. 1
figure 1

The \(\texttt {life}^{\text{x}}\) library encompasses essential features and a framework for numerically solving Finite Element problems. \(\texttt {life}^{\text{x}}\)-ep is a publicly released package designed for cardiac electrophysiology simulations based on \(\texttt {life}^{\text{x}}\). Left picture icons sourced from https://fontawesome.com/license

Fig. 2
figure 2

Multiscale cardiac electrophysiology model. From the largest to the smallest spatial scale: a organ; b tissue; c cell; d membrane

Fig. 3
figure 3

a Computational domain \(\Omega _{} \subset {\mathbb {R}}^3\) in an idealized slab of ventricular myocardium, partitioned in three subdomains (\({\overline{\Omega }}_{} = {\overline{\Omega }}_{1} \cup {\overline{\Omega }}_{2} \cup {\overline{\Omega }}_{3}\)), where \(\Gamma _{12}\) and \(\Gamma _{23}\) denote the interfaces among the corresponding subdomains; (b) Representation of the orthonormal triplet for the myofiber orientations \({{\textbf{f}}_0}\) (fiber, in red), \({{\textbf{s}}_0}\) (sheet, in blu) and \({{\textbf{n}}_0}\) (sheet-normal, in green) in an idealized slab of ventricular tissue

Fig. 4
figure 4

Domains and meshes used in the numerical examples. Domains in the top row ac are composed of a single subdomain, while for domains in the middle row d–f, zoomed on the bottom gi, colors are used to differentiate the subdomains

Fig. 5
figure 5

Fiber field computed using Laplace-Dirichlet Rule-Based Methods [15, 17] visualized as streamlines: a Slab tissue; b Idealized left atrium; c Idealized left ventricle; d Ventricular slab; e Realistic left atrium; f Realistic left ventricle. The Laplace solution \(\phi\) is the transmural function where \(\phi =0\) on the endocardium and \(\phi =1\) on the epicardium

Fig. 6
figure 6

Snapshots of the transmembrane potential for all test cases: a Slab tissue; b Idealized left atrium; c Idealized left ventricle; d Ventricular slab; e Realistic left atrium; f Realistic left ventricle

Fig. 7
figure 7

Activation maps computed for all the test cases: a Slab tissue; b Idealized left atrium; c Idealized left ventricle; d Ventricular slab; e Realistic left atrium; f Realistic left ventricle

Fig. 8
figure 8

Activation times evaluated along the cuboid diagonal line in the N-version benchmark problem [84] for all the numerical solutions performed with \(\texttt {life}^{\text{x}}\)-ep at different refinements in space and time. Red lines=Hexahedral simulations; Blue lines=Tetrahedral simulations

Fig. 9
figure 9

a Activation times evaluted along the cuboid diagonal and in the points P1, P9 and P8 in the N-version benchmark problem [84]. b Comparison of the \(\texttt {life}^{\text{x}}\)-ep numerical solutions (with Red line=Hexahedral mesh and Blue line=Tetrahedral mesh) with respect to the other codes partecipating to the benchmark problem [84]

Fig. 10
figure 10

Computational time (left) and parallel speedup (right) against the number of cores for the strong scalability test. Dashed lines indicate the ideal linear scaling

Fig. 11
figure 11

Total memory occupation (left) and memory occupation per core with three differently refined meshes

We run a strong scalability test, varying the number of parallel processes used in the computation and measuring the wall time necessary for the solution of the ionic models, the assembly of the monodomain system and the solution of the linear system. The test was run on the CINECA GALILEO100 supercomputer.

The wall times, plotted in Fig. 10 for the different steps of the solver, scales linearly up to over 1000 cores, confirming the results of [13] on the scalability properties of the \(\texttt {life}^{\text{x}}\) core components. We report in Table 1 a breakdown of the computational cost of the different sections for the simulation with 960 cores. Most of the computational time is spent for the solution of the linear system (6). The matrix of the linear system is assembled only once (since it is the same at every time step), and the right-hand-side is efficiently recomputed at every time step by means of matrix–vector products, resulting in a very small computational cost for the system assembly phase. Moreover, since the ionic model is solved independently at each degree of freedom, the associated computational cost scales almost perfectly, and becomes very small if a sufficiently large number of parallel processes is employed.

Table 1 Summary of the computational costs for the strong scalability test, using 960 parallel cores. For each section, we report the total wall time, the wall time for each time step and the wall time relative to the total. Sections are sorted in descending order of cost

Memory occupation

We consider the same setting used for the scalability test, with three differently refined meshes of 772497, 6038305 and 47744577 nodes (corresponding to a mesh size of approximately \({1.5 \times 10^{-4}}{\hbox {m}}\), \({7.3 \times 10^{-5}}{\hbox {m}}\) and \({3.6 \times 10^{-5}}{\hbox {m}}\), respectively). We measure the peak memory occupation of the simulation, varying the number of parallel processes. Results are reported in Fig. 11. We observe that the total memory occupation increases as the number of processes increases. This is expected due to distributed memory parallelization, that requires some duplicated information across processes. However, as shown in Fig. 11 (right), this overhead grows more slowly than the number of processes itself. Therefore, increasing the number of computing cores, which generally implies increasing the number of available RAM, does not lead to out-of-memory issues.

Output and visualization

Two different types of output file formats are available in the \(\texttt {life}^{\text{x}}\)-ep release: HDF5 and csv. Both of them can be enabled and configured in the Output subsections:

figure ac

The HDF5 output is available in the following subsections:

  • Electrophysiology/Output,

  • Electrophysiology/Activation time,

  • Fiber generation/Output.

This generates an XDMF file named output_filename.xdmf (which links to a corresponding HDF5 output file output_filename.h5). These files can be visualized using ParaView,Footnote 41 an open-source multi-platform data analysis and visualization application, see e.g., Fig. 5,  6 and  7. The HDF5 format ensures that the output can be easily post-processed, not only for visualization purposes but also as input for more advanced computational pipelines.

The Comma-Separated Values (csv) format consists of delimited text files where values are separated by commas, with each line representing a specific data record. The csv files can be found in different subsections:

  • Ionic model parameters/Output,

  • Ionic model parameters/0D Output,

  • Electrophysiology/Output.

These csv output files can be conveniently used to plot electrophysiology variable (min, max and pointwise) values over time in the computed numerical simulation.

Conclusions

In this work, we introduced \(\texttt {life}^{\text{x}}\)-ep, a robust and advanced software specifically designed for simulating the electrophysiology activity of the cardiac muscle. With the goal of addressing the computational challenges associated with cardiac simulations, \(\texttt {life}^{\text{x}}\)-ep provides efficient numerical methods while maintaining precision and accuracy. \(\texttt {life}^{\text{x}}\)-ep incorporates a numerical solver for the monodomain equation coupled with both phenomenological and second-generation ionic models, namely Aliev-Panfilov [42], Bueno-Orovio [43], TTP06 [44], and CRN [45]. These models are discretized in time using the BDF scheme and in space using the FE method of orders 1 and 2 on tetrahedral meshes, and of arbitrary degree on hexahedral meshes, thus providing a comprehensive framework for modeling the electrical activity of the heart under both normal and pathological conditions.

Leveraging the capabilities of \(\texttt {life}^{\text{x}}\), \(\texttt {life}^{\text{x}}\)-ep provides users with a user-friendly and flexible interface, facilitated by self-documenting parameter files for easy simulation setup. For enhanced accessibility, \(\texttt {life}^{\text{x}}\)-ep is distributed in an AppImage binary format, rendering it universally compatible with any recent x86-64 Linux system. Researchers from diverse backgrounds, such as medicine and bio-engineering, can readily access and utilize \(\texttt {life}^{\text{x}}\)-ep for in-silico simulations. The underlying principles and structure of \(\texttt {life}^{\text{x}}\)-ep can be readily understood thanks to the comprehensive technical and mathematical documentation.

As unique and distinctive features, \(\texttt {life}^{\text{x}}\)-ep provides two options for prescribing myocardial fibers. Users can either import them from a file or generate them online by exploiting the LDRBMs presented in [17] and implemented in the previous release \(\texttt {life}^{\text{x}}\)-fiber [15]. Moreover, it supports spatial heterogeneity in the choice of both models and physical coefficients, easily configurable through a convenient parameter file, without the need to access and modify the source code.

\(\texttt {life}^{\text{x}}\)-ep benefits from its high-performance computing capabilities, achieving ideal parallel speedup on thousands of cores. The accuracy and reliability of \(\texttt {life}^{\text{x}}\)-ep have been verified through a benchmark for computational electrophysiology, ensuring the validity of its results. Furthermore, a range of idealized and realistic cardiac simulations in both normal and pathological settings highlights its capabilities and versatility in capturing complex cardiac dynamics and its potential for patient-specific studies. \(\texttt {life}^{\text{x}}\)-ep offers the capability to facilitate the simulation of pathological scenarios, allowing the creation of scars, grey zones, and an arbitrary number of conduction “levels”. This potential impact in the study of pathologies can prove to be highly valuable.

In conclusion, \(\texttt {life}^{\text{x}}\)-ep provides to the scientific community a comprehensive, high-performance, and user-friendly software for conducting in-silico cardiac electrophysiology simulations.

In the future, efforts will be made to further improving the accuracy and efficiency of \(\texttt {life}^{\text{x}}\)-ep. One possible approach is the adoption of lookup tables instead of repeatedly evaluating the expensive functions in the ionic current term, which could provide a speed-up, although with a potential impact on accuracy [86]. Inexact solvers, such as those based on domain decomposition methods, show promising optimality and scalability properties [87]. Recent studies have also highlighted the advantages of employing high-order discretization schemes for accurately capturing the intricate electrical wavefront propagation observed in cardiac electrophysiology [19]. These schemes not only offer improved accuracy but also enable the implementation of efficient matrix-free solvers that require minimal memory usage [19]. This opens up the possibility of leveraging GPU architectures for accelerated computations [88]. Additionally, the use of higher-order or adaptive time stepping schemes and \(hp\)-adaptive FE holds the potential to achieve greater accuracy in simulations while optimizing computational efficiency [89]. These developments aim to further enhance the capabilities of \(\texttt {life}^{\text{x}}\)-ep and expand its range of features.

Availability and requirements

Project name: \(\texttt {life}^{\text{x}}\)-ep

Project home page: https://lifex.gitlab.io/

Operating system(s): Linux (x86-64)

Programming language: C++

Other requirements: glibc version 2.28 or higher

License: CC BY-NC-ND 4.0

Any restrictions to use by non-academics: no additional restriction.

Availability of data and materials

All input data, meshes and the binary executable of \(\texttt {life}^{\text{x}}\)-ep can be found at https://doi.org/10.5281/zenodo.8085266.

Notes

  1. https://www.numericor.at/

  2. https://appimage.org/

  3. https://www.debian.org/releases/

  4. https://docs.appimage.org/introduction/concepts.html

  5. https://www.gnu.org/software/libc/

  6. https://www.kernel.org/doc/html/latest/filesystems/fuse.html

  7. https://docs.appimage.org/user-guide/troubleshooting/fuse.html

  8. http://creativecommons.org/licenses/by-nc-nd/4.0/

  9. https://lifex.gitlab.io/

  10. https://www.dealii.org/

  11. https://trilinos.github.io/

  12. https://www.boost.org/

  13. https://vtk.org/

  14. https://github.com/coin-or/ADOL-C

  15. https://github.com/opencollab/arpack-ng

  16. https://www.netlib.org/blacs/

  17. https://eigen.tuxfamily.org/

  18. https://www.fftw.org/

  19. https://www.gnu.org/software/glpk/

  20. https://www.hdfgroup.org/solutions/hdf5/

  21. https://www.llnl.gov/casc/hypre/

  22. http://glaros.dtc.umn.edu/gkhome/metis/metis/overview

  23. https://www.mpich.org/

  24. http://mumps.enseeiht.fr/index.php?page=home

  25. https://www.unidata.ucar.edu/software/netcdf/

  26. https://www.openblas.net/

  27. https://www.mcs.anl.gov/petsc/

  28. http://glaros.dtc.umn.edu/gkhome/metis/parmetis/overview

  29. https://www.netlib.org/scalapack/

  30. https://gitlab.inria.fr/scotch/scotch

  31. https://people.engr.tamu.edu/davis/suitesparse.html

  32. https://portal.nersc.gov/project/sparse/superlu/

  33. https://oneapi-src.github.io/oneTBB/

  34. https://www.p4est.org/

  35. https://gmsh.info/doc/texinfo/gmsh.html

  36. https://cubit.sandia.gov/

  37. https://zenodo.org/record/5801337

  38. https://kcl.figshare.com/articles/dataset/A_Virtual_Cohort_of_Twenty-four_Left-ventricular_Models_of_Ischemic_Cardiomyopathy_Patients/16473903

  39. https://polimi365-my.sharepoint.com/:v:/g/personal/10594253_polimi_it/EcRfIyP2ya5EsD6IixPi_4ABCaZfJViTwwgeDpKBzCsquw?e=VUMMP9

  40. https://polimi365-my.sharepoint.com/:v:/g/personal/10594253_polimi_it/EUOytZLaj1VGnsxoZxuZcb4BXSqKaGG63oUZPN6EPVt67g?e=77zmAP

  41. https://www.paraview.org

Abbreviations

FE:

Finite Element

BDF:

Backward differentiation formula

TTP06:

ten Tusscher-Panfilov 2006

APf:

Aliev-Panfilov

BO:

Bueno-Orovio

CRN:

Courtemanche-Ramirez-Nattel

AF:

Atrial Fibrillation

VT:

Ventricular Tachycardia

CRT:

Cardiac Resynchronization Therapy

MRI:

Magnetic Resonance Imaging

LGE-MRI:

Late Gadolinium Enhancement-Magnetic Resonance Imaging

AP:

Action Potential

CV:

Conduction Velocity

VT:

Ventricular Tachycardia

LBBB:

Left Bundle Branch Block

ECG:

Electrocardiogram

BSPM:

Body Surface Potential Mapping

CT:

Computed Tomography

ODE:

Ordinary Differential Equation

PDE:

Partial Differential Equation

ICI:

Ionic Current Interpolation

RL:

Rossi-Lassila

LDRBM:

Laplace–Dirichlet Rule-Based Method

IIR:

Imaging Itensity Ratio

MPI:

Message Passing Interface

References

  1. Tortora GJ, Derrickson BH. Principles of anatomy and physiology. London: Wiley; 2008.

    Google Scholar 

  2. Katz AM. Physiology of the heart. Wilkins: Lippincott Williams; 2010.

    Google Scholar 

  3. Klabunde R. Cardiovascular physiology concepts. Wilkins: Lippincott Williams; 2011.

    Google Scholar 

  4. Harrington RA, Narula J, Eapen ZJ. Hurst’s the Heart. MacGraw-Hill 2011.

  5. Quarteroni A, Dede’ L, Manzoni A, Vergara C. Mathematical modelling of the human cardiovascular system: data, numerical approximation. Clinical Applications: Cambridge University Press; 2019.

    Book  Google Scholar 

  6. Niederer SA, Lumens J, Trayanova NA. Computational models in cardiology. Nat Rev Cardiol. 2019;16:100–11.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. ...Corral-Acero J, Margara F, Marciniak M, Rodero C, Loncaric F, Feng Y, Gilbert A, Fernandes JF, Bukhari HA, Wajdan A, Martinez MV, Santos MS, Shamohammdi M, Luo H, Westphal P, Leeson P, DiAchille P, Gurev V, Mayr M, Geris L, Pathmanathan P, Morrison T, Cornelussen R, Prinzen F, Delhaas T, Doltra A, Sitges M, Vigmond EJ, Zacur E, Grau V, Rodriguez B, Remme EW, Niederer S, Mortier P, McLeod K, Potse M, Pueyo E, Bueno-Orovio A, Lamata P. The “Digital Twin” to enable the vision of precision cardiology. Eur Heart J. 2020;41(48):4556–64.

    Article  PubMed  PubMed Central  Google Scholar 

  8. Peirlinck M, Costabal FS, Yao J, Guccione JM, Tripathy S, Wang YS, Ozturk D, Segars P, Morrison TM, Levine S. Precision medicine in human heart modeling. Biomech Model Mechanobiol. 2021;20:803–31.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Arevalo HJ, Vadakkumpadan F, Guallar E, Jebb A, Malamas P, Wu KC, Trayanova NA. Arrhythmia risk stratification of patients after myocardial infarction using personalized heart models. Nat Commun. 2016;7:11437.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Prakosa A, Arevalo HJ, Deng D, Boyle PM, Nikolov PP, Ashikaga H, Blauer JJE, Ghafoori E, Park CJ, Blake RC, Han FT, MacLeod RS, Halperin HR, Callans DJ, Ranjan R, Chrispin J, Nazarian S, Trayanova NA. Personalized virtual-heart technology for guiding the ablation of infarct-related ventricular tachycardia. Nat Biomed Eng. 2018;2:732–40.

    Article  PubMed  PubMed Central  Google Scholar 

  11. Strocchi M, Lee AWC, Neic A, Bouyssier J, Gillette K, Plank G, Elliott MK, Gould J, Behar JM, Sidhu B, Mehta V, Bishop MJ, Vigmond EJ, Rinaldi CA, Niederer SA: His-bundle and left bundle pacing with optimized atrioventricular delay achieve superior electrical synchrony over endocardial and epicardial pacing in left bundle branch block patients. Heart Rhythm 1922-1929 (2020)

  12. Campos FO, Neic A, Mendonca Costa C, Whitaker J, O’Neill M, Razavi R, Rinaldi CA, Scherr D, Niederer SA, Plank G, Bishop MJ. An automated near-real time computational method for induction and treatment of scar-related ventricular tachycardias. Med Image Anal. 2022;80: 102483.

    Article  PubMed  PubMed Central  Google Scholar 

  13. Africa PC. lifex: a flexible, high performance library for the numerical solution of complex finite element problems. SoftwareX. 2022;20: 101252.

    Article  Google Scholar 

  14. Arndt D, Bangerth W, Blais B, Fehling M, Gassmöller R, Heister T, Heltai L, Köcher U, Kronbichler M, Maier M, Munch P, Pelteret JP, Proell S, Konrad S, Turcksin B, Wells D, Zhang J. The dealII library, version 9.3. J Numer Math 2021; 29(3), 171–186.

  15. Africa PC, Piersanti R, Fedele M, Dede’ L, Quarteroni A. lifex-fiber: an open tool for myofibers generation in cardiac computational models. BMC Bioinformatics. 2023;24(1):143.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Africa PC, Fumagalli I, Bucelli M, Zingaro A, Fedele M, Dede’ L, Quarteroni A. lifex-cfd: an open-source computational fluid dynamics solver for cardiovascular applications. arXiv preprint arXiv:2304.12032v3 2023.

  17. Piersanti R, Africa PC, Fedele M, Vergara C, Dede’ L, Corno AF, Quarteroni A. Modeling cardiac muscle fibers in ventricular and atrial electrophysiology simulations. Comput Methods Appl Mech Eng. 2021;373: 113468.

    Article  Google Scholar 

  18. Pagani S, Dede’ L, Frontera A, Salvador M, Limite LR, Manzoni A, Lipartiti F, Tsitsinakis G, Hadjis A, Della Bella P, Quarteroni A. A computational study of the electrophysiological substrate in patients suffering from atrial fibrillation. Front Physiol 12,2021.

  19. Africa PC, Salvador M, Gervasio P, Dede’ L, Quarteroni A. A matrix-free high-order solver for the numerical solution of cardiac electrophysiology. J Comput Phys. 2023;478: 111984.

    Article  Google Scholar 

  20. Salvador M, Regazzoni F, Pagani S, Dede’ L, Trayanova N, Quarteroni A. The role of mechano-electric feedbacks and hemodynamic coupling in scar-related ventricular tachycardia. Comput Biol Med 2022; 105203.

  21. Cicci L, Fresca S, Manzoni A. Deep-hyromnet: a deep learning-based operator approximation for hyper-reduction of nonlinear parametrized pdes. J Sci Comput. 2022;93(2):57.

    Article  Google Scholar 

  22. Cicci L, Fresca S, Manzoni A, Quarteroni A. Efficient approximation of cardiac mechanics through reduced order modeling with deep learning-based operator approximation. arXiv preprint arXiv:2202.03904 2022.

  23. Cicci L, Fresca S, Pagani S, Manzoni A, Quarteroni A. Projection-based reduced order models for parameterized nonlinear time-dependent problems arising in cardiac mechanics. Math Eng. 2022;5(2):1–38.

    Article  Google Scholar 

  24. Fedele M, Piersanti R, Regazzoni F, Salvador M, Africa PC, Bucelli M, Zingaro A, Quarteroni A, et al. A comprehensive and biophysically detailed computational model of the whole human heart electromechanics. Comput Methods Appl Mech Eng. 2023;410: 115983.

    Article  Google Scholar 

  25. Salvador M, Fedele M, Africa PC, Sung E, Prakosa A, Chrispin J, Trayanova N, Quarteroni A. Electromechanical modeling of human ventricles with ischemic cardiomyopathy: numerical simulations in sinus rhythm and under arrhythmia. Comput Biol Med. 2021;136: 104674.

    Article  PubMed  Google Scholar 

  26. Regazzoni F, Salvador M, Africa PC, Fedele M, Dede’ L, Quarteroni A. A cardiac electromechanical model coupled with a lumped-parameter model for closed-loop blood circulation. J Comput Phys 2022; 111083.

  27. Piersanti R, Regazzoni F, Salvador M, Corno AF, Vergara C, Quarteroni A. 3D–0D closed-loop model for the simulation of cardiac biventricular electromechanics. Comput Methods Appl Mech Eng. 2022;391: 114607.

    Article  Google Scholar 

  28. Corti M, Dede’ L, Zingaro A, Quarteroni A. Impact of atrial fibrillation on left atrium haemodynamics: a computational fluid dynamics study. Comput Biol Med 2022,106143.

  29. Zingaro A, Fumagalli I, Dede’ L, Fedele M, Africa PC, Corno AF, Quarteroni AM. A geometric multiscale model for the numerical simulation of blood flow in the human left heart. Discrete Contin Dyn Syst S. 2022;15(8):2391–427.

    Article  Google Scholar 

  30. Zingaro, A., Bucelli, M., Piersanti, R., Regazzoni, F., Dede’, L., Quarteroni, A.: An electromechanics-driven fluid dynamics model for the simulation of the whole human heart. arXiv preprint arXiv:2301.02148 (2023)

  31. Fumagalli, I., Vitullo, P., Vergara, C., Fedele, M., Corno, A.F., Ippolito, S., Scrofani, R., Quarteroni, A.: Image-based computational hemodynamics analysis of systolic obstruction in hypertrophic cardiomyopathy. Front Physiol 2437 (2022)

  32. Marcinno’, F., Zingaro, A., Fumagalli, I., Dede’, L., Vergara, C.: A computational study of blood flow dynamics in the pulmonary arteries. Vietnam J Math 1–23 (2022)

  33. Bennati, L., Vergara, C., Giambruno, V., Fumagalli, I., Corno, A.F., Quarteroni, A., Puppini, G., Luciani, G.B.: An image-based computational fluid dynamics study of mitral regurgitation in presence of prolapse. Cardiovasc Eng Technol 1–19 (2023)

  34. Bennati, L., Giambruno, V., Renzi, F., Di Nicola, V., Maffeis, C., Puppini, G., Luciani, G.B., Vergara, C.: Turbulence and blood washout in presence of mitral regurgitation: a computational fluid-dynamics study in the complete left heart. bioRxiv, 2023–03 (2023)

  35. Zingaro, A., Bucelli, M., Fumagalli, I., Dede’, L., Quarteroni, A.: Modeling isovolumetric phases in cardiac flows by an augmented resistive immersed implicit surface method. arXiv preprint arXiv:2208.09435 (2022)

  36. Bucelli M, Dede’ L, Quarteroni A, Vergara C. Partitioned and monolithic algorithms for the numerical solution of cardiac fluid-structure interaction. Commun Comput Phys. 2022;32(5):1217–56.

    Article  Google Scholar 

  37. Bucelli M, Zingaro A, Africa PC, Fumagalli I, Dede’ L, Quarteroni AM. A mathematical model that integrates cardiac electrophysiology, mechanics and fluid dynamics: application to the human left heart. Int J Numer Methods Biomed Eng. 2023;39(3):3678.

    Article  Google Scholar 

  38. Bucelli, M., Gabriel, M.G., Quarteroni, A., Gigante, G., Vergara, C.: A stable loosely-coupled scheme for cardiac electro-fluid-structure interaction. J Comput Phys 112326 (2023)

  39. Di Gregorio, S., Vergara, C., Pelagi, G.M., Baggiano, A., Zunino, P., Guglielmo, M., Fusini, L., Muscogiuri, G., Rossi, A., Rabbat, M.G., et al.: Prediction of myocardial blood flow under stress conditions by means of a computational model. Eurp J Nuclear Med Mol Imaging, 1–12 (2022)

  40. Zingaro A, Vergara C, Dede’ L, Regazzoni F, Quarteroni A. A comprehensive mathematical model for cardiac perfusion. Sci Rep. 2023;13:14220.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Franzone PC, Pavarino LF, Scacchi S. Mathematical cardiac electrophysiology. Cham: Springer; 2014.

    Book  Google Scholar 

  42. Aliev RR, Panfilov AV. A simple two-variable model of cardiac excitation. Chaos, Solitons & Fractals. 1996;7(3):293–301.

    Article  Google Scholar 

  43. Bueno-Orovio A, Cherry EM, Fenton FH. Minimal model for human ventricular action potentials in tissue. J Theor Biol. 2008;253(3):544–60.

    Article  PubMed  Google Scholar 

  44. ten Tusscher KH, Panfilov AV. Alternans and spiral breakup in a human ventricular tissue model. Am J Physiol Heart Circ Physiol. 2006;291:1088–100.

    Article  Google Scholar 

  45. Courtemanche M, Ramirez RJ, Nattel S. Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model. Am J Physiol Heart Circ Physiol. 1998;275(1):301–21.

    Article  Google Scholar 

  46. Frontera A, Pagani S, Limite LR, Hadjis A, Manzoni A, Dede’ L, Quarteroni A, Della Bella P. Outer loop and isthmus in ventricular tachycardia circuits: characteristics and implications. Heart Rhythm. 2020;17(10):1719–28.

    Article  PubMed  Google Scholar 

  47. Frontera A, Pagani S, Limite LR, Peirone A, Fioravanti F, Enache B, Silva JC, Vlachos K, Meyer C, Montesano G, Manzoni A, Dede’ L, Quarteroni A, Laṭcu DG, Rossi P, Della Bella P. Slow conduction corridors and pivot sites characterize the electrical remodeling in atrial fibrillation. JACC Clin Electrophysiol. 2022;8(5):561–77.

    Article  PubMed  Google Scholar 

  48. Zahid S, Cochet H, Boyle PM, Schwarz EL, Whyte KN, Vigmond EJ, Dubois R, Hocini M, Haïssaguerre M, Jaïs P, et al. Patient-derived models link re-entrant driver localization in atrial fibrillation to fibrosis spatial pattern. Cardiovasc Res. 2016;110(3):443–54.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Arevalo HJ, Vadakkumpadan F, Guallar E, Jebb A, Malamas P, Wu KC, Trayanova NA. Arrhythmia risk stratification of patients after myocardial infarction using personalized heart models. Nat Commun. 2016;7(1):11437.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Quarteroni A, Sacco R, Saleri F. Numerical mathematics. Cham: Springer; 2010.

    Google Scholar 

  51. Quarteroni A. Numerical models for differential problems. Cham: Springer; 2017.

    Book  Google Scholar 

  52. Pathmanathan P, Mirams GR, Southern J, Whiteley JP. The significant effect of the choice of ionic current integration method in cardiac electro-physiological simulations. Int J Numer Methods Biomed Eng. 2011;27(11):1751–70.

    Article  Google Scholar 

  53. Rossi S, Lassila T, Ruiz-Baier R, Sequeira A, Quarteroni A. Thermodynamically consistent orthotropic activation model capturing ventricular systolic wall thickening in cardiac electromechanics. Eurp J Mech-A/Solids. 2014;48:129–42.

    Article  Google Scholar 

  54. Hodgkin AL, Huxley AF. A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol. 1952;117(4):500.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Noble D. A modification of the hodgkin–huxley equations applicable to purkinje fibre action and pacemaker potentials. J Physiol. 1962;160(2):317.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Lines GT, Buist ML, Grottum P, Pullan AJ, Sundnes J, Tveito A. Mathematical models and numerical methods for the forward problem in cardiac electrophysiology. Comput Vis Sci. 2003;5:215–39.

    Article  Google Scholar 

  57. Clayton RH, Aboelkassem Y, Cantwell CD, Corrado C, Delhaas T, Huberts W, Lei CL, Ni H, Panfilov AV, Roney C, dos Santos RW. An audit of uncertainty in multi-scale cardiac electrophysiology models. Philos Trans R Soc Math Phys Eng Sci. 2020;378(2173):20190335.

    Google Scholar 

  58. Plank G, Loewe A, Neic A, Augustin C, Huang Y-L, Gsell MAF. The openCARP simulation environment for cardiac electrophysiology. Comput Methods Progr Biomed. 2021;208: 106223.

    Article  Google Scholar 

  59. Lloyd C, Lawson J, Hunter P, Nielsen P. The cellml model repository. Bioinformatics. 2008;24:2122–3.

    Article  CAS  PubMed  Google Scholar 

  60. Mirams GR, Arthurs CJ, Bernabeu MO, Bordas R, Cooper J, Corrias A, Davit Y, Dunn S, Fletcher AG, Harvey DG, Marsh ME, Osborne JM, Pathmanathan P, Pitt-Francis J, Southern J, Zemzemi N, Gavaghan DJ. Chaste: an open source c++ library for computational physiology and biology. PLoS Comput Biol. 2013;9(3):1002970.

    Article  Google Scholar 

  61. Richards DF, Glosli JN, Draeger EW, Mirin AA, Chan B, Fattebert J-L, Krauss WD, Oppelstrup T, Butler CJ, Gunnels JA, et al. Towards real-time simulation of cardiac electrophysiology in a human heart at high resolution. Comput Methods Biomech Biomed Engin. 2013;16(7):802–5.

    Article  PubMed  Google Scholar 

  62. Vázquez M, Houzeaux G, Koric S, Artigues A, Aguado-Sierra J, Arís R, Mira D, Calmet H, Cucchietti F, Owen H, Taha A, Burness ED, Cela JM, Valero M. Alya: multiphysics engineering simulation toward exascale. J Comput Sci. 2016;14:15–27.

    Article  Google Scholar 

  63. Baillargeon B, Rebelo N, Fox DD, Taylor RL, Kuhl E. The living heart project: a robust and integrative simulator for human heart function. Eur J Mech A Solids. 2014;48:38–47.

    Article  PubMed  PubMed Central  Google Scholar 

  64. Levine S, Battisti T, Butz B, D’Souza K, Costabal F, Peirlinck M: Dassault Systèmes’ Living Heart Project, pp. 245–259 (2022)

  65. Vigmond E, Dos Santos RW, Prassl A, Deo M, Plank G. Solvers for the cardiac bidomain equations. Prog Biophys Mol Biol. 2008;96(1–3):3–18.

    Article  CAS  PubMed  Google Scholar 

  66. Courtemanche M, Ramirez RJ, Nattel S. Ionic targets for drug therapy and atrial fibrillation-induced electrical remodeling: insights from a mathematical model. Cardiovasc Res. 1999;42(2):477–89.

    Article  CAS  PubMed  Google Scholar 

  67. Neic A, Campos FO, Prassl AJ, Niederer SA, Bishop MJ, Vigmond EJ, Plank G. Efficient computation of electrograms and ecgs in human whole heart simulations using a reaction-eikonal model. J Comput Phys. 2017;346:191–211.

    Article  PubMed  PubMed Central  Google Scholar 

  68. Gillette K, Gsell MAF, Prassl AJ, Karabelas E, Reiter U, Reiter G, Grandits T, Payer C, Štern D, Urschler M, Bayer JD, Augustin CM, Neic A, Pock T, Vigmond EJ, Plank G. A framework for the generation of digital twins of cardiac electrophysiology from clinical 12-leads ECGs. Med Image Anal. 2021;71: 102080.

    Article  PubMed  Google Scholar 

  69. Loewe A, Poremba E, Oesterlein T, Luik A, Schmitt C, Seemann G, Dössel O. Patient-specific identification of atrial flutter vulnerability-a computational approach to reveal latent reentry pathways. Front Physiol 2019;9.

  70. Azzolin L, Eichenlaub M, Nagel C, Nairn D, Sanchez J, Unger L, Dössel O, Jadidi A, Loewe A. Personalized ablation vs. conventional ablation strategies to terminate atrial fibrillation and prevent recurrence. EP Europace 2022.

  71. Gillette K, Gsell MAF, Bouyssier J, Prassl AJ, Neic A, Vigmond EJ, Plank G. Automated framework for the inclusion of a his-purkinje system in cardiac digital twins of ventricular electrophysiology. Ann Biomed Eng. 2021;49:3143–53.

    Article  PubMed  PubMed Central  Google Scholar 

  72. Sung E, Prakosa A, Aronis KN, Zhou S, Zimmerman SL, Tandri H, Nazarian S, Berger RD, Chrispin J, Trayanova NA. Personalized digital-heart technology for ventricular tachycardia ablation targeting in hearts with infiltrating adiposity. Circ Arrhythmia Electrophysiol. 2020;13(12):8912.

    Article  Google Scholar 

  73. Roney CH, Sim I, Yu J, Beach M, Mehta A, Solis-Lemus JA, Kotadia I, Whitaker J, Corrado C, Razeghi O, Vigmond E, Narayan SM, O’Neill M, Williams SE, Niederer SA. Predicting atrial fibrillation recurrence by combining population data and virtual cohorts of patient-specific left atrial models. Circ Arrhythmia Electrophysiol. 2022;15(2): 010253.

    Article  Google Scholar 

  74. Margara F, Wang ZJ, Levrero-Florencio F, Santiago A, Vázquez M, Bueno-Orovio A, Rodriguez B. In-silico human electro-mechanical ventricular modelling and simulation for drug-induced pro-arrhythmia and inotropic risk assessment. Prog Biophys Mol Biol. 2021;159:58–74.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Peirlinck M, Sahli Costabal F, Kuhl E. Sex differences in drug-induced arrhythmogenesis. Front Physiol. 2021;12: 708435.

    Article  PubMed  PubMed Central  Google Scholar 

  76. González-Martín P, Sacco F, Butakoff C, Doste R, Bederián C, Gutierrez Espinosa de los Monteros LK, Houzeaux G, Iaizzo PA, Iles TL, Vázquez M, Aguado-Sierra J. Ventricular anatomical complexity and gender differences impact predictions from computational models. PLoS ONE 2023;18(2), 0263639.

  77. Trilinos project website. https://trilinos.github.io 2023.

  78. Schäling B. The boost C++ libraries. Boris Schäling 2011.

  79. Schroeder W, Martin KM, Lorensen WE. The Visualization Toolkit: An Object-Oriented Approach to 3D Graphics. USA: Prentice-Hall Inc; 2006.

    Google Scholar 

  80. Fedele M, Quarteroni A. Polygonal surface processing and mesh generation tools for the numerical simulation of the cardiac function. Int J Numer Methods Biomed Eng. 2021;37(4):3435.

    Article  Google Scholar 

  81. Antiga L, Steinman DA. The vascular modeling toolkit, 2008.

  82. Roney CH, Sim I, Yu J, Beach M, Mehta A, Alonso Solis-Lemus J, Kotadia I, Whitaker J, Corrado C, Razeghi O, et al. Predicting atrial fibrillation recurrence by combining population data and virtual cohorts of patient-specific left atrial models. Circ Arrhythmia Electrophysiol. 2022;15(2): 010253.

    Article  Google Scholar 

  83. Costa CM, Neic A, Kerfoot E, Gillette K, Porter B, Sieniewicz B, Gould J, Sidhu B, Chen Z, Elliott M, Mehta V, Plank G, Rinaldi C, Bishop M, Niederer S. A virtual cohort of twenty-four left-ventricular models of ischemic cardiomyopathy patients. King’s College London 2020.

  84. Niederer SA, Kerfoot E, Benson AP, Bernabeu MO, Bernus O, Bradley C, Cherry EM, Clayton R, Fenton FH, Garny A, et al. Verification of cardiac tissue electrophysiology simulators using an n-version benchmark. Philos Trans R Soc Math Phys Eng Sci. 2011;369(1954):4331–51.

    Google Scholar 

  85. Piersanti R, Regazzoni F, Salvador M, Corno AF, Dede’ L, Vergara C, Quarteroni A. 3D–0D closed-loop model for the simulation of cardiac biventricular electromechanics. Comput Methods Appl Mech Eng. 2022;391: 114607.

    Article  Google Scholar 

  86. Cooper J, Spiteri RJ, Mirams GR. Cellular cardiac electrophysiology modeling with chaste and cellml. Front Physiol. 2015;5:511.

    Article  PubMed  PubMed Central  Google Scholar 

  87. Zampini S. Inexact bddc methods for the cardiac bidomain model. In: domain decomposition methods in science and engineering XXI, pp. 247–255. Springer. 2014.

  88. Del Corso G, Verzicco R, Viola F. A fast computational model for the electrophysiology of the whole human heart. J Comput Phys. 2022;457: 111084.

    Article  Google Scholar 

  89. Chamakuri N, Kügler P. Parallel space-time adaptive numerical simulation of 3d cardiac electrophysiology. Appl Numer Math. 2022;173:295–307.

    Article  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No 740132, iHEART - An Integrated Heart Model for the simulation of the cardiac function, P.I. Prof. A. Quarteroni). P.C.A, R.P., F.R., M.B., M.F., S.P. and L.D. are members of the INdAM research group GNCS.

Author information

Authors and Affiliations

Authors

Contributions

PCA: conceptualization, methodology, software (development and maintenance), formal analysis, supervision, writing (original draft). RP: conceptualization, methodology, software (development and simulation), formal analysis, writing (original draft). FR: conceptualization, methodology, software (development and simulation), formal analysis, writing (original draft). MB: conceptualization, methodology, software (development, simulation, testing and maintenance), formal analysis, writing (original draft). MS: conceptualization, methodology, software (development and simulation), formal analysis, writing (original draft). MF: conceptualization, methodology, software (development), formal analysis, writing (review and editing). SP: conceptualization, methodology, software (development and simulation), formal analysis, writing (original draft). LD: project administration, writing (review and editing), supervision. AQ: funding acquisition, project administration, writing (review and editing), supervision.

Corresponding author

Correspondence to Roberto Piersanti.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Conflict of interest

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Africa, P.C., Piersanti, R., Regazzoni, F. et al. lifex-ep: a robust and efficient software for cardiac electrophysiology simulations. BMC Bioinformatics 24, 389 (2023). https://doi.org/10.1186/s12859-023-05513-8

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12859-023-05513-8

Keywords

Mathematics Subject Classification