Reconstructing Complex Cardiac Excitation Waves From Incomplete Data Using Echo State Networks and Convolutional Autoencoders

The mechanical contraction of the pumping heart is driven by electrical excitation waves running across the heart muscle due to the excitable electrophysiology of heart cells. With cardiac arrhythmias these waves turn into stable or chaotic spiral waves (also called rotors) whose observation in the heart is very challenging. While mechanical motion can be measured in 3D using ultrasound, electrical activity can (so far) not be measured directly within the muscle and with limited resolution on the heart surface, only. To bridge the gap between measurable and not measurable quantities we use two approaches from machine learning, echo state networks and convolutional autoencoders, to solve two relevant data modelling tasks in cardiac dynamics: Recovering excitation patterns from noisy, blurred or undersampled observations and reconstructing complex electrical excitation waves from mechanical deformation. For the synthetic data sets used to evaluate both methods we obtained satisfying solutions with echo state networks and good results with convolutional autoencoders, both clearly indicating that the data reconstruction tasks can in principle be solved by means of machine learning.


INTRODUCTION
Cardiac arrhythmias, such as ventricular or atrial fibrillation, are electro-mechanical dysfunctions of the heart that are associated with complex, chaotic spatio-temporal excitation waves within the heart muscle resulting in incoherent mechanical contraction and a significant loss of pump function [1][2][3]. Ventricular fibrillation (VF) is the most common deadly manifestation of a cardiac arrhythmia and requires immediate defibrillation using high-energy electric shocks. Atrial fibrillation (AF) is the most common form of a cardiac arrhythmia, affecting 33 million patients worldwide [62]. While not immediately life-threatening, AF is considered to be responsible for 15% of strokes if left untreated [63,64]. The structural substrate and functional mechanisms that underlie the onset and perpetuation of VF and AF are not fully understood. It is generally agreed that imaging of the cardiac electrical and mechanical function is key to an improved mechanistic understanding of cardiac disease and the development of novel diagnosis and therapy. This has motivated the development of non-invasive and invasive electrophysiological measurement and imaging modalities. Electrical activity of the heart can (so far) noninvasively be measured on its surface, only. Direct measurements can be made in vivo inside the heart using socalled basket catheters with typically 64 electrodes or in ex-vivo experiments, where an extracted heart in a Langendorff perfusion set-up is kept beating and the cell membrane voltage on the epicardial surface is made visible using fluorescent dyes (a method also known as optical mapping) [4]. A method for indirect observation of electrical excitation waves is ECG imaging where an array of EEG-electrodes is placed on the body surface and an (ill-posed) inverse problem is solved to estimate the potential on the surface of the heart. Mechanical contraction and deformation of the heart tissue can be studied in full 3D using ultrasound, in 2D using real-time MRT [5] or (using optical mapping) by motion tracking in Langendorff experiments.
The reconstruction of patterns of action potential wave propagation in cardiac tissue from ultrasound has been introduced by Otani et al. [6,7]. They proposed to use ultrasound to visualize the propagation of these waves through the mechanical deformations they induce and to reconstruct action potential-induced active stress from the deformation. Provost et al. [8] introduced electromechanical wave imaging to map the mechanical deformation of cardiac tissue at high temporal and spatial resolutions. The observed deformations resulting from the electrical activation were found to be closely correlated with electrical activation sequences. The cardiac excitation-contractioncoupling (ECC) [9] has also been studied in optical mapping experiments in Langendorff-perfused isolated hearts [10][11][12]. Using electromechanical optical mapping [12], it was shown that during ventricular tachyarrhythmias electrical rotors introduce corresponding rotating mechanical waves. These co-existing electro-mechanical rotors were observed on the epicardial surface of isolated Langendorff-perfused intact pig and rabbit hearts using optical mapping [13]. Using high-resolution ultrasound, these mechanical rotors were also observed inside the ventricular wall during ventricular tachycardia and fibrillation [13].
All these measurement modalities are limited, in particular those suitable for in vivo applications. Measurements with basket catheters are effectively undersampling the spatio-temporal wave pattern. Inverse ECGs suffer from ill-posedness and require regularization that may lead to loss of spatial resolution and blurring. Limited spatial resolution is also an issue with ultrasound measurements, but they are currently the only way to "look inside" the heart, albeit measuring only mechanical motion. Electrical excitation waves inside the heart muscle are so far not accessible by any measurement modality available.
These limitations motivated the search for algorithms to reconstruct electro-mechanical wave dynamics in cardiac tissue from measurable quantities. Berg et al. [14] devised synchronization-based system identification of extended excitable media, in which model parameters are estimated by minimizing the synchronization error. Using this approach, Lebert and Christoph [15] demonstrated that electro-mechanic wave dynamics of excitable-deformable media can be recovered from a limited set of observables using a synchronization-based data assimilation approach. Hoffman et al. reconstructed electrical wave dynamics using ensemble Kalman filters [16,17]. In another approach, it was shown that echo state networks [18] and deep convolutional neural networks [19,20] provide excellent cross estimation results for different variables of a mathematical model describing complex electrical excitation waves during cardiac arrhythmias. Following this approach, Christoph and Lebert [21] demonstrated the reconstruction of electrical excitation and active stress from deformation using a simulated deformable excitable medium. To continue this research and to address the general challenge of missing or impaired observations we consider in this article two tasks: (i) recovering electrical excitation patterns from noisy, blurred or undersampled observations and (ii) reconstructing electrical excitation waves from mechanical deformation. To solve the corresponding data processing and cross-prediction tasks two machine learning methods are employed and evaluated: echo state networks and convolutional autoencoders. Both algorithms are applied to synthetical data generated by prototypical models for electrophysiology and electromechanical coupling.

METHODS
In this section we will first introduce in Sections 2.1 and 2.2 the mathematical models describing cardiac dynamics which were used to generate the example data for the two tasks to be solved: (i) recovering electrical wave pattern from impaired observations and (ii) crosspredicting electrical excitation from mechanical deformation. Then in Section 2.3 both machine learning methods used for solving these tasks, echo state networks (Section 2.3.1) and convolutional autoencoders (Section 2.3.2), will be briefly introduced.

Recovering Complex Spatio-Temporal Wave Patterns From Impaired Observations
For motivating, illustrating, and evaluating the employed methods for dealing with incomplete or distorted observations we shall use spatio-temporal time series generated with the Bueno-Orovio-Cherry-Fenton (BOCF) model [22] describing complex electrical excitation patterns in the heart during cardiac arrhythmias. The BOCF model is a set of partial differential equations (PDEs) with four variables and will be introduced in Section 2.1.1. In Section 2.1.2 a formal description of the data recovery tasks will be given.  [22].
Cardiac dynamics is controlled by electrical excitation waves triggering mechanical contractions of the heart. In the case of cardiac arrhythmias like lethal ventricular fibrillation, wave break-up and complex chaotic wave patterns occur resulting in significantly reduced pump performance of the heart. From the broad range of mathematical models describing this spatio-temporal dynamics [23] we chose the Bueno-Orovio-Cherry-Fenton (BOCF) model [22] to generate spatio-temporal time series that are used as a benchmark to validate our approaches for reconstructing complex wave patterns in excitable media from incomplete data. The BOCF model consists of four system variables whose evolution is given by four (partial) differential equations The variable u represents the continuum limit representation of the membrane voltage of cardiac cells and the variables v, w, and s are gating variables controlling ionic transmembrane currents J si , J fi and J so given by the equations Here H(·) denotes the Heaviside function and the currents depend on the following seven voltage controlled variables (3)  For simulating the dynamics we used the set of parameters given in Table 1 for which the BOCF model was found [22] to exhibit excitation wave dynamics similar to the Ten Tusscher-Noble-Noble-Panfilov (TNNP) model [24] describing human heart tissue.
Typical snapshots of the four variables during a chaotic evolution are shown in Figure 1. The spatio-temporal chaotic dynamics of this system is actually transient chaos whose lifetime grows exponentially with system size [25,26]. To obtain chaotic dynamics with a sufficiently long lifetime the system has been simulated on a domain of 512 × 512 grid points with a grid constant of Δx 1.0 space units and a diffusion constant D 0.2. Furthermore, an explicit Euler stepping in time with Δt 0.1, a 5 point approximation of the Laplace operator, and no-flux boundary conditions were used for solving the PDEs.

Reconstruction Tasks
Experimental measurements of the dynamics of a system of interest often allow only the observation of some state variables (e.g., the membrane voltage) and may provide only incomplete or distorted information about the measured observable. Typical limitations are (additive) measurement noise and low-spatial resolution (due to the experimental conditions and/or the available hardware). Formally, measurements impaired due to noise, blurring or undersampling can be described as follows: Let X n ∈ R r×c be the measured data (here: snapshots of the field u) where r and c specify the two spatial dimensions. Each sample X n with n 1, . . . , N corresponds to a true system output X ′ n ∈ R r′×c′ that is assumed to be known only during the training phase in terms of a training set D {Z 1 (X 1 , X ′ 1 ), . . . , Z N (X N , X ′ N )}. Note that with coarse graining r ≤ r ′ and c ≤ c ′ . The task is to predict the true system output X ′ from impaired observations X which belong to one of the following three cases: 1. Noisy data: To add noise each element of X ′ is replaced with probability p by 0 or 1 drawn from a Bernoulli distribution B(0.5) (note that in our case X ′ is given by the variable u of the BOCF model which has a range of [0, 1]). To simulate different levels of noise different probabilities p 0.1, 0.2, . . . , 0.9 are used to generate noisy data sets {X n }. In the following p is called the noise level. 2. Blurred data: Date with reduced spatial resolution are obtained as Fourier low-pass filtered data X F −1 (P m (F (X ′ ))) where F and F −1 denote the Fourier transform and its inverse, respectively, and P m is a projection where frequencies outside a radius m ∈ [2, 4, 8, . . . , 18] (Manhattan distance) centered at frequency zero are set to zero. 3. Undersampled data: To generate undersampled date X ′ is down-sampled R r′×c′ → R r×c with r < r ′ and c < c ′ by accessing every 2 i -th value of X ′ , where i ∈ [1,7]. Figure 2 shows examples of the three types of impaired observations.

Predicting Electrical Excitation From Mechanical Contraction
To learn the relation between mechanical deformation and electrical excitation inverse modelling data were generated by a conceptual electro-mechanical model consisting of an Aliev-Panfilov model describing the electrical activity and a driven mass-spring-system [15].

Aliev-Panfilov Model
Specifically developed to mimic cardiac action potentials in the myocardium, the Aliev-Panfilov model is a modification of the FitzHugh-Nagumo model, which reproduces the characteristic shape of electric pulses occurring in the heart [27]. It is given by a set of two differential equations, FIGURE 3 | Two dimensional mass-spring damper system with one active spring modelling fibre orientation (red) and one passive spring (gray), the centre of mass x cm , the four points of attachment q i to the structural springs and the orientation parameter η.
Frontiers in Applied Mathematics and Statistics | www.frontiersin.org March 2021 | Volume 6 | Article 616584 in which u and v are the normalized membrane voltage and the recovery variable, respectively, and a, b and k are model parameters. The term ∇(D · ∇u) accounts for the diffusion, in which the tensor D can be used to model anisotropies in the myocardial tissue. In addition, the term ϵ(u, v) is introduced to adjust the shape of the restitution curve by modulating the parameters μ 1 and μ 2 . The computational advantage of the Aliev-Panfilov model lies in its simplicity over other ion-flowbased models which allows shorter runtimes and combined with the elastomechanical model, keeps computational costs fairly reasonable. For this reason, the Aliev-Panfilov model was chosen for generating synthetic data from complex chaotic electromechnical wave dynamics.
Within the heart muscle, the myocardium, cells contract upon electrical excitation through a passing action potential. At this point it is important to note that muscle fibre contracts along its principal orientation which has to be considered during the implementation of the mechanical part of the simulation. To couple the mechanical contraction of the muscle fibre to electrical excitation of a cell, as an extension to the Aliev-Panfilov model the active stress T a was introduced by Nash and Panfilov [28] which leads to contraction in the principal orientation of the muscle fibre. The change of the active stress is described by where k T controls the strength of the build-up of active stress. The term ϵ T (u) regulates the influence of u on T a for large u. In our simulations we use a smooth function introduced by Göktepe and Kuhl [29] given by Here, ξ T controls the steepness of the transition between ϵ ∞ and ϵ T,0 and u 0 denotes the potential threshold for the activation of the active stress, with ϵ ∞ < ϵ T,0 to achieve a physiological time course [30].

Mass-Spring Damper System
The elasto-mechanical properties of the cardiac muscle fibre were implemented using a modified two-dimensional mass-spring damper system [31]. For the current study the mass-spring system was implemented in two dimensions because this allows shorter runtimes and primarily serves as a proof-ofprinciple for the evaluated reconstruction approach. In its two-dimensional form this mechanical model might correspond best to a cut-out of the atrium's wall, since there the muscle tissue is less than 4 mm thick. However, this massspring system can easily be expanded to three dimensions (see [15]).  March 2021 | Volume 6 | Article 616584 FIGURE 7 | Autoencoder architecture used for reconstruction from undersampled observations. Each block is a set of layers. The layer labeling is the same as in Figure 6. Visualized with Net2Vis [55].   Placed on a regular lattice, one mechanical cell is made up of four particles x i at the corners connected by structural springs and two sets of orthogonal springs connecting the centre of mass x → cm to each side of the cell (see Figure 3). The springs in the middle of the cell are called axial springs, of which one is made to be active (red). Here it is important to point out that one cell in the electrical model corresponds to one cell in the mechanical mass-spring system. For setting the fibre orientation through the active axial spring, the orientation parameter η ∈ [0, 1] has been introduced, with which the four points of attachment q i can be computed easily. This parameter can be set individually for each cell, so that various fibre orientations can be modelled.
Using x → Here l j,0 , k j and l a,0 , k a denote the resting lengths and spring constants of the passive and active spring, respectively. From Eq. (9) it can be seen that, upon a rise in active stress T a from Eq. (6), the active spring contracts and an inward force is generated. The parameter c a represents a scaling factor to modulate the influence of the active stress. Through the orientation parameter the forces from the active and passive spring can be redistributed to the corresponding particles at the corners.
For example for q → 1 , the force on x 0 would be f and on x 1 it would amount to f → In addition, the mechanical grid is held together by structural forces between the corner particles, which can be computed using with l ij being the resting length between particle x i and x j . Finally, with all the above forces acting on particle x i with mass m i , its motion is determined according to with the sum {a}{j}{ij} f k → of all relevant springs pulling or pushing the particle. The damping constant ν sets the strength of the damping to increase the stability of the mechanical system as a whole.
The area of each cell was calculated with a simple formula for a general quadrilateral using the positions of its four corners. As a measure of contraction, the relative change of area has been used. The numerical algorithm for solving the full set of electro-mechanical ODEs is summarized in the Appendix.

Reconstruction Task
The inverse modelling data are generated by forward modelling M : u1ΔA using the output of Equations (4) and (13). The task is to train an ESN or CAE to approximate M −1 : ΔA1u. To fulfill this task we use the membrane voltages and the local deformations at all r × c grid points sampled at times t n . The training data set D {Z 1 (X 1 , X ′ 1 ), . . . , Z N (X N , X ′ N )} thus consists of snapshots X n ∈ R r×c and X ′ n ∈ R r×c of the relative mechanical deformation ΔA(t n ) and the membrane voltage u(t n ), respectively, and we aim at approximating M −1 : X n 1X ′ n with r, c 100.  Table 3).

Machine Learning Methods
In this section we will introduce the two machine learning approaches, echo state networks (ESN) [32] and convolutional autoencoders (CAE) [33], that will be applied to solve the reconstruction tasks defined in Section 2.1.2.

Echo State Network
Echo state networks have been introduced in 2001 by Jaeger [32] as a simplified type of recurrent neural network, in which the weights describing the strength of the connections within the network are fixed. In its general composition an ESN subdivides into three sections [32], as illustrated in Figure 4. First of all, there is the input layer into which the input signal u → n ∈ R Nu and a constant bias b in are fed. Secondly, the intermediate reservoir consists of N nonlinear units and its state is given by s → n ∈ R N . And lastly, the output layer provides the output signal y → n ∈ R N y . Here, n denotes the discrete time steps n 1, . . . , T.
The concatenated bias-input vector [b in ; u → n ] is fed into the reservoir through the input matrix W in ∈ R N×(1+N u ) . Inside the reservoir connections are given by the weight matrix W ∈ R N×N , where N is the reservoir size. Together with the input matrix it is possible to determine the state of the reservoir at time n through the update rule in which [·; ·] denotes a concatenated vector. The input bias b in , as well as the later introduced output bias b out were both set to 1 in the following. The parameter α ∈ (0, 1] in Eq. (14) represents the leaking rate which controls how much of a neuron's activation is carried over to the next time step and can be used as a parameter to enhance predictions. As for the transfer function f in (·) we use tanh(·) and the network dynamics used has no feedback loop. Only the weights W out providing the output signal are adapted during the training process by minimizing the cost function [34] C(W out ) n y →true where Tr denotes the trace of a matrix and λ controls the impact of the regularization term that prevents overfitting [35]. The final output matrix is given by the minimum of the cost function at W out YX T (XX T + λ1) − 1 where X and Y are matrices whose columns are given by the vectors x → n and y →true n , respectively. Both matrices W in and W, are initialized with random values from the interval [−0.5, 0.5]. Since in experiments it turned out that more diverse dynamics could be modelled using networks in which only a small percentage ϵ of weights inside the reservoir remained non-zero [32], the weight matrix W is made sparse with only a portion ϵ of its values remaining non-zero. Furthermore, it is scaled by a factor ρ |μmax| where μ max denotes here the largest eigenvalue of W and ρ is a hyperparameter for optimizing the performance (by ensuring the so-called echo state property [36]). To reduce the probability of drawing an dysfunctional set of matrix entries the randomly generated matrices W in and W were selected from four different realisations. To optimize the performance of the ESN five hyperparameters (N, ϵ, ρ, α and λ) are tuned. Reservoir computing using ESNs for predicting chaotic dynamics has already been demonstrated in 2004 by Jaeger and Haas [37]. Since then many studies appeared analyzing and optimizing this approach (see, for example [38][39][40][41][42][43][44], and references cited therein). In particular, it has been pointed out how reservoir computing exploits generalized synchronization of uni-directionally coupled systems [45,46].
Recently, applications of ESNs to spatio-temporal time series have been presented [18,47] employing many networks operating in parallel at different spatial locations based on the concept of (reconstructed) local states [48]. In particular, using this mode of reservoir computing it was possible to perform a cross-prediction between the four different variables of the BOCF model [18]. Therefore, for the current task of reconstructing data from impaired observations we build on the previous ESN design and modelling procedure. For each pixel an ESN is trained receiving input from neighboring pixels, only, representing the local state at the location of the reference pixel as illustrated in Figure 5. This design introduces two new hyperparameters σ and Δσ to the default ESN, where σ is the size of the stencil to define the local state and Δσ specifies the spatial distance of adjacent pixels included in the local state. Optimal values for all hyperparameters are determined by a grid search.

Convolutional Autoencoder
A convolutional autoencoder [33] is a special architecture of a feed forward network (FFN) with convolutional layers similar to convolutional neural networks (CNNs) [49]. Generally a CAE learns a representation of the training set D with the purpose of dimensionality reduction. For each pair Z i (X i , X ′ i ) ∈ D the CAE is trained to perform a nonlinear transformation from the input representation of X i to the output representation of X ′ i . Like CNNs a CAE is a partially locally connected feed forward network, which is typically composed of the following layers: Convolutional layers: Convolution of the input by a kernel sliding over the input. The number of rows and columns of the kernel are hyper-parameters, in this work they are set to be 3 × 3. Batch normalization layer: Normalization of the activations of the previous layer during training and for each batch. Batch normalization allows the use of higher learning rates, being computationally more efficient, and also acts as a regularizer [50]. Leaky ReLU [51] layer: Leaky version of a rectified linear unit (ReLU) [52], such that: Max pooling layer: Sample-based operation for discretization based on a kernel that slides over the input like the convolutional operator but only the maximum value of the kernel is passed to the next layer. Width and height of the kernel are hyperparameters (in this work 2 × 2). In contrast to the convolutional layer a pooling layer is not trainable. Dropout layer: Regularization method to prevent overfitting where during training some weights are set randomly to zero [53]. In this work the probability of setting the weights to zero is 0.05.
The eponymous part of the CAEs are the convolutional layers, a convolution of A (a ij ) ∈ K n×n with a kernel F (f ij ) ∈ K k×k , where k < n, is given by: with x, y ∈ 1, . . . , n. If i − x− or j − y exceeds the range of A zeropadding is applied [54]. In this work two architectures are used. The first one employed to reconstruct the data from noisy, blurred and inverse modelling data is illustrated in Figure 6. The architecture is the same for the tasks in Section 2.1.2 and Section 2.2.3 but the sizes of X and X ′ are different, with X, X ′ ∈ R 512×512 and X, X ′ ∈ R 100×100 , respectively. Due to the smaller input size, the data for the inverse modelling reconstruction is transformed into a latent space with the size of 25 × 25. The second architecture is sketched in Figure 7 and deals with the undersampled data reconstruction.

RESULTS
In the following both machine learning methods will be applied to two tasks: (i) Reconstructing electrical excitation waves from noisy blurred and under sampled data (Section 3.1) and (ii) Predicting electrical excitation from mechanical contraction (Section 3.2).

Recovering Complex Spatio-Temporal Wave Patterns From Impaired Observations
To benchmark both reconstruction methods, using ESNs and CAEs, we use time series generated by the BOCF model introduced in Section 2.1.1. The same data were used for both methods, consisting of 5,002 samples in the training data set, 2,501 samples in the validation data set, and 2,497 samples in the test data set. The sampling time of all time series equalled 10Δt 1. We considered nine cases of noisy data (with different noise levels), ten cases of (differently) blurred data and seven examples of (spatially) undersampled time series.
For the implemention of the ESN we used the software package easyesn [56]. To determine the optimal ESN hyperparameters a grid search is performed as described in [18] using the training and validation subsets of the data. This search consists of two stages: first, for each combination of the local states' hyperparameters σ and Δσ as listed in Table 2 a grid search is performed to find the optimal five hyperparameters of the ESN resulting in 37 sets of optimal hyperparameters. To make these grid searches more feasible, they were performed just for a single input patch (area covered by the stencil, see Figure 5) in the spatial center of the training set and thus not using the full spatial data.
In the second stage, for each of the 37 sets of optimal hyperparameters determined before (for each combination of σ and Δσ), an ESN is trained on a larger subset of the training data and not just on a single patch. Ideally, this step should be performed on the entire spatial domain of the training set, however, as we did not notice significant differences in the results when the ESNs were trained on a spatial subset of size 250 × 250 to speed up the training process. Following the same methodology as in [18], for each pixel from this spatial subset a single ESN is trained and then the obtained output matrices W out of these ESNs are averaged over all pixels. Compared to the procedure used in [18], the handling of boundary values has been changed. As for boundary pixels fewer adjacent pixels exist than for those inside, the creation of local states is obstructed, and boundary pixels require special treatment. In our previous work [18] individual ESNs have been trained for the boundary pixels using local states of lower dimensionality. In the following we use an alternative approach based on padding the boundary pixels by mirroring their values (motivated by the no-flux boundary conditions used). In this way, local states can be formally defined for boundary pixels in the same way as for inner pixels.
Next, the different optimal ESNs obtained for different stencils (σ, Δσ) were evaluated by comparing their performance on the validation subset. In this way optimal values for σ and Δσ were selected by choosing the combination (σ, Δσ) with the lowest ℓ 2 difference between the prediction and ground truth on the validation set. This process yields an ESN whose  Figure 9B boxes are horizontally shifted for better visibility. A tabular overview of the values can be found in Table 3 hyperparameters and weights are optimized to yield minimal ℓ 2 error. Finally, without training the network again on the entire training set, the optimal ESN found before is used to perform the prediction on the entire test set. As a pre-processing step, both the input and target data of the training, validation and test set are rescaled with min-max scaling, where the minimal and maximal value are determined over all pixels of the training set. The CAE was trained using the ADAM optimizer [57], implemented with TensorFlow [58] in version 2.3, with early stopping when the validation loss has not improved at least by 10 − 6 for 20 epochs. The learning rate was reduced by a factor of 0.2 when the loss metric stopped improving at least by 10 − 5 for ten epochs. Dropout was set to be 0.05 in all cases. As loss function the mean absolute error (MAE) was chosen: where N is the number of elements in the data set, X i the network output, X ′ i desired ground-truth and |.| stands for the absolute values. Figure 8 shows snapshots of the noisy input data, the corresponding ground truth, the outputs provided by the CAE and the ESN, respectively, and the absolute values of their prediction errors with respect to the ground truth. The evolution of the loss function during the training epochs is shown in Figure 9A. In all cases the error decreases and the training converges, but the duration of the training depends on the complexity of the case.

Noisy Data
Figures 9B shows a comparison of the performance of the CAE and the ESN for noisy data with nine different noise levels p 0.1, . . . , 0.9. While the mean absolute error of the CAE remains below 0.02, the reconstruction error of the ESN increases from 0.06 for p 0.1 to 0.18 for p 0.9, the associated ESN hyperparameter can be found in Table 4.

Blurred Data
To evaluate the performance of CAE and ESN for recovering full resolution (ground truth) data from blurred observations we consider nine cases where the radius m of Fourier low-pass filtering ranges from m 2 to m 18 (in steps of 2). Figure 10 shows snapshots of reconstructions of the u-variable of the BOCF model using CAE and ESN with filter parameters m 20 (A-F), m 14 (G-L), and m 8 (M-R). Similar to Figure 8 the errors are largest at fronts of the excitation waves, but in contrast to noisy images the performances of CAE and ESN differ not much for blurred data. This observation is also confirmed by a systematic comparison of the mean absolute errors of both methods for different manhatten distances m given in Figure 11A. The errors decrease with m because the larger m the less blurred are the input data of the CAE or ESN (for hyperparameter see Table 5). Figure 11B shows the evolution of the loss function during training of the CAE.  Figure 12 shows examples of data reconstructed from undersampled data. In total we considered seven cases of undersampling by 2 i pixels, where i ranges from i 1 to i 7. For i 1 input images have a resolution of X ∈ R 256×256 and for i 7, X ∈ R 4×4 . In all cases the desired output (ground truth) X ′ has a size of 512 × 512 pixels. The used hyperparameter for the ESN can be found in Table 6). Figure 13B shows the evolution of the corresponding loss function.

Predicting Electrical Excitation From Mechanical Contraction
Echo state networks as well as convolutional autoencoders have been trained with time series generated by the electromechanical model introduced in Section 2.2 to predict the membrane voltage u(x, t) Eq. (4) from the local contraction ΔA(x, t) given in Eq. (13). The sampling time of all time series equalled 6Δt 0.48. Since for periodically rotating spiral waves this cross-prediction task is quite straightforward we focus here on the much more demanding case of an example exhibiting spatio-temporal chaos (corresponding to atrial or ventricular fibrillation). Figure 14 shows at three instants of time snapshots of the observed contraction ΔA (first column), the ground truth of the voltage u (second column), the prediction of the CAE u CAE and its absolute error |u − u CAE | (third and fourth column) and the prediction of the ESN u ESN and the corresponding absolute error |u − u ESN | (columns five and six, respectively). Both, the ESN and the CAE were trained and tested with the same spatio-temporal time series (lengths of training, validation and test sets are 15,000, 2000 and 2000 samples, respectively). Hyperparameters of the ESN are N 600, α 0.5, ρ 1.1, ϵ 0.05, λ 5 · 10 − 3 , σ 7 and Δσ 1. To make the ESN more robust normally distributed noise with zero mean and a variance of 10 − 4 was added to the arguments of the activation function. As illustrated in Figure 14 both networks can solve the inverse problem and reconstruct the electrical potential field u from Eq. (4). However, the reconstruction of the CAE is more precise, which is particularly noticeable at the edges of the reconstructed electrical potential field. Considering the median of the MAE over the entire test data the ESN  approach achieves an error of 0.1963 ± 0.0260 while the median error of the CAE equals 0.0164 ± 0.0028.

CONCLUSION
Using synthetic data generated with conceptual models describing complex cardiac dynamics we have demonstrated possible applications of machine learning methods to complete and enhance experimental observations. It was shown that echo state networks as well as convolution autoencoders provide promising results, where the latter turned out to be the method of choice in terms of more faithful reconstructions. At this point, however, we would like to stress, that we didn't try to fully optimize the algorithms employed. One could, for example, increase the size of the ESNs used or extend and refine the grid search of hyperparameters. Also with the CAE several options exist to improve the performance even more. Instead of the MAE in the loss function one could use an adaptive robust loss function [59] or the Jensen-Shannon divergence [19]. The weights of the CAE could be optimized with a stochastic gradient descend approach instead of the ADAM algorithm [57]. But we expect such modifications would show only minor improvements (if at all) [60,61]. A comparison of the computing times with ESNs and CAEs is unfortunately not immediately possible. For the CAE the training time depends on the convergence, as illustrated in Figures 9A, 11B and 13B. In contrast to this, for ESNs training times depend strongly on the search for optical hyperparameters, especially the size of the reservoir, and this search is strongly dependent on the search space size and the number of parameters. For the task where electrical excitation is predicted from mechanical contraction our computations took 3,382 s on two Intel Xeon CPU E5-4620. While the CAE simulations have been run on GPUs the training and application of the ESN was on CPUs which makes a direct comparison difficult. Furthermore, the runtime of the programmes used is highly dependent on the libraries used and how well they have been adapted to special system architectures. In general we would estimate that in this work the effort to train and search hyperparameters for an ESN was a bit less demanding, in the sense of computational resources, compared to the training of the CAE. However, we would not consider the difference big enough to be an advantage for the ESN approach. Once trained both approaches need comparable execution times when applied to new data and executed on a CPU. In future work, using more realistic numerical simulations (and experimental data) such an optimization should be performed to achieve the best possible result for the intended medical application. Since here we used only data from conceptual models we refrained from fully optimizing the machine learning methods applied. The fact that already a straight-forward application of known algorithms and architectures provided very good results for the considered reconstruction tasks is very promising and encourages to address in future work extended tasks (including other variables, like calcium concentration, mechanical stress and strain, etc.) and reconstruction tasks with more realistic synthetic data (from 3D models, for example) combined with experimental measurements. Frontiers in Applied Mathematics and Statistics | www.frontiersin.org