Abstract
Predicting the electrical behavior of the heart, from the cellular scale to the tissue level, relies on the formulation and 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 (RB) 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 inputoutput maps that we intend to reconstruct numerically. To overcome this difficulty, in this paper we propose a new, nonlinear approach which exploits deep learning (DL) algorithms to obtain accurate and efficient ROMs, whose dimensionality matches the number of system parameters. Our DL approach combines deep feedforward neural networks (NNs) and convolutional autoencoders (AEs). We show that the proposed DLROM framework can efficiently provide solutions to parametrized electrophysiology problems, thus enabling multiscenario analysis in pathological cases. We investigate three challenging test cases in cardiac electrophysiology and prove that DLROM outperforms classical projectionbased 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 highfidelity, 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. [4, 5, 3]. 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. Multiquery 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 intersubject 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 projectionbased reduced order models (ROMs) built, e.g., through the reduced basis (RB) method [14], yields inefficient ROMs when dealing with nonlinear timedependent parametrized PDEODE 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 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 hyperreduction stage to reduce the solution of the ODE system to a few, selected mesh nodes [15].
To overcome the limitations of projectionbased ROMs, we propose a new, nonintrusive ROM technique based on deep learning (DL) algorithms, which we refer to as DLROM. Combining in a suitable way a convolutional autoencoder (AE) and a deep feedforward neural network (DFNN), the DLROM 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 DLROM technique has already been presented in [16], albeit on simpler – yet challenging – test cases.
The proposed DLROM technique is a combination of a datadriven 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 nonintrusive 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 DLROM 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 DLROM 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 patientspecific 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 projectionbased ROMs, but still relying on a linear trial manifold built, e.g., through POD. In [21, 22, 23] the construction of ROMs for nonlinear, timedependent problems has been replaced by the evaluation of ANNbased 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 projectionbased ROMs. We demonstrate that our DLROM provides accurate results by constructing ROMs with extremely lowdimension 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.
Materials and methods
Cardiac electrophysiology
Muscle contraction and relaxation drive the pump function of the heart. In particular, tissue contraction is triggered by electrical signals selfgenerated 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 socalled monodomain model for the transmembrane potential with a phenomenological model for the ionic currents – involving a single gating variable – in a domain representing, e.g., a portion of the myocardium, results in the following nonlinear timedependent system
(1) 
Here denotes a rescaled time^{1}^{1}1Dimensional times and potential [27] are given by and . The transmembrane potential ranges from the resting state of mV to the excited state of mV.,
denotes the outward directed unit vector normal to the boundary
of , whereas is an applied current representing, e.g., the initial activation of the tissue. The nonlinear diffusionreaction equation for is twoways coupled with the ODE system, which must be in principle solved at any point ; indeed, the reaction term and the function depend on both and . The most common choices for the two functions and in order to efficiently reproduce the actionpotential are, e.g., the FitzHughNagumo [28, 29], the AlievPanfilov [27, 30] or the Mitchell and Schaeffer models [31]. The diffusivity tensor
usually depends on the fiberssheet 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(2) 
where and are the conductivities in the fibers and the transversal directions.
When a simple phenomenological ionic model is considered, such as the FitzHughNagumo or the AlievPanfilov (AP) model, the ionic current takes the form of a cubic nonlinear function of 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 AlievPanfilov model, which consists in taking
(3)  
The parameters , , , , , are related to the cell. Here represents an oscillation threshold, whereas the weighting factor was introduced in [27] to tune the restitution curve to experimental observations by adjusting the parameters and ; see, e.g., [32, 1, 2, 3] for a detailed review. In the remaining part of the paper, we denote by a parameter vector listing all the input parameters characterizing physical (and, possibly, geometrical) properties we might be interested to vary; is a subset of , denoting the parameter space. Relevant physical situations are those in which input parameters affect the diffusivity matrix (through the conduction velocities) and the applied current ; previous analyses focused instead on the gating variable dynamics (through ) and the ionic current , see [15].
Projectionbased ROMs
From an algebraic standpoint, the spatial discretization of system (1) through the Galerkinfinite element (FE) approximation [33] yields the following nonlinear dynamical system for , , representing our full order model (FOM):
(4) 
Here is a matrix arising from the diffusion operator (thus including the conductivity tensor , which can vary within the myocardium due to fiber orientation and conditions, such as the possible presence of ischemic regions); is the mass matrix; are vectors arising from the nonlinear terms; is a vector collecting the applied currents; finally, are the initial data, possibly depending on . The dimension is related to the dimension of the FE space and, ultimately, depends on the size of the computational mesh used to discretize the domain . Note that the system of ODEs arises from the collocation of the ODE (1) at the nodes used for the numerical integration.
The intrinsic dimension of the solution manifold
(5) 
obtained by solving (4) when varies in , is usually much smaller than and, under suitable conditions, is at most , where 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 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
collecting a set of FOM snapshotsthis is a set of solutions obtained for selected input parameters at (a subset, possibly, of) the time instants in which is partitioned for the sake of time discretization. The most common choice is to set where .
When using a projectionbased ROM, the approximation of is sought as a linear superimposition of modes, under the form
(6) 
thus yielding a linear ROM, in which the columns of the matrix form an orthonormal basis of a space , an dimensional subspace of . In the case of POD, provides the best rank approximation of in the Frobenius norm, that is, are the first (left) singular vectors of corresponding to the largest singular values of , such that the projection error is smaller than a desired tolerance . To meet this requirement, it is sufficient to choose as the smallest integer such that
i.e., the energy retained by the last POD modes is equal or smaller than .
The approximation of is given instead by its restriction
to a (possibly, small) subset of degrees of freedom, where , at which the nonlinear term
is interpolated exploiting a problemdependent basis, spanned by the columns of a matrix
, which is built according to a suitable hyperreduction strategy; see, e.g., [15] for further details. Here denotes a matrix formed by the columns of the identity matrix corresponding to the selected degrees of freedom.A GalerkinPOD ROM for system (1) is then obtained by (i) first, substituting equation (6) into equation 4 and projecting it onto ; then, (ii) solving the system of ODEs at selected degrees of freedom, thus yielding the following nonlinear dynamical system for and the selected components of :
(7) 
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 hyperreduction stage (exploiting, e.g., the discrete empirical interpolation method (DEIM)) must be performed in order to be able to evaluate any  or dependent quantities appearing in (7), that is, without relying on dimensional arrays. Moreover, whenever the solution of the differential problem features coherent structures that propagate over time, such as steep wavefronts, the dimension of the projectionbased 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
clusters, according to suitable clustering (or unsupervised learning) algorithms
[15].Deep learningbased reduced order modeling (DLROM)
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 DLROM technique aims at learning both the nonlinear trial manifold (corresponding to the matrix 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 and the number of training and testing parameter instances, respectively; the ROM dimension is again denoted by . 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
(8) 
where
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 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, (a task which we refer to as reduced trial manifold learning) we employ the decoder function of a convolutional autoencoder^{2}^{2}2The AE is a particular type of neural network aiming at learning the identity function
(10) 
where 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 collecting all the corresponding weights and biases.
As a matter of fact, the approximation provided by the DLROM technique is defined as
(11) 
The encoder function of the convolutional AE can then be exploited to map the FOM solution associated to onto a lowdimensional representation
(12) 
denotes the encoder function, depending upon a vector of parameters.
Computing the DLROM approximation of , for any possible and , corresponds to the testing stage of a DFNN and of the decoder function of a convolutional AE; this does not require the evaluation of the encoder function. We remark that our DLROM strategy overcomes the three major computational bottlenecks implied by the use of projectionbased ROMs, since:

the dimension of the DLROM can be kept extremely small;

the time resolution required by the DLROM can be chosen to be larger than the one required by the numerical solution of dynamical systems in cardiac electrophysiology;

the DLROM can be queried at any desired time instant, without requiring the solution of a dynamical system until that time;

the DLROM does not require to account for the dynamics of the gating variables, thus avoiding any hyperreduction 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 socalled 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 , after the snapshot matrix has been formed:
(13) 
where and
(14) 
with . The perexampleloss function (14) combines the reconstruction error (that is, the error between the FOM solution and the DLROM approximation) and the error between the intrinsic coordinates and the output of the encoder.
The architecture of DLROM 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 DLROM.
We highlight that the DLROM 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 , for each . Numerical results are extremely accurate, the mean relative error is indeed below (see, e.g., Test 2), even if the causality intrinsic to the parabolic nature of the diffusionreaction equation providing the monodomain model is broken when computing the DLROM approximation. Moreover, the map features an extremely low dimension, in the most favorable scenario equal to . From a computational perspective, remarkable gains and simplifications can be obtained against a linear ROM, since (i) no hyperreduction is required to enhance the evaluation of any  or dependent quantity, and (ii) even more interestingly, there is no need to evaluate the dynamics of the recovery variable if one is only interested in the electrical potential.
Results and discussion
We now assess the computational performances of the proposed DLROM strategy on three relevant test cases in cardiac electrophysiology. Our choice of the numerical tests is aimed at highlighting the performance of our DLROM 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 . Moreover, we perform crossvalidation by splitting the data in training and validation and following a proportion 8:2 and we implement an earlystopping regularization technique to reduce overfitting [34].
To evaluate the performance of the DLROM, we use the loss function (14) and on an error indicator defined as
(15) 
Neural networks required by our DLROM 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: Twodimensional slab with ischemic region
We consider the computation of the transmembrane potential in a square slab of cardiac tissue in presence of an ischemic (nonconductive) 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 nonconductive 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 , defined in (2), with , where the function is given by
(16)  
In this case, parameters are considered, representing the coordinates of the center of the scar, belong to the parameter space . Moreover, cm^{2}, , the transversal and longitudinal conductivities are cm^{2}/ms and cm^{2}/ms, respectively, and , meaning that the tissue fibers are parallel to the axis. Similarly, the ionic current in (1) is replaced by . The applied current takes the form
where mA, cm^{2} and ms. The parameters appearing in (3) are set to , , , , , and , see [39]. The equations have been discretized in space through linear finite elements by considering grid points. For the time discretization and the treatment of nonlinear terms, we use a onestep, semiimplicit, first order scheme (see [15] for further details) by considering a time step over with ms.
For the training phase, we uniformly sample time instances over and consider trainingparameter instances, with ,
. The maximum number of epochs is set equal to
, the batch size is and, regarding the earlystopping criterion, we stop the training if the loss function does not decrease in 500 epochs. For the testing phase, testingparameter instances , , have been considered.In Figs 2 and 3 we show the FOM and the DLROM solutions, the latter obtained with for the testingparameter instance cm at and 356 ms, respectively, together with the relative error , for , defined as
(17) 
While (15) is a synthetic indicator, the quantity defined in (17) is instead a function of the space independent variable. In 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 DLROM solution is able to capture the diseased (nonconductive) nature of this portion of tissue.
In Fig 5 we show the action potentials (APs) computed at the six points reported in Fig 4. The DLROM is able to provide an accurate reconstruction of the AP at almost all points; the maximum error is associated to the point , the closest one to the center of the scar, for ms. However, even in this case, the DLROM technique is able to capture the difference, in terms of AP values, between the diseased and the healthy tissue.
The AP variability across the parameter space characterizing both the FOM and the DLROM solutions is shown in Fig 6. Still with a DLROM dimension , we report the APs for cm, with , evaluated at cm. The DLROM is able to capture such variability over ; moreover, the larger , the smaller the distance between the point and the scar, with their proximity impacting on the shape and the values of the AP. In particular, for , the point falls into the grey zone.
By using the DLROM technique and setting the dimension of the nonlinear trial manifold equal to the dimension of the solution manifold, i.e. , we obtain an error indicator (15) of . In order to assess the computational efficiency of DLROM, we compare it with the PODGalerkin ROM relying on local reduced bases; we report in Table 1 the maximum and minimum number of basis functions, among all the clusters, required by the PODGalerkin ROM [14, 15] to achieve the same accuracy.
250  219  200  193 
107  35  26 
Maximum and minimum dimensions of the local reduced bases (that is, linear trial manifolds) built by the PODGalerkin ROM for different numbers of clusters.
In Fig 7 we compare the CPU time required to solve the FOM (through linear finite elements) over the time interval , with the one needed by DLROM with , and the PODGalerkin ROM with local reduced bases, at testing time, by varying the FOM dimension . Here, with testing time we refer, both for the DLROM and the PODGalerkin ROM, to the time needed to query the ROM over the whole interval , by using for each technique the proper time resolution, for a given testingparameter instance. Since the DLROM solution can be queried at a given time without requiring any solution of a dynamical system to recover the former time instances, the DLROM can employ larger time windows compared to the time steps required by the solution of the FOM and PODGalerkin ROM dynamical systems for the cases at hand. This fact also has a positive impact on the data used during the training phase^{3}^{3}3Indeed, in order to build the snapshot matrix, we uniformly sample time instances of the FOM solution over time steps; for each training parameter instance, only of 4000 snapshots are retained from the FOM solution in the DLROM case, against snapshots in the PODGalerkin ROM case.. The speedup obtained, for each value of considered, is reported in Table 1. Both the DLROM and the PODGalerkin 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 PODGalerkin ROM during the testing phase are remarkably higher than the ones required by a DLROM with .
Both the DLROM and the PODGalerkin ROM depend on the FOM dimension . In the case of DLROM, the dependency on at testing time, for a fixed value of , 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 PODGalerkin ROM on the FOM dimension also impacts on the dimension of the local linear trial manifolds: in general, by increasing the dimension of each local linear subspace also increases. Referring to Fig 7 and Table 2, the CPU time required by the DLROM at testing time scales linearly with , instead the one required by the PODGalerkin ROM scales linearly with . In particular, even for the larger FOM dimension considered ( for this test case), our DLROM is 19 times faster than the PODGalerkin ROM. We are not able to run simulations for , 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 DLROM technique, practice indicates that the CPU time for DLROM is smaller than the one for the PODGalerkin ROM for small values of , in other words only with very large values of the PODGalerkin ROM outperforms the DLROM strategy. Indeed, a linear fitting of the DLROM and the PODGalerkin ROM CPU times^{4}^{4}4 and for this test case represent FOM dimensions corresponding to mesh sizes needed to solve, by means of linear finite elements, the problem on a 3D slab geometry both for physiological and pathological electrophysiology in the case a ten TusscherPanfilov ionic model [40] is used. This latter would indeed require smaller values of compared to the AlievPanfilov model, due to the shape of the AP. See, e.g., [41, 42] for further details. in Fig 7 highlights that for and , DLROM could be almost 10 and 5 times, respectively, faster than the PODGalerkin ROM for the same values of . Note that the results of this section have been obtained by employing the DLROM on a single CPU, an architecture which is not favorable to neural networks^{5}^{5}5Indeed, all tests are performed on a node (20 Intel^{®} Xeon^{®} E52640 v4 2.4GHz cores), using 5 cores, of our inhouse HPC cluster.. Further improvements are expected when employing our DLROM on a GPU for a given testingparameter instance.
FOM vs. DLROM  472  536  539  412 

FOM vs. PODGalerkin ROM  3  6  12  22 
DLROM and PODGalerkin ROM vs. FOM speedup by varying . The DLROM speedup is remarkably higher than the one obtained by using the PODGalerkin ROM.
In Figs 8 and 9 we show the feature maps of the first convolutional layer of the encoder function , for , in the DLROM neural network when the FOM solution for the testingparameter instances cm and cm at
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 testingparameter instances, we observe the translation equivariance 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 highlevel features are still related both to the ischemic region and the electrical activation pattern.We highlight that since the DLROM solution can be evaluated at any desired time instance without solving any dynamical system, the resulting computational time entailed by the DLROM at testing time are drastically reduced compared to the ones required by the FOM or the PODGalerkin ROM to compute solutions at a particular time instance. In Fig 10 we show the DLROM, FOM and PODGalerkin ROM CPU time needed to compute the approximated solution at , for 1, 10, 100 and 400 ms averaged over the testing set and with . We perform the training phase of the PODGalerkin ROM over the original time interval ms and we report the results for , the number of clusters for which the smallest computational time is obtained. The DLROM CPU time to compute does not vary over and, by choosing , the DLROM speedups are equal to and with respect to the FOM and the PODGalerkin ROM, with , computational times.
Test 2: Twodimensional slab with figure of eight reentry
The most recognized cellular mechanisms sustaining atrial tachycardia is reentry [43]. The particular kind of reentry we deal with in this test case is called figure of eight reentry, and can be obtained by solving equations (1). To induce the reentry, we apply a classical S1S2 protocol [44, 3]. In particular, we consider a square slab of cardiac tissue and apply an initial stimulus at the bottom edge of the domain, i.e.
(18) 
where , ms and ms.
A second stimulus under the form
(19) 
with , ms and ms, is then applied. The parameter , consisting in the coordinate of the center of the second circular stimulus, ranges in the parameter space cm. This choice has been made to obtain a reentry elicited and sustained until ms.
We restrict our study to the time interval [95, 175] ms, i.e. we do not consider the first time instances in which the reentry has not arisen yet, being them equal over . The time step is . We consider grid points, implying a mesh size 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 axis and the conductivities in the longitudinal and transversal directions to the fibers are cm^{2}/ms and cm^{2}/ms, respectively. The parameters appearing in (3) are set to , , , , , and , 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 semiimplicit scheme, over time instances. Moreover, we consider
trainingparameter instances uniformly distributed in the parameter space and
testingparameter instances, each of them corresponding to the midpoint of two consecutive trainingparameter instances. The maximum number of epochs is set equal to , the batch size is , due to the high GPU memory occupation of each sample. Regarding the earlystopping criterion, we stop the training if the loss does not decrease in 1000 epochs.In Fig 11 we show the FOM solution and the DLROM one obtained by setting the reduced dimension to , for the testingparameter instance cm, at ms and ms, together with the relative error , for , defined as
(20) 
The trend of the relative error (20) over time, for the selected testingparameter instance cm, is depicted in Fig 12; we highlight that the error is always smaller than 1%. In particular, in Fig 12
we show the mean (over the domain), the median, and the first and third quartile of the relative error, as well as its minimum. The interquartile range (IQR) shows that the distribution of the error is almost uniform over time, and that the maximum error is associated to the first time instant – this latter being the time instant at which the solution is most different over
.In Fig 13 we show the FOM and the DLROM solutions, the latter obtained by setting , for the last time instance, i.e. at ms, for cm and cm, in order to point out the variability of the solution over the parameter space and the ability of DLROM to capture it.
We now compare the computational times required by the FOM, the PODGalerkin ROM (for different values of ) and the DLROM, keeping for all the same degree of accuracy achieved by DLROM, i.e. , and running the code on the hardware each implementation is optimized for – a CPU for the FOM and the PODGalerkin ROM, a GPU^{6}^{6}6Indeed, at each layer of a neural network thousands of identical computations must be performed. The most suitable hardware architectures to carry out this kind of operations are GPUs because (i) they have more computational units (cores) and (ii) they have a higher bandwidth to retrieve from memory. Moreover, in applications requiring image processing, as CNNs, the graphics specific capabilities can be further exploited to speed up calculations. for the DLROM. In Table 3 we report the CPU time needed to compute the FOM solution, and the PODGalerkin ROM (at the testing phase), both on a full 64 GB node (20 Intel^{®} Xeon^{®} E52640 v4 2.4GHz cores), and the GPU time required by the DLROM 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 PODGalerkin ROM) at testing time, by fixing its dimension to , on an Nvidia GeForce GTX 1070 8 GB GPU. For the sake of completeness, we also report the computational time required by the DLROM when employing a single CPU node. It is evident that a PODGalerkin ROM, built employing a global reduced basis , is not amenable to a complex and challenging pathological cardiac electrophysiology problem like the figure of eight reentry. Using a nonlinear approach, for which the solution manifold is approximated through a piecewise linear trial manifold (as in the case of or local reduced bases) reduces the online computational time. However, the DLROM 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 PODGalerkin ROM with (or ) local reduced bases.
time [s]  FOM/ROM dimensions  

FOM (CPU)  382  
DLROM (CPU/GPU)  15/1.2  
PODGalerkin ROM (CPU)  103  
PODGalerkin ROM (CPU)  70  
PODGalerkin ROM (CPU)  33 
PODGalerkin ROM and DLROM computational times along with FOM and reduced trial manifold(s) dimensions. DLROM provides a more efficient ROM with respect the PODGalerkin ROMs.
In Fig 14 we show the trend of the error indicator (15) over the testing set versus the CPU computational time both for the DLROM and the PODGalerkin ROM at testing phase. Slight improvements of the performance of DLROM, in terms of accuracy, are obtained for a small increase of the DLROM dimension , coherently with our previous findings reported in [16]. Indeed, the DLROM is able, also in this case, to accurately represent the solution manifold by a reduced nonlinear trial manifold of dimension ; for the case at hand, we report the results for (very close to the intrinsic dimension of the problem, and much smaller than the PODGalerkin ROM dimension), providing slightly smaller values of the error indicator (15) than in the case . Regarding instead the PODGalerkin ROM, in Fig 14 we report results obtained for different tolerances , , , , . In the cases and we only report the results related to the smallest POD tolerances, which indeed allow us to meet the prescribed accuracy on the approximation of the gating variable, which would otherwise impact dramatically on the overall accuracy of the PODGalerkin ROM. Moreover, we do not consider more than local reduced bases in order not to generate too small local linear subspaces. As shown in Fig 14, the proposed DLROM outperforms the PODGalerkin ROM in terms of both efficiency and accuracy.
In Fig 15 we show the solutions obtained through the PODGalerkin ROM with (top) and (bottom) local reduced bases, along with the relative error defined in (20), for the testingparameter instance cm at ms. In both cases, we have considered the setting yielding the most efficient PODGalerkin ROM approximation, which require about 30 (40, respectively) seconds to be evaluated. By comparing Fig 15 and Fig 11 (bottom), we observe that the DLROM outperforms the PODGalerkin ROM in terms of accuracy.
In Fig 16 we show the action potentials obtained through the FOM, the DLROM and the PODGalerkin ROM (with local reduced bases), for the testingparameter instance cm, and evaluated at cm and cm. These two points are close to the left core of the figure of eight reentry, where a shorter action potential duration, and lower values of AP due to the meandering of the cores, are observed. The AP dynamics at those points is accurately captured by the DLROM, while the PODGalerkin ROM might fail in this respect.
Test 3: Threedimensional ventricle geometry
We finally consider the solution of the coupled system (1) in a threedimensional left ventricle (LV) geometry, obtained from the 3D Human Heart Model provided by Zygote [46]. Here, we consider a single () parameter, given by the longitudinal conductivity in the fibers direction. The conductivity tensor takes the form
(21) 
where mm^{2}/ms; 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
where ms, mA, , mm^{2}, mm.
In order to build the snapshot matrix , we solve problem (1) completed with the conductivity tensor (21) by means of linear finite elements, on a mesh made by vertices, and a semiimplicit scheme in time over a uniform partition of with ms and time step . We uniformly sample time instances in
and we zeropadded
[34] the snapshot matrix to reshape each column in a 2D square matrix. The parameter space is provided by mm^{2}/ms; here we consider trainingparameter instances and testingparameter instances computed as in Test 2. In this case, the maximum number of epochs is set to , the batch size is and the training is stopped if the loss does not decrease over 4000 epochs.In Fig 18 we report the FOM solution for two testingparameter instances, i.e. mm^{2}/ms and mm^{2}/ms, at 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 .
In Fig 1920 we report the FOM and DLROM solutions, the latter with , at ms and ms, for two testingparameter instances, mm^{2}/ms and mm^{2}/ms. The DLROM approximation is essentially as accurate as the FOM solution.
Also for this test case, it is possible to build a reduced nonlinear trial manifold of dimension very close to the intrinsic one – – as long as the maximum number of epochs is increased; the choice is obtained as the best tradeoff between accuracy and efficiency of the DLROM approximation in this case.
The DLROM approximation can also replace the FOM solution when evaluating outputs of interest. For instance, in Fig 21 and 22 we show the FOM and DLROM activation maps, the latter obtained by choosing as DLROM dimension. Given the electric potential , the (unipolar) activation map at a point is evaluated as the minimum time at which the AP peak reaches , that is,
Here we compare the activation maps and obtained through the FOM and the DLROM, respectively, by evaluating the maximum of the relative error
over the mesh points; in the case , the maximum relative error is equal to .
In Fig 23 (left) we report the action potentials obtained with the FOM and the DLROM, this latter with , computed at point mm for the testingparameter instance mm^{2}/ms. Moreover, we also report the best approximation of the FOM action potential over a POD space of same dimension for the sake of comparison. Clearly, in dimension the DLROM approximation is much more accurate than the POD best approximation; to reach the same accuracy (about , measured through the error indicator (15)) achieved by the DLROM with , POD modes would be required.
In Fig 23 (right) we highlight instead the improvements, in terms of efficiency, enabled by the use of the DLROM technique. In particular, we point out the CPU time required to solve the FOM for a testing parameter instance, and the one required by DLROM (of dimension ) at testing time, by using the time resolutioin each solution computation requires and by varying the FOM dimension on a 6core platform^{7}^{7}7Numerical tests have been performed on a MacBook Pro Intel Core i7 6core with 16 GB RAM. the FOM solution with degrees of freedom requires about 40 minutes to be computed, against 57 seconds required by the DLROM approximation, thus implying a speedup almost equal to 41 times.
Conclusion
In this work we have proposed a new efficient reduced order model obtained using deep learning algorithms to boost the solution of parametrized problems in cardiac electrophysiology. Numerical results show that the resulting DLROM technique, formerly introduced in [16], allows one to accurately capture complex wave propagation processes, both in physiological and pathological scenarios.
The proposed DLROM technique provides ROMs that are orders of magnitude more efficient than the ones provided by common linear (projectionbased) ROMs, built for instance through a PODGalerkin reduced basis method, for a prescribed level of accuracy. Through the use of DLROM, it is possible to overcome the main computational bottlenecks shown by PODGalerkin ROMs, when addressing parametrized problems in cardiac electrophysiology. The most critical points related to (i) the linear superimposition of modes which linear ROMs are based on; (ii) the need to account for the gating variables when solving the reduced dynamics, even if not required; and (iii) the necessity to use (very often, expensive) hyperreduction techniques to deal with terms that depend nonlinearly on either the transmembrane potential or the input parameters, are all addressed by the DLROM technique, which finally yields more efficient and accurate approximation than PODGalerkin ROMs. Moreover, larger time resolutions can be employed when using a DLROM, compared to the ones required by the numerical solution of a dynamical systems through a FOM or a PODGalerkin ROM. Indeed, the DLROM approximation can be queried at any desired time, without requiring to solve a dynamical system until that time, thus drastically decreasing the computational time required to compute the approximated solution at any given time.
Outputs of clinical interest, such as activation maps and action potentials, can be more efficiently evaluated by the DLROM technique than by a FOM built through the finite element method, while maintaining a high level of accuracy. This work is a proofofconcept of the DLROM technique ability to investigate intra and inter subjects variability, towards performing multiscenario analyses in real time and, ultimately, supporting decisions in clinical practice. In this respect, the use of DLROM techniques can foster assimilation of clinical data with physicsdriven computational models.
Supporting information
SI Code.
Code and data. The code used in this work can be downloaded from: https://github.com/stefaniafresca/DLROM. The training and testing datasets will be made available upon request to the authors.
SI Appendix.
DLROM neural network architecture.
Here we report the configuration of the DLROM neural network used for our numerical tests. We employ a 12layers DFNN equipped with 50 neurons per hidden layer and
neurons in the output layer, where corresponds to the dimension of the reduced nonlinear trial manifold. The architectures of the encoder and decoder functions are instead reported in Table 4 and 5.Acknowledgments
The authors have been partially supported by the ERC Advanced Grant iHEART, “An integrated heart model for the simulation of the cardiac function”, 20172022, P.I. A. Quarteroni (ERC2016AdG, project ID: 740132). Moreover, the authors acknowledge Dr. S. Pagani (MOX, Politecnico di Milano) for the kind help in the FOM software development, and Dr. M. Fedele for providing us the computational mesh of the Zygote Solid 3D heart model.
References
 1. Quarteroni A, Manzoni A, Vergara C. The cardiovascular system: Mathematical modeling, numerical algorithms, clinical applications. Acta Numerica. 2017;26:365–590.
 2. Quarteroni A, Dedè L, Manzoni A, Vergara C. Mathematical modelling of the human cardiovascular system: Data, numerical approximation, clinical applications. Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press; 2019.
 3. Colli Franzone P, Pavarino LF, Scacchi S. Mathematical cardiac electrophysiology. vol. 13 of MS&A. Springer; 2014.
 4. Sundnes J, Lines GT, Cai X, Nielsen BF, Mardal KA, Tveito A. Computing the electrical activity in the heart. vol. 1. Springer Science & Business Media; 2007.
 5. Colli Franzone P, Pavarino LF. A parallel solver for reaction–diffusion systems in computational electrocardiology. Mathematical Models and Methods in Applied Sciences. 2004;14(06):883–911.
 6. Mirams GR, Pathmanathan P, Gray RA, Challenor P, Clayton RH. Uncertainty and variability in computational and mathematical models of cardiac physiology. The Journal of Physiology. 2016;594.23:6833–6847.
 7. Johnstone RH, Chang ETY, Bardenet R, de Boer TP, Gavaghan DJ, Pathmanathan P, et al. Uncertainty and variability in models of the cardiac action potential: Can we build trustworthy models? Journal of Molecular and Cellular Cardiology. 2016;96:49–62.
 8. Hurtado DE, Castro S, Madrid P. Uncertainty quantification of two models of cardiac electromechanics. International Journal for Numerical Methods in Biomedical Engineering. 2017; 33(12):e2894.

9.
Dhamala J, Arevalo HJ, Sapp J, Horácek BM, Wu KC, Trayanova NA, et al.
Quantifying the uncertainty in model parameters using Gaussian processbased Markov chain Monte Carlo in cardiac electrophysiology.
Medical Image Analysis. 2018;48:43–57.  10. Quaglino A, Pezzuto S, Koutsourelakis PS, Auricchio A, Krause R. Fast uncertainty quantification of activation sequences in patientspecific cardiac electrophysiology meeting clinical time constraints. International Journal for Numerical Methods in Biomedical Engineering. 2018;34(7):e2985.
 11. Johnston BM, Coveney S, Chang ET, Johnston PR, Clayton RH. Quantifying the effect of uncertainty in input parameters in a simplified bidomain model of partial thickness ischaemia. Medical & Biological Engineering & Computing. 2018;56(5):761–780.
 12. Pathmanathan P, Cordeiro JM, Gray RA. Comprehensive uncertainty quantification and sensitivity analysis for cardiac action potential models. Frontiers in Physiology. 2019;10.
 13. LevreroFlorencio F, Margara F, Zacur E, BuenoOrovio A, Wang Z, Santiago A, et al. Sensitivity analysis of a stronglycoupled humanbased electromechanical cardiac model: Effect of mechanical parameters on physiologically relevant biomarkers. Computer Methods in Applied Mechanics and Engineering. 2020;361:112762.
 14. Quarteroni A, Manzoni A, Negri F. Reduced basis methods for partial differential equations: An introduction. vol. 92. Springer; 2016.
 15. Pagani S, Manzoni A, Quarteroni A. Numerical approximation of parametrized problems in cardiac electrophysiology by a local reduced basis method. Computer Methods in Applied Mechanics and Engineering. 2018;340:530–558.
 16. Fresca S, Dedé L, Manzoni A. A comprehensive deep learningbased approach to reduced order modeling of nonlinear timedependent parametrized PDEs. arXiv preprint arXiv:200104001. 2020.
 17. Guo M, Hesthaven JS. Reduced order modeling for nonlinear structural analysis using Gaussian process regression. Computer Methods in Applied Mechanics and Engineering. 2018;341:807–826.
 18. Guo M, Hesthaven JS. Datadriven reduced order modeling for timedependent problems. Computer Methods in Applied Mechanics and Engineering. 2019;345:75–99.
 19. Hesthaven J, Ubbiali S. Nonintrusive reduced order modeling of nonlinear problems using neural networks. Journal of Computational Physics. 2018;363:55–78.
 20. San O, Maulik R. Neural network closures for nonlinear model order reduction. Advances in Computational Mathematics. 2018;44(6):1717–1750.

21.
Kani JN, Elsheikh AH.
DRRNN: A deep residual recurrent neural network for model reduction.
arXiv preprint arXiv:170900939. 2017.  22. Mohan A, Gaitonde DV. A deep learning based approach to reduced order modeling for turbulent flow control using LSTM neural networks. arXiv preprint arXiv:18040926. 2018.
 23. Wan Z, Vlachas P, Koumoutsakos P, Sapsis T. Dataassisted reducedorder modeling of extreme events in complex dynamical systems. PLoS one. 2018;13.
 24. González FJ, Balajewicz M. Deep convolutional recurrent autoencoders for learning lowdimensional feature dynamics of fluid systems. arXiv preprint arXiv:180801346. 2018.
 25. Lee K, Carlberg K. Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders. arXiv preprint arXiv:181208373. 2018.
 26. Klabunde R. Cardiovascular Physiology Concepts. Lippincott Williams & Wilkins; 2011.
 27. Aliev RR, Panfilov AV. A simple twovariable model of cardiac excitation. Chaos Solitons Fractals. 1996;7(3):293–301.
 28. FitzHugh R. Impulses and physiological states in theoretical models of nerve membrane. Biophysical Journal. 1961;1(6):445–466.
 29. Nagumo J, Arimoto S, Yoshizawa S. An active pulse transmission line simulating nerve axon. Proceedings of the IRE. 1962;50(10):2061–2070.
 30. Nash MP, Panfilov AV. Electromechanical model of excitable tissue to study reentrant cardiac arrhythmias. Progress in Biophysics and Molecular Biology. 2004;85:501–522.
 31. Mitchell CC, Schaeffer DG. A twocurrent model for the dynamics of cardiac membrane. Bulletin of Mathematical Biology. 2003;65(5):767–793.
 32. Clayton RH, Bernus O, Cherry EM, Dierckx H, Fenton FH, Mirabella L, et al. Models of cardiac tissue electrophysiology: Progress, challenges and open questions. Progress in Biophysics and Molecular Biology. 2011;104(1):22–48.
 33. Quarteroni A, Valli A. Numerical approximation of partial differential equations. vol. 23. Springer; 1994.
 34. Goodfellow I, Bengio Y, Courville A. Deep Learning. MIT Press; 2016. Available from: http://www.deeplearningbook.org.
 35. LeCun Y, Bottou L, Bengio Y, Haffner P. Gradient based learning applied to document recognition. Proceedings of the IEEE. 1998; p. 533–536.
 36. Hinton GE, Zemel RS. Autoencoders, minimum description length, and Helmholtz free energy. Proceedings of the 6th International Conference on Neural Information Processing Systems (NIPS’1993). 1994;.
 37. Kingma DP, Ba J. Adam: A method for stochastic optimization. International Conference on Learning Representations (ICLR); 2015.

38.
Abadi M, Barham P, Chen J, Chen Z, Davis A, Dean J, et al.
TensorFlow: A system for largescale machine
Comments
There are no comments yet.