POD-Enhanced Deep Learning-Based Reduced Order Models for the Real-Time Simulation of Cardiac Electrophysiology in the Left Atrium

The numerical simulation of multiple scenarios easily becomes computationally prohibitive for cardiac electrophysiology (EP) problems if relying on usual high-fidelity, full order models (FOMs). Likewise, the use of traditional reduced order models (ROMs) for parametrized PDEs to speed up the solution of the aforementioned problems can be problematic. This is primarily due to the strong variability characterizing the solution set and to the nonlinear nature of the input-output maps that we intend to reconstruct numerically. To enhance ROM efficiency, we proposed a new generation of non-intrusive, nonlinear ROMs, based on deep learning (DL) algorithms, such as convolutional, feedforward, and autoencoder neural networks. In the proposed DL-ROM, both the nonlinear solution manifold and the nonlinear reduced dynamics used to model the system evolution on that manifold can be learnt in a non-intrusive way thanks to DL algorithms trained on a set of FOM snapshots. DL-ROMs were shown to be able to accurately capture complex front propagation processes, both in physiological and pathological cardiac EP, very rapidly once neural networks were trained, however, at the expense of huge training costs. In this study, we show that performing a prior dimensionality reduction on FOM snapshots through randomized proper orthogonal decomposition (POD) enables to speed up training times and to decrease networks complexity. Accuracy and efficiency of this strategy, which we refer to as POD-DL-ROM, are assessed in the context of cardiac EP on an idealized left atrium (LA) geometry and considering snapshots arising from a NURBS (non-uniform rational B-splines)-based isogeometric analysis (IGA) discretization. Once the ROMs have been trained, POD-DL-ROMs can efficiently solve both physiological and pathological cardiac EP problems, for any new scenario, in real-time, even in extremely challenging contexts such as those featuring circuit re-entries, that are among the factors triggering cardiac arrhythmias.

1. Introduction

Computational cardiac electrophysiology (EP) is built upon mathematical and numerical models that aim at simulating both physiological and pathological heart rhythm, such as, e.g., ventricular tachycardia and atrial fibrillation (see, e.g., Vigmond et al., 2002, 2008; Niederer et al., 2009, 2011; Trayanova, 2011; Prakosa et al., 2018; Strocchi et al., 2020). Simulating the electrical behavior of the heart, from the cellular scale to the tissue level, requires the numerical approximation of coupled nonlinear dynamical systems, such as, e.g. the Bidomain equations (see, e.g., Colli Franzone et al., 2005, 2006), coupled with suitable ionic models, such as the FitzHugh-Nagumo (FitzHugh, 1961; Nagumo et al., 1962), the Aliev-Panfilov (Aliev and Panfilov, 1996; Nash and Panfilov, 2004), the Roger-McCulloch (Rogers and McCulloch, 1994), the ten Tusscher-Panfilov (ten Tusscher and Panfilov, 2006), or the Mitchell and Schaeffer models (Mitchell and Schaeffer, 2003). Multiple solutions of these systems, corresponding to different model inputs parameters and data, such as, e.g., electrical conductivities, ionic model parameters, and applied currents, need to be computed to evaluate outputs of clinical interest, such as activation maps (ACs) and action potential (AP) duration. All these instances can be cast either in multi-query or real-time contexts. In the former case, the input-output map is repetitively evaluated in order to perform multi-scenario analysis, to deal with uncertainties and with inter- and intra-subject variability and to consider specific pathological scenarios; in the latter, outputs of interest must be computed in a very limited amount of time, in view of a possible integration in the clinical setting. Performing the numerical approximation of cardiac EP problems in multi-query context or solving them in real-time is in general out of reach for high-fidelity techniques or full order models (FOMs), such as the finite element (FE) method (Quarteroni and Valli, 1994) or isogeometric analysis (IGA) (Cottrell et al., 2009). To enhance their computational efficiency, multi-query and real-time problems may benefit from suitable surrogate models that can be built according to different strategies (see, e.g., Niederer et al., 2020 for a recent review). In particular, reduced order modeling techniques, can potentially provide more accurate approximations than data fitting techniques such as, e.g., data-driven emulators built through polynomial chaos expansions or gaussian processes. Moreover, they yield more significant computational savings than low-fidelity models (such as, e.g., FOMs built on coarser meshes) by replacing the FOM by a reduced order model (ROM), featuring a much lower dimension, yet capable to express the physical features of the problem at hand.

Cardiac EP problems are extremely challenging for traditional ROMs. Indeed, the latters tend to be inaccurate and/or computationally inefficient. This is primarily due to the high variability characterizing the solution manifold (with respect to the problem parameters), as well as to the nonlinear nature of the input-output maps that are more frequently approximated. Indeed, cardiac EP models feature coherent structures that propagate over time. In particular, as soon as re-entries, the most recognized cellular mechanisms sustaining atrial tachycardia and atrial fibrillation (Nattel, 2002) are considered, and wavefronts show abnormal activation patterns. These systems can hardly be reduced to lower dimensional problems by traditional ROMs for parametrized problems such as, e.g., the reduced basis (RB) method (Quarteroni et al., 2016). The most advanced example of efficient and accurate ROM in cardiac EP can be found in Pagani et al. (2018), where a local POD-Galerkin ROM has been proposed to handle physiological cardiac EP described in terms of the simpler Monodomain equation. However, to the best of our knowledge, no attempt to construct a comprehensive and systematic ROM framework to efficiently deal with parameter-dependent Bidomain equations involving pathological scenarios, such as re-entries, has been made yet.

Recently, we have introduced a new class of non-intrusive—since just a collection of FOM snapshots is required—nonlinear ROM techniques based on deep learning (DL) algorithms, named DL-ROMs, for the construction of efficient ROMs for parameter-dependent PDEs; in particular, we have focused so far on the Monodomain equations for cardiac EP, both in physiological and pathological scenarios (Fresca et al., 2020), as well as on several other nonlinear time-dependent parametrized problems, see (Fresca et al., 2021). DL-ROMs proved to be computationally efficient during the testing stage, that is, for any new scenario unseen during the training stage, but they might imply overwhelming training costs (and times) when the FOM dimension becomes moderately large. POD-enhanced DL-ROMs, first introduced and analyzed in Fresca et al. (2020), also enable fast training stages, improving on the weakest aspect—however, taking advantage of the key properties—of DL-ROMs.

So far, limited attempts have been made to solve, by means of DL algorithms, problems featuring traveling waves or front propagation processes. For example, recurrent and convolutional deep neural networks have been employed to predict the propagation of surface waves in Fotiadis et al. (2020). Regarding the cardiac EP context, in Cantwell et al. (2019) machine learning (ML) techniques have been considered for time prediction or parameters estimation, and EP-nets have been proposed in Ayed et al. (2019) and Kashtanova et al. (2021) to replace numerical integration of PDEs. On the other hand, deep neural networks have been extensively exploited to address several issues in computational fluid dynamics (see, e.g., Kutz, 2017; Bhatnagar et al., 2019; Ströfer et al., 2019; Brunton et al., 2020; Thuerey et al., 2020; Eichinger et al., 2021; Fresca and Manzoni, 2021b).

In this study, we show that POD-DL-ROMs can handle parametrized problems in cardiac EP effectively and provide fast and accurate solutions to EP problems set on realistic geometries. In particular, the performance of POD-DL-ROMs is assessed on cardiac EP on a left atrium (LA) surface geometry, both in physiological and pathological scenarios. These problems are challenging for traditional ROMs, due to (i) the presence of steep wavefronts, (ii) the complex activation patterns associated with pathological scenarios, (iii) the high FOM dimension, and (iv) the geometrical complexity. POD-DL-ROMs yield accurate and extremely efficient numerical approximations, irrespectively of the concurrence of these challenging features. This is particularly useful in view of the evaluation of patient-specific features to enable the integration of computational methods in current clinical practice; indeed, outputs of clinical interest, such as ACs, APs, electrograms, and ablation targets, can be more efficiently evaluated by the POD-DL-ROMs than by a FOM, while maintaining a high level of accuracy. The numerical tests carried out in this study represent a proof-of-concept of the POD DL-ROM technique ability to investigate intra- and inter-subjects variability toward performing multi-scenario analyses in real-time and, ultimately, supporting decisions in clinical practice.

To build our ROMs, we rely on a FOM obtained by means of an IGA spatial discretization. This choice is motivated by the suitability of high order polynomials, with high order global continuity, to control and limit numerical dispersions and, thus, to accurately capture wavefronts (Dedè et al., 2015; Pegolotti et al., 2019) and the smoothness in the representation of the computational domain (Cottrell et al., 2009). These relevant features have been exploited to address cardiac EP problems in Patelli et al. (2017), Pegolotti et al. (2019), and Bucelli et al. (2021). It is also worthy to highlight that, so far, only few works provide a combination of IGA-based FOMs and reduced order modeling techniques. IGA POD-Galerkin ROMs have been first applied to potential flows (Manzoni et al., 2015) and shell structural problems (Rinaldi, 2015), then to linear parabolic PDEs (Zhu et al., 2017a) and time-dependent parameterized acoustic wave equations (Zhu et al., 2017b); see also (Salmoiraghi et al., 2016; Garotta et al., 2020).

The structure of this study is as follows. In section 2, we introduce the FOM used to approximate the problem at hand and the POD-DL-ROM technique. The numerical assessment of this latter is carried out in section 3 on three different test cases; a more in-depth discussion is reported in section 4.

2. Materials and Methods

This section provides an overview of the mathematical and numerical models describing cardiac EP, including the reduced order modeling technique we employ to achieve computational efficiency in the solution of the Bidomain equations.

2.1. Mathematical Models for Cardiac Electrophysiology

The electrical activation of the heart, which drives its contraction, is the result of two processes (Klabunde, 2011; Colli Franzone et al., 2014): the generation of ionic currents through the cellular membrane producing a local AP, at the microscopic scale, and the propagation of the AP from cell to cell in the form of a transmembrane potential, at the macroscopic scale. The latter process can be described by means of PDEs, suitably coupled with systems of ODEs accounting for the former (Quarteroni et al., 2017, 2019). To model the propagation of the electrical signal in the heart, we may consider the so-called Bidomain equations (Geselowitz and Miller III, 1983; Colli Franzone et al., 2014) in a domain Ω ⊂ ℝ^d, with d = 2, 3, representing a portion of the myocardium, considered as a continuum composed of two interpenetrating domains, the intracellular and the extracellular spaces. Each point x ∈ Ω is associated with the intracellular potential u_i, the extracellular potential u_e, and the transmembrane potential u = u_i−u_e. Coupling the parabolic-elliptic formulation of the Bidomain model for the transmembrane potential u = u(x, t) and the extracellular potential u_e = u_e(x, t) with a phenomenological¹ model for the ionic currents—involving a single gating variable w = w(x, t)—results in the following nonlinear time-dependent system:

$\begin{array}{ll}\{\begin{array}{ll}\frac{\partial u}{\partial t}-\text{div}(D_i\nabla u)-\text{div}(D_i\nabla u_e)+I_{ion}(u,w)=I_{app}^i\hfill & (x,t)\in \Omega \times (0,T),\hfill \\ -\text{div}(D_i\nabla u)-\text{div}((D_i+D_e)\nabla u_e)=I_{app}^i+I_{app}^e\hfill & (x,t)\in \Omega \times (0,T),\hfill \\ \frac{\partial w}{\partial t}+g(u,w)=0\hfill & (x,t)\in \Omega \times (0,T),\hfill \\ D_i\nabla (u+u_e)·n=0\hfill & (x,t)\in \partial \Omega \times (0,T),\hfill \\ (D_i+D_e)\nabla u_e·n+D_i\nabla u·n=0\hfill & (x,t)\in \partial \Omega \times (0,T),\hfill \\ u(x,0)=u_0,w(x,0)=w_0\hfill & x\in \Omega .\hfill \end{array}& (1)\end{array}$

Here, t and u denote a rescaled and dimensionless time and trasmembrane potential, depending on the ionic model considered², n denotes the outward directed unit vector normal to the boundary ∂Ω of Ω, whereas $I_{app}^i=I_{app}^i(x,t)$ and $I_{app}^e=I_{app}^e(x,t)$ are the intracellular and the extracellular applied currents representing, e.g., the initial activation of the tissue. The parabolic nonlinear diffusion-reaction equation for u is two-way coupled with the ODE system; this latter must be solved, in principle, at any point x ∈ Ω. Indeed, both I_ion and g depend on u and w, and the most common choices to efficiently reproduce the AP are, e.g., the FitzHugh-Nagumo (FitzHugh, 1961; Nagumo et al., 1962), the Aliev-Panfilov (Aliev and Panfilov, 1996; Nash and Panfilov, 2004), the Roger-McCulloch (Rogers and McCulloch, 1994), or the Mitchell-Schaeffer models (Mitchell and Schaeffer, 2003). The diffusivity tensors D_i, D_e usually depend on the fibers-sheet structure of the tissue, affecting directional conduction velocities and direction. In particular, by assuming an axisymmetric distribution of the fibers, the intracellular and extracellular conductivity tensors take the form

$\begin{array}{l}{D}_i(x)={\sigma }_t^{i}I+({\sigma }_l^i-{\sigma }_t^i)f_0(x)\otimes f_0(x),\\ D_e(x)={\sigma }_t^{e}I+({\sigma }_l^e-{\sigma }_t^e)f_0(x)\otimes f_0(x),& (2)\end{array}$

where ${\sigma }_l^i,{\sigma }_l^e$ and ${\sigma }_t^i,{\sigma }_t^e$ are the electrical conductivities in the fibers and the transversal directions, for the intracellular and extracellular conductivity tensors. A simplified model is given by the Monodomain equation (Colli Franzone et al., 2014), written only in terms of the transmembrane potential u.

For most of the basic phenomenological ionic models, such as the FitzHugh-Nagumo, the Aliev-Panfilov (A-P) (Aliev and Panfilov, 1996) or the Roger-McCulloch (R-M) (Rogers and McCulloch, 1994) model, the ionic current takes the form of a cubic nonlinear function of u and a single (dimensionless) gating variable plays the role of a recovery function, allowing to model cell refractoriness. In this study, we focus on the simple phenomenological A-P and R-M ionic models in order to lessen the computational costs associated with the approximation of Equation (1) through a FOM. The A-P model consists in taking

$\begin{array}{l}{I}_{ion}(u,w)=Ku(u-a)(u-1)+uw,\\ g(u,w)=({ϵ}_0+\frac{c_{1}w}{c_2+u})(-w-Ku(u-b-1)),& (3)\end{array}$

where the parameters K, a, b, ε₀, c₁, c₂ are related to the cell. Here, a represents an oscillation threshold, the weighting factor ${\epsilon }_0+\frac{c_{1}w}{c_2+u}$ was introduced in Aliev and Panfilov (1996) to tune the restitution curve to experimental observations by adjusting the parameters c₁ and c₂, whereas K and b are coefficients set according to Aliev and Panfilov (1996); see, e.g., (Clayton et al., 2011; Colli Franzone et al., 2014) for a detailed review. For the R-M ionic model, we rely on the following variant provided in Rogers and McCulloch (1994)

$\begin{array}{l}{I}_{ion}(u,w)=Gu(1-\frac{u}{u_{th}})(1-\frac{u}{u_p})+{\eta }_{1}uw,\\ g(u,w)={\eta }_2(\frac{u}{u_p}-{\eta }_{3}w),& (4)\end{array}$

where G, η₁, η₂, η₃ are positive coefficients, v_th is a threshold potential, and v_p is the peak potential.

The coupled system (Equation 1) depends on several parameters representing either functional or geometric data such as, e.g., material properties, initial and boundary conditions, or the shape of the domain. In the remaining part of the study, we denote by $\mu \in {P}\subset {ℝ}^{n_{\mu }}$ a parameter vector listing all the n_μ input parameters characterizing physical (and, possibly, geometrical) properties; ${P}$ is a subset of ${ℝ}^{n_{\mu }}$, denoting the parameter space. Relevant physical situations are those in which input parameters affect the diffusivity matrix D (through the conduction velocities) and the applied current I_app; for previous analyses focused instead on the gating variable dynamics (through g) and the ionic current I_ion in the case of the Monodomain equation (see, e.g., Pagani et al., 2018).

Regarding the spatial discretization of the system (Equation 1), we consider NURBS-based IGA on surfaces (e.g., the LA), in the framework of Galerkin methods (Quarteroni, 2017). Here, the same NURBS basis functions are employed both to define the computational domain and to construct the finite-dimensional space in which the numerical solution of the PDE is sought (Cottrell et al., 2009). Globally high order continuous polynomials have proved to control and limit numerical dispersion (Dedè et al., 2015), which may lead to artificial fractionated potential fronts, when dealing with the sharp but smooth fronts arising in cardiac EP. To correctly describe cardiac EP, capturing propagating fronts and their velocity is essential. The use of NURBS basis functions with high polynomial degree (say, p) and global high order continuity (say, ${{C}}^{p-1}$) is beneficial, in terms of both accuracy and efficiency, to deal with Monodomain/Bidomain equations since they limit dispersion effects typical of traveling wave phenomena (Patelli et al., 2017; Pegolotti et al., 2019). Moreover, NURBS basis functions also allow a smooth representation of the computational domain starting from medical images, compared to methods exploiting polyhedral elements, as it usually happens when dealing with finite element approximations (Cottrell et al., 2009). In particular, we employ a two-dimensional NURBS surface of the LA built starting from B-spline basis functions of degree p = 2. For further details on the construction of the LA computational domain, we refer to Patelli et al. (2017). The smoothness of the computational domain, together with the regularity of NURBS basis functions, makes IGA well-suited for surface problems requiring high order polynomials.

2.2. Proper Orthogonal Decomposition-Enhanced Deep Learning-Based Reduced Order Models (POD-DL-ROMs)

From an algebraic standpoint, the spatial discretization of the system (Equation 1) through a NURBS-based IGA approximation yields the following nonlinear dynamical system for u_h = u_h(t, μ), u_e,h = u_e,h(t, μ) and w_h = w_h(t, μ), representing our FOM:

$\begin{array}{ll}\{\begin{array}{ll}M(\mu )\frac{\partial u_h}{\partial t}+A_i(\mu )u_h+A_i(\mu )u_{e,h}=I_{app}^i(t;\mu )\hfill & t\in (0,T),\hfill \\ +I_{ion}(t,u_h,w_h;\mu )\hfill & \hfill \\ A_i(\mu )u_h+A(\mu )u_{e,h}=I_{app}^i(t;\mu )+I_{app}^e(t;\mu )\hfill & t\in (0,T),\hfill \\ \frac{\partial w_h}{\partial t}=g(u_h,w_h;\mu )\hfill & t\in (0,T),\hfill \\ u_h(0)=u_0(\mu ),w_h(0)=w_0(\mu ),\hfill & \hfill \end{array}& (5)\end{array}$

where u_h, u_e,h, and $w_h\in {ℝ}^{N_h}$, being the dimension N_h related to the dimension of the NURBS space, and $\mu \in {P}\subset {ℝ}^{n_{\mu }}$. In the remaining part of this study, we consider as initial data u₀(μ) = 0 and w₀(μ) = 0. A detailed derivation of the FOM (Equation 5) is reported in the Supplementary Material.

Solving (Equation 5) is computationally demanding and far beyond the possibility to provide solutions or compute outputs of interest in real-time applications. Indeed, the propagation of the electrical signal is characterized by the fast dynamics of very steep fronts, thus requiring very fine space and time discretizations (Colli Franzone and Pavarino, 2004; Sundnes et al., 2006). This is even more true if such a coupled system must be solved for several parameters instances, that is, in a multi-query context, in order to investigate different scenarios or intra- and inter-subject variability. ROM techniques replace the FOM (Equation 5) by a model featuring a much lower complexity but still able to retain the physical features of the problem at hand. Traditional projection-based ROMs built, e.g., through the RB method (Quarteroni et al., 2016), yields inefficient ROMs when dealing with nonlinear time-dependent parametrized PDE-ODE system as the one arising from cardiac EP (Fresca et al., 2020). To overcome the limitation of traditional projection-based ROMs, we have recently proposed in Fresca et al. (2021) a strategy to construct, in a non-intrusive/data-driven way (indeed neither access or solution to the governing equations are required), DL-based ROMs (DL-ROMs) for nonlinear time-dependent parametrized problems, exploiting deep neural networks (Goodfellow et al., 2016) as a main building block, and a set of FOM snapshots. A first attempt to solve, by means of DL-ROMs, parametrized benchmark test cases in cardiac EP described by the Monodomain equations, has been carried out in Fresca et al. (2020). Although extremely efficient at testing (i.e., online) time, when evaluating the problem solution for any new testing-parameter instance, DL-ROMs require an expensive training (i.e., offline) stage, because of the extremely large number of network parameters to be estimated. POD-DL-ROMs provide a possible enhancement of DL-ROMs, which avoids expensive training stages, by (i) performing a prior dimensionality reduction through proper orthogonal decomposition (POD), and (ii) using a multi-fidelity pretraining stage, where different physical models can be efficiently combined, as recently shown in Fresca and Manzoni (2021a). In particular, through the use of randomized POD, the POD-DL-ROM training phase is extremely fast, especially if compared to the training stage of DL-ROMs. For example, in Fresca and Manzoni (2021a), where we consider the solution of the parametrized Monodomain equation in a square slab of cardiac tissue on a FOM dimension N_h = 4096, the use of the POD-enhanced DL-ROM reduces the GPU training time from 15 h to 24 min, while preserving extremely efficient testing times.

Tailored on the applications at hand, the goal of POD-DL-ROMs is to approximate the map (t, μ) ↦ u_h(t, μ), where t ∈ (0, T) denotes time, $\mu \in {P}\subset {ℝ}^{n_{\mu }}$ a vector of input parameters, and $u_h(t,\mu )\in {ℝ}^{N_h}$ the trasmembrane potential solution of Equation (5). This may be achieved without taking into account, and then expensively solving, the dynamics of the extracellular potential u_e,h(t, μ) and the gating variable w_h(t, μ) in the construction of the ROM. More precisely, we build a nonlinear ROM to approximate $V^T{u}_h(t;\mu )\approx {\tilde{u}}_N(t;\mu )$ by

$\begin{array}{ll}{\tilde{u}}_N(t;\mu )={\text{Ψ}}_N(u_n(t;\mu )),& (6)\end{array}$

where ${\text{Ψ}}_N:{ℝ}^n\to {ℝ}^N$, Ψ_N:s_n ↦ Ψ_N(s_n), n ≪ N, is a nonlinear, differentiable function and $V\in {ℝ}^{N_h\times N}$ is the rPOD basis matrix of a N-dimensional subspace of ${ℝ}^{N_h}$. In particular, the columns of V form an orthonormal basis of dimension N, computed by means of randomized SVD (rSVD) (Halko et al., 2011). In this way, the manifold ${{S}}_N=\left\{V^T{u}_h(t;\mu )|t\in [0,T)\text{and}\mu \in {P}\subset {ℝ}^{n_{\mu }}\right\}\subset {ℝ}^N$ is approximated by the n-dimensional reduced nonlinear trial manifold

$\begin{array}{l}{\tilde{S}}_n=\{{\Psi }_N(u_n(t;\mu ))|u_n(t;\mu )\in {ℝ}^n,\\ t\in [0,T)\text{ and }\mu \in \subset {ℝ}^{n_{\mu }}\}\subset {ℝ}^N,& (7)\end{array}$

where ${\tilde{u}}_N:[0,T)\times {P}\to {\tilde{{S}}}_n$. The function $u_n:[0,T)\times {P}\to {ℝ}^n$ denotes the minimal coordinates of ${\tilde{u}}_N$ on the nonlinear trial manifold ${\tilde{{S}}}_n$. Our goal is to set-up a ROM whose dimension n is as close as possible to the intrinsic dimension n_μ + 1 (time plays the role of an additional coordinate) of the solution manifold ${{S}}_h$, i.e. n ≥ n_μ + 1, to correctly capture the degrees of freedom of the set ${{S}}_N$ by containing its size (Lee and Carlberg, 2020). To model the relationship between each pair (t, μ) ↦ u_n(t, μ), and to describe the reduced dynamics on the reduced nonlinear trial manifold ${\tilde{{S}}}_n$, we consider a nonlinear map under the form

$\begin{array}{ll}{u}_n(t;\mu )={\Phi }_n(t,\mu ),& (8)\end{array}$

where ${\Phi }_n:[0,T)\times {ℝ}^{n_{\mu }}\to {ℝ}^n$ is a differentiable, nonlinear function. As for DL-ROMs (see e.g., Fresca et al., 2021), both the reduced dynamics and the reduced nonlinear manifold where the ROM solution is sought (or trial manifold) must be learnt. In particular,

• reduced dynamics learning: We aim at learning the dynamics of the set ${{S}}_N$ on the nonlinear trial manifold ${\tilde{{S}}}_n$ in terms of minimal coordinates, by means of a deep feedforward neural network (DFNN). Indeed, we set the function Φ_n in Equation (8) equal to

$\begin{array}{l}{\Phi }_n(t;\mu ,{\theta }_{DF})={\varphi }_n^{DF}(t;\mu ,{\theta }_{DF}),\end{array}$

where θ_DF denotes the vector of parameters of the DFNN, collecting all the corresponding weights and biases of each layer of the DFNN;

• nonlinear trial manifold learning: We employ the decoder function of a convolutional autoencoder (AE), that is, we define the function in Equation (6) as

$\begin{array}{l}{\text{Ψ}}_N(u_n(t;\mu ,{\theta }_{DF});{\theta }_D)=f_N^D(u_n(t;\mu ,{\theta }_{DF});{\theta }_D),\end{array}$

where $f_N^D$ depends on the vector θ_D of parameters of the convolutional/dense layers of the decoder.

By combining the two previous stages, the POD-DL-ROM approximation ${\tilde{u}}_N$ finally takes the form

$\begin{array}{ll}{\tilde{u}}_N(t;\mu ,{\theta }_{DF},{\theta }_D)=f_N^D({\varphi }_n^{DF}(t;\mu ,{\theta }_{DF});{\theta }_D).& (9)\end{array}$

The encoder function of the convolutional AE can then be exploited to map the intrinsic coordinates $V^T{u}_h$ associated to (t, μ) onto a low-dimensional representation

$\begin{array}{l}{\tilde{u}}_n(t;\mu ,{\theta }_E)=f_n^E(V^T{u}_h(t;\mu );{\theta }_E),\end{array}$

where $f_n^E$ denotes the encoder function, depending upon a vector θ_E of parameters. The architecture of the POD-DL-ROM neural network, employed at training time, is the one shown in Figure 1. At testing time we can discard the encoder function.

FIGURE 1

Figure 1. Starting from the FOM solution u_h(t; μ), the intrinsic coordinates $V^T{u}_h(t;\mu )$ are computed through rSVD; their approximation ${\tilde{u}}_N(t;\mu )$ is provided by the neural network as output, so that the reconstructed solution ${\tilde{u}}_h(t;\mu )$ is recovered through the rPOD basis matrix. In particular, the intrinsic coordinates $V^T{u}_h(t;\mu )$ are provided as input to block (A), which returns as output ${\tilde{u}}_n(t;\mu )$. The same parameter instance (t; μ) enters block (B), which provides as output u_n(t; μ), and the error between the low-dimensional vectors is accumulated. The minimal coordinates u_n(t; μ) are given as input to block (C), which returns the approximated intrinsic coordinates ${\tilde{u}}_N(t;\mu )$. Then, the reconstruction error is computed.

Computing the POD-DL-ROM approximation (Equation 9) thus consists of solving the optimization problem

$\begin{array}{ll}{\displaystyle \underset{\theta }{min}}{J}(\theta )={\displaystyle \underset{\theta }{min}}\frac{1}{N_s}{\displaystyle \sum _{i=1}^{N_{train}}}{\displaystyle \sum _{k=1}^{N_t}}{L}(t^k,{\mu }_i;\theta ),& (10)\end{array}$

where the per-example loss function is given by

$\begin{array}{l}{L}(t^k,{\mu }_i;\theta )=\frac{{\omega }_h}{2}\left|\right|V^T{u}_h(t^k;{\mu }_i)-{\tilde{u}}_N(t^k;{\mu }_i,{\theta }_{DF},{\theta }_D)|{|}^2\\ +\frac{1-{\omega }_h}{2}\left|\right|{\tilde{u}}_n(t^k;{\mu }_i,{\theta }_E)-u_n(t^k;{\mu }_i,{\theta }_{DF})|{|}^2,& (11)\end{array}$

N_train and N_test are the number of training- and testing-parameter instances, respectively, N_t is the number of time instances, N_s = N_trainN_t, and ω_h ∈ [0, 1]. The POD-DL-ROM approximation of the FOM solution ${\tilde{u}}_h(t;\mu )\approx u_h(t;\mu )$ is then recovered by means of the rPOD basis matrix as.

$\begin{array}{l}{\tilde{u}}_h(t;\mu ,{\theta }_{DF},{\theta }_D)=V{\tilde{u}}_N(t;\mu ,{\theta }_{DF},{\theta }_D).\end{array}$

3. Results

In this section, we apply the POD-DL-ROM technique to relevant problems in cardiac EP, both in physiological and pathological scenarios, solved on a rectangular slab and a left atrium surface geometry. Dealing with realistic geometries, large-scale problems, i.e., high FOM dimensions N_h, and pathological scenarios, such as re-entries, show the feasibility of POD-DL-ROM to be integrated in to the clinical practice in order to compute outputs of interest, e.g., ACs, action potential durations, electrograms, and location of cores of rotors. To evaluate the performance of POD-DL-ROM, we rely on the loss function (Equation 11) and on:

• the error indicator ϵ_rel ∈ ℝ given by

$\begin{array}{l}{ϵ}_{rel}={ϵ}_{rel}(u_h,{\tilde{u}}_h)\\ =\frac{1}{N_{test}}{\displaystyle \sum _{i=1}^{N_{test}}}\left(\frac{\sqrt{{\displaystyle {\sum }_{k=1}^{N_t}}||u_h^k({\mu }_{test,i})-{\tilde{u}}_h^k({\mu }_{test,i})|{|}^2}}{\sqrt{{\displaystyle {\sum }_{k=1}^{N_t}}||u_h^k({\mu }_{test,i})|{|}^2}}\right);& (12)\end{array}$

• the relative error ${ϵ}_k\in {ℝ}^{N_h}$, for k = 1, …, N_t, defined as

$\begin{array}{ll}{ϵ}_k={ϵ}_k(u_h,{\tilde{u}}_h)=\frac{|u_h^k({\mu }_{test})-{\tilde{u}}_h^k({\mu }_{test})|}{\sqrt{\frac{1}{N_t}{\displaystyle {\sum }_{k=1}^{N_t}}||u_h^k({\mu }_{test})|{|}^2}}.& (13)\end{array}$

While Equation (12) is a synthetic indicator, the quantity defined in Equation (13) is instead a spatially distributed function.

The configuration of the POD DL-ROM neural network, together with the values of the hyperparameters not reported in this study, used for our numerical tests is the same provided as in Fresca and Manzoni (2021a). The FOM simulations are carried out on a MacBook Pro Intel Core i7 6-core with 16 GB RAM CPU, while the POD-DL-ROM training and testing phases on a Tesla V100 32GB GPU.

3.1. Test 1: Slab and Left Atrium Surface Geometry by Varying Conductivities

3.1.1. Test 1.1: Slab Geometry

We consider the Bidomain Equation (1) coupled with the R-M ionic model (Equation 4) in a two-dimensional rectangular slab of cardiac tissue Ω = (0, 10)cm × (0, 2)cm. In order to characterize the bidomain nature of the tissue, we focus on the reconstruction of both the transmembrane and the extracellular potentials. To achieve this task, the intrinsic coordinates of the two field variables, i.e. v_h(t) and u_e,h(t), are stacked together, thus forming a tensor with d = 2 input channels, which represent the actual input (output) of the POD-DL-ROM neural network. The parameter (n_μ = 1) consists of the electrical extracellular conductivity in the longitudinal direction to the fibers, i.e., the conductivity tensor D_e(x; μ) takes the form

$\begin{array}{l}{D}_e(x;\mu )={\sigma }_t^{e}I+(\mu -{\sigma }_t^e)f_0\otimes f_0,\end{array}$

where $f_0={(1,0)}^T$ and the parameter space is ${P}=1.5·\left[10^{-3},10^{-2}\right]{\Omega }^{-1}{\text{cm}}^{-1}$. The remaining intracellular and extracellular conductivities are set equal to ${\sigma }_l^i=2.3\times 10^{-3}{\Omega }^{-1}{\text{cm}}^{-1}$, ${\sigma }_t^i=2.4\times 10^{-4}{\Omega }^{-1}{\text{cm}}^{-1}$, and ${\sigma }_t^e=1\times 10^{-3}{\Omega }^{-1}{\text{cm}}^{-1}$, respectively. The parameters of the R-M ionic model are given by u_th = 13 mV, v_p = 100 mV, G = 1.5 ms⁻¹, η₁ = 4.4 ms⁻¹, ${\eta }_2=1.2\times 10^{-2}$, and η₃ = 1 (see, e.g., Gerardo-Giorda, 2007). We provide snapshots computed by means of P3/C2 NURBS basis functions, where N_h = 165 × 35 = 5705, with n_el = 5120 mesh elements. Time integration is performed over the interval (0, T), with T = 150 ms and a time-step Δt = 0.05 ms, through a BDF of order 2. The intracellular applied current takes the form

$\begin{array}{ll}{I}_{app}^i(x,t)=C{1}_{{\Omega }_{app}}(x)1_{[t^i,t^f]}(t),& (14)\end{array}$

where C = 100 mA, Ω_app = {x ∈ Ω:x ≤ 0.2}, tⁱ = 0 ms, and t^f = 1 ms.

For the training phase, we consider N_t = 1500 time instances in the interval (0, T) and N_train = 11 training-parameter instances uniformly distributed in the parameter space. For the testing phase, N_test = 10 testing-parameter instances have been considered, each of them corresponding to the midpoint of two consecutive training-parameter instances. The maximum number of epochs is N_epochs = 20, 000, the batch size is N_b = 40, and regarding the early-stopping criterion, we stop the training if the loss function does not decrease along 1,000 epochs. In Figure 2, we report both the FOM and the POD-DL-ROM solution, the latter with n = 2 and N = 256, together with the resulting relative error, both for the transmembrane and the extracellular potentials, for the testing-parameter instance ${\mu }_{test}=0.0143{\Omega }^{-1}{\text{cm}}^{-1}$ at t = 150 ms.

FIGURE 2

Figure 2. Test 1.1: FOM solution (top) and POD-DL-ROM one (center), with n = 2 and N = 256, along with the relative error ϵ_k for u_h (bottom-left) and u_e,h (bottom-right), for the testing-parameter instance ${\mu }_{test}=0.0143{\Omega }^{-1}{\text{cm}}^{-1}$ at t = 150 ms.

The trend of the mean (with respect to the spatial coordinates) of the relative error ϵ_k over time, for the selected testing-parameter instance ${\mu }_{test}=0.0143{\Omega }^{-1}{\text{cm}}^{-1}$, is shown in Figure 3, for both the trasmembrane (left) and the extracellular (right) potentials. We highlight that the errors are, on average, always smaller than 0.15%. The distribution of the errors is almost uniform over time; indeed, due to the fact that the snapshots associated with different time instances are treated as independent by the POD-DL-ROM, errors do not accumulate over time. In this manner, neither instability issues nor specific error patterns are found. In particular, the error related to the extracellular potential is higher than the one associated with the transmembrane potential. As a matter of fact, the former can be more difficult to approximate than the latter, because of the different ranges the extracellular potential can vary in for different parameters and time instances—and on the other hand, the transmembrane potential always takes values in the [0, 100] mV range.

FIGURE 3

Figure 3. Test 1.1: Relative error ϵ_k, averaged with respect to the spatial coordinates, for the transmembrane (A) and the extracellular (B) potentials, for the testing-parameter instance ${\mu }_{test}=0.0143{\Omega }^{-1}{\text{cm}}^{-1}$, over time.

In Table 1, we report the GPU POD-DL-ROM total training and validation times, together with the testing time, and the DL-ROM total time; the time needed to assemble the snapshot matrix S is not included. Using a POD-DL-ROM, that is, employing a prior dimensionality reduction through rSVD, drastically accelerates the training stage. We point out that in this test case, in contrast with the following ones, we did not perform any sampling in time, considering all the time instances provided by the IGA solver.

TABLE 1

Table 1. Test 1.1: POD-DL-ROM and DL-ROM computational times.

3.1.2. Test 1.2: Left Atrium Surface Geometry

We now consider the solution of the Bidomain (Equation 1) coupled with the A-P ionic model (Equation 3) on an idealized LA surface geometry. We are interested in the reconstruction of both the transmembrane and the extracellular potentials as in the previous test. The direction of the cardiac fibers is determined by following the same strategy adopted in Rossi et al. (2014) and Patelli et al. (2017), where a vector field directed as the gradient of the solution of a Laplace problem defined on the atrial surface is assigned to the LA. The resulting distribution of fibers on the atrial surface is displayed in the Supplementary Material.

System Equation (1) has been first discretized in space by means of P2 NURBS basis functions, with a global C1 continuity, yielding a FOM dimension equal to N_h = 61, 732. Time integration over (0, T), with T = 200 ms, has been performed introducing a time-step Δt = 0.2 ms. Provided the position of the Bachmann bundle $\bar{x}={(\bar{x},ȳ,\bar{z})}^T={(-1.51,0.1,-1.71)}^T$ cm, the intracellular applied stimulus is given by

$\begin{array}{l}{I}_{app}^i(x,t)=C{1}_{{\Omega }_{app}}(x)1_{[t^i,t^f]}(t),\end{array}$

with C = 1 mA, ${\Omega }_{app}=\left\{x\in \Omega :{(x-\bar{x})}^2+{(y-ȳ)}^2+{(z-\bar{z})}^2\le {(0.5)}^2\right\}$, tⁱ = 0 ms and t^f = 5 ms.

The parameter (n_μ = 1) consists of the electrical intracellular conductivity in the longitudinal direction to the fibers, i.e., the conductivity tensor D_i(x; μ) takes the form

$\begin{array}{l}{D}_i(x;\mu )={\sigma }_t^{i}I+(\mu -{\sigma }_t^i)f_0(x)\otimes f_0(x),\end{array}$

where the parameter space is ${P}=3.1·\left[10^{-4},10^{-2}\right]{\Omega }^{-1}{\text{cm}}^{-1}$. The remaining intracellular and extracellular conductivities are equal to ${\sigma }_t^i=2\times 10^{-2}{\Omega }^{-1}{\text{cm}}^{-1}$, ${\sigma }_l^e=1.3\times 10^{-4}{\Omega }^{-1}{\text{cm}}^{-1}$, and ${\sigma }_t^e=2\times 10^{-3}{\Omega }^{-1}{\text{cm}}^{-1}$. The parameters of the A-P ionic model (Equation 3) are given by K = 8, a = 0.1, ϵ₀ = 0.01, b = 0.1, c₂ = 0.3, and c₁ = 0.05 (ten Tusscher, 2004).

For the training phase, we uniformly sample N_t = 500 time instances in the interval (0, T) and consider N_train = 11 training-parameter instances, uniformly distributed over ${P}$. For the testing phase, N_test = 10 testing-parameter instances have been considered, each of them corresponding to the midpoint of two consecutive training-parameter instances. The maximum number of epochs is set to N_epochs = 20, 000, the batch size is N_b = 40, and regarding the early-stopping criterion, we stop the training if the loss function does not decrease along 1,000 epochs. In Figures 4, 5, we report the FOM transmembrane and extracellular potentials and their POD-DL-ROM approximation, obtained by selecting n = 2 and N = 256, for the testing-parameter instance ${\mu }_{test}=0.0295{\Omega }^{-1}{\text{cm}}^{-1}$ at t = 52.8 ms and t = 112 ms.

FIGURE 4

Figure 4. Test 1.2: FOM solution (left) and POD-DL-ROM one (right), with n = 2 and N = 256, for the testing-parameter instance ${\mu }_{test}=0.0295{\Omega }^{-1}{\text{cm}}^{-1}$ at t = 52.8 ms (A) and t = 112 ms (B).

FIGURE 5

Figure 5. Test 1.2: FOM solution (left) and POD-DL-ROM one (right), with n = 2 and N = 256, for the testing-parameter instance ${\mu }_{test}=0.0295{\Omega }^{-1}{\text{cm}}^{-1}$ at t = 52.8 ms (A) and t = 112 ms (B).

In Figure 6, we show the FOM and POD-DL-ROM APs and extracellular potentials evaluated at a point x^*, with n = 2 and N = 256, for the testing-parameter instance μ_test = 0.0295 Ω⁻¹cm⁻¹. Despite the POD-DL-ROM solution is affected by some tiny oscillations related to the truncated rPOD modes, it is able to capture the shape and the different phases of the electrical propagation over time.

FIGURE 6

Figure 6. Test 1.2: FOM and POD-DL-ROM, with n = 2 and N = 256, APs for the testing-parameter instance μ_test = 0.0295 Ω⁻¹cm⁻¹.

In Table 2, we report the CPU computational time needed to solve the FOM by means of NURBS-based IGA and the GPU POD-DL-ROM total training and validation times, together with the testing time. Also in this case, the time needed to assemble the snapshot matrix S is not included. We notice that, using a POD-DL-ROM, we achieve the possibility to solve the problem in several different scenarios, during the testing stage, in real-time, since the final time T = 0.2 s coincides with the computational time entailed by the evaluation of the POD-DL-ROM.

TABLE 2

Table 2. Test 1.2: FOM and POD-DL-ROM computational times.

3.2. Test 2: Left Atrium Surface Geometry by Varying Stimulation Site

We still consider the LA surface geometry and the direction of cardiac fibers as in test 1.2 and deal with the Bidomain (Equation 1) coupled with the R-M model (Equation 4), thus selecting a different ionic model than the one of the previous example. The equations have been discretized in space by means of P2 NURBS basis functions, with a global C1 continuity, yielding a FOM dimension equal to N_h = 154, 036; time integration has been performed over the interval (0, T), with T = 200 ms and a time-step Δt = 0.1 ms. Here, we consider n_μ = 3 parameters, consisting of the coordinates of the center of an intracellular applied current, and belonging to the subregion highlighted in Figure 7—and the portion of the domain affected by the corresponding stimulus is highlighted, too. The intracellular applied current is thus defined by setting C = 100 mA and

$\begin{array}{l}{I}_{app}^i(x,t)=C{1}_{{\Omega }_{app}(\mu )}(x)1_{[t^i,t^f]}(t),\end{array}$

with ${\Omega }_{app}(\mu )=\left\{x\in \Omega :{(x-{\mu }_1)}^2+{(y-{\mu }_2)}^2+{(z-{\mu }_3)}^2\le {(0.5)}^2\right\}$, tⁱ = 0 ms, and t^f = 5 ms.

FIGURE 7

Figure 7. Test 2: Parameter space (dark magenta region) and portion of domain affected by the stimulus (light magenta region).

We set the rPOD dimension equal to N = 256 and the dimension n of the POD-DL-ROM approximation equal to n = n_μ + 1 = 4. For the training phase, we uniformly sample N_t = 200 time instances in the interval (0, T) and consider N_train = 18 training-parameter instances randomly sampled from the parameter space. For the testing phase, N_test = 14 randomly sampled testing-parameter instances have been considered. The maximum number of epochs is N_epochs = 40, 000, the batch size is N_b = 40, the starting learning rate is η = 2 · 10⁻⁴, and regarding the early-stopping criterion, we stop the training if the loss function does not decrease along 2,000 epochs.

We remark that the POD-DL-ROM approximation to the FOM solution is also efficient in computing several outputs of interest. We compare, for instance, the ACs obtained through the FOM and by POD-DL-ROM. Given the transmembrane potential u = u(x, t; μ), the (unipolar) AC at a point x ∈ Ω is evaluated as the minimum time which the AP peak reaches the point x at,

$\begin{array}{l}AC(x;\mu )=arg{\displaystyle \underset{t\in (0,T)}{min}}\left(u(x,t;\mu )={\displaystyle \underset{t\in (0,T)}{max}}u(x,t;\mu )\right).\end{array}$

In Figure 8, we compare the FOM and the POD-DL-ROM outputs, together with the associated relative error ϵ_k, for the testing-parameter instances μ_test = (1.7168, −0.353198, −1.70097) cm and μ_test = (1.43862, −0.803806, −1.43678) cm. We highlight the strong variability of the solution over the parameter space, shown by the different shape of the contour lines in Figures 8A,B, and the ability of the POD-DL-ROM to capture it accurately.

FIGURE 8

Figure 8. Test 3: FOM (left) and POD-DL-ROM (center), with n = 4 and N = 256, ACs and relative error ϵ_k (right), for the testing-parameter instances μ_test = (1.7168, −0.353198, −1.70097) cm (A) and μ_test = (1.43862, −0.803806, −1.43678) cm (B).

Finally, in Table 3 we report the FOM CPU computational time and the POD-DL-ROM GPU training and testing times; the time needed to assemble the snapshot matrix S is not included. Solving the FOM, for a single testing-parameter instance, requires 10 h, with respect to the POD-DL-ROM total training and validation time, which is equal to 5 h. POD-DL-ROM also proves to be extremely efficient at testing time, since it provides, once again, accurate results in almost real-time.

TABLE 3

Table 3. Test 2: FOM and POD-DL-ROM computational times.

3.3. Test 3: Figure of Eight Re-entry on Left Atrium Surface Geometry

We finally investigate the generation of the figure of eight re-entries on the left atrium surface geometry as a consequence of a S1-S2 electrical stimulation protocol, to highlight the ability of the POD-DL-ROM the technique of solving cardiac EP problems in a more challenging pathological scenario as well. The set-up of the FOM is the one provided in the Supplementary Material, except for the final time equal to T = 500 ms. Here, we consider n_μ = 3 parameters, consisting of the coordinates of the center of the S2 intracellular applied currents, which can vary in the three-dimensional region highlighted in Figure 10 (left). The choice of the parameter space is motivated by the fact that ectopic complexes usually arise in correspondence of pulmonary veins (PVs). We first apply a physiological stimulus (S1) in correspondence of the posterior septum and then a second stimulus (S2) acting on ${\Omega }_2(\mu )=\left\{x\in \Omega :{(x-{\mu }_1)}^2+{(y-{\mu }_2)}^2+{(z-{\mu }_3)}^2\le {(0.5)}^2\right\}$, which takes the form

$\begin{array}{l}{I}_{app}^{i,2}(x,t)=C{1}_{{\Omega }_2(\mu )}(x)1_{[t_2^i,t_2^f]}(t),\end{array}$

with C = 100 mA, $t_2^i=210$ ms, and $t_2^f=215$ ms.

This test case represents a proof-of-concept of the strategy used in the clinical practice to identify possible re-entrant circuits, part of which may be latent, by conducting a virtual multi-site delivery of electrical stimuli from a number of possible atrial locations (Arevalo et al., 2016; Boyle et al., 2018; Prakosa et al., 2018).

As before, we set the rPOD dimension equal to N = 256, and the dimension n of the POD-DL-ROM approximation equal to n = n_μ + 1 = 4. We consider N_t = 1, 000 time instances in the interval (300, 500) ms and randomly sample N_train = 15 training-parameter and N_test = 5 testing-parameter instances from the parameter space. The maximum number of epochs is N_epochs = 30, 000, the batch size is N_b = 40, and regarding the early-stopping criterion, we stop the training if the loss function does not decrease along 2,000 epochs. Choosing the rPOD dimension equal to N = 256 yields, over the testing set, a projection error indicator ${ϵ}_{rel}(u_h,V{V}^T{u}_h)=6.8\times 10^{-2}$ and the projection relative error ${ϵ}_k(u_h,V{V}^T{u}_h)$ shown in the Supplementary Material. This value can be used as the lower bound of the reconstruction error indicator, being smaller than the previous values over the testing set.

In Figure 9A, we compare the FOM and POD-DL-ROM solutions, the latter with n = 4 and N = 256, together with ϵ_k, for the testing-parameter instance μ_test = (0.2508, 0.7932, 1.66) cm at t = 316.4 ms. The error indicator ${ϵ}_{rel}(u_h,{\tilde{u}}_h)$ is equal to 7.06 ×10⁻², meaning that the projection error provides an upper bound to the error ${ϵ}_{rel}(u_h,{\tilde{u}}_h)$ over the testing set. However, the POD-DL-ROM is able to completely capture the location and the shape of the re-entry, and the moving front; the error is indeed related to the reconstruction of the steep fronts. Obtained results are thus satisfying, keeping into account the extreme complexity of the problem at hand. Then, we investigate the impact of a higher value for the rPOD dimension, setting it equal to N = 1, 024. In this case, the projection error indicator ${ϵ}_{rel}(u_h,V{V}^T{u}_h)$ is equal to 2.84 ×10⁻² and the error indicator (Equation 12) becomes ${ϵ}_{rel}=5.4\times 10^{-2}$. In Figure 9B, we report the FOM solution and the POD-DL-ROM approximation, obtained with n = 4 and N = 1, 024, together with the relative error (Equation 13), for the testing-parameter instance μ_test = (0.2508, 0.7932, 1.66) cm at t = 316.4 ms. By comparing Figures 9A,B, we can note how the use of a larger N leads to only slightly more accurate results.

FIGURE 9

Figure 9. Test 3: FOM (left) and POD-DL-ROM (center) solutions, the latter obtained with n = 4 and N = 256 (A), and n = 4 and N = 1, 024 (B), together with ϵ_k (right), for the testing-parameter instance μ_test = (0.2508, 0.7932, 1.66) cm at t = 316.4 ms.

In Table 4, we report the FOM CPU computational time and the POD-DL-ROM GPU total, i.e., training and validation time, and testing times, and the total number of epochs n_e, by varying N. As expected, both the training and the testing times are larger for N = 1, 024 than for N = 256, since the number of parameters of the neural network is higher in the former case. We highlight that, if we do not take into account the time needed to assemble the snapshot matrix, the time required to train the POD-DL-ROM over the parameter space, for N = 256, is smaller than performing a FOM simulation for a single parameter instance. We remark that we started from a learning rate equal to η = 2 · 10⁻⁴ for N = 256 and η = 10⁻⁴ for N = 1, 024, the latter resulting in a longer total training and validation time; indeed, in this case training stops because of the maximum number of epochs achieved, however, yielding a higher accuracy. At testing time, both the networks show to be extremely efficient.

TABLE 4

Table 4. Test 3: FOM and POD-DL-ROM computational times.

As done in Fresca et al. (2020), we increase the complexity of the problem by enlarging the dimension of the parameter space, thus considering both re-entry and non re-entry dynamics. We randomly sample N_train = 20+20 = 40 training-parameter and N_test = 10+10 = 20 testing-parameter instances from the parameter space. We set the rPOD dimension equal to N = 1, 024. In this case, the projection error indicator value is ${ϵ}_{rel}(u_h,V{V}^T{u}_h)=4.34\times 10^{-2}$, while the reconstruction error is ${ϵ}_{rel}(u_h,{\tilde{u}}_h)=7.7\times 10^{-2}$. We set the maximum number of epochs N_epochs to 30,000—by increasing this value it is possible to achieve a reconstruction error equal to the projection one. The parameter space is the one shown in Figure 10 (right).

FIGURE 10

Figure 10. Test 3: Possible sites of S2 stimulus applications in the case of re-entry dynamics (magenta region) (A) and including both re-entry and non-re-entry dynamics (magenta region) (B). The coordinates of the points belonging to the highlighted region are the input parameters.

In Figure 11, we report the FOM and POD-DL-ROM solutions, with n = 4 and N = 1, 024, along with ϵ_k, for the testing-parameter instances μ_test = (0.3162, 0.8638, 0.6864) cm and μ_test = (0.2508, 0.7932, 0.8895) cm at t = 300.8 ms. The POD-DL-ROM is then able to reproduce the main features of the dynamics of the solution, and the error is mainly associated with the truncated POD modes.

FIGURE 11

Figure 11. Test 3: FOM (left) and POD-DL-ROM (center) solutions, the latter obtained with n = 4 and N = 1, 024, together with ϵ_k (right), for the testing-parameter instances μ_test = (0.3162, 0.8638, 0.6864) cm (A) and μ_test = (0.2508, 0.7932, 0.8895) cm (B) at t = 300.8 ms.

4. Discussion

The cardiac EP problems addressed in this paper fit into both (i) a multi-query context, since repetitive evaluations of the input-output map are required in order to perform multi-scenario analysis, in order to deal with inter- and intra-subject variability and to consider specific pathological scenarios, and a (ii) real-time context, due to the need, in a clinical setting, to compute outputs of interest in a very limited amount of time. Performing the numerical approximation of cardiac EP problems in these contexts, by means of traditional FOMs, such as the FE method or NURBS-based IGA, is prohibitive because of the huge computational costs associated to the solution of the equations. Indeed, small time-step sizes must be selected to ensure stability; small mesh sizes are required in order to capture the steep fronts and preserve accuracy.

We have taken advantage of a recently proposed technique (Fresca and Manzoni, 2021a) to build low-dimensional ROMs by exploiting DL algorithms. This strategy allows us to overcome typical computational bottlenecks shown by classical, linear projection-based ROM techniques (such as POD-Galerkin ROMs) when dealing with problems featuring coherent structures propagating over time. The DL-ROM technique allows to approximate the solution manifold of a given parametrized nonlinear, time-dependent PDE by means of a low-dimensional, nonlinear trial manifold, and the nonlinear dynamics of the generalized coordinates on such reduced manifold, as a function of the time coordinate and the parameters. Both the nonlinear trial manifold and the reduced dynamics are learnt in a non-intrusive way, thus avoiding to query the FOM high-dimensional arrays. The solution manifold is learnt by means of the decoder function of a convolutional AE neural network; the reduced dynamics is approximated through a DFNN and the encoder function of the convolutional AE. Through the use of the DL-ROM, it is possible to boost the solution of parametrized problems in cardiac EP remarkably, thus overcoming the main computational bottlenecks that affect POD-Galerkin ROMs in this context (Fresca et al., 2020). A key aspect in the setting of DL-ROMs concerns their computational efficiency during the offline (or training) stage, which is also related with the curse of dimensionality entailed by the (possibly, extremely large) dimension of the FOM. This gain, which makes the offline training stage dramatically faster, hinges upon (i) a preliminary dimensionality reduction in the FOM snapshots, by means of rPOD, and (ii) a suitable multi-fidelity pretraining stage, exploiting snapshots computed through different low-fidelity models to initialize the parameters of the neural networks in a sequential procedure.

So far, only few works have focused on the solution, by means of DL algorithms, of problems featuring traveling waves or front propagation processes in the cardiac EP context. For example, in Court and Kunisch (2021) the ionic model is designed to exploit an artificial neural network, in order to identify the nonlinearity in the Monodomain model from given data, yet without providing information about neither the spatial distribution of the electrical signal in the heart, nor the whole range of time and spatial scales of the transmembrane potential. The reconstruction of ACs by means of a physics-informed neural network (PINN) trained by minimizing the residual associated with the Eikonal equation is addressed by Sahli Costabal et al. (2020); several techniques based on ML algorithms are reviewed in Cantwell et al. (2019), for the sake of addressing either classification or estimation problems, such as, e.g., prediction making from the contact electrogram. Finally, neural networks are used for the numerical integration of the Monodomain equation coupled with the Mitchell-Schaeffer ionic model, assessing their performance on two-dimensional benchmarks, in Ayed et al. (2019) and Kashtanova et al. (2021).

In this study, we assessed the performance of the POD DL-ROM technique when applied to the solution of cardiac EP problems on a left atrium geometry, in both physiological and pathological scenarios, by showing its ability in providing an accurate and efficient ROM, which multi-query and real-time problems may rely on. Indeed, POD-DL-ROMs enable to explore the parameter space, thus accounting for different scenarios, and it not only provide real-time solutions to parametrized cardiac EP problems at the testing stage—being able to match the intrinsic dimension of the problems investigated—but can also be trained very efficiently. Moreover, we point out that it is also possible to include more complex ionic models in the FOM, for a more accurate description of the electrical activity of the heart at the microscopic level, without affecting the computational times of the POD-DL-ROM. Indeed, due to the non-intrusive nature of this technique, the dynamics of the gating variables is not taken explicitly into account by the networks in order to compute the electrical potential. In the same way, the choice of a particular model of fibers and the definition of the conductivity tensor (possibly accounting for the presence of ischemic, non-conductive regions as in Fresca et al., 2020; Kashtanova et al., 2021), are considered by the neural network only through the effects they produce on the FOM snapshots. The accuracy and the efficiency obtained by the POD-DL-ROM approximations make them amenable, in the clinical setting, to replace high-fidelity, FOM solvers, for the computation of quantities of interest, such as ACs and APs.

Finally, we highlight that a possible pitfall of the proposed methodology is represented by the amount/quality of training data: If too few (or low-quality) snapshots are considered, further operations like (i) increasing the number of parameters of the network, (ii) training the network for a larger number of epochs, or (iii) generating more data by means of data augmentation techniques can be required. A relevant issue is also related to the generalization properties of the network outside the parameter range and/or the time interval where snapshots are sampled. Ensuring good approximation properties when interested in long-time scenarios, even in presence of almost periodic regimes, without more specific network architectures, is an open issue our efforts are focusing on; however, this represents a general aspect shared by several ROM techniques.

To the best of our knowledge, this study represents the first attempt of reducing the computational complexity associated with the reconstruction of both the transmembrane and the extracellular potentials and re-entry problems, this virtually opening a new path toward the model personalization in real-time, even when dealing with extremely challenging, and computationally involved, settings. We remark that the performance of the POD-DL-ROM technique evaluated on new, unseen scenarios with respect to the ones used during the training phase of the network, thus virtually allowing to compute, during interventions, outputs related to subject-specific data such as, e.g., ACs o voltage maps, in real-time. The possibility to perform real-time numerical simulations, in cardiac EP, can be seen as the first step toward the translation of computational methods into the clinical practice enabling a cooperation for supporting decisions, quantifying risks related to cardiac pathologies, planning therapies, and interventions.

Data Availability Statement

Publicly available datasets were analyzed in this study. The code used in the numerical tests is freely available at https://github.com/stefaniafresca/POD-DL-ROM.

Author Contributions

SF, AM, LD, and AQ: conceptualization. SF: formal analysis, investigation, software, validation, and visualization. AQ: funding acquisition. SF and AM: methodology and writing—original draft. AM, LD, and AQ: supervision and writing—review and editing. All authors contributed to the article and approved the submitted version.

Funding

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 740132).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2021.679076/full#supplementary-material

Footnotes

1. ^We can distringuish among phenomenological, first generation, and second generation ionic models (Sundnes et al., 2006; Colli Franzone et al., 2014). Compared to phenomenological models, first, and second-generation models attempt to include also a description of the cell mechanisms, with the latter also including sub-cellular processes. In the following, we will focus only on phenomenological models.

2. ^In the A-P ionic model, dimensional times and potentials are given by $\tilde{t}=12.9t\left[ms\right]$ and ũ = (100u − 80)[mV]. The transmembrane potential ranges from the resting state of −80 mV to the excited state of +20 mV.

References

Aliev, R. R., and Panfilov, A. V. (1996). A simple two-variable model of cardiac excitation. Chaos Solitons Fractals 7, 293–301. doi: 10.1016/0960-0779(95)00089-5

CrossRef Full Text | Google Scholar

Arevalo, H., Vadakkumpadan, F., Guallar, E., Jebb, A., Malamas, P., Wu, K., et al. (2016). Arrhythmia risk stratification of patients after myocardial infarction using personalized heart models. Nat. Commun. 7:11437. doi: 10.1038/ncomms11437

PubMed Abstract | CrossRef Full Text | Google Scholar

Ayed, I., Cedilnik, N., Gallinari, P., and Sermesant, M. (2019). “Ep-net: learning cardiac electrophysiology models for physiology-based constraints in data-driven predictions,” in Functional Imaging and Modeling of the Heart, eds Y. Coudiére, V. Ozenne, E. Vigmond, and N. Zemzemi (Cham: Springer).

Google Scholar

Bhatnagar, S., Afshar, Y., Pan, S., Duraisamy, K., and Kaushik, S. (2019). Prediction of aerodynamic flow fields using convolutional neural networks. Comp. Mech. 64, 525–545. doi: 10.1007/s00466-019-01740-0

CrossRef Full Text | Google Scholar

Boyle, P., Hakim, J. B., Zahid, S., Franceschi, W. H., Murphy, M. J., Prakosa, A., et al. (2018). The fibrotic substrate in persistent atrial fibrillation patients: comparison between predictions from computational modeling and measurements from focal impulse and rotor mapping. Front. Physiol. 9:1151. doi: 10.3389/fphys.2018.01151

PubMed Abstract | CrossRef Full Text | Google Scholar

Brunton, S., Noack, B., and Koumoutsakos, P. (2020). Machine learning for fluid mechanics. Ann. Rev. Fluid Mech. 52, 477–508. doi: 10.1146/annurev-fluid-010719-060214

CrossRef Full Text | Google Scholar

Bucelli, M., Salvador, M., and Dedè, L. (2021). Multipatch isogeometric analysis for electrophysiology: simulation in a human heart. Comput. Methods Appl. Mech. Eng. 376:113666. doi: 10.1016/j.cma.2021.113666

CrossRef Full Text | Google Scholar

Cantwell, C. D., Mohamied, Y., Tzortzis, K. N., Garasto, S., Houston, C., Chowdhury, R. A., et al. (2019). Rethinking multiscale cardiac electrophysiology with machine learning and predictive modelling. Comput. Biol. Med. 104, 339–351. doi: 10.1016/j.compbiomed.2018.10.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Clayton, R., Bernus, O., Cherry, E., Dierckx, H., Fenton, F., Mirabella, L., et al. (2011). Models of cardiac tissue electrophysiology: progress, challenges and open questions. Progr. Biophys. Mol. Biol. 104, 22–48. doi: 10.1016/j.pbiomolbio.2010.05.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Colli Franzone, P., Pavarino, L., and Taccardi, B. (2005). Simulating patterns of excitation, repolarization and action potential duration with cardiac bidomain and monodomain models. Math. Biosci. 197, 35–66. doi: 10.1016/j.mbs.2005.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Colli Franzone, P., and Pavarino, L. F. (2004). A parallel solver for reaction-diffusion systems in computational electrocardiology. Math. Models Methods Appl. Sci. 14, 883–911. doi: 10.1142/S0218202504003489

CrossRef Full Text | Google Scholar

Colli Franzone, P., Pavarino, L. F., and Savaré, G. (2006). “Computational electrocardiology: mathematical and numerical modeling,” in Complex Systems in Biomedicine. eds A. Quarteroni, L. Formaggia, A. Veneziani, Milan: Springer, doi: 10.1007/88-470-0396-2_6

CrossRef Full Text | Google Scholar

Colli Franzone, P., Pavarino, L. F., and Scacchi, S. (2014). Mathematical Cardiac Electrophysiology, Vol. 13 of Modeling, Simulation &Applications. Cham: Springer International Publishing.

Google Scholar

Cottrell, J. A., Hughes, T. J. R., and Bazilevs, Y. (2009). Isogeometric Analysis: Toward Integration of CAD and FEA, 1st Edn. Chichester: Wiley Publishing.

Google Scholar

Court, S., and Kunisch, K. (2021). Design of the monodomain model by artificial neural networks. arXiv:2107.03136.

Google Scholar

Dedè, L., Jäggli, C., and Quarteroni, A. (2015). Isogeometric numerical dispersion analysis for two-dimensional elastic wave propagation. Comput Methods Appl. Mech. Eng. 284, 320–348. doi: 10.1016/j.cma.2014.09.013

CrossRef Full Text | Google Scholar

Eichinger, M., Heinlein, A., and Klawonn, A. (2021). “Stationary flow predictions using convolutional neural networks,” in Numerical Mathematics and Advanced Applications ENUMATH 2019, Vol. 139 of Lecture Notes in Computational Science and Engineering, eds F. Vermolen and C. Vuik (Cham: Springer), 541–549.

Google Scholar

FitzHugh, R. (1961). Impulses and physiological states in theoretical models of nerve membrane. Biophys. J. 1, 445–466. doi: 10.1016/S0006-3495(61)86902-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Fotiadis, S., Pignatelli, E., Lino Valencia, M., Cantwell, C. D., Storkey, A. J., and Bharath, A. A. (2020). Comparing recurrent and convolutional neural networks for predicting wave propagation. arXiv:2002.08981.

Google Scholar

Fresca, S., Dedè, L., and Manzoni, A. (2021). A comprehensive deep learning-based approach to reduced order modeling of nonlinear time-dependent parametrized PDEs. J. Sci. Comput. 87:61. doi: 10.1007/s10915-021-01462-7

CrossRef Full Text | Google Scholar

Fresca, S., and Manzoni, A. (2021a). POD DL-ROM: enhancing deep learning-based reduced order models for nonlinear parametrized PDEs by proper orthogonal decomposition. arXiv:2101.11845.

Google Scholar

Fresca, S., and Manzoni, A. (2021b). Real-time simulation of parameter-dependent fluid flows through deep learning-based reduced order models. Fluids 6:259. doi: 10.3390/fluids6070259

CrossRef Full Text | Google Scholar

Fresca, S., Manzoni, A., Dedè, L., and Quarteroni, A. (2020). Deep learning-based reduced order models in cardiac electrophysiology. PLoS ONE 15:e0239416. doi: 10.1371/journal.pone.0239416

PubMed Abstract | CrossRef Full Text | Google Scholar

Garotta, F., Demo, N., Tezzele, M., Carraturo, M., Reali, A., and Rozza, G. (2020). “Reduced order isogeometric analysis approach for pdes in parametrized domains,” in Quantification of Uncertainty: Improving Efficiency and Technology: QUIET Selected Contributions, eds M. D'Elia, M. Gunzburger, and G. Rozza (Cham: Springer International Publishing), 153–170.

Google Scholar

Gerardo-Giorda, L. (2007). Modeling and numerical simulation of action potential patterns in human atrial tissues. Preprint hal-00132706.

Google Scholar

Geselowitz, D. B., and Miller, I. I. I. W. T. (1983). A bidomain model for anisotropic cardiac muscle. Ann. Biomed. Eng. 11, 191–206. doi: 10.1007/BF02363286

PubMed Abstract | CrossRef Full Text | Google Scholar

Goodfellow, I., Bengio, Y., and Courville, A. (2016). Deep Learning. Cambridge, MA: MIT Press.

Google Scholar

Halko, N., Martinsson, P., and Tropp, J. A. (2011). Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. SIAM Rev. 53, 217–288. doi: 10.1137/090771806

CrossRef Full Text | Google Scholar

Kashtanova, V., Ayed, I., Cedilnik, N., Gallinari, P., and Sermesant, M. (2021). “Ep-net 2.0: out-of-domain generalisation for deep learning models of cardiac electrophysiology,” in Functional Imaging and Modeling of the Heart, eds D. B. Ennis, L. E. Perotti, and V. Y. Wang (Cham: Springer International Publishing), 482–492.

Google Scholar

Klabunde, R. (2011). Cardiovascular Physiology Concepts. Philadelphia, PA: Wolters Kluwer Health/Lippincott Williams &Wilkins.

Google Scholar

Kutz, N. J. (2017). Deep learning in fluid dynamics. J. Fluid Mech. 814, 1–4. doi: 10.1017/jfm.2016.803

CrossRef Full Text | Google Scholar

Lee, K., and Carlberg, K. (2020). Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders. J. Comput. Phys. 404:108973. doi: 10.1016/j.jcp.2019.108973

CrossRef Full Text | Google Scholar

Manzoni, A., Salmoiraghi, F., and Heltai, L. (2015). Reduced basis isogeometric methods (RB-IGA) for the real-time simulation of potential flows about parametrized NACA airfoils. Comput. Methods Appl. Mech. Eng. 284, 1147–1180. doi: 10.1016/j.cma.2014.11.037

CrossRef Full Text | Google Scholar

Mitchell, C. C., and Schaeffer, D. G. (2003). A two-current model for the dynamics of cardiac membrane. Bull. Math. Biol. 65, 767–793. doi: 10.1016/S0092-8240(03)00041-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagumo, J., Arimoto, S., and Yoshizawa, S. (1962). An active pulse transmission line simulating nerve axon. Proc. IRE 50, 2061–2070. doi: 10.1109/JRPROC.1962.288235

CrossRef Full Text | Google Scholar

Nash, M. P., and Panfilov, A. V. (2004). Electromechanical model of excitable tissue to study reentrant cardiac arrhythmias. Prog. Biophys. Mol. Biol. 85, 501–522. doi: 10.1016/j.pbiomolbio.2004.01.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Nattel, S. (2002). New ideas about atrial fibrillation 50 years on. Nature 415, 219–226. doi: 10.1038/415219a

PubMed Abstract | CrossRef Full Text | Google Scholar

Niederer, S., Mitchell, L., Smith, N., and Plank, G. (2011). Simulating human cardiac electrophysiology on clinical time-scales. Front. Physiol. 2:14. doi: 10.3389/fphys.2011.00014

PubMed Abstract | CrossRef Full Text | Google Scholar

Niederer, S. A., Aboelkassem, Y., Cantwell, C. D., Corrado, C., Coveney, S., Cherry, E. M., et al. (2020). Creation and application of virtual patient cohorts of heart models. Philos. Trans. R. Soc. Math. Phys. Eng. Sci. 378:20190558. doi: 10.1098/rsta.2019.0558

PubMed Abstract | CrossRef Full Text | Google Scholar

Niederer, S. A., Fink, M., Noble, D., and Smith, N. P. (2009). A meta-analysis of cardiac electrophysiology computational models. Exp. Physiol. 94, 486–495. doi: 10.1113/expphysiol.2008.044610

PubMed Abstract | CrossRef Full Text | Google Scholar

Pagani, S., Manzoni, A., and Quarteroni, A. (2018). Numerical approximation of parametrized problems in cardiac electrophysiology by a local reduced basis method. Comput. Methods Appl. Mech. Eng. 340, 530–558. doi: 10.1016/j.cma.2018.06.003

CrossRef Full Text | Google Scholar

Patelli, A. S., Dedè, L., Lassila, T., Bartezzaghi, A., and Quarteroni, A. (2017). Isogeometric approximation of cardiac electrophysiology models on surfaces: an accuracy study with application to the human left atrium. Comput. Methods Appl. Mech. Eng. 317, 248–273. doi: 10.1016/j.cma.2016.12.022

CrossRef Full Text | Google Scholar

Pegolotti, L., Dedè, L., and Quarteroni, A. (2019). Isogeometric analysis of the electrophysiology in the human heart: numerical simulation of the bidomain equations on the atria. Comput. Methods Appl. Mech. Eng. 343, 52–73. doi: 10.1016/j.cma.2018.08.032

CrossRef Full Text | Google Scholar

Prakosa, A., Arevalo, H. J., Deng, D., Boyle, P. M., Nikolov, P. P., Ashikaga, H., et al. (2018). Personalized virtual-heart technology for guiding the ablation of infarct-related ventricular tachycardia. Nat. Biomed. Eng. 2, 732–740. doi: 10.1038/s41551-018-0282-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Quarteroni, A. (2017). Numerical Models for Differential Problems, 3rd Edn. Cham: Springer International Publishing.

Google Scholar

Quarteroni, A., Dedè, L., Manzoni, A., and Vergara, C. (2019). Mathematical Modelling of the Human Cardiovascular System: Data, Numerical Approximation, Clinical Applications.. Cambridge, UK: Cambridge University Press.

Google Scholar

Quarteroni, A., Manzoni, A., and Negri, F. (2016). Reduced Basis Methods for Partial Differential Equations: An Introduction. Cham: Springer International Publishing.

Google Scholar

Quarteroni, A., Manzoni, A., and Vergara, C. (2017). The cardiovascular system: Mathematical modeling, numerical algorithms, clinical applications. Acta Numerica 26, 365–590. doi: 10.1017/S0962492917000046

PubMed Abstract | CrossRef Full Text | Google Scholar

Quarteroni, A., and Valli, A. (1994). Numerical Approximation of Partial Differential Equations. Berlin Heidelberg: Springer-Verlag.

PubMed Abstract | Google Scholar

Rinaldi, M. (2015). Reduced Basis Method for Isogeometric Analysis: Application to Structural Problems. Master thesis. Politecnico di Milano.

Google Scholar

Rogers, J. M., and McCulloch, A. D. (1994). A collocation-galerkin finite element model of cardiac action potential propagation. IEEE Trans. Biomed. Eng. 41, 743–757. doi: 10.1109/10.310090

PubMed Abstract | CrossRef Full Text | Google Scholar

Rossi, S., Lassila, T., Ruiz Baier, R., Sequeira, A., and Quarteroni, A. (2014). Thermodynamically consistent orthotropic activation model capturing ventricular systolic wall thickening in cardiac electromechanics. Eur. J. Mech. 48, 129–142. doi: 10.1016/j.euromechsol.2013.10.009

CrossRef Full Text | Google Scholar

Sahli Costabal, F., Yang, Y., Perdikaris, P., Hurtado, D. E., and Kuhl, E. (2020). Physics-informed neural networks for cardiac activation mapping. Front. Phys. 8:42. doi: 10.3389/fphy.2020.00042

CrossRef Full Text | Google Scholar

Salmoiraghi, F., Ballarin, F., Heltai, L., and Rozza, G. (2016). Isogeometric analysis-based reduced order modelling for incompressible linear viscous flows in parametrized shapes. Adv. Model. Simul. Eng. Sci. 3:21. doi: 10.1186/s40323-016-0076-6

CrossRef Full Text | Google Scholar

Strocchi, M., Augustin, C. M., Gsell, M. A. F., Karabelas, E., Neic, A., Gillette, K., et al. (2020). A publicly available virtual cohort of four-chamber heart meshes for cardiac electro-mechanics simulations. PLoS ONE 15:e0235145. doi: 10.1371/journal.pone.0235145

PubMed Abstract | CrossRef Full Text | Google Scholar

Ströfer, C., Wu, J.-L., Xiao, H., and Paterson, E. (2019). Data-driven, physics-based feature extraction from fluid flow fields using convolutional neural networks. Commun. Comput. Phys. 25, 625–650. doi: 10.4208/cicp.OA-2018-0035

CrossRef Full Text | Google Scholar

Sundnes, J., Lines, G. T., Cai, X., Nielsen, B. F., Mardal, K., and Tveito, A. (2006). Computing the Electrical Activity in the Heart. Berlin; Heidelberg: Springer-Verlag.

PubMed Abstract | Google Scholar

ten Tusscher, K. (2004). Spiral Wave Dynamics and Ventricular Arrhythmias. (Ph.D Thesis), Department of Theoretical Biology and Bioinformatics, Utrecht University.

Google Scholar

ten Tusscher, K. H. W. J., and Panfilov, A. V. (2006). Alternans and spiral breakup in a human ventricular tissue model. Am. J. Physiol. Heart Circ. Physiol. 291, 1088–1100. doi: 10.1152/ajpheart.00109.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Thuerey, N., Weißenow, K., Prantl, L., and Hu, X. (2020). Deep learning methods for reynolds-averaged Navier-Stokes simulations of airfoil flows. AIAA J. 58, 25–36. doi: 10.2514/1.J058291

CrossRef Full Text | Google Scholar

Trayanova, N. A. (2011). Whole-heart modeling applications to cardiac electrophysiology and electromechanics. Circ. Res. 108, 113–128. doi: 10.1161/CIRCRESAHA.110.223610

PubMed Abstract | CrossRef Full Text | Google Scholar

Vigmond, E., Dos Santos, R. W., Prassl, A., Deo, M., and Plank, G. (2008). Solvers for the cardiac bidomain equations. Prog. Biophys. Mol. Biol. 96, 3–18. doi: 10.1016/j.pbiomolbio.2007.07.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Vigmond, E. J., Aguel, F., and Trayanova, N. A. (2002). Computational techniques for solving the bidomain equations in three dimensions. IEEE Trans. Biomed. Eng. 49, 1260–1269. doi: 10.1109/TBME.2002.804597

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, S., Dedè, L., and Quarteroni, A. (2017a). Isogeometric analysis and proper orthogonal decomposition for parabolic problems. Numer. Math. 135, 333–370. doi: 10.1007/s00211-016-0802-5

CrossRef Full Text | Google Scholar

Zhu, S., Dedè, L., and Quarteroni, A. (2017b). Isogeometric analysis and proper orthogonal decomposition for the acoustic wave equation. ESAIM: Math. Mod. Numer. Anal. 51, 1197–1221. doi: 10.1051/m2an/2016056

CrossRef Full Text | Google Scholar