Deep learning-based reduced order models in cardiac electrophysiology

Predicting the electrical behavior of the heart, from the cellular scale to the tissue level, relies on the numerical approximation of coupled nonlinear dynamical systems. These systems describe the cardiac action potential, that is the polarization/depolarization cycle occurring at every heart beat that models the time evolution of the electrical potential across the cell membrane, as well as a set of ionic variables. Multiple solutions of these systems, corresponding to different model inputs, are required to evaluate outputs of clinical interest, such as activation maps and action potential duration. More importantly, these models feature coherent structures that propagate over time, such as wavefronts. These systems can hardly be reduced to lower dimensional problems by conventional reduced order models (ROMs) such as, e.g., the reduced basis method. This is primarily due to the low regularity of the solution manifold (with respect to the problem parameters), as well as to the nonlinear nature of the input-output maps that we intend to reconstruct numerically. To overcome this difficulty, in this paper we propose a new, nonlinear approach relying on deep learning (DL) algorithms—such as deep feedforward neural networks and convolutional autoencoders—to obtain accurate and efficient ROMs, whose dimensionality matches the number of system parameters. We show that the proposed DL-ROM framework can efficiently provide solutions to parametrized electrophysiology problems, thus enabling multi-scenario analysis in pathological cases. We investigate four challenging test cases in cardiac electrophysiology, thus demonstrating that DL-ROM outperforms classical projection-based ROMs.


Introduction
The electrical activation of the heart is the main responsible of its contraction, is the result of two processes: at the microscopic scale, the generation of ionic currents through the cellular membrane producing a local action potential; and at the macroscopic scale, the propagation of the action potential from cell to cell in the form of a transmembrane potential [1][2][3].This latter process can be described by means of partial differential equations (PDEs), suitably coupled with systems of ordinary differential equations (ODEs) modeling the ionic currents in the cells.
Solving this system using a high-fidelity, full order model (FOM) such as, e.g., the finite element (FE) method, is computationally demanding.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.[3][4][5].Using a FOM may quickly become unaffordable if such a coupled system must be solved for several values of parameters representing either functional or geometric data such as, e.g., material properties, initial and boundary conditions, or the shape of the domain.Multi-query analysis is relevant in a variety of situations: when analysing multiple scenarios, when dealing with sensitivity analysis and uncertainty quantification (UQ) problems in order to account for inter-subject variability [6][7][8], for parameter estimation or data assimilation, in which some unknown (or unaccessible) quantities characterizing the mathematical model must be inferred from a set of measurements [9][10][11][12][13].
Conventional projection-based reduced order models (ROMs) built, e.g., through the reduced basis (RB) method [14], yields inefficient ROMs when dealing with nonlinear time-dependent parametrized PDE-ODE system as the one arising from cardiac electrophysiology.The three major computational bottlenecks shown by such kind of ROMs for cardiac electrophysiology are due the fact that: -the linear superimposition of modes, on which they are based, would cause the dimension of the ROM to be excessively large to guarantee an acceptable accuracy; -evaluating the ROM requires the solution of a dynamical system, which might be unstable unless the size of time step ∆t is very small; -the ROM must also account for the dynamics of the gating variables, even when aiming at computing just the electrical potential.This fact entails an extremely intrusive and costly hyper-reduction stage to reduce the solution of the ODE system to a few, selected mesh nodes [15].
To overcome the limitations of projection-based ROMs, we propose a new, nonintrusive ROM technique based on deep learning (DL) algorithms, which we refer to as DL-ROM.Combining in a suitable way a convolutional autoencoder (AE) and a deep feedforward neural network (DFNN), the DL-ROM technique enables the construction of an efficient ROM, whose dimension is as close as possible to the number of parameters upon which the solution of the differential problem depends.A preliminary numerical assessment of our DL-ROM technique has already been presented in [16], albeit on simpler -yet challenging -test cases.
The proposed DL-ROM technique is a combination of a data-driven with a physics based model approach.Indeed, it exploits snapshots taken from a set of FOM solutions (for selected parameter values and time instances) and deep neural network architectures to learn, in a non-intrusive way, both (i) the nonlinear trial manifold where the ROM solution is sought, and (ii) the nonlinear reduced dynamics.In a linear ROM built, e.g., thorugh proper orthogonal decomposition (POD), the former quantity is nothing but a set of basis functions, while the latter task corresponds to the projection stage in the subspace spanned by these basis functions.Here, our goal is to show that DL-ROM can be effectively used to handle parametrized problems in cardiac electrophysiology, accounting for both physiological and pathological conditions, in order to provide fast and accurate solutions.The proposed DL-ROM is computationally efficient during the testing stage, that is for any new scenario unseen during the training stage.This is particularly useful in view of the evaluation of patient-specific features to enable the integration of computational methods in current clinical platforms.
DL techniques for parametrized PDEs have previously been proposed in other contexts.In [17][18][19][20] feedforward neural networks have been employed to model the reduced dynamics in a less intrusive way, that is, avoiding the costs entailed by projection-based ROMs, but still relying on a linear trial manifold built, e.g., through POD.In [21][22][23] the construction of ROMs for nonlinear, time-dependent problems has been replaced by the evaluation of ANN-based regression models.In [24,25] the reduced trial manifold where the approximation is sought has been modeled through ANNs thus avoiding the linear superimposition of POD modes, on a minimum residual formulation to derive the ROM [25], or without considering an explicit parameter dependence in the differential problem that is considered [24].In all these works, coupled problems have never been considered.Moreover, very often DL techniques have been exploited to address problems which require only a moderate dimension of projection-based ROMs.We demonstrate that our DL-ROM provides accurate results by constructing ROMs with extremely low-dimension in prototypical test cases.These tests exhibit all the relevant physical features which make the numerical approximation of parametrized problems in cardiac electrophysiology a challenging task.

Cardiac electrophysiology
Muscle contraction and relaxation drive the pump function of the heart.In particular, tissue contraction is triggered by electrical signals self-generated in the heart and propagated through the myocardium thanks to the excitability of the cardiac cells, the cardiomyocites [3,26].When suitably stimulated, cardiomyocites produce a variation of the potential across the cellular membrane, called transmembrane potential.Its evolution in time is usually referred to as action potential, involving a polarization and a depolarization in the early stage of every heart beat.The action potential is generated by several ion channels (e.g., calcium, sodium, potassium) that open and close, and by the resulting ionic currents crossing the membrane.For instance, coupling the so-called monodomain model for the transmembrane potential u = u(x, t) with a phenomenological model for the ionic currents -involving a single gating variable w = w(x, t) -in a domain Ω representing, e.g., a portion of the myocardium, results in the following nonlinear time-dependent system Here t denotes a rescaled time1 , n denotes the outward directed unit vector normal to the boundary ∂Ω of Ω, whereas I app is an applied current representing, e.g., the initial activation of the tissue.The nonlinear diffusion-reaction equation for u is two-ways coupled with the ODE system, which must be in principle solved at any point x ∈ Ω; indeed, the reaction term I ion and the function g depend on both u and w.The most common choices for the two functions I ion and g in order to efficiently reproduce the action-potential are, e.g., the FitzHugh-Nagumo [28,29], the Aliev-Panfilov [27,30] or the Mitchell and Schaeffer models [31].The diffusivity tensor D usually depends on the fibers-sheet structure of the tissue, affecting directional conduction velocities and directions.In particular, by assuming an axisymmetric distribution of the fibers, the conductivity tensor takes the form where σ l and σ t are the conductivities in the fibers and the transversal directions.When a simple phenomenological ionic model is considered, such as the FitzHugh-Nagumo or the Aliev-Panfilov (A-P) 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 refractariness of cells.In this paper, we focus on the Aliev-Panfilov model, which consists in taking The parameters K, a, b, ε 0 , c 1 , c 2 are related to the cell.Here a represents an oscillation threshold, whereas the weighting factor ε 0 + c1w c2+u was introduced in [27] to tune the restitution curve to experimental observations by adjusting the parameters c 1 and c 2 ; see, e.g., [1][2][3]32] for a detailed review.In the remaining part of the paper, we denote by µ ∈ P ⊂ R nµ a parameter vector listing all the n µ input parameters characterizing physical (and, possibly, geometrical) properties we might be interested to vary; P is a subset of R nµ , 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 ; previous analyses focused instead on the gating variable dynamics (through g) and the ionic current I ion , see [15].

Projection-based ROMs
From an algebraic standpoint, the spatial discretization of system (1) through the Galerkin-finite element (FE) approximation [33] yields the following nonlinear dynamical system for u = u(t; µ), w = w(t; µ), representing our full order model (FOM): Here A(µ) ∈ R N ×N is a matrix arising from the diffusion operator (thus including the conductivity tensor D(µ) = D(x; µ), which can vary within the myocardium due to fiber orientation and conditions, such as the possible presence of ischemic regions); M(µ) ∈ R N ×N is the mass matrix; I ion , g ∈ R N are vectors arising from the nonlinear terms; I app ∈ R N is a vector collecting the applied currents; finally, u 0 , w 0 ∈ R N are the initial data, possibly depending on µ.The dimension N is related to the dimension of the FE space and, ultimately, depends on the size h > 0 of the computational mesh used to discretize the domain Ω.Note that the system of ODEs arises from the collocation of the ODE (1) 2 at the nodes used for the numerical integration.The intrinsic dimension of the solution manifold obtained by solving (4) when (t; µ) varies in [0, T ) × P, is usually much smaller than N and, under suitable conditions, is at most n µ + 1 ≪ N , where n µ is the number of parameters -in this respect, the time independent variable plays the role of a parameter.For this reason, ROMs attempt at approximating S by introducing a suitable trial manifold of lower dimension.The most popular approach is proper orthogonal decomposition (POD), which exploits a linear trial manifold built through the singular value decomposition of a matrix S ∈ R N ×Ns collecting a set of FOM snapshots this is a set of solutions obtained for N train selected input parameters at (a subset, possibly, of) the time instants {t k } Nt k=1 in which (0, T ) is partitioned for the sake of time discretization.The most common choice is to set t k = k∆t where ∆t = T /(N t − 1).
When using a projection-based ROM, the approximation of u(t; µ) is sought as a linear superimposition of modes, under the form thus yielding a linear ROM, in which the columns of the matrix V = [ζ 1 , . . ., ζ n ] ∈ R N ×n form an orthonormal basis of a space V n , an n-dimensional subspace of R N .In the case of POD, V n provides the best n-rank approximation of S in the Frobenius norm, that is, ζ 1 , . . ., ζ n are the first n (left) singular vectors of S corresponding to the n largest singular values σ 1 , . . ., σ n of S, such that the projection error is smaller than a desired tolerance ε P OD .To meet this requirement, it is sufficient to choose n as the smallest integer such that i.e., the energy retained by the last N s − n POD modes is equal or smaller than ε 2 P OD .The approximation of w is given instead by its restriction w(t; µ) ≈ Pw m (t; µ), to a (possibly, small) subset of m degrees of freedom, where m ≪ n, at which the nonlinear term I ion is interpolated exploiting a problem-dependent basis, spanned by the columns of a matrix Φ ∈ R N ×m , which is built according to a suitable hyperreduction strategy; see, e.g., [15] for further details.Here P = [e 1 , . . ., e m ] ∈ R N ×m denotes a matrix formed by the columns of the N × N identity matrix corresponding to the m selected degrees of freedom.
A Galerkin-POD ROM for system (1) is then obtained by (i) first, substituting equation ( 6) into equation 4 and projecting it onto V n ; then, (ii) solving the system of ODEs at m selected degrees of freedom, thus yielding the following nonlinear dynamical system for u n = u n (t; µ) and the selected components P T w = P T w(t; µ) of w: This strategy is the essence of the reduced basis (RB) method for nonlinear timedependent parametrized PDEs.However, using (7) as an approximation to (4) is known to suffer from several problems.First of all, an extensive hyper-reduction stage (exploiting, e.g., the discrete empirical interpolation method (DEIM)) must be performed in order to be able to evaluate any µ-or u-dependent quantities appearing in (7), that is, without relying on N -dimensional arrays.Moreover, whenever the solution of the differential problem features coherent structures that propagate over time, such as steep wavefronts, the dimension n of the projection-based ROM (7) might easily become very large, due to the basic linearity assumption, by which the solution is given by a linear superimposition of POD modes, thus severely degrading the computational efficiency of the ROM.A possible way to overcome this bottleneck is to rely on local reduced bases, built through POD after the set of snapshots has been split into N c > 1 clusters, according to suitable clustering (or unsupervised learning) algorithms [15].

Deep learning-based reduced order modeling (DL-ROM)
To overcome the limitations of linear ROMs, we consider a new, nonlinear ROM technique based on deep learning models.First introduced in [16] and assessed on onedimensional benchmark problems, the DL-ROM technique aims at learning both the nonlinear trial manifold (corresponding to the matrix V in the case of a linear ROM) in which we seek the solution to the parametrized system (1) and the nonlinear reduced dynamics (corresponding to the projection stage in a linear ROM).This method is not intrusive; it relies on DL algorithms trained on a set of FOM solutions obtained for different parameter values.
We denote by N train and N test the number of training and testing parameter instances, respectively; the ROM dimension is again denoted by n ≪ N .In order to describe the system dynamics on a suitable reduced nonlinear trial manifold (a task which we refer to as reduced dynamics learning), the intrinsic coordinates of the ROM approximation are defined as where φ DF n (•; •, θ DF ) : R (nµ+1) → R n is a deep feedforward neural network (DFNN), consisting in the subsequent composition of a nonlinear activation function, applied to a linear transformation of the input, multiple times [34].Here θ DF denotes the vector of parameters of the DFNN, collecting all the corresponding weights and biases of each layer of the DFNN.
Regarding instead the description of the reduced nonlinear trial manifold, approximating the solution one, S ≈ S (a task which we refer to as reduced trial manifold learning) we employ the decoder function of a convolutional autoencoder 2 (AE) [35,36].More precisely, S takes the form where f D (•; θ D ) : R n → R N consists in the decoder function of a convolutional AE.This latter results from the composition of several layers (some of which are convolutional), depending upon a vector θ D collecting all the corresponding weights and biases.
As a matter of fact, the approximation ũ(t; µ) ≈ u(t; µ) provided by the DL-ROM technique is defined as The encoder function of the convolutional AE can then be exploited to map the FOM solution associated to (t, µ) onto a low-dimensional representation f E n (•; θ E ) : R N → R n denotes the encoder function, depending upon a vector θ E of parameters.
Computing the DL-ROM approximation of u(t; µ test ), for any possible t ∈ (0, T ) and µ test ∈ P, corresponds to the testing stage of a DFNN and of the decoder function of a 2 The AE is a particular type of neural network aiming at learning the identity function . It is composed by two main parts: convolutional AE; this does not require the evaluation of the encoder function.We remark that our DL-ROM strategy overcomes the three major computational bottlenecks implied by the use of projection-based ROMs, since: -the dimension of the DL-ROM can be kept extremely small; -the time resolution required by the DL-ROM can be chosen to be larger than the one required by the numerical solution of dynamical systems in cardiac electrophysiology; -the DL-ROM can be queried at any desired time instant, without requiring the solution of a dynamical system until that time; -the DL-ROM does not require to account for the dynamics of the gating variables, thus avoiding any hyper-reduction stage.This advantage, already visible when employing a single gating variable as in the test cases addressed later in this paper, might become even more effective when dealing with more realistic ionic models (the so-called I and II generation models), when dozens of additional variables in the system of ODEs must be accounted for [3].
The training stage consists in solving the following optimization problem, in the variable θ = (θ E , θ DF , θ D ), after the snapshot matrix S has been formed: where N s = N train N t and with ω h ∈ [0, 1].The per-example loss function (14) combines the reconstruction error (that is, the error between the FOM solution and the DL-ROM approximation) and the error between the intrinsic coordinates and the output of the encoder.The architecture of DL-ROM is the one shown in Fig 1 .The encoder function is used only during the training and validation steps; it is instead discarded during the testing phase.See [16] for further algorithmic details about the training and the testing algorithms required to build and evaluate a DL-ROM.
We highlight that the DL-ROM technique does not require to solve a (reduced) nonlinear dynamical system for the reduced degrees of freedom as in (7); rather, it evaluates a nonlinear map for any given couple (t, µ test ), for each t ∈ (0, T ).Numerical results are extremely accurate, the mean relative error is indeed below 1% (see, e.g., Test 2), even if the causality intrinsic to the parabolic nature of the diffusion-reaction equation providing the monodomain model is broken when computing the DL-ROM approximation.Moreover, the map features an extremely low dimension, in the most favorable scenario equal to n µ + 1.From a computational perspective, remarkable gains and simplifications can be obtained against a linear ROM, since (i) no hyper-reduction is required to enhance the evaluation of any µ-or u-dependent quantity, and (ii) even more interestingly, there is no need to evaluate the dynamics of the recovery variable w if one is only interested in the electrical potential.

Results and discussion
We now assess the computational performances of the proposed DL-ROM strategy on three relevant test cases in cardiac electrophysiology.Our choice of the numerical tests is aimed at highlighting the performance of our DL-ROM method in challenging electrophysiology problems, namely pathological cases in portion of cardiac tissues or physiological scenarios on realistic left ventricle geometries.
The architecture used to perform all the numerical tests is the one reported in the SI Appendix.To solve the optimization problem ( 13)-( 14) we use the ADAM algorithm [37] with a starting learning rate equal to η = 10 −4 .Moreover, we perform crossvalidation by splitting the data in training and validation and following a proportion 8:2 and we implement an early-stopping regularization technique to reduce overfitting [34].
To evaluate the performance of the DL-ROM, we use the loss function (14) and on an error indicator defined as Neural networks required by our DL-ROM technique have been implemented by means of the Tensorflow deep learning framework [38]; numerical simulations have been carried out on a workstation equipped with an Nvidia GeForce GTX 1070 8 GB GPU.
Test 1: Two-dimensional slab with ischemic region We consider the computation of the transmembrane potential in a square slab Ω = (0, 10 cm) 2 of cardiac tissue in presence of an ischemic (non-conductive) region.The ischemic region may act as anatomical driver of cardiac arrhythmias like tachycardias and fibrillations.The system we want to solve is a slight modification of equations (1), accounting for the presence of a non-conductive region which affects both the conductivity tensor and the ionic current term.The ischemic portion of the domain is modeled by replacing the conductivity tensor D(x), defined in (2), with D(x; µ) = σ(x, µ)D(x), where the function σ(x, µ) is given by In this case, n µ = 2 parameters are considered, representing the coordinates of the center of the scar, belong to the parameter space P = [3.5, 6.5 cm] 2 .Moreover, α = 7 cm 2 , σ 0 = 10 −4 , the transversal and longitudinal conductivities are σ t = 12.9 • 0.1 cm 2 /ms and σ l = 12.9 • 0.2 cm 2 /ms, respectively, and f 0 = (1, 0) T , meaning that the tissue fibers are parallel to the x−axis.Similarly, the ionic current I ion (u, w) in ( 1) is replaced by Īion (u, w; µ) = ρ(x; µ)I ion (u, w).The applied current takes the form where C = 100 mA, β = 0.02 cm 2 and t = 2 ms.The parameters appearing in (3) are set to K = 8, a = 0.01, b = 0.15, ε 0 = 0.002, c 1 = 0.2, and c 2 = 0.3, see [39].The equations have been discretized in space through linear finite elements by considering N = 64 × 64 = 4096 grid points.For the time discretization and the treatment of nonlinear terms, we use a one-step, semi-implicit, first order scheme (see [15] for further details) by considering a time step ∆t = 0.1/12.9over (0, T ) with T = 400 ms.
In Figs 2 and 3 we show the FOM and the DL-ROM solutions, the latter obtained with n = 3 for the testing-parameter instance µ test = (6.25,6.25) cm at t = 100 and 356 ms, respectively, together with the relative error ǫ k ∈ R N , for k = 1, . . ., N t , defined as While ( 15) is a synthetic indicator, the quantity defined in ( 17) is instead a function of the space independent variable.In   AP at almost all points; the maximum error is associated to the point P 3 , the closest one to the center of the scar, for t ≥ 200 ms.However, even in this case, the DL-ROM technique is able to capture the difference, in terms of AP values, between the diseased and the healthy tissue.with their proximity impacting on the shape and the values of the AP.In particular, for µ test = 6.25, the point P falls into the grey zone.By using the DL-ROM technique and setting the dimension of the nonlinear trial manifold equal to the dimension of the solution manifold, i.e. n = 3, we obtain an error indicator (15) of ǫ rel = 2.01 •10 −2 .In order to assess the computational efficiency of DL-ROM, we compare it with the POD-Galerkin ROM relying on N c local reduced bases; we report in Table 1 the maximum and minimum number of basis functions, among all the clusters, required by the POD-Galerkin ROM [14,15] to achieve the same accuracy.
In the whole interval (0, T ), by using for each technique the proper time resolution, for a given testing-parameter instance.Since the DL-ROM solution can be queried at a given time without requiring any solution of a dynamical system to recover the former time instances, the DL-ROM can employ larger time windows compared to the time steps required by the solution of the FOM and POD-Galerkin ROM dynamical systems for the cases at hand.This fact also has a positive impact on the data used during the training phase 3 .The speed-up obtained, for each value of N considered, is reported in Table 1.Both the DL-ROM and the POD-Galerkin ROM allow us to decrease the computational costs associated to the computation of the FOM solution for a testingparameter instance.However, for a desired level of accuracy, CPU times required by the POD-Galerkin ROM during the testing phase are remarkably higher than the ones required by a DL-ROM with n = 3.  2, the CPU time required by the DL-ROM at testing time scales linearly with N , instead the one required by the POD-Galerkin ROM scales linearly with √ N .In particular, even for the larger FOM dimension considered (N = 16384 for this test case), our DL-ROM is 19 times faster than the POD-Galerkin ROM.We are not able to run simulations for N > 16384, because of the limitation of the computing resources we have at our disposal.Despite the trend in Fig 7 is apparently not favorable for the DL-ROM technique, practice indicates that the CPU time for DL-ROM is smaller than the one for the POD-Galerkin ROM for small values of N , in other words only with very large values of N the POD-Galerkin ROM outperforms the DL-ROM strategy.Indeed, a linear fitting of the DL-ROM and the POD-Galerkin ROM CPU times4 in Fig 7 highlights that for N = 65536 and N = 262144, DL-ROM could be almost 10 and 5 times, respectively, faster than the POD-Galerkin ROM for the same values of N .Note that the results of this section have been obtained by employing the DL-ROM on a single CPU, an architecture which is not favorable to neural networks 5 .Further improvements are expected when employing our DL-ROM on a GPU for a given testing-parameter instance.DL-ROM and POD-Galerkin ROM vs. FOM speed-up by varying N .The DL-ROM speed-up is remarkably higher than the one obtained by using the POD-Galerkin ROM.
In Figs 8 and 9 we show the feature maps of the first convolutional layer of the encoder function σ , for k = 1, . . ., 8, in the DL-ROM neural network when the FOM solution for the testing-parameter instances µ test = (3.75,3.75) cm and µ test = (6.25,6.25) cm at t = 0.2 ms, are provided as inputs.At this stage, the feature maps retain most of the information present in the FOM solution.Moreover, by considering the two testing-parameter instances, we observe the translation equi-variance property [34] that convolutional layers hold when applied to the part of cardiac tissue corresponding to the scar.Moving to deeper layers, feature maps become increasingly abstract, and less visually interpretable; however, the extracted high-level features are still related both to the ischemic region and the electrical activation pattern.
We highlight that since the DL-ROM solution can be evaluated at any desired time instance without solving any dynamical system, the resulting computational time entailed by the DL-ROM at testing time are drastically reduced compared to the ones required by the FOM or the POD-Galerkin ROM to compute solutions at a particular time instance.In Fig 10 we show the DL-ROM, FOM and POD-Galerkin ROM CPU time needed to compute the approximated solution at t, for t = 1, 10, 100 and 400 ms averaged over the testing set and with N = 4096.We perform the training phase of the POD-Galerkin ROM over the original time interval (0, T ) ms and we report the results for N c = 6, the number of clusters for which the smallest computational time is obtained.The DL-ROM CPU time to compute ũ( t; µ test ) does not vary over t and,  The most recognized cellular mechanisms sustaining atrial tachycardia is re-entry [43].The particular kind of re-entry we deal with in this test case is called figure of eight re-entry, and can be obtained by solving equations (1).To induce the re-entry, we apply a classical S1-S2 protocol [3,44].In particular, we consider a square slab of cardiac tissue Ω = (0, 2 cm) 2 and apply an initial stimulus at the bottom edge of the domain, i.e.
We restrict our study to the time interval [95, 175] ms, i.e. we do not consider the first time instances in which the re-entry has not arisen yet, being them equal over P. The time step is ∆t = 0.2/12.9.We consider N = 256 × 256 = 65536 grid points, implying a mesh size h = 0.0784 mm; this mesh size is recognized to correclty solve the tiny transition front developing during depolarization of the tissue, as highlighted in [41,42].The fibers are parallel to the x-axis and the conductivities in the longitudinal and transversal directions to the fibers are σ l = 2 × 10 −3 cm 2 /ms and σ t = 3.1 × 10 −4 cm 2 /ms, respectively.The parameters appearing in (3) are set to K = 8, a = 0.1, b = 0.1, ε 0 = 0.01, c 1 = 0.14, and c 2 = 0.3, see [45].
The snapshot matrix is built by solving problem (1), completed with the applied currents 18 and 19, by means of linear finite elements and a semi-implicit scheme, over N t = 400 time instances.Moreover, we consider N train = 13 training-parameter instances uniformly distributed in the parameter space and N test = 12 testing-parameter instances, each of them corresponding to the midpoint of two consecutive trainingparameter instances.The maximum number of epochs is set equal to N epochs = 6000, the batch size is N b = 3, due to the high GPU memory occupation of each sample.Regarding the early-stopping criterion, we stop the training if the loss does not decrease in 1000 epochs. In The trend of the relative error (20) over time, for the selected testing-parameter instance µ test = 0.9625 cm, is depicted in Fig 12; we highlight that the error is always   We now compare the computational times required by the FOM, the POD-Galerkin ROM (for different values of N c ) and the DL-ROM, keeping for all the same degree of accuracy achieved by DL-ROM, i.e. ǫ rel = 7.87 × 10 −3 , and running the code on the hardware each implementation is optimized for -a CPU for the FOM and the POD-Galerkin ROM, a GPU6 for the DL-ROM.In Table 3 we report the CPU time needed to compute the FOM solution, and the POD-Galerkin ROM (at the testing phase), both on a full 64 GB node (20 Intel ® Xeon ® E5-2640 v4 2.4GHz cores), and the GPU time required by the DL-ROM to compute 875 time instances (the same number of time instants considered in the solution of the dynamical systems associated to the FOM and the POD-Galerkin ROM) at testing time, by fixing its dimension to n = 5, on an Nvidia GeForce GTX 1070 8 GB GPU.For the sake of completeness, we also report the computational time required by the DL-ROM when employing a single CPU node.It is evident that a POD-Galerkin ROM, built employing a global reduced basis (N c = 1), is not amenable to a complex and challenging pathological cardiac electrophysiology problem like the figure of eight re-entry.Using a nonlinear approach, for which the solution manifold is approximated through a piecewise linear trial manifold (as in the case of N c = 2 or N c = 4 local reduced bases) reduces the online computational time.However, the DL-ROM still confirms to provide a more efficient ROM, almost 5 (or 2) times faster on the CPU, and 39 (or 19) faster on the GPU, than the POD-Galerkin ROM with N c = 2 (or N c = 4) local reduced bases. In

Test 3: Three-dimensional ventricle geometry
We finally consider the solution of the coupled system (1) in a three-dimensional left ventricle (LV) geometry, obtained from the 3D Human Heart Model provided by Zygote [46].Here, we consider a single (n µ = 1) parameter, given by the longitudinal conductivity in the fibers direction.The conductivity tensor takes the form where σ t = 12.9 • 0.02 mm 2 /ms; f 0 is determined at each mesh point through a rulebased approach, by solving a suitable Laplace problem [47].The resulting fibers field is reported in Fig 17 .The applied current is defined as   In order to build the snapshot matrix S, we solve problem (1) completed with the conductivity tensor ( 21) by means of linear finite elements, on a mesh made by N = 16365 vertices, and a semi-implicit scheme in time over a uniform partition of (0, T ) with T = 300 ms and time step ∆t = 0.1/12.9.We uniformly sample N t = 1000 time instances in (0, T ) and we zero-padded [34] the snapshot matrix to reshape each column in a 2D square matrix.The parameter space is provided by P = 12.9 • [0.04, 0.4] mm 2 /ms; here we consider N train = 25 training-parameter instances and N test = 24 testing-parameter instances computed as in Test 2. In this case, the maximum number of epochs is set to N epochs = 30000, the batch size is N b = 40 and the training is stopped if the loss does not decrease over 4000 epochs.
In Fig 18 we report the FOM solution for two testing-parameter instances, i.e. µ = 12.9 • 0.0739 mm 2 /ms and µ = 12.9 • 0.1991 mm 2 /ms, at t = 276 ms, to show the variability of the FOM solution over the parameter space.As expected, front propagation is faster for larger values of the parameter µ.Also for this test case, it is possible to build a reduced nonlinear trial manifold of dimension very close to the intrinsic one -n µ + 1 = 2 -as long as the maximum number of epochs N epochs is increased; the choice n = 10 is obtained as the best trade-off between accuracy and efficiency of the DL-ROM approximation in this case.
The DL-ROM approximation can also replace the FOM solution when evaluating outputs of interest.For instance, in     FOM, DL-ROM and optimal-POD APs for the testing-parameter instance µ test = 12.9 • 0.31 mm 2 /ms (left).For the same n, the DL-ROM is able to provide more accurate results than the optimal-POD.CPU time required to solve the FOM and by DL-ROM at testing time with n = 10 vs N (right).The DL-ROM is able to provide a speed-up equal to 41.
n and n ≪ N , mapping the high-dimensional input x onto a low-dimensional code xn; • the decoder function f D (•; θ D ) : xn → x = f D (xn; θ D ), where f D (•; θ D ) : R n → R N , mapping the low-dimensional code xn to an approximation of the original high-dimensional input x.

Fig 1 .
Fig 1. DL-ROM architecture.DL-ROM architecture used during the training phase.In the red box, the DL-ROM to be queried for any new selected couple (t, µ) during the testing phase.The FOM solution u(t; µ) is provided as input to block (A) which outputs ũn (t; µ).The same parameter instance associated to the FOM, i.e. (t; µ), enters block (B) which provides as output u n (t; µ) and the error between the low-dimensional vectors (dashed green box) is accumulated.The intrinsic coordinates u n (t; µ) are given as input to block (C) returning the ROM approximation ũ(t; µ).Then the reconstruction error (dashed black box) is computed.
Fig 2 the tissue is depolarized except for the region occupied by scar and surrounding it, which is clearly characterized by a slower conduction.In Fig 3 the tissue is starting to repolarize and even if the shape of the ischemic region is not sharply reproduced, the DL-ROM solution is able to capture the diseased (non-conductive) nature of this portion of tissue.In Fig 5 we show the action potentials (APs) computed at the six points P 1 , . . ., P 6 reported in Fig 4. The DL-ROM is able to provide an accurate reconstruction of the

Fig 2 .
Fig 2. Test 1: comparison between FOM and DL-ROM solutions for a testing-parameter instance.FOM solution (left), DL-ROM solution with n = 3 (center) and relative error ǫ k (right) for the testing-parameter instance µ test = (6.25,6.25) cm at t = 100 ms.The maximum of the relative error ǫ k is 10 −3 and it is associated to the diseased tissue.

Fig 3 .
Fig 3. Test 1: comparison between FOM and DL-ROM solutions for a testing-parameter instance.FOM solution (left), DL-ROM solution with n = 3 (center) and relative error ǫ k (right) for the testing-parameter instance µ test = (6.25,6.25) cm at t = 356 ms.The maximum of the relative error ǫ k is 10 −3 and it is associated to the diseased tissue.

Fig 4 .
Fig 4. Test 1: location of points P i .FOM solution evaluated for µ test = (6.25,6.25) cm at t = 400 ms together with the points P 1 , . . ., P 6 .The AP variability across the parameter space characterizing both the FOM and the DL-ROM solutions is shown in Fig 6.Still with a DL-ROM dimension n = 3, we report the APs for µ test = (µ test , µ test ) cm, with µ test = 3.75, 4.25, 4.75, 5.25, 5.75, 6.25, evaluated at P = (7.46,6.51) cm.The DL-ROM is able to capture such variability over P; moreover, the larger µ test , the smaller the distance between the point P and the scar,

Fig 5 .
Fig 5. Test 1: comparison between the FOM and DL-ROM APs at P i .APs evaluated for µ test = (6.25,6.25) cm at points P 1 , . . ., P 6 .The DL-ROM, with n = 3, is able to to sharply reconstruct the AP in almost all the points and the main features are captured also for the points close to the scar.

Fig 6 .
Fig 6.Test 1: variability of the FOM and DL-ROM solutions over the parameter space.FOM (right) and DL-ROM (left) AP variability over P at P = (7.46,6.51) cm.The DL-ROM sharply reconstructs the FOM variability over P.
Fig 7 we compare the CPU time required to solve the FOM (through linear finite elements) over the time interval (0, T ), with the one needed by DL-ROM with n = 3, and the POD-Galerkin ROM with N c = 6 local reduced bases, at testing time, by varying the FOM dimension N .Here, with testing time we refer, both for the DL-ROM and the POD-Galerkin ROM, to the time needed to query the ROM over Table 1.Test 1: dimensions of the POD-Galerkin ROM linear trial manifolds by varying the number of clusters.N c = 1 N c = 2 N c = 4 N c dimensions of the local reduced bases (that is, linear trial manifolds) built by the POD-Galerkin ROM for different numbers N c of clusters.

Fig 7 .
Fig 7. Test 1: FOM, DL-ROM and POD-Galerkin ROM CPU computational times.CPU time required to solve the FOM, by DL-ROM at testing time with n = 3 and by the POD-Galerkin ROM at testing time with N c = 6 vs. N .The DL-ROM provides the smallest testing computational time for each N considered.Both the DL-ROM and the POD-Galerkin ROM depend on the FOM dimension N .In the case of DL-ROM, the dependency on N at testing time, for a fixed value of ∆t, is due to the presence of the decoder function; indeed, the process of learning the reduced dynamics (and so the dimension of the nonlinear trial manifold) does not depend on the FOM dimension.On the other hand, the dependence of the POD-Galerkin ROM on the FOM dimension also impacts on the dimension of the local linear trial manifolds: in general, by increasing N the dimension of each local linear subspace

Fig 8 .
Fig 8. Test 1: activations of the first convolutional layer of the encoder function for a testing-parameter instance.Feature maps of the first convolutional layer of the encoder function in the DL-ROM neural network for the testing-parameter instance µ test = (3.75,3.75) cm at t = 0.2 ms.

Fig 9 .
Fig 9. Test 1: activations of the first convolutional layer of the encoder function for a testing-parameter instance.Activations of the first convolutional layer of the encoder function in the DL-ROM neural network for the testing-parameter instance µ test = (6.25,6.25) cm at t = 0.2 ms.

Fig 10 .
Fig 10.Test 1: FOM, POD-Galerkin ROM and DL-ROM CPU computational times.FOM, POD-Galerkin ROM and DL-ROM CPU computational times to compute ũ( t; µ test ) vs. t averaged over the testing set.Thanks to the fact that the DL-ROM can be queried at any time istance it is extremely efficient in computing ũ( t; µ test ) with respect to both the FOM and the POD-Galerkin ROM.
Fig 11 we show the FOM solution and the DL-ROM one obtained by setting the reduced dimension to n = 5, for the testing-parameter instance µ test = 0.9625 cm, at t = 141.2ms and t = 157.2ms, together with the relative error ǫ s k ∈ R N , for k = 1, . . ., N t , defined as

Fig 11 .
Fig 11.Test 2: comparison between FOM and DL-ROM solutions for a testing-parameter instance.FOM solution (left), DL-ROM one (center) with n = 5, and relative error ǫ s k (right) at t = 141.2ms (top) and t = 157.2ms (bottom), for the testing-parameter instance µ test = 0.9625 cm.The relative error ǫ s k is below the 2% both for t = 141.2ms and t = 157.2ms, the maximum value of the error being associated to very few points of the domain.

Fig 12 .
Fig 12. Test 2: trend of the relative error over time.Relative error ǫ s k vs. t with n = 5 for the testing-parameter instance µ test = 0.9625 cm (the red band indicates the IQR).The distribution of the error maintains uniform over time.

Fig 13 .
Fig 13.Test 2: comparison between FOM and DL-ROM solutions for different testing-parameter instances.FOM solution (left), DL-ROM one (center) with n = 5, and relative error ǫ s k (right) at t = 175 ms, for the testing-parameter instance µ = 0.8125 cm (top) and µ = 1.0625 cm (bottom).The relative error ǫ s k is below the 2.8% both for µ = 0.8125 cm and µ = 1.0625 cm, the maximum value of the error being associated to very few points of the domain.
Fig 14 we show the trend of the error indicator (15) over the testing set versus the CPU computational time both for the DL-ROM and the POD-Galerkin ROM at testing phase.Slight improvements of the performance of DL-ROM, in terms of accuracy, are In both cases, we have considered the setting yielding the most efficient POD-Galerkin ROM approximation, which require about 30 (40, respectively) seconds to be evaluated.By comparing Fig 15 and Fig 11 (bottom), we observe that the DL-ROM outperforms the POD-Galerkin ROM in terms of accuracy.

Fig 19 .
Fig 19.Test 3: comparison between FOM and DL-ROM solutions for a testing-parameter instance at different time instances.FOM solution (left) and DL-ROM one (right), with n = 10, at t = 42.1 ms (top) and t = 276 ms (bottom), for the testing-parameter instance µ test = 12.9 • 0.1435 mm 2 /ms.
Fig 21 and 22 we show the FOM and DL-ROM activation maps, the latter obtained by choosing n = 10 as DL-ROM dimension.Given the electric potential u = u(x, t; µ), the (unipolar) activation map at a point x ∈ Ω is

Fig 20 .
Fig 20.Test 3: comparison between FOM and DL-ROM solutions for a testing-parameter instance at different time instances.FOM solution (left) and DL-ROM one (right), with n = 10, at t = 42.1 ms (top) and t = 276 ms (bottom), for the testing-parameter instance µ test = 12.9 • 0.3243 mm 2 /ms.

Fig 21 .
Fig 21.Test 3: FOM activation map.FOM activation map for the testing-parameter instance µ test = 12.9 • 0.31 mm 2 /ms.In Fig 23 (left) we report the action potentials obtained with the FOM and the DL-ROM, this latter with n = 20, computed at point P = [36.56,1329.59,28.82] mm for the

Fig 23 .
Fig 23.Test 3: FOM, DL-ROM and optimal-POD APs for a testing-parameter instance.FOM and DL-ROM CPU computational times.FOM, DL-ROM and optimal-POD APs for the testing-parameter instance µ test = 12.9 • 0.31 mm 2 /ms (left).For the same n, the DL-ROM is able to provide more accurate results than the optimal-POD.CPU time required to solve the FOM and by DL-ROM at testing time with n = 10 vs N (right).The DL-ROM is able to provide a speed-up equal to 41.