Deep learning-based reduced order models in cardiac electrophysiology

06/02/2020 ∙ by Stefania Fresca, et al. ∙ Politecnico di Milano 0

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 input-output 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 DL-ROM framework can efficiently provide solutions to parametrized electrophysiology problems, thus enabling multi-scenario analysis in pathological cases. We investigate three challenging test cases in cardiac electrophysiology and prove that DL-ROM outperforms classical projection-based ROMs.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 10

page 14

page 16

page 17

page 19

page 21

page 22

page 23

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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 input-output 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 DL-ROM framework can efficiently provide solutions to parametrized electrophysiology problems, thus enabling multi-scenario analysis in pathological cases. We investigate three challenging test cases in cardiac electrophysiology and prove that DL-ROM outperforms classical projection-based ROMs.

Introduction

The electrical activation of the heart is the main responsible of its contraction, is the result of two processes: at the microscopic scale, the generation of ionic currents through the cellular membrane producing a local action potential; and at the macroscopic scale, the propagation of the action potential from cell to cell in the form of a transmembrane potential [1, 2, 3]

. This latter process can be described by means of partial differential equations (PDEs), suitably coupled with systems of ordinary differential equations (ODEs) modeling the ionic currents in the cells.

Solving this system using a high-fidelity, full order model (FOM) such as, e.g., the finite element (FE) method, is computationally demanding. Indeed, the propagation of the electrical signal is characterized by the fast dynamics of very steep fronts, thus requiring very fine space and time discretizations. [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. Multi-query analysis is relevant in a variety of situations: when analysing multiple scenarios, when dealing with sensitivity analysis and uncertainty quantification (UQ) problems in order to account for inter-subject variability [6, 7, 8]

, for parameter estimation or data assimilation, in which some unknown (or unaccessible) quantities characterizing the mathematical model must be inferred from a set of measurements

[9, 10, 11, 12, 13].

Conventional projection-based reduced order models (ROMs) built, e.g., through the reduced basis (RB) method [14], yields inefficient ROMs when dealing with nonlinear time-dependent parametrized PDE-ODE system as the one arising from cardiac electrophysiology. The three major computational bottlenecks shown by such kind of ROMs for cardiac electrophysiology are due the fact that:

  • the linear superimposition of modes, on which they are based, would cause the dimension of the ROM to be excessively large to guarantee an acceptable accuracy;

  • evaluating the ROM requires the solution of a dynamical system, which might be unstable unless the size of time step is very small;

  • the ROM must also account for the dynamics of the gating variables, even when aiming at computing just the electrical potential. This fact entails an extremely intrusive and costly hyper-reduction stage to reduce the solution of the ODE system to a few, selected mesh nodes [15].

To overcome the limitations of projection-based ROMs, we propose a new, non-intrusive ROM technique based on deep learning (DL) algorithms, which we refer to as DL-ROM. Combining in a suitable way a convolutional autoencoder (AE) and a deep feedforward neural network (DFNN), the DL-ROM technique enables the construction of an efficient ROM, whose dimension is as close as possible to the number of parameters upon which the solution of the differential problem depends. A preliminary numerical assessment of our DL-ROM technique has already been presented in [16], albeit on simpler – yet challenging – test cases.

The proposed DL-ROM technique is a combination of a data-driven with a physics based model approach. Indeed, it exploits snapshots taken from a set of FOM solutions (for selected parameter values and time instances) and deep neural network architectures to learn, in a non-intrusive way, both (i) the nonlinear trial manifold where the ROM solution is sought, and (ii) the nonlinear reduced dynamics. In a linear ROM built, e.g., thorugh proper orthogonal decomposition (POD), the former quantity is nothing but a set of basis functions, while the latter task corresponds to the projection stage in the subspace spanned by these basis functions. Here, our goal is to show that DL-ROM can be effectively used to handle parametrized problems in cardiac electrophysiology, accounting for both physiological and pathological conditions, in order to provide fast and accurate solutions. The proposed DL-ROM is computationally efficient during the testing stage, that is for any new scenario unseen during the training stage. This is particularly useful in view of the evaluation of patient-specific features to enable the integration of computational methods in current clinical platforms.

DL techniques for parametrized PDEs have previously been proposed in other contexts. In [17, 18, 19, 20] feedforward neural networks have been employed to model the reduced dynamics in a less intrusive way, that is, avoiding the costs entailed by projection-based ROMs, but still relying on a linear trial manifold built, e.g., through POD. In [21, 22, 23] the construction of ROMs for nonlinear, time-dependent problems has been replaced by the evaluation of ANN-based regression models. In [24, 25] the reduced trial manifold where the approximation is sought has been modeled through ANNs thus avoiding the linear superimposition of POD modes, on a minimum residual formulation to derive the ROM [25], or without considering an explicit parameter dependence in the differential problem that is considered [24]. In all these works, coupled problems have never been considered. Moreover, very often DL techniques have been exploited to address problems which require only a moderate dimension of projection-based ROMs. We demonstrate that our DL-ROM provides accurate results by constructing ROMs with extremely low-dimension in prototypical test cases. These tests exhibit all the relevant physical features which make the numerical approximation of parametrized problems in cardiac electrophysiology a challenging task.

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 self-generated in the heart and propagated through the myocardium thanks to the excitability of the cardiac cells, the cardiomyocites [3, 26]. When suitably stimulated, cardiomyocites produce a variation of the potential across the cellular membrane, called transmembrane potential. Its evolution in time is usually referred to as action potential, involving a polarization and a depolarization in the early stage of every heart beat. The action potential is generated by several ion channels (e.g., calcium, sodium, potassium) that open and close, and by the resulting ionic currents crossing the membrane. For instance, coupling the so-called monodomain model for the transmembrane potential 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 time-dependent system

(1)

Here denotes a rescaled time111Dimensional 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 diffusion-reaction equation for is two-ways 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 action-potential are, e.g., the FitzHugh-Nagumo [28, 29], the Aliev-Panfilov [27, 30] or the Mitchell and Schaeffer models [31]

. The diffusivity tensor

usually depends on the fibers-sheet structure of the tissue, affecting directional conduction velocities and directions. In particular, by assuming an axisymmetric distribution of the fibers, the conductivity tensor takes the form

(2)

where and are the conductivities in the fibers and the transversal directions.

When a simple phenomenological ionic model is considered, such as the FitzHugh-Nagumo or the Aliev-Panfilov (A-P) model, the ionic current takes the form of a cubic nonlinear function of and a single (dimensionless) gating variable plays the role of a recovery function, allowing to model refractariness of cells. In this paper, we focus on the Aliev-Panfilov model, which consists in taking

(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].

Projection-based ROMs

From an algebraic standpoint, the spatial discretization of system (1) through the Galerkin-finite element (FE) approximation [33] yields the following nonlinear dynamical system for , , 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 snapshots

this 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 projection-based 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 problem-dependent basis, spanned by the columns of a matrix

, which is built according to a suitable hyper-reduction 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 Galerkin-POD 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 time-dependent parametrized PDEs. However, using (7) as an approximation to (4) is known to suffer from several problems. First of all, an extensive hyper-reduction stage (exploiting, e.g., the discrete empirical interpolation method (DEIM)) must be performed in order to be able to evaluate any - or -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 projection-based ROM (7) might easily become very large, due to the basic linearity assumption, by which the solution is given by a linear superimposition of POD modes, thus severely degrading the computational efficiency of the ROM. A possible way to overcome this bottleneck is to rely on local reduced bases, built through POD after the set of snapshots has been split into

clusters, according to suitable clustering (or unsupervised learning) algorithms

[15].

Deep learning-based reduced order modeling (DL-ROM)

To overcome the limitations of linear ROMs, we consider a new, nonlinear ROM technique based on deep learning models. First introduced in [16] and assessed on one-dimensional benchmark problems, the DL-ROM 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 autoencoder222The AE is a particular type of neural network aiming at learning the identity function

(9)
. It is composed by two main parts: the encoder function , where and , mapping the high-dimensional input onto a low-dimensional code ; the decoder function , where , mapping the low-dimensional code to an approximation of the original high-dimensional input . (AE) [35, 36]. More precisely, takes the form

(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 DL-ROM 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 low-dimensional representation

(12)

denotes the encoder function, depending upon a vector of parameters.

Computing the DL-ROM 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 DL-ROM strategy overcomes the three major computational bottlenecks implied by the use of projection-based ROMs, since:

  • the dimension of the DL-ROM can be kept extremely small;

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

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

  • the DL-ROM does not require to account for the dynamics of the gating variables, thus avoiding any hyper-reduction stage. This advantage, already visible when employing a single gating variable as in the test cases addressed later in this paper, might become even more effective when dealing with more realistic ionic models (the so-called I and II generation models), when dozens of additional variables in the system of ODEs must be accounted for [3].

The training stage consists in solving the following optimization problem, in the variable , after the snapshot matrix has been formed:

(13)

where and

(14)

with . The per-exampleloss function (14) combines the reconstruction error (that is, the error between the FOM solution and the DL-ROM approximation) and the error between the intrinsic coordinates and the output of the encoder.

The architecture of DL-ROM is the one shown in Fig 1. The encoder function is used only during the training and validation steps; it is instead discarded during the testing phase. See [16] for further algorithmic details about the training and the testing algorithms required to build and evaluate a DL-ROM.

Fig 1: DL-ROM architecture. DL-ROM architecture used during the training phase. In the red box, the DL-ROM to be queried for any new selected couple during the testing phase. The FOM solution is provided as input to block (A) which outputs . The same parameter instance associated to the FOM, i.e. , enters block (B) which provides as output and the error between the low-dimensional vectors (dashed green box) is accumulated. The intrinsic coordinates are given as input to block (C) returning the ROM approximation . Then the reconstruction error (dashed black box) is computed.

We highlight that the DL-ROM technique does not require to solve a (reduced) nonlinear dynamical system for the reduced degrees of freedom as in (7); rather, it evaluates a nonlinear map for any given couple , 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 diffusion-reaction equation providing the monodomain model is broken when computing the DL-ROM approximation. Moreover, the map features an extremely low dimension, in the most favorable scenario equal to . From a computational perspective, remarkable gains and simplifications can be obtained against a linear ROM, since (i) no hyper-reduction is required to enhance the evaluation of any - or -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 DL-ROM strategy on three relevant test cases in cardiac electrophysiology. Our choice of the numerical tests is aimed at highlighting the performance of our DL-ROM method in challenging electrophysiology problems, namely pathological cases in portion of cardiac tissues or physiological scenarios on realistic left ventricle geometries.

The architecture used to perform all the numerical tests is the one reported in the SI Appendix.. To solve the optimization problem (13)-(14) we use the ADAM algorithm [37] with a starting learning rate equal to . Moreover, we perform cross-validation by splitting the data in training and validation and following a proportion 8:2 and we implement an early-stopping regularization technique to reduce overfitting [34].

To evaluate the performance of the DL-ROM, we use the loss function (14) and on an error indicator defined as

(15)

Neural networks required by our DL-ROM technique have been implemented by means of the Tensorflow deep learning framework [38]; numerical simulations have been carried out on a workstation equipped with an Nvidia GeForce GTX 1070 8 GB GPU.

Test 1: Two-dimensional slab with ischemic region

We consider the computation of the transmembrane potential in a square slab of cardiac tissue in presence of an ischemic (non-conductive) region. The ischemic region may act as anatomical driver of cardiac arrhythmias like tachycardias and fibrillations. The system we want to solve is a slight modification of equations (1), accounting for the presence of a non-conductive region which affects both the conductivity tensor and the ionic current term. The ischemic portion of the domain is modeled by replacing the conductivity tensor , 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, cm2, , the transversal and longitudinal conductivities are cm2/ms and cm2/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, cm2 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 one-step, semi-implicit, 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 training-parameter instances, with ,

. The maximum number of epochs is set equal to

, the batch size is and, regarding the early-stopping criterion, we stop the training if the loss function does not decrease in 500 epochs. For the testing phase, testing-parameter instances , , have been considered.

In Figs 2 and 3 we show the FOM and the DL-ROM solutions, the latter obtained with for the testing-parameter 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 DL-ROM solution is able to capture the diseased (non-conductive) nature of this portion of tissue.

Fig 2: Test 1: comparison between FOM and DL-ROM solutions for a testing-parameter instance. FOM solution (left), DL-ROM solution with (center) and relative error (right) for the testing-parameter instance cm at ms. The maximum of the relative error is and it is associated to the diseased tissue.
Fig 3: Test 1: comparison between FOM and DL-ROM solutions for a testing-parameter instance. FOM solution (left), DL-ROM solution with (center) and relative error (right) for the testing-parameter instance cm at ms. The maximum of the relative error is and it is associated to the diseased tissue.

In Fig 5 we show the action potentials (APs) computed at the six points reported in Fig 4. The DL-ROM 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 DL-ROM technique is able to capture the difference, in terms of AP values, between the diseased and the healthy tissue.

Fig 4: Test 1: location of points . FOM solution evaluated for cm at ms together with the points .
Fig 5: Test 1: comparison between the FOM and DL-ROM APs at . APs evaluated for cm at points . The DL-ROM, with , is able to to sharply reconstruct the AP in almost all the points and the main features are captured also for the points close to the scar.

The AP variability across the parameter space characterizing both the FOM and the DL-ROM solutions is shown in Fig 6. Still with a DL-ROM dimension , we report the APs for cm, with , evaluated at cm. The DL-ROM 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.

Fig 6: Test 1: variability of the FOM and DL-ROM solutions over the parameter space. FOM (right) and DL-ROM (left) AP variability over at cm. The DL-ROM sharply reconstructs the FOM variability over .

By using the DL-ROM technique and setting the dimension of the nonlinear trial manifold equal to the dimension of the solution manifold, i.e. , we obtain an error indicator (15) of . In order to assess the computational efficiency of DL-ROM, we compare it with the POD-Galerkin 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 POD-Galerkin 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 POD-Galerkin ROM for different numbers of clusters.

Table 1: Test 1: dimensions of the POD-Galerkin ROM linear trial manifolds by varying the number 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 DL-ROM with , and the POD-Galerkin ROM with local reduced bases, at testing time, by varying the FOM dimension . Here, with testing time we refer, both for the DL-ROM and the POD-Galerkin ROM, to the time needed to query the ROM over the whole interval , by using for each technique the proper time resolution, for a given testing-parameter instance. Since the DL-ROM solution can be queried at a given time without requiring any solution of a dynamical system to recover the former time instances, the DL-ROM can employ larger time windows compared to the time steps required by the solution of the FOM and POD-Galerkin ROM dynamical systems for the cases at hand. This fact also has a positive impact on the data used during the training phase333Indeed, 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 DL-ROM case, against snapshots in the POD-Galerkin ROM case.. The speed-up obtained, for each value of considered, is reported in Table 1. Both the DL-ROM and the POD-Galerkin ROM allow us to decrease the computational costs associated to the computation of the FOM solution for a testing-parameter instance. However, for a desired level of accuracy, CPU times required by the POD-Galerkin ROM during the testing phase are remarkably higher than the ones required by a DL-ROM with .

Fig 7: Test 1: FOM, DL-ROM and POD-Galerkin ROM CPU computational times. CPU time required to solve the FOM, by DL-ROM at testing time with and by the POD-Galerkin ROM at testing time with vs. . The DL-ROM provides the smallest testing computational time for each considered.

Both the DL-ROM and the POD-Galerkin ROM depend on the FOM dimension . In the case of DL-ROM, 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 POD-Galerkin 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 DL-ROM at testing time scales linearly with , instead the one required by the POD-Galerkin ROM scales linearly with . In particular, even for the larger FOM dimension considered ( for this test case), our DL-ROM is 19 times faster than the POD-Galerkin 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 DL-ROM technique, practice indicates that the CPU time for DL-ROM is smaller than the one for the POD-Galerkin ROM for small values of , in other words only with very large values of the POD-Galerkin ROM outperforms the DL-ROM strategy. Indeed, a linear fitting of the DL-ROM and the POD-Galerkin ROM CPU times444 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 Tusscher-Panfilov ionic model [40] is used. This latter would indeed require smaller values of compared to the Aliev-Panfilov model, due to the shape of the AP. See, e.g., [41, 42] for further details. in Fig 7 highlights that for and , DL-ROM could be almost 10 and 5 times, respectively, faster than the POD-Galerkin ROM for the same values of . Note that the results of this section have been obtained by employing the DL-ROM on a single CPU, an architecture which is not favorable to neural networks555Indeed, all tests are performed on a node (20 Intel® Xeon® E5-2640 v4 2.4GHz cores), using 5 cores, of our in-house HPC cluster.. Further improvements are expected when employing our DL-ROM on a GPU for a given testing-parameter instance.

FOM vs. DL-ROM 472 536 539 412
FOM vs. POD-Galerkin ROM 3 6 12 22

DL-ROM and POD-Galerkin ROM vs. FOM speed-up by varying . The DL-ROM speed-up is remarkably higher than the one obtained by using the POD-Galerkin ROM.

Table 2: Test 1: DL-ROM and POD-Galerkin ROM vs. FOM speed-up.

In Figs 8 and 9 we show the feature maps of the first convolutional layer of the encoder function , for , in the DL-ROM neural network when the FOM solution for the testing-parameter 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 testing-parameter instances, we observe the translation equi-variance property

[34] that convolutional layers hold when applied to the part of cardiac tissue corresponding to the scar. Moving to deeper layers, feature maps become increasingly abstract, and less visually interpretable; however, the extracted high-level features are still related both to the ischemic region and the electrical activation pattern.

Fig 8: Test 1: activations of the first convolutional layer of the encoder function for a testing-parameter instance. Feature maps of the first convolutional layer of the encoder function in the DL-ROM neural network for the testing-parameter instance cm at ms.
Fig 9: Test 1: activations of the first convolutional layer of the encoder function for a testing-parameter instance. Activations of the first convolutional layer of the encoder function in the DL-ROM neural network for the testing-parameter instance cm at ms.

We highlight that since the DL-ROM solution can be evaluated at any desired time instance without solving any dynamical system, the resulting computational time entailed by the DL-ROM at testing time are drastically reduced compared to the ones required by the FOM or the POD-Galerkin ROM to compute solutions at a particular time instance. In Fig 10 we show the DL-ROM, FOM and POD-Galerkin ROM CPU time needed to compute the approximated solution at , for 1, 10, 100 and 400 ms averaged over the testing set and with . We perform the training phase of the POD-Galerkin 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 DL-ROM CPU time to compute does not vary over and, by choosing , the DL-ROM speed-ups are equal to and with respect to the FOM and the POD-Galerkin ROM, with , computational times.

Fig 10: Test 1: FOM, POD-Galerkin ROM and DL-ROM CPU computational times. FOM, POD-Galerkin ROM and DL-ROM CPU computational times to compute vs. averaged over the testing set. Thanks to the fact that the DL-ROM can be queried at any time istance it is extremely efficient in computing with respect to both the FOM and the POD-Galerkin ROM.

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 apply a classical S1-S2 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 re-entry 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 re-entry 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 cm2/ms and cm2/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 semi-implicit scheme, over time instances. Moreover, we consider

training-parameter instances uniformly distributed in the parameter space and

testing-parameter instances, each of them corresponding to the midpoint of two consecutive training-parameter 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 early-stopping criterion, we stop the training if the loss does not decrease in 1000 epochs.

In Fig 11 we show the FOM solution and the DL-ROM one obtained by setting the reduced dimension to , for the testing-parameter instance cm, at ms and ms, together with the relative error , for , defined as

(20)
Fig 11: Test 2: comparison between FOM and DL-ROM solutions for a testing-parameter instance. FOM solution (left), DL-ROM one (center) with , and relative error (right) at ms (top) and ms (bottom), for the testing-parameter instance cm. The relative error is below the 2 both for ms and ms, the maximum value of the error being associated to very few points of the domain.

The trend of the relative error (20) over time, for the selected testing-parameter 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

.

Fig 12: Test 2: trend of the relative error over time. Relative error vs. with for the testing-parameter instance cm (the red band indicates the IQR). The distribution of the error maintains uniform over time.

In Fig 13 we show the FOM and the DL-ROM 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 DL-ROM to capture it.

Fig 13: Test 2: comparison between FOM and DL-ROM solutions for different testing-parameter instances. FOM solution (left), DL-ROM one (center) with , and relative error (right) at ms, for the testing-parameter instance cm (top) and cm (bottom). The relative error is below the 2.8 both for cm and cm, the maximum value of the error being associated to very few points of the domain.

We now compare the computational times required by the FOM, the POD-Galerkin ROM (for different values of ) and the DL-ROM, keeping for all the same degree of accuracy achieved by DL-ROM, i.e. , and running the code on the hardware each implementation is optimized for – a CPU for the FOM and the POD-Galerkin ROM, a GPU666Indeed, 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 DL-ROM. In Table 3 we report the CPU time needed to compute the FOM solution, and the POD-Galerkin ROM (at the testing phase), both on a full 64 GB node (20 Intel® Xeon® E5-2640 v4 2.4GHz cores), and the GPU time required by the DL-ROM to compute 875 time instances (the same number of time instants considered in the solution of the dynamical systems associated to the FOM and the POD-Galerkin ROM) at testing time, by fixing its dimension to , on an Nvidia GeForce GTX 1070 8 GB GPU. For the sake of completeness, we also report the computational time required by the DL-ROM when employing a single CPU node. It is evident that a POD-Galerkin ROM, built employing a global reduced basis , is not amenable to a complex and challenging pathological cardiac electrophysiology problem like the figure of eight re-entry. Using a nonlinear approach, for which the solution manifold is approximated through a piecewise linear trial manifold (as in the case of or local reduced bases) reduces the online computational time. However, the DL-ROM still confirms to provide a more efficient ROM, almost 5 (or 2) times faster on the CPU, and 39 (or 19) faster on the GPU, than the POD-Galerkin ROM with (or ) local reduced bases.

time [s] FOM/ROM dimensions
FOM (CPU) 382
DL-ROM (CPU/GPU) 15/1.2
POD-Galerkin ROM (CPU) 103
POD-Galerkin ROM (CPU) 70
POD-Galerkin ROM (CPU) 33

POD-Galerkin ROM and DL-ROM computational times along with FOM and reduced trial manifold(s) dimensions. DL-ROM provides a more efficient ROM with respect the POD-Galerkin ROMs.

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

In Fig 14 we show the trend of the error indicator (15) over the testing set versus the CPU computational time both for the DL-ROM and the POD-Galerkin ROM at testing phase. Slight improvements of the performance of DL-ROM, in terms of accuracy, are obtained for a small increase of the DL-ROM dimension , coherently with our previous findings reported in [16]. Indeed, the DL-ROM 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 POD-Galerkin ROM dimension), providing slightly smaller values of the error indicator (15) than in the case . Regarding instead the POD-Galerkin 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 POD-Galerkin 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 DL-ROM outperforms the POD-Galerkin ROM in terms of both efficiency and accuracy.

Fig 14: Test 2: trend of the error indicator versus the CPU testing computational time. Error indicator vs. CPU testing computational time for different values of and . The DL-ROM outperforms the POD-Galerkin ROM in terms of both efficiency and accuracy.

In Fig 15 we show the solutions obtained through the POD-Galerkin ROM with (top) and (bottom) local reduced bases, along with the relative error defined in (20), for the testing-parameter instance cm at ms. In both cases, we have considered the setting yielding the most efficient POD-Galerkin ROM approximation, which require about 30 (40, respectively) seconds to be evaluated. By comparing Fig 15 and Fig 11 (bottom), we observe that the DL-ROM outperforms the POD-Galerkin ROM in terms of accuracy.

Fig 15: Test 2: POD-Galerkin ROM solutions for different testing-parameter instances. POD-Galerkin ROM solution (left) and relative error (right) for (top) and (bottom) at ms, for cm.

In Fig 16 we show the action potentials obtained through the FOM, the DL-ROM and the POD-Galerkin ROM (with local reduced bases), for the testing-parameter instance cm, and evaluated at cm and cm. These two points are close to the left core of the figure of eight re-entry, 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 DL-ROM, while the POD-Galerkin ROM might fail in this respect.

Fig 16: Test 2: FOM, POD-Galerkin ROM and DL-ROM APs at and . AP obtained through the FOM, the DL-ROM and the POD-Galerkin ROM with , for the testing-parameter instance cm, at cm and cm. The POD-Galerkin ROM approximations are obtained by imposing a POD tolerance and , resulting in error indicator (15) values equal to and , respectively.

Test 3: Three-dimensional ventricle geometry

We finally consider the solution of the coupled system (1) in a three-dimensional left ventricle (LV) geometry, obtained from the 3D Human Heart Model provided by Zygote [46]. Here, we consider a single () parameter, given by the longitudinal conductivity in the fibers direction. The conductivity tensor takes the form

(21)

where mm2/ms; is determined at each mesh point through a rule-based approach, by solving a suitable Laplace problem [47]. The resulting fibers field is reported in Fig 17.

Fig 17: Test 3: fibers distribution. Fibers field on the Zygote LV geometry.

The applied current is defined as

where ms, mA, , mm2, 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 semi-implicit scheme in time over a uniform partition of with ms and time step . We uniformly sample time instances in

and we zero-padded

[34] the snapshot matrix to reshape each column in a 2D square matrix. The parameter space is provided by mm2/ms; here we consider training-parameter instances and testing-parameter 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 testing-parameter instances, i.e. mm2/ms and mm2/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 .

Fig 18: Test 3: FOM solutions for different testing-parameter instances. FOM solutions for mm2/ms (left) and mm2/ms (right) at ms. By increasing the value of the wavefront propagates faster.

In Fig 19-20 we report the FOM and DL-ROM solutions, the latter with , at ms and ms, for two testing-parameter instances, mm2/ms and mm2/ms. The DL-ROM approximation is essentially as accurate as the FOM solution.

Fig 19: Test 3: comparison between FOM and DL-ROM solutions for a testing-parameter instance at different time instances. FOM solution (left) and DL-ROM one (right), with , at ms (top) and ms (bottom), for the testing-parameter instance mm2/ms.
Fig 20: Test 3: comparison between FOM and DL-ROM solutions for a testing-parameter instance at different time instances. FOM solution (left) and DL-ROM one (right), with , at ms (top) and ms (bottom), for the testing-parameter instance mm2/ms.

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 trade-off between accuracy and efficiency of the DL-ROM approximation in this case.

The DL-ROM approximation can also replace the FOM solution when evaluating outputs of interest. For instance, in Fig 21 and 22 we show the FOM and DL-ROM activation maps, the latter obtained by choosing as DL-ROM 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 DL-ROM, respectively, by evaluating the maximum of the relative error

over the mesh points; in the case , the maximum relative error is equal to .

Fig 21: Test 3: FOM activation map. FOM activation map for the testing-parameter instance mm2/ms.
Fig 22: Test 3. DL-ROM activation map. DL-ROM activation map for the testing-parameter instance mm2/ms with .

In Fig 23 (left) we report the action potentials obtained with the FOM and the DL-ROM, this latter with , computed at point mm for the testing-parameter instance mm2/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 DL-ROM 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 DL-ROM 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 DL-ROM technique. In particular, we point out the CPU time required to solve the FOM for a testing parameter instance, and the one required by DL-ROM (of dimension ) at testing time, by using the time resolutioin each solution computation requires and by varying the FOM dimension on a 6-core platform777Numerical tests have been performed on a MacBook Pro Intel Core i7 6-core with 16 GB RAM. the FOM solution with degrees of freedom requires about 40 minutes to be computed, against 57 seconds required by the DL-ROM approximation, thus implying a speed-up almost equal to 41 times.

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

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 DL-ROM technique, formerly introduced 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 more efficient than the ones provided by common linear (projection-based) ROMs, built for instance through a POD-Galerkin reduced basis method, for a prescribed level of accuracy. Through the use of DL-ROM, it is possible to overcome the main computational bottlenecks shown by POD-Galerkin 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) hyper-reduction techniques to deal with terms that depend nonlinearly on either the transmembrane potential or the input parameters, are all addressed by the DL-ROM technique, which finally yields more efficient and accurate approximation than POD-Galerkin ROMs. Moreover, larger time resolutions can be employed when using a DL-ROM, compared to the ones required by the numerical solution of a dynamical systems through a FOM or a POD-Galerkin ROM. Indeed, the DL-ROM 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 DL-ROM technique than by a FOM built through the finite element method, while maintaining a high level of accuracy. This work is a proof-of-concept of the DL-ROM technique ability to investigate intra- and inter- subjects variability, towards performing multi-scenario analyses in real time and, ultimately, supporting decisions in clinical practice. In this respect, the use of DL-ROM techniques can 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 will be made available upon request to the authors.

SI Appendix.

DL-ROM neural network architecture.

Here we report the configuration of the DL-ROM neural network used for our numerical tests. We employ a 12-layers 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.

Layer Input Output Kernel size of filters Stride Padding dimension dimension 1 [5, 5] 8 1 SAME 2 [5, 5] 16 2 SAME 3 [5, 5] 32 2 SAME 4 [5, 5] 64 2 SAME 5 256 6 256

Table 4: Attributes of convolutional and dense layers in the encoder .

Layer Input Output Kernel size of filters Stride Padding dimension dimension 1 256 2 256 3 [5, 5] 64 2 SAME 4 [5, 5] 32 2 SAME 5 [5, 5] 16 2 SAME 6 [5, 5] 1 1 SAME

Table 5: Attributes of dense and transposed convolutional layers in the decoder .

Acknowledgments

The authors have been partially supported by the ERC Advanced Grant iHEART, “An integrated heart model for the simulation of the cardiac function”, 2017-2022, 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 process-based 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 patient-specific 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. Levrero-Florencio F, Margara F, Zacur E, Bueno-Orovio A, Wang Z, Santiago A, et al. Sensitivity analysis of a strongly-coupled human-based 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 learning-based approach to reduced order modeling of nonlinear time-dependent 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. 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 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.

    DR-RNN: 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. Data-assisted reduced-order modeling of extreme events in complex dynamical systems. PLoS one. 2018;13.
  •  24. González FJ, Balajewicz M. Deep convolutional recurrent autoencoders for learning low-dimensional 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 two-variable 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 two-current 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 large-scale machine