Deep learning-based reduced order models in cardiac electrophysiology
Stefania Fresca, Andrea Manzoni, Luca Dedè, Alfio Quarteroni
aa r X i v : . [ phy s i c s . c o m p - ph ] J un Deep learning-based reduced order models in cardiacelectrophysiology
Stefania Fresca , Andrea Manzoni , Luca Ded´e , Alfio Quarteroni , MOX - Dipartimento di Matematica, Politecnico di Milano, Milano, Italy Mathematics Institute, ´Ecole Polytechnique F´ed´erale de Lausanne, Lausanne, Switzerland* [email protected]
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 dynamicalsystems. These systems describe the cardiac action potential, that is the polariza-tion/depolarization cycle occurring at every heart beat that models the time evolutionof 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 requiredto evaluate outputs of clinical interest, such as activation maps and action potential du-ration. More importantly, these models feature coherent structures that propagate overtime, such as wavefronts. These systems can hardly be reduced to lower dimensionalproblems 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 (withrespect to the problem parameters) as well as to the nonlinear nature of the input-outputmaps that we intend to reconstruct numerically. To overcome this difficulty, in this pa-per we propose a new, nonlinear approach which exploits deep learning (DL) algorithmsto obtain accurate and efficient ROMs, whose dimensionality matches the number of sys-tem parameters. Our DL approach combines deep feedforward neural networks (NNs)and convolutional autoencoders (AEs). We show that the proposed DL-ROM frame-work can efficiently provide solutions to parametrized electrophysiology problems, thusenabling multi-scenario analysis in pathological cases. We investigate three challengingtest cases in cardiac electrophysiology and prove that DL-ROM outperforms classicalprojection-based ROMs.
Introduction
The electrical activation of the heart is the main responsible of its contraction, is the re-sult of two processes: at the microscopic scale, the generation of ionic currents throughthe 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 transmembranepotential [1–3]. This latter process can be described by means of partial differential equa-tions (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., thefinite element (FE) method, is computationally demanding. Indeed, the propagation ofthe electrical signal is characterized by the fast dynamics of very steep fronts, thus requir-ing very fine space and time discretizations. [3–5]. Using a FOM may quickly becomeunaffordable if such a coupled system must be solved for several values of parametersJune 5, 2020 1/28epresenting either functional or geometric data such as, e.g., material properties, initialand boundary conditions, or the shape of the domain. Multi-query analysis is relevantin a variety of situations: when analysing multiple scenarios, when dealing with sen-sitivity analysis and uncertainty quantification (UQ) problems in order to account forinter-subject variability [6–8], for parameter estimation or data assimilation, in whichsome unknown (or unaccessible) quantities characterizing the mathematical model mustbe inferred from a set of measurements [9–13].Conventional projection-based reduced order models (ROMs) built, e.g., through thereduced basis (RB) method [14], yields inefficient ROMs when dealing with nonlineartime-dependent parametrized PDE-ODE system as the one arising from cardiac electro-physiology. The three major computational bottlenecks shown by such kind of ROMsfor cardiac electrophysiology are due the fact that:- the linear superimposition of modes, on which they are based, would cause thedimension of the ROM to be excessively large to guarantee an acceptable accuracy;- evaluating the ROM requires the solution of a dynamical system, which might beunstable unless the size of time step ∆ t is very small;- the ROM must also account for the dynamics of the gating variables, even whenaiming at computing just the electrical potential. This fact entails an extremelyintrusive and costly hyper-reduction stage to reduce the solution of the ODEsystem to a few, selected mesh nodes [15].To overcome the limitations of projection-based ROMs, we propose a new, non-intrusive ROM technique based on deep learning (DL) algorithms, which we refer to asDL-ROM. Combining in a suitable way a convolutional autoencoder (AE) and a deepfeedforward neural network (DFNN), the DL-ROM technique enables the constructionof an efficient ROM, whose dimension is as close as possible to the number of parametersupon which the solution of the differential problem depends. A preliminary numericalassessment of our DL-ROM technique has already been presented in [16], albeit onsimpler – yet challenging – test cases.The proposed DL-ROM technique is a combination of a data-driven with a physicsbased model approach. Indeed, it exploits snapshots taken from a set of FOM solutions(for selected parameter values and time instances) and deep neural network architecturesto learn, in a non-intrusive way, both (i) the nonlinear trial manifold where the ROMsolution 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 aset of basis functions, while the latter task corresponds to the projection stage in thesubspace spanned by these basis functions. Here, our goal is to show that DL-ROMcan be effectively used to handle parametrized problems in cardiac electrophysiology,accounting for both physiological and pathological conditions, in order to provide fastand accurate solutions. The proposed DL-ROM is computationally efficient during thetesting stage, that is for any new scenario unseen during the training stage. This isparticularly useful in view of the evaluation of patient-specific features to enable theintegration of computational methods in current clinical platforms.DL techniques for parametrized PDEs have previously been proposed in other con-texts. In [17–20] feedforward neural networks have been employed to model the reduceddynamics in a less intrusive way, that is, avoiding the costs entailed by projection-basedROMs, but still relying on a linear trial manifold built, e.g., through POD. In [21–23]the construction of ROMs for nonlinear, time-dependent problems has been replaced bythe evaluation of ANN-based regression models. In [24, 25] the reduced trial manifoldwhere the approximation is sought has been modeled through ANNs thus avoiding thelinear superimposition of POD modes, on a minimum residual formulation to deriveJune 5, 2020 2/28he ROM [25], or without considering an explicit parameter dependence in the differ-ential problem that is considered [24]. In all these works, coupled problems have neverbeen considered. Moreover, very often DL techniques have been exploited to addressproblems which require only a moderate dimension of projection-based ROMs. Wedemonstrate that our DL-ROM provides accurate results by constructing ROMs withextremely low-dimension in prototypical test cases. These tests exhibit all the relevantphysical features which make the numerical approximation of parametrized problems incardiac electrophysiology a challenging task. Materials and methods
Cardiac electrophysiology
Muscle contraction and relaxation drive the pump function of the heart. In particu-lar, tissue contraction is triggered by electrical signals self-generated in the heart andpropagated through the myocardium thanks to the excitability of the cardiac cells, thecardiomyocites [3, 26]. When suitably stimulated, cardiomyocites produce a variationof the potential across the cellular membrane, called transmembrane potential . Its evo-lution in time is usually referred to as action potential , involving a polarization and adepolarization in the early stage of every heart beat. The action potential is generatedby several ion channels (e.g., calcium, sodium, potassium) that open and close, and bythe resulting ionic currents crossing the membrane. For instance, coupling the so-calledmonodomain model for the transmembrane potential u = u ( x , t ) with a phenomenolog-ical model for the ionic currents – involving a single gating variable w = w ( x , t ) – ina domain Ω representing, e.g., a portion of the myocardium, results in the followingnonlinear time-dependent system ∂u∂t − div( D ∇ u ) + I ion ( u, w ) = I app ( x , t ) ( x , t ) ∈ Ω × (0 , T ) ,∂w∂t + g ( u, w ) = 0 ( x , t ) ∈ Ω × (0 , T ) , ∇ u · n = 0 ( x , t ) ∈ ∂ Ω × (0 , T ) ,u ( x ,
0) = 0 , w ( x ,
0) = 0 x ∈ Ω . (1)Here t denotes a rescaled time , n denotes the outward directed unit vector normal tothe boundary ∂ Ω of Ω, whereas I app is an applied current representing, e.g., the initialactivation of the tissue. The nonlinear diffusion-reaction equation for u is two-wayscoupled 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 mostcommon choices for the two functions I ion and g in order to efficiently reproduce theaction-potential are, e.g., the FitzHugh-Nagumo [28, 29], the Aliev-Panfilov [27, 30] orthe Mitchell and Schaeffer models [31]. The diffusivity tensor D usually depends onthe fibers-sheet structure of the tissue, affecting directional conduction velocities anddirections. In particular, by assuming an axisymmetric distribution of the fibers, theconductivity tensor takes the form D ( x ) = σ t I + ( σ l − σ t ) f ⊗ f , (2)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 Dimensional times and potential [27] are given by ˜ t [ ms ] = 12 . t and ˜ u [ mV ] = 100 u −
80. Thetransmembrane potential ranges from the resting state of −
80 mV to the excited state of +20 mV.
June 5, 2020 3/28onlinear function of u and a single (dimensionless) gating variable plays the role of arecovery function, allowing to model refractariness of cells. In this paper, we focus onthe Aliev-Panfilov model, which consists in taking I ion ( u, w ) = Ku ( u − a )( u −
1) + uw,g ( u, w ) = (cid:16) ǫ + c wc + u (cid:17) ( − w − Ku ( u − b − . (3)The parameters K , a , b , ε , c , c are related to the cell. Here a represents an oscillationthreshold, whereas the weighting factor ε + c wc + u was introduced in [27] to tune therestitution curve to experimental observations by adjusting the parameters c and c ;see, e.g., [1–3, 32] for a detailed review. In the remaining part of the paper, we denoteby µ ∈ P ⊂ R n µ a parameter vector listing all the n µ input parameters characterizingphysical (and, possibly, geometrical) properties we might be interested to vary; P is asubset of R n µ , denoting the parameter space. Relevant physical situations are thosein which input parameters affect the diffusivity matrix D (through the conductionvelocities) and the applied current I app ; previous analyses focused instead on the gatingvariable 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 theGalerkin-finite element (FE) approximation [33] yields the following nonlinear dynam-ical system for u = u ( t ; µ ), w = w ( t ; µ ), representing our full order model (FOM): M ( µ ) ∂ u ∂t = A ( µ ) u + I ion ( t, u , w ; µ ) + I app ( t ; µ ) , t ∈ (0 , T ) ,∂ w ∂t ( t ; µ ) = g ( t, u , w ; µ ) , t ∈ (0 , T ) , u (0) = u , w (0) = w . (4)Here A ( µ ) ∈ R N × N is a matrix arising from the diffusion operator (thus includingthe conductivity tensor D ( µ ) = D ( x ; µ ), which can vary within the myocardium dueto 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 nonlinearterms; I app ∈ R N is a vector collecting the applied currents; finally, u , w ∈ R N are theinitial data, possibly depending on µ . The dimension N is related to the dimension ofthe FE space and, ultimately, depends on the size h > at the nodes used for the numerical integration.The intrinsic dimension of the solution manifold S = { u ( t ; µ ) | t ∈ [0 , T ) and µ ∈ P ⊂ R n µ } ⊂ R N , (5)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 ofparameters – in this respect, the time independent variable plays the role of a parame-ter. For this reason, ROMs attempt at approximating S by introducing a suitable trialmanifold of lower dimension. The most popular approach is proper orthogonal decom-position (POD), which exploits a linear trial manifold built through the singular valuedecomposition of a matrix S ∈ R N × N s collecting a set of FOM snapshots S = (cid:2) u ( t ; µ ) | . . . | u ( t N t ; µ ) | . . . | u ( t ; µ N train ) | . . . | u ( t N t ; µ N train ) (cid:3) ;this is a set of solutions obtained for N train selected input parameters at (a subset,possibly, of) the time instants { t k } N t k =1 in which (0 , T ) is partitioned for the sake of timediscretization. The most common choice is to set t k = k ∆ t where ∆ t = T / ( N t − u ( t ; µ ) is sought as alinear superimposition of modes, under the form u ( t ; µ ) ≈ Vu n ( t ; µ ) , (6)thus yielding a linear ROM, in which the columns of the matrix V = [ ζ , . . . , ζ n ] ∈ R N × n form an orthonormal basis of a space V n , an n -dimensional subspace of R N . Inthe case of POD, V n provides the best n -rank approximation of S in the Frobeniusnorm, that is, ζ , . . . , ζ n are the first n (left) singular vectors of S corresponding to the n largest singular values σ , . . . , σ n of S , such that the projection error is smaller thana desired tolerance ε P OD . To meet this requirement, it is sufficient to choose n as thesmallest integer such that P Ni =1 σ i P N s i =1 σ i > − ε P OD , i.e., the energy retained by the last N s − n POD modes is equal or smaller than ε 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 thenonlinear term I ion is interpolated exploiting a problem-dependent basis, spanned bythe columns of a matrix Φ ∈ R N × m , which is built according to a suitable hyper-reduction strategy; see, e.g., [15] for further details. Here P = [ e , . . . , e m ] ∈ R N × m denotes a matrix formed by the columns of the N × N identity matrix correspondingto the m selected degrees of freedom.A Galerkin-POD ROM for system (1) is then obtained by (i) first, substitutingequation (6) into equation 4 and projecting it onto V n ; then, (ii) solving the system ofODEs at m selected degrees of freedom, thus yielding the following nonlinear dynamicalsystem for u n = u n ( t ; µ ) and the selected components P T w = P T w ( t ; µ ) of w : V T M ( µ ) V ∂ u n ∂t + V T A ( µ ) V T u n + V T Φ ( P T Φ ) − I ion ( t, P T Vu n , P T w ; µ ) − V T I app ( t ; µ ) = , t ∈ (0 , T ) , P T ∂ w ∂t + g ( t, P T Vu n , P T w ; µ ) = , t ∈ (0 , T ) , u n (0) = V T u , P T w (0) = P T w . (7)This strategy is the essence of the reduced basis (RB) method for nonlinear time-dependent parametrized PDEs. However, using (7) as an approximation to (4) is knownto suffer from several problems. First of all, an extensive hyper-reduction stage (exploit-ing, e.g., the discrete empirical interpolation method (DEIM)) must be performed inorder to be able to evaluate any µ - or u -dependent quantities appearing in (7), thatis, without relying on N -dimensional arrays. Moreover, whenever the solution of thedifferential problem features coherent structures that propagate over time, such as steepwavefronts, the dimension n of the projection-based ROM (7) might easily become verylarge, due to the basic linearity assumption, by which the solution is given by a linearsuperimposition of POD modes, thus severely degrading the computational efficiencyof the ROM. A possible way to overcome this bottleneck is to rely on local reducedbases, built through POD after the set of snapshots has been split into N c > eep learning-based reduced order modeling (DL-ROM) To overcome the limitations of linear ROMs, we consider a new, nonlinear ROM tech-nique based on deep learning models. First introduced in [16] and assessed on one-dimensional benchmark problems, the DL-ROM technique aims at learning both thenonlinear 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 reduceddynamics (corresponding to the projection stage in a linear ROM). This method is notintrusive; it relies on DL algorithms trained on a set of FOM solutions obtained fordifferent parameter values.We denote by N train and N test the number of training and testing parameter in-stances, respectively; the ROM dimension is again denoted by n ≪ N . In order todescribe the system dynamics on a suitable reduced nonlinear trial manifold (a taskwhich we refer to as reduced dynamics learning ), the intrinsic coordinates of the ROMapproximation are defined as u n ( t ; µ , θ DF ) = φ DFn ( t ; µ , θ DF ) , (8)where φ DFn ( · ; · , θ DF ) : R ( n µ +1) → R n is a deep feedforward neural network (DFNN),consisting in the subsequent composition of a nonlinear activation function, applied toa linear transformation of the input, multiple times [34]. Here θ DF denotes the vectorof parameters of the DFNN, collecting all the corresponding weights and biases of eachlayer of the DFNN.Regarding instead the description of the reduced nonlinear trial manifold, approx-imating the solution one, ˜ S ≈ S (a task which we refer to as reduced trial manifoldlearning ) we employ the decoder function of a convolutional autoencoder (AE) [35, 36].More precisely, ˜ S takes the form˜ S = { f D ( u n ( t ; µ , θ DF ); θ D ) | u n ( t ; µ , θ DF ) ∈ R n , t ∈ [0 , T ) and µ ∈ P ⊂ R n µ } (10)where f D ( · ; θ D ) : R n → R N consists in the decoder function of a convolutional AE. Thislatter 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 ˜u ( t ; µ ) ≈ u ( t ; µ ) provided by the DL-ROMtechnique is defined as ˜u ( t ; µ , θ DF , θ D ) = f D ( φ DFn ( t ; µ , θ DF ); θ D ) . (11)The encoder function of the convolutional AE can then be exploited to map the FOMsolution associated to ( t, µ ) onto a low-dimensional representation ˜u n ( t ; µ , θ E ) = f En ( u ( t ; µ ); θ E ); (12) f En ( · ; θ E ) : R N → R n denotes the encoder function, depending upon a vector θ E ofparameters.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 The AE is a particular type of neural network aiming at learning the identity function f AE ( · ; θ E , θ D ) : x ˜x with ˜x ≃ x . (9). It is composed by two main parts: • the encoder function f En ( · ; θ E ) : x ˜x n = f En ( x ; θ E ), where f En ( · ; θ E ) : R N → R n and n ≪ N ,mapping the high-dimensional input x onto a low-dimensional code ˜x n ; • the decoder function f D ( · ; θ D ) : ˜x n ˜x = f D ( ˜x n ; θ D ), where f D ( · ; θ D ) : R n → R N , mappingthe low-dimensional code ˜ x n to an approximation of the original high-dimensional input ˜x . June 5, 2020 6/28onvolutional AE; this does not require the evaluation of the encoder function. We re-mark that our DL-ROM strategy overcomes the three major computational bottlenecksimplied 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 theone required by the numerical solution of dynamical systems in cardiac electro-physiology;- the DL-ROM can be queried at any desired time instant, without requiring thesolution 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 whenemploying 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 inthe system of ODEs must be accounted for [3].The training stage consists in solving the following optimization problem, in thevariable θ = ( θ E , θ DF , θ D ), after the snapshot matrix S has been formed:min θ J ( θ ) = min θ N s N train X i =1 N t X k =1 L ( t k , µ i ; θ ) , (13)where N s = N train N t and L ( t k , µ i ; θ ) = ω h k u ( t k ; µ i ) − ˜u ( t k ; µ i , θ DF , θ D ) k + 1 − ω h k ˜ u n ( t k ; µ i , θ E ) − u n ( t k ; µ i , θ DF ) k , (14)with ω h ∈ [0 , per-example loss function (14) combines the reconstruction error(that is, the error between the FOM solution and the DL-ROM approximation) and theerror 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 isused only during the training and validation steps; it is instead discarded during thetesting phase. See [16] for further algorithmic details about the training and the testingalgorithms 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, itevaluates a nonlinear map for any given couple ( t, µ test ), for each t ∈ (0 , T ). Numericalresults 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-reactionequation providing the monodomain model is broken when computing the DL-ROMapproximation. Moreover, the map features an extremely low dimension, in the mostfavorable scenario equal to n µ + 1. From a computational perspective, remarkable gainsand simplifications can be obtained against a linear ROM, since (i) no hyper-reductionis required to enhance the evaluation of any µ - or u -dependent quantity, and (ii) evenmore interestingly, there is no need to evaluate the dynamics of the recovery variable w if one is only interested in the electrical potential.June 5, 2020 7/28 ig 1. DL-ROM architecture. DL-ROM architecture used during the trainingphase. 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 ˜ u 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 thelow-dimensional vectors (dashed green box) is accumulated. The intrinsic coordinates u n ( t ; µ ) are given as input to block (C) returning the ROM approximation ˜u ( t ; µ ).Then the reconstruction error (dashed black box) is computed. Results and discussion
We now assess the computational performances of the proposed DL-ROM strategy onthree relevant test cases in cardiac electrophysiology. Our choice of the numerical testsis aimed at highlighting the performance of our DL-ROM method in challenging elec-trophysiology problems, namely pathological cases in portion of cardiac tissues or phys-iological scenarios on realistic left ventricle geometries.The architecture used to perform all the numerical tests is the one reported in theSI Appendix. To solve the optimization problem (13)-(14) we use the ADAM algorithm[37] with a starting learning rate equal to η = 10 − . Moreover, we perform cross-validation by splitting the data in training and validation and following a proportion 8:2and 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 onan error indicator defined as ǫ rel = 1 N test N test X i =1 qP N t k =1 || u k ( µ test,i ) − ˜u k ( µ test,i ) || qP N t k =1 || u k ( µ test,i ) || . (15)Neural networks required by our DL-ROM technique have been implemented bymeans of the Tensorflow deep learning framework [38]; numerical simulations havebeen carried out on a workstation equipped with an Nvidia GeForce GTX 1070 8 GBGPU.June 5, 2020 8/28 est 1: Two-dimensional slab with ischemic region
We consider the computation of the transmembrane potential in a square slab Ω =(0 ,
10 cm) of cardiac tissue in presence of an ischemic (non-conductive) region. Theischemic region may act as anatomical driver of cardiac arrhythmias like tachycardiasand 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 conduc-tivity tensor and the ionic current term. The ischemic portion of the domain is modeledby replacing the conductivity tensor D ( x ), defined in (2), with ¯D ( x ; µ ) = σ ( x , µ ) D ( x ),where the function σ ( x , µ ) is given by σ ( x ; µ ) = ρ ( x ; µ ) + σ (1 − ρ ( x ; µ )) ,ρ ( x ; µ ) = 1 − exp (cid:18) − ( x − µ ) + ( x − µ ) α (cid:19) . (16)In this case, n µ = 2 parameters are considered, representing the coordinates of the centerof the scar, belong to the parameter space P = [3 . , . . Moreover, α = 7 cm , σ = 10 − , the transversal and longitudinal conductivities are σ t = 12 . · . /msand σ l = 12 . · . /ms, respectively, and f = (1 , T , meaning that the tissue fibersare parallel to the x − axis. Similarly, the ionic current I ion ( u, w ) in (1) is replaced by¯ I ion ( u, w ; µ ) = ρ ( x ; µ ) I ion ( u, w ). The applied current takes the form I app ( x , t ) = C exp (cid:18) − || x || β (cid:19) [0 , ¯ t ] (˜ t ) , where C = 100 mA, β = 0 .
02 cm and ¯ t = 2 ms. The parameters appearing in (3) areset to K = 8, a = 0 . b = 0 . ε = 0 . c = 0 .
2, and c = 0 .
3, see [39]. Theequations have been discretized in space through linear finite elements by considering N = 64 ×
64 = 4096 grid points. For the time discretization and the treatment ofnonlinear terms, we use a one-step, semi-implicit, first order scheme (see [15] for furtherdetails) by considering a time step ∆ t = 0 . / . , T ) with T = 400 ms.For the training phase, we uniformly sample N t = 1000 time instances over (0 , T ) andconsider N train = 49 training-parameter instances, with µ train = (3 . i . , . j . i, j = 0 , . . . ,
6. The maximum number of epochs is set equal to N epochs = 10000, thebatch size is N b = 40 and, regarding the early-stopping criterion, we stop the trainingif the loss function does not decrease in 500 epochs. For the testing phase, N test = 36testing-parameter instances µ test = (3 .
75 + i . , .
75 + j . i, j = 0 , . . . ,
5, have beenconsidered.In Figs 2 and 3 we show the FOM and the DL-ROM solutions, the latter obtainedwith n = 3 for the testing-parameter instance µ test = (6 . , .
25) cm at ˜ t = 100 and 356ms, respectively, together with the relative error ǫ k ∈ R N , for k = 1 , . . . , N t , defined as ǫ k = | u k ( µ test ) − ˜u k ( µ test ) | q N t P N t k =1 || u k ( µ test ) || . (17)While (15) is a synthetic indicator, the quantity defined in (17) is instead a functionof the space independent variable. In Fig 2 the tissue is depolarized except for theregion occupied by scar and surrounding it, which is clearly characterized by a slowerconduction. In Fig 3 the tissue is starting to repolarize and even if the shape of theischemic region is not sharply reproduced, the DL-ROM solution is able to capture thediseased (non-conductive) nature of this portion of tissue.In Fig 5 we show the action potentials (APs) computed at the six points P , . . . , P reported in Fig 4. The DL-ROM is able to provide an accurate reconstruction of theJune 5, 2020 9/28 ig 2. Test 1: comparison between FOM and DL-ROM solutions for atesting-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) cm at ˜ t = 100 ms. The maximum of the relative error ǫ k is 10 − and it is associated to the diseased tissue. Fig 3. Test 1: comparison between FOM and DL-ROM solutions for atesting-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) cm at ˜ t = 356 ms. The maximum of the relative error ǫ k is 10 − and it is associated to the diseased tissue.AP at almost all points; the maximum error is associated to the point P , the closestone to the center of the scar, for ˜ t ≥
200 ms. However, even in this case, the DL-ROMtechnique is able to capture the difference, in terms of AP values, between the diseasedand the healthy tissue.
Fig 4. Test 1: location of points P i . FOM solution evaluated for µ test = (6 . , .
25) cm at ˜ t = 400 ms together with the points P , . . . , P .The AP variability across the parameter space characterizing both the FOM andthe DL-ROM solutions is shown in Fig 6. Still with a DL-ROM dimension n = 3, wereport the APs for µ test = ( µ test , µ test ) cm, with µ test = 3 . , . , . , . , . , . P = (7 . , .
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,June 5, 2020 10/28 ig 5. Test 1: comparison between the FOM and DL-ROM APs at P i . APsevaluated for µ test = (6 . , .
25) cm at points P , . . . , P . The DL-ROM, with n = 3,is able to to sharply reconstruct the AP in almost all the points and the main featuresare captured also for the points close to the scar.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 . Fig 6. Test 1: variability of the FOM and DL-ROM solutions over theparameter space.
FOM (right) and DL-ROM (left) AP variability over P at P = (7 . , .
51) cm. The DL-ROM sharply reconstructs the FOM variability over P .By using the DL-ROM technique and setting the dimension of the nonlinear trialmanifold equal to the dimension of the solution manifold, i.e. n = 3, we obtain an errorindicator (15) of ǫ rel = 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 allthe clusters, required by the POD-Galerkin ROM [14, 15] to achieve the same accuracy.In Fig 7 we compare the CPU time required to solve the FOM (through linearfinite 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 theDL-ROM and the POD-Galerkin ROM, to the time needed to query the ROM overJune 5, 2020 11/28 able 1. Test 1: dimensions of the POD-Galerkin ROM linear trialmanifolds by varying the number of clusters. N c = 1 N c = 2 N c = 4 N c = 6250 219 200 193107 35 26Maximum and minimum dimensions of the local reduced bases (that is, linear trialmanifolds) built by the POD-Galerkin ROM for different numbers N c of clusters.the whole interval (0 , T ), by using for each technique the proper time resolution, for agiven testing-parameter instance. Since the DL-ROM solution can be queried at a giventime without requiring any solution of a dynamical system to recover the former timeinstances, the DL-ROM can employ larger time windows compared to the time stepsrequired by the solution of the FOM and POD-Galerkin ROM dynamical systems forthe cases at hand. This fact also has a positive impact on the data used during thetraining phase . The speed-up obtained, for each value of N considered, is reportedin Table 1. Both the DL-ROM and the POD-Galerkin ROM allow us to decrease thecomputational costs associated to the computation of the FOM solution for a testing-parameter instance. However, for a desired level of accuracy, CPU times required bythe POD-Galerkin ROM during the testing phase are remarkably higher than the onesrequired by a DL-ROM with n = 3. Fig 7. Test 1: FOM, DL-ROM and POD-Galerkin ROM CPUcomputational times.
CPU time required to solve the FOM, by DL-ROM at testingtime 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 valueof ∆ t , is due to the presence of the decoder function; indeed, the process of learningthe reduced dynamics (and so the dimension of the nonlinear trial manifold) does notdepend 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 lineartrial manifolds: in general, by increasing N the dimension of each local linear subspace Indeed, in order to build the snapshot matrix, we uniformly sample N t time instances of the FOMsolution over T/ ∆ t = 4000 time steps; for each training parameter instance, only 25% of 4000 snapshotsare retained from the FOM solution in the DL-ROM case, against 4000 snapshots in the POD-GalerkinROM case. June 5, 2020 12/28lso increases. Referring to Fig 7 and Table 2, the CPU time required by the DL-ROM attesting time scales linearly with N , instead the one required by the POD-Galerkin ROMscales 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-GalerkinROM. We are not able to run simulations for N > small values of N , inother words only with very large values of N the POD-Galerkin ROM outperforms theDL-ROM strategy. Indeed, a linear fitting of the DL-ROM and the POD-Galerkin ROMCPU times in Fig 7 highlights that for N = 65536 and N = 262144, DL-ROM could bealmost 10 and 5 times, respectively, faster than the POD-Galerkin ROM for the samevalues of N . Note that the results of this section have been obtained by employing theDL-ROM on a single CPU, an architecture which is not favorable to neural networks .Further improvements are expected when employing our DL-ROM on a GPU for a giventesting-parameter instance. Table 2. Test 1: DL-ROM and POD-Galerkin ROM vs. FOM speed-up. N = 256 N = 1024 N = 4096 N = 16384FOM vs. DL-ROM 472 536 539 412FOM vs. POD-Galerkin ROM 3 6 12 22DL-ROM and POD-Galerkin ROM vs. FOM speed-up by varying N . The DL-ROMspeed-up is remarkably higher than the one obtained by using the POD-GalerkinROM.In Figs 8 and 9 we show the feature maps of the first convolutional layer of the en-coder function σ ( W k ∗ u ( µ test )+ b k ), for k = 1 , . . . ,
8, in the DL-ROM neural networkwhen the FOM solution for the testing-parameter instances µ test = (3 . , .
75) cm and µ test = (6 . , .
25) cm at t = 0 . solving any dynamical system, the resulting computational time en-tailed by the DL-ROM at testing time are drastically reduced compared to the onesrequired by the FOM or the POD-Galerkin ROM to compute solutions at a particulartime instance. In Fig 10 we show the DL-ROM, FOM and POD-Galerkin ROM CPUtime needed to compute the approximated solution at ¯ t , for ¯ t = 1, 10, 100 and 400ms averaged over the testing set and with N = 4096. We perform the training phaseof the POD-Galerkin ROM over the original time interval (0 , T ) ms and we report theresults for N c = 6, the number of clusters for which the smallest computational timeis obtained. The DL-ROM CPU time to compute ˜ u (¯ t ; µ test ) does not vary over ¯ t and, N = 65536 and N = 262144 for this test case represent FOM dimensions corresponding to meshsizes h needed to solve, by means of linear finite elements, the problem on a 3D slab geometry both forphysiological and pathological electrophysiology in the case a ten Tusscher-Panfilov ionic model [40] isused. This latter would indeed require smaller values of h compared to the Aliev-Panfilov model, dueto the shape of the AP. See, e.g., [41, 42] for further details. Indeed, all tests are performed on a node (20 Intel ® Xeon ® E5-2640 v4 2.4GHz cores), using 5cores, of our in-house HPC cluster.
June 5, 2020 13/28 ig 8. Test 1: activations of the first convolutional layer of the encoderfunction for a testing-parameter instance.
Feature maps of the firstconvolutional layer of the encoder function in the DL-ROM neural network for thetesting-parameter instance µ test = (3 . , .
75) cm at ˜ t = 0 . Fig 9. Test 1: activations of the first convolutional layer of the encoderfunction for a testing-parameter instance.
Activations of the first convolutionallayer of the encoder function in the DL-ROM neural network for the testing-parameterinstance µ test = (6 . , .
25) cm at ˜ t = 0 . t = T , the DL-ROM speed-ups are equal to 7 . × and 6 . × withrespect to the FOM and the POD-Galerkin ROM, with N c = 6, computational times. Test 2: Two-dimensional slab with figure of eight re-entry
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 applya classical S1-S2 protocol [3, 44]. In particular, we consider a square slab of cardiactissue Ω = (0 , and apply an initial stimulus at the bottom edge of the domain,i.e. I app ( x , t ) = Ω ( x ) [ t i ,t f ] (˜ t ) , (18)where Ω = { x ∈ Ω : y ≤ . } , t i = 0 ms and t f = 5 ms.A second stimulus under the form I app ( x , t ; µ ) = Ω ( µ ) ( x ) [ t i ,t f ] (˜ t ) , (19)June 5, 2020 14/28 ig 10. Test 1: FOM, POD-Galerkin ROM and DL-ROM CPUcomputational times. FOM, POD-Galerkin ROM and DL-ROM CPUcomputational times to compute ˜ u (¯ t ; µ test ) vs. ¯ t averaged over the testing set. Thanksto the fact that the DL-ROM can be queried at any time istance it is extremelyefficient in computing ˜ u (¯ t ; µ test ) with respect to both the FOM and thePOD-Galerkin ROM.with Ω ( µ ) = { x ∈ Ω : ( x − + ( y − µ ) ≤ (0 . } , t i = 70 ms and t f = 75 ms, is thenapplied. The parameter µ , consisting in the y -coordinate of the center of the secondcircular stimulus, ranges in the parameter space P = [0 . , .
1] cm. This choice has beenmade to obtain a re-entry elicited and sustained until T = 175 ms.We restrict our study to the time interval [95, 175] ms, i.e. we do not consider thefirst time instances in which the re-entry has not arisen yet, being them equal over P .The time step is ∆ t = 0 . / .
9. We consider N = 256 ×
256 = 65536 grid points,implying a mesh size h = 0 . x -axis and the conductivities in the longitudinaland transversal directions to the fibers are σ l = 2 × − cm /ms and σ t = 3 . × − cm /ms, respectively. The parameters appearing in (3) are set to K = 8, a = 0 . b = 0 . ε = 0 . c = 0 .
14, and c = 0 .
3, see [45].The snapshot matrix is built by solving problem (1), completed with the appliedcurrents 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-parameterinstances uniformly distributed in the parameter space and N test = 12 testing-parameterinstances, each of them corresponding to the midpoint of two consecutive training-parameter 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 decreasein 1000 epochs.In Fig 11 we show the FOM solution and the DL-ROM one obtained by settingthe reduced dimension to n = 5, for the testing-parameter instance µ test = 0 . t = 141 . t = 157 . ǫ sk ∈ R N , for k = 1 , . . . , N t , defined as ǫ sk = | u k ( µ test ) − ˜u k ( µ test ) |k u k ( µ test ) k . (20)The trend of the relative error (20) over time, for the selected testing-parameterinstance µ test = 0 . ig 11. Test 2: comparison between FOM and DL-ROM solutions for atesting-parameter instance. FOM solution (left), DL-ROM one (center) with n = 5, and relative error ǫ sk (right) at ˜ t = 141 . t = 157 . µ test = 0 . ǫ sk is belowthe 2% both for ˜ t = 141 . t = 157 . P . Fig 12. Test 2: trend of the relative error over time.
Relative error ǫ sk vs. ˜ t with n = 5 for the testing-parameter instance µ test = 0 . n = 5, for the last time instance, i.e. at ˜ t = 175 ms, for µ test = 0 . µ test = 1 . ig 13. Test 2: comparison between FOM and DL-ROM solutions fordifferent testing-parameter instances. FOM solution (left), DL-ROM one(center) with n = 5, and relative error ǫ sk (right) at ˜ t = 175 ms, for thetesting-parameter instance µ = 0 . µ = 1 . ǫ sk is below the 2.8% both for µ = 0 . µ = 1 . N c ) and the DL-ROM, keeping for all the same degree ofaccuracy achieved by DL-ROM, i.e. ǫ rel = 7 . × − , and running the code on thehardware each implementation is optimized for – a CPU for the FOM and the POD-Galerkin ROM, a GPU for the DL-ROM. In Table 3 we report the CPU time needed tocompute the FOM solution, and the POD-Galerkin ROM (at the testing phase), bothon a full 64 GB node (20 Intel ® Xeon ® E5-2640 v4 2.4GHz cores), and the GPU timerequired by the DL-ROM to compute 875 time instances (the same number of timeinstants considered in the solution of the dynamical systems associated to the FOMand the POD-Galerkin ROM) at testing time, by fixing its dimension to n = 5, on anNvidia GeForce GTX 1070 8 GB GPU. For the sake of completeness, we also report thecomputational time required by the DL-ROM when employing a single CPU node. Itis 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 electrophysiologyproblem like the figure of eight re-entry. Using a nonlinear approach, for which thesolution manifold is approximated through a piecewise linear trial manifold (as in thecase 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-GalerkinROM with N c = 2 (or N c = 4) local reduced bases.In Fig 14 we show the trend of the error indicator (15) over the testing set versus theCPU computational time both for the DL-ROM and the POD-Galerkin ROM at testingphase. Slight improvements of the performance of DL-ROM, in terms of accuracy, are Indeed, 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) theyhave 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 canbe further exploited to speed up calculations.
June 5, 2020 17/28 able 3. Test 2: FOM, POD-Galerkin ROM and DL-ROM computationaltimes. time [s] FOM/ROM dimensionsFOM (CPU) 382 N = 65536DL-ROM (CPU/GPU) 15/1.2 n = 5POD-Galerkin ROM N c = 1 (CPU) 103 n = 1538POD-Galerkin ROM N c = 2 (CPU) 70 n = 1158 , N c = 4 (CPU) 33 n = 435 , , , n , coherently with our previousfindings reported in [16]. Indeed, the DL-ROM is able, also in this case, to accuratelyrepresent the solution manifold by a reduced nonlinear trial manifold of dimension n µ + 1 = 2; for the case at hand, we report the results for n = 5 (very close to theintrinsic dimension n µ +1 = 2 of the problem, and much smaller than the POD-GalerkinROM dimension), providing slightly smaller values of the error indicator (15) than inthe case n = 2. Regarding instead the POD-Galerkin ROM, in Fig 14 we report resultsobtained for different tolerances ε P OD = 10 − , 5 · − , 10 − , 5 · − , 10 − . In the cases N c = 2 and N c = 4 we only report the results related to the smallest POD tolerances,which indeed allow us to meet the prescribed accuracy on the approximation of thegating variable, which would otherwise impact dramatically on the overall accuracy ofthe POD-Galerkin ROM. Moreover, we do not consider more than N c = 4 local reducedbases in order not to generate too small local linear subspaces. As shown in Fig 14, theproposed DL-ROM outperforms the POD-Galerkin ROM in terms of both efficiencyand accuracy. Fig 14. Test 2: trend of the error indicator versus the CPU testingcomputational time.
Error indicator ǫ rel vs. CPU testing computational time fordifferent values of N c and ε P OD . The DL-ROM outperforms the POD-Galerkin ROMin terms of both efficiency and accuracy.In Fig 15 we show the solutions obtained through the POD-Galerkin ROM with N c = 2 (top) and N c = 4 (bottom) local reduced bases, along with the relative errordefined in (20), for the testing-parameter instance µ test = 0 . t = 157 . Fig 15. Test 2: POD-Galerkin ROM solutions for differenttesting-parameter instances.
POD-Galerkin ROM solution (left) and relative error ǫ sk (right) for N c = 2 (top) and N c = 4 (bottom) at ˜ t = 157 . µ test = 0 . N c = 4 local reduced bases), for the testing-parameterinstance µ test = 0 . P = (0 . , .
11) cm and P = (0 . , . Test 3: Three-dimensional ventricle geometry
We finally consider the solution of the coupled system (1) in a three-dimensional leftventricle (LV) geometry, obtained from the 3D Human Heart Model provided by Zy-gote [46]. Here, we consider a single ( n µ = 1) parameter, given by the longitudinalconductivity in the fibers direction. The conductivity tensor takes the form D ( x ; µ ) = σ t I + ( µ − σ t ) f ⊗ f , (21)where σ t = 12 . · .
02 mm /ms; f is determined at each mesh point through a rule-based approach, by solving a suitable Laplace problem [47]. The resulting fibers field isreported in Fig 17. The applied current is defined as I app ( x , t ) = C (2 π ) / α exp (cid:18) − || x − ¯x || β (cid:19) [0 , ¯ t ] (˜ t ) , where ¯ t = 2 ms, C = 1000 mA, α = 50, β = 50 mm , ¯x = [44 . , . , . T mm.June 5, 2020 19/28 ig 16. Test 2: FOM, POD-Galerkin ROM and DL-ROM APs at P andP . AP obtained through the FOM, the DL-ROM and the POD-Galerkin ROM with N c = 4, for the testing-parameter instance µ test = 0 . P = (0 . , .
11) cmand P = (0 . , .
03) cm. The POD-Galerkin ROM approximations are obtained byimposing a POD tolerance ε P OD = 10 − and 10 − , resulting in error indicator (15)values equal to 5 . × − and 7 . × − , respectively. Fig 17. Test 3: fibers distribution.
Fibers field on the Zygote LV geometry.In order to build the snapshot matrix S , we solve problem (1) completed with theconductivity 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 . / .
9. We uniformly sample N t = 1000time instances in (0 , T ) and we zero-padded [34] the snapshot matrix to reshape eachcolumn in a 2D square matrix. The parameter space is provided by P = 12 . · [0 . , . /ms; here we consider N train = 25 training-parameter instances and N test = 24testing-parameter instances computed as in Test 2. In this case, the maximum numberof epochs is set to N epochs = 30000, the batch size is N b = 40 and the training is stoppedif the loss does not decrease over 4000 epochs.In Fig 18 we report the FOM solution for two testing-parameter instances, i.e. µ = 12 . · . /ms and µ = 12 . · . /ms, at ˜ t = 276 ms, to showthe variability of the FOM solution over the parameter space. As expected, front prop-agation is faster for larger values of the parameter µ .June 5, 2020 20/28 ig 18. Test 3: FOM solutions for different testing-parameter instances. FOM solutions for µ = 12 . · . /ms (left) and µ = 12 . · . /ms(right) at ˜ t = 276 ms. By increasing the value of mu the wavefront propagates faster.In Fig 19-20 we report the FOM and DL-ROM solutions, the latter with n = 10, at˜ t = 42 . t = 222 . µ test = 12 . · . /ms and µ test = 12 . · . /ms. The DL-ROM approximation is essentiallyas accurate as the FOM solution. Fig 19. Test 3: comparison between FOM and DL-ROM solutions for atesting-parameter instance at different time instances.
FOM solution (left)and DL-ROM one (right), with n = 10, at ˜ t = 42 . t = 276 ms (bottom),for the testing-parameter instance µ test = 12 . · . /ms.Also for this test case, it is possible to build a reduced nonlinear trial manifold ofdimension very close to the intrinsic one – n µ + 1 = 2 – as long as the maximum numberof epochs N epochs is increased; the choice n = 10 is obtained as the best trade-offbetween accuracy and efficiency of the DL-ROM approximation in this case.The DL-ROM approximation can also replace the FOM solution when evaluatingoutputs of interest. For instance, in Fig 21 and 22 we show the FOM and DL-ROMactivation maps, the latter obtained by choosing n = 10 as DL-ROM dimension. Giventhe electric potential u = u ( x , t ; µ ), the (unipolar) activation map at a point x ∈ Ω isJune 5, 2020 21/28 ig 20. Test 3: comparison between FOM and DL-ROM solutions for atesting-parameter instance at different time instances.
FOM solution (left)and DL-ROM one (right), with n = 10, at ˜ t = 42 . t = 276 ms (bottom),for the testing-parameter instance µ test = 12 . · . /ms.evaluated as the minimum time at which the AP peak reaches x , that is, AC ( x ; µ ) = arg min t ∈ (0 ,T ) (cid:18) u ( x , t ; µ ) = max t ∈ (0 ,T ) u ( x , t ; µ ) (cid:19) . Here we compare the activation maps AC F OM and AC DL − ROM obtained through theFOM and the DL-ROM, respectively, by evaluating the maximum of the relative error ǫ AC ( x ; µ ) = | AC F OM ( x ; µ ) − AC DL − ROM ( x ; µ ) || AC F OM ( x ; µ ) | over the N mesh points; in the case µ = µ test = 12 . · .
31, the maximum relative erroris equal to 4 . × − . Fig 21. Test 3: FOM activation map.
FOM activation map for thetesting-parameter instance µ test = 12 . · .
31 mm /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 . , . , .
82] mm for theJune 5, 2020 22/28 ig 22. Test 3. DL-ROM activation map.
DL-ROM activation map for thetesting-parameter instance µ test = 12 . · .
31 mm /ms with n = 20.testing-parameter instance µ test = 12 . · .
31 mm /ms. Moreover, we also report the bestapproximation of the FOM action potential over a POD space of same dimension n = 20for the sake of comparison. Clearly, in dimension n = 20 the DL-ROM approximationis much more accurate than the POD best approximation; to reach the same accuracy(about ǫ rel = 5 . × − , measured through the error indicator (15)) achieved by theDL-ROM with n = 20, n = 120 POD modes would be required.In Fig 23 (right) we highlight instead the improvements, in terms of efficiency, en-abled by the use of the DL-ROM technique. In particular, we point out the CPU timerequired to solve the FOM for a testing parameter instance, and the one required byDL-ROM (of dimension n = 10) at testing time, by using the time resolutioin each solu-tion computation requires and by varying the FOM dimension N on a 6-core platform the FOM solution with N = 16365 degrees of freedom requires about 40 minutes to becomputed, against 57 seconds required by the DL-ROM approximation, thus implyinga speed-up almost equal to 41 times. Fig 23. Test 3: FOM, DL-ROM and optimal-POD APs for atesting-parameter instance. FOM and DL-ROM CPU computational times.
FOM, DL-ROM and optimal-POD APs for the testing-parameter instance µ test = 12 . · .
31 mm /ms (left). For the same n , the DL-ROM is able to providemore accurate results than the optimal-POD. CPU time required to solve the FOMand by DL-ROM at testing time with n = 10 vs N (right). The DL-ROM is able toprovide a speed-up equal to 41. Numerical tests have been performed on a MacBook Pro Intel Core i7 6-core with 16 GB RAM.
June 5, 2020 23/28 onclusion
In this work we have proposed a new efficient reduced order model obtained usingdeep learning algorithms to boost the solution of parametrized problems in cardiac elec-trophysiology. Numerical results show that the resulting DL-ROM technique, formerlyintroduced in [16], allows one to accurately capture complex wave propagation processes,both in physiological and pathological scenarios.The proposed DL-ROM technique provides ROMs that are orders of magnitude moreefficient than the ones provided by common linear (projection-based) ROMs, built forinstance through a POD-Galerkin reduced basis method, for a prescribed level of accu-racy. Through the use of DL-ROM, it is possible to overcome the main computationalbottlenecks shown by POD-Galerkin ROMs, when addressing parametrized problems incardiac electrophysiology. The most critical points related to (i) the linear superimposi-tion of modes which linear ROMs are based on; (ii) the need to account for the gatingvariables when solving the reduced dynamics, even if not required; and (iii) the neces-sity to use (very often, expensive) hyper-reduction techniques to deal with terms thatdepend nonlinearly on either the transmembrane potential or the input parameters, areall addressed by the DL-ROM technique, which finally yields more efficient and accu-rate approximation than POD-Galerkin ROMs. Moreover, larger time resolutions canbe employed when using a DL-ROM, compared to the ones required by the numericalsolution of a dynamical systems through a FOM or a POD-Galerkin ROM. Indeed, theDL-ROM approximation can be queried at any desired time, without requiring to solvea dynamical system until that time, thus drastically decreasing the computational timerequired to compute the approximated solution at any given time.Outputs of clinical interest, such as activation maps and action potentials, can bemore efficiently evaluated by the DL-ROM technique than by a FOM built throughthe finite element method, while maintaining a high level of accuracy. This work is aproof-of-concept of the DL-ROM technique ability to investigate intra- and inter- sub-jects variability, towards performing multi-scenario analyses in real time and, ultimately,supporting decisions in clinical practice. In this respect, the use of DL-ROM techniquescan foster assimilation of clinical data with physics-driven computational models.
Supporting information
SI Code. Code and data.
The code used in this work can be downloaded from: https://github.com/stefaniafresca/DL-ROM . The training and testing datasets willbe made available upon request to the authors.
SI Appendix. DL-ROM neural network architecture.
Here we report the con-figuration of the DL-ROM neural network used for our numerical tests. We employ a12-layers DFNN equipped with 50 neurons per hidden layer and n neurons in the out-put layer, where n corresponds to the dimension of the reduced nonlinear trial manifold.The architectures of the encoder and decoder functions are instead reported in Table 4and 5. Acknowledgments
The authors have been partially supported by the ERC Advanced Grant iHEART, “Anintegrated heart model for the simulation of the cardiac function”, 2017-2022, P.I. A.Quarteroni (ERC2016AdG, project ID: 740132). Moreover, the authors acknowledgeDr. S. Pagani (MOX, Politecnico di Milano) for the kind help in the FOM softwareJune 5, 2020 24/28 ayer Input Output Kernel size N n Table 4.
Attributes of convolutional and dense layers in the encoder f En . Layer Input Output Kernel size n N Table 5.
Attributes of dense and transposed convolutional layers in the decoder f D .development, and Dr. M. Fedele for providing us the computational mesh of the ZygoteSolid 3D heart model. References
1. Quarteroni A, Manzoni A, Vergara C. The cardiovascular system: Mathe-matical modeling, numerical algorithms, clinical applications. Acta Numerica.2017;26:365–590.2. Quarteroni A, Ded`e L, Manzoni A, Vergara C. Mathematical modelling of the hu-man cardiovascular system: Data, numerical approximation, clinical applications.Cambridge Monographs on Applied and Computational Mathematics. CambridgeUniversity 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 theelectrical activity in the heart. vol. 1. Springer Science & Business Media; 2007.5. Colli Franzone P, Pavarino LF. A parallel solver for reaction–diffusion systems incomputational electrocardiology. Mathematical Models and Methods in AppliedSciences. 2004;14(06):883–911.6. Mirams GR, Pathmanathan P, Gray RA, Challenor P, Clayton RH. Uncertaintyand 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, Path-manathan P, et al. Uncertainty and variability in models of the cardiac actionpotential: Can we build trustworthy models? Journal of Molecular and CellularCardiology. 2016;96:49–62.June 5, 2020 25/28. Hurtado DE, Castro S, Madrid P. Uncertainty quantification of two modelsof cardiac electromechanics. International Journal for Numerical Methods inBiomedical Engineering. 2017; 33(12):e2894.9. Dhamala J, Arevalo HJ, Sapp J, Hor´acek BM, Wu KC, Trayanova NA, et al.Quantifying the uncertainty in model parameters using Gaussian process-basedMarkov 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 uncer-tainty quantification of activation sequences in patient-specific cardiac electro-physiology meeting clinical time constraints. International Journal for NumericalMethods in Biomedical Engineering. 2018;34(7):e2985.11. Johnston BM, Coveney S, Chang ET, Johnston PR, Clayton RH. Quantifyingthe effect of uncertainty in input parameters in a simplified bidomain modelof partial thickness ischaemia. Medical & Biological Engineering & Computing.2018;56(5):761–780.12. Pathmanathan P, Cordeiro JM, Gray RA. Comprehensive uncertainty quantifi-cation and sensitivity analysis for cardiac action potential models. Frontiers inPhysiology. 2019;10.13. Levrero-Florencio F, Margara F, Zacur E, Bueno-Orovio A, Wang Z, SantiagoA, et al. Sensitivity analysis of a strongly-coupled human-based electrome-chanical cardiac model: Effect of mechanical parameters on physiologically rel-evant biomarkers. Computer Methods in Applied Mechanics and Engineering.2020;361:112762.14. Quarteroni A, Manzoni A, Negri F. Reduced basis methods for partial differentialequations: An introduction. vol. 92. Springer; 2016.15. Pagani S, Manzoni A, Quarteroni A. Numerical approximation of parametrizedproblems in cardiac electrophysiology by a local reduced basis method. ComputerMethods in Applied Mechanics and Engineering. 2018;340:530–558.16. Fresca S, Ded´e L, Manzoni A. A comprehensive deep learning-based approach toreduced order modeling of nonlinear time-dependent parametrized PDEs. arXivpreprint arXiv:200104001. 2020.17. Guo M, Hesthaven JS. Reduced order modeling for nonlinear structural analysisusing Gaussian process regression. Computer Methods in Applied Mechanics andEngineering. 2018;341:807–826.18. Guo M, Hesthaven JS. Data-driven reduced order modeling for time-dependent problems. Computer Methods in Applied Mechanics and Engineering.2019;345:75–99.19. Hesthaven J, Ubbiali S. Non-intrusive reduced order modeling of nonlinear prob-lems 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. DR-RNN: A deep residual recurrent neural network formodel reduction. arXiv preprint arXiv:170900939. 2017.June 5, 2020 26/282. Mohan A, Gaitonde DV. A deep learning based approach to reduced order mod-eling for turbulent flow control using LSTM neural networks. arXiv preprintarXiv:18040926. 2018.23. Wan Z, Vlachas P, Koumoutsakos P, Sapsis T. Data-assisted reduced-order mod-eling of extreme events in complex dynamical systems. PLoS one. 2018;13.24. Gonz´alez FJ, Balajewicz M. Deep convolutional recurrent autoencoders forlearning low-dimensional feature dynamics of fluid systems. arXiv preprintarXiv:180801346. 2018.25. Lee K, Carlberg K. Model reduction of dynamical systems on nonlinear manifoldsusing 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 two-variable model of cardiac excitation. ChaosSolitons Fractals. 1996;7(3):293–301.28. FitzHugh R. Impulses and physiological states in theoretical models of nervemembrane. Biophysical Journal. 1961;1(6):445–466.29. Nagumo J, Arimoto S, Yoshizawa S. An active pulse transmission line simulatingnerve axon. Proceedings of the IRE. 1962;50(10):2061–2070.30. Nash MP, Panfilov AV. Electromechanical model of excitable tissue to studyreentrant cardiac arrhythmias. Progress in Biophysics and Molecular Biology.2004;85:501–522.31. Mitchell CC, Schaeffer DG. A two-current model for the dynamics of cardiacmembrane. 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 ques-tions. 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. Availablefrom: .35. LeCun Y, Bottou L, Bengio Y, Haffner P. Gradient based learning applied todocument recognition. Proceedings of the IEEE. 1998; p. 533–536.36. Hinton GE, Zemel RS. Autoencoders, minimum description length, andHelmholtz free energy. Proceedings of the 6th International Conference on NeuralInformation Processing Systems (NIPS’1993). 1994;.37. Kingma DP, Ba J. Adam: A method for stochastic optimization. InternationalConference on Learning Representations (ICLR); 2015.38. Abadi M, Barham P, Chen J, Chen Z, Davis A, Dean J, et al. Tensor-Flow: A system for large-scale machine learning; 2016. Available from: