mathematics Article Detection of Abnormal Cardiac Response Patterns in Cardiac Tissue Using Deep Learning Xavier Marimon 1,2,3, * , Sara Traserra 4 , Marcel Jiménez 4 , Andrés Ospina 5 and Raúl Benítez 1,3 1 Automatic Control Department, Universitat Politècnica de Catalunya (UPC-BarcelonaTECH), 08034 Barcelona, Spain 2 Bioengineering Institute of Technology, Universitat Internacional de Catalunya (UIC), 08195 Barcelona, Spain 3 Institut de Recerca Sant Joan de Déu (IRSJD), 08950 Barcelona, Spain 4 Faculty of Veterinary Medicine, Universitat Autònoma de Barcelona, 08193 Barcelona, Spain 5 Computer Science Department, Universitat Oberta de Catalunya (UOC), 08018 Barcelona, Spain * Correspondence: xavier.marimon@upc.edu or xmarimon@uic.es Abstract: This study reports a method for the detection of mechanical signaling anomalies in cardiac tissue through the use of deep learning and the design of two anomaly detectors. In contrast to anomaly classifiers, anomaly detectors allow accurate identification of the time position of the anomaly. The first detector used a recurrent neural network (RNN) of long short-term memory (LSTM) type, while the second used an autoencoder. Mechanical contraction data present several challanges, including high presence of noise due to the biological variability in the contraction response, noise introduced by the data acquisition chain and a wide variety of anomalies. Therefore, we present a robust deep-learning-based anomaly detection framework that addresses these main issues, which are difficult to address with standard unsupervised learning techniques. For the time series recording, an experimental model was designed in which signals of cardiac mechanical contraction (right and left atria) of a CD-1 mouse could be acquired in an automatic organ bath, reproducing the physiological conditions. In order to train the anomaly detection models and validate Citation: Marimon, X.; Traserra, S.; their performance, a database of synthetic signals was designed (n = 800 signals), including a wide Jiménez, M.; Ospina, A.; Benítez, R. range of anomalous events observed in the experimental recordings. The detector based on the LSTM Detection of Abnormal Cardiac neural network was the most accurate. The performance of this detector was assessed by means of Response Patterns in Cardiac Tissue experimental mechanical recordings of cardiac tissue of the right and left atria. Using Deep Learning. Mathematics 2022, 10, 2786. https://doi.org/ Keywords: deep learning; autoencoder; cardiac tissue; electrophysiology; electrostimulation; anomaly 10.3390/math10152786 detection; recurrent neural network; long short-term memory; CD-1 mouse model Academic Editor: Zhaobin Wang MSC: 92; 68T07; 92C55 Received: 6 July 2022 Accepted: 3 August 2022 Published: 5 August 2022 Publisher’s Note: MDPI stays neutral 1. Introduction with regard to jurisdictional claims in Anomaly detection is a broad field of research focused on the detection of abnormal published maps and institutional affil- patterns within a given dataset. Anomaly detectors have been traditionally used in the iations. industrial world for predictive maintenance to identify faults from sensor signals placed on motors and machinery [1]. Anomaly detection is also important in other fields, such as healthcare, structural health monitoring, security, surveillance and fraud detection, where Copyright: © 2022 by the authors. anomalies convey actionable and critical information. Apart from the use of temporal Licensee MDPI, Basel, Switzerland. signals, video sequences are also commonly used—for instance, in security surveillance This article is an open access article based on the detection of unusual movements in video surveillance footage [2]. Time series distributed under the terms and anomalies can be classified into two categories: point anomalies and structural anoma- conditions of the Creative Commons lies [3]. Point anomalies have historically been the most widely studied and correspond to Attribution (CC BY) license (https:// sample values of the signal statistically different from other signal values [4]. In contrast, creativecommons.org/licenses/by/ structural anomalies are defined as patterns presenting a morphology that differs from the 4.0/). behavior observed in the signal [5]. Mathematics 2022, 10, 2786. https://doi.org/10.3390/math10152786 https://www.mdpi.com/journal/mathematics Mathematics 2022, 10, 2786 2 of 21 Currently, machine learning plays an important role in the study of structural anoma- lies of time series [6]. When correctly labeled signals are available in normal and abnormal classes, the features in the time domain (statistical or morphological features) and in the frequency domain (spectral features) can be examined in order to better identify the fea- tures that distinguish between regular and abnormal conditions. These attributes can then be used to develop an anomaly classification algorithm based on labeled data through supervised learning. When labeled data are not available, the standard approach is to use unsupervised learning algorithms such as clustering to identify different behavior in the recordings. The resulting clusters can then be identified as normal and abnormal situations by an expert. This approximation allows for classifying correctly whether it is a normal or abnormal signal but does not give temporal information on the position in which the anomaly responsible for a deviation with respect to the normal signals has taken place. In the contraction of the heart chambers, the sinus node, which acts as a pacemaker, is responsible for generating an electrical stimulus that activates the two atria. Then, the electrical stimulus propagates to the atrioventricular node through the Purkinje fibers to contract the left and right ventricles [7]. A cardiac anomaly emerges when the conduction pathway generated at the sinus node is interrupted [8]. Therefore, it is necessary to monitor the electrical and mechanical activity of the heart to effectively identify indicative patterns of heart disease to ensure an early diagnosis. For this reason, in the field of translational cardiology, it is relevant to study the behavior of the heart at various scales and under different domains: electrical, mechanical, fluid and chemical (molecular). Different authors have used anomaly detectors in the electrical conduction system from the study of electrocardiogram (ECG) signals [9–24]. However, little attention has been paid to assess anomaly detection in the mechanical contractile signals at the tissue level. Therefore, the present study is focused on the detection of mechanical anomalies in cardiac tissue. Previous studies have used traditional machine learning methods to detect anomalies in electrical conduction from an ECG signal [9]. Unsupervised learning algorithms used for the detection of anomalies in the ECG signal include the dynamic time warping distance (DTW) [10,11], the Euclidean distance from an adaptive window discord discovery (AWDD) comparing the heartbeats to each other [12,13] and the k-means clustering algorithm [14]. In relation to supervised learning algorithms, the most commonly used are k-nearest neighbor (k- NN) [15], linear discriminant analysis (LDA) [16], quadratic discriminant analysis (QDA) [17], support vector machine (SVM) [18–22] and neural networks (NNs) [23,24]. The rapid development of graphics processing units (GPUs) versus central processing units (CPUs) has shown that deep learning can address artificial intelligence problems more quickly and reliably. Moreover, compared to traditional machine learning algorithms, they do not require cardiology experts to define the most relevant features to be extracted for identifying an anomaly in a given signal as deep learning networks are also able to extract features automatically. Thus, the characteristics extraction phase and the classification phases are integrated into the same deep learning network. The drawback, however, is that deep learning models require large training datasets to achieve high classification accuracy and are difficult to interpret by humans. Regarding deep learning algorithms, the most frequently used for the detection of ECG signal anomalies is long short-term memory (LSTM) [25]. Convolutional neural networks (CNNs) [25–29] can also be used after the conversion of the signal into a spectrogram image using wavelet transform or short-time Fourier transform [30]. The good results obtained in two-dimensional convolutional neural networks (2D CNNs) in computer vision applica- tions, such as AlexNet, VGG16 or ResNet networks, together with the emergence of the large datasets CIFAR100 and ImageNet have prompted the signal processing community to have a renewed interest in 1D CNNs [31]. In the detection of anomalies in health sciences, there is not always a sufficiently large pathological dataset to design an artificial intelligence algorithm to classify abnormalities based on the features extracted from normal and pathological data. Therefore, the present investigation focuses on studying anomalies by the construction of a model of the behavior Mathematics 2022, 10, 2786 3 of 21 under normal conditions. This approach is well suited for biomedical data as it requires few or no abnormal data. Thus, this article presents a framework that uses the one-class classification approach to construct a deep learning model based on normal data—i.e., a model trained with data under normal conditions. The input of the model is a signal to be analyzed, and the output is the reconstruction of the input signal based on the trained model. The resulting model will present a small reconstruction error when the input is a normal signal and a large reconstruction error when the input is a signal under abnormal conditions. The reconstruction error, defined as the difference between the input and its reconstruction, is indicative of the presence of an anomaly. The structure of this study is described next. First, Section 2 describes the biological preparation of atrial tissue samples (Section 2.1), the experimental setup for recording the mechanical contraction of cardiac tissue (Section 2.2), the creation of the synthetic dataset based on anomalous contractions observed in the experimental recordings (Section 2.3) and the theory and methodology used for the detection of anomalies with deep learning (Section 2.4). Next, in Section 3, the correct operation of the two anomaly detectors (LSTM and autoencoder) is validated with synthetic data (Section 3.1), and then, different application examples with real data are shown using the detector with the highest accuracy (LSTM detector). Finally, Section 4 is reserved for the critical discussion of the results. 2. Materials and Methods 2.1. Preparation of Atrial Tissue Samples For the present study, tissue samples obtained from the left and right atria of outbred albino CD-1 (Cluster of Differentiation 1) mice bred at the animal facility center of the Universitat Autònoma de Barcelona (UAB) were used. The protocol was approved by the Ethics Committee of the Autonomous University of Barcelona and complied with the Helsinki declaration ethical principles of the World Medical Association (WMA). 2.1.1. Mouse Preparation CD-1 mice aged between 15 and 18 weeks, weighing between 20 and 40 g, were selected. Housing was under controlled conditions: constant temperature (22 ± 2 ◦ C) and humidity (55 ± 10%), a 12-h light/dark cycle and ad libitum access to water and food. The animal was anesthetized by intraperitoneal injections of anesthetics, as shown in Table 1. Table 1. Composition of anesthetic preparations. Injection Composition Via 0.03 mL sodium heparin 50% First Intraperitoneal 0.3 mL sodium chloride heparin 0.9% 1 mg/kg medetomidine Second Intraperitoneal 75 mg/kg ketamine After intraperitoneal injection, the correct effect of anesthesia was confirmed by abol- ishment of the palpebral and swallowing reflexes. Animals were sacrificed by atlanto- occipital dislocation. 2.1.2. Atria Isolation Animals were placed in a supine decubitus position, and two thoracic incisions were made to expose the heart and lungs. The heart was extracted from the thoracic cavity and placed on a Sylgard-coated (Dow Corning. Midland, MI, USA) Petri dish containing carboxygenated Krebs solution (Table 2). Using a microscope, the left and right atria were identified and carefully dissected from the remaining cardiac tissue. Once atria had been isolated and prior to embedding in the bath, they were individually sutured to the force transductor using 3-0 silk sutures. placed on a Sylgard-coated (Dow Corning. Midland, MI, USA) Petri dish containing car- boxygenated Krebs solution (Table 2). Table 2. Composition of the carboxygenated Krebs–Ringer solution. Mathematics 2022, 10, 2786 Solution Composition Bubbles pH Temperature 4 of 21 115.48 mM NaCl 4.61 mM KCl Table 1.16 mM MgSO4 2. Composition of the carboxygenated Krebs–Ringer solution. Carboxygenated 5% CO 2 21.90 mM NaHCO3 7.4 37.5 °C Krebs 95% O2 Solution 1.14 mM NaHComposition 2PO4 Bubbles pH Temperature 2.50 mM CaCl 115.48 mM NaCl 2 10.10 mM4.61 mM KCl glucose 1.16 mM MgSO4 5% CO2 Carboxygenated Krebs 21.90 mM NaHCO3 7.4 37.5 ◦ C Using a microscope, the left and right atria were 95%identified O2 and carefully dissected 1.14 mM NaH2 PO4 from the remaining cardiac tissue. Once 2.50 mM atria had been isolated and prior to embedding CaCl2 in the bath, they were individually sutured 10.10 mM to the force transductor using 3-0 silk sutures. glucose 2.2. Mechanical Recordings 2.2. Mechanical Recordings The The mechanical mechanicalactivity activityofofthe CD-1 the CD-1atria waswas atria assessed using assessed an automatic using organorgan an automatic bath system (Letica Scientific Instruments, Barcelona, Spain). This organ bath system bath system (Letica Scientific Instruments, Barcelona, Spain). This organ bath system is is com- posed of four composed independent of four independent20 20mLmLchambers chambers filled filledwith withcarboxygenated carboxygenatedKrebs Krebs solution solution (5% CO 2 and 95% O2), as shown in Figure 1. (5% CO and 95% O ), as shown in Figure 1. 2 2 Figure 1. Figure Experimental setup 1. Experimental setup for for the the recording recording of of mechanical mechanical contractions contractions of of cardiac cardiac tissue. tissue. The temperature The temperaturewas wasregulated regulatedusing usinga Letica a Letica 13206 13206 thermostat thermostat (Letica (Letica Scientific Scientific In- Instruments, Barcelona, Spain), keeping a constant temperature of 37.5 ◦ C. Each atrium, struments, Barcelona, Spain), keeping a constant temperature of 37.5 °C. Each atrium, left left or right, was placed individually inside a chamber and sutured to the isometric force or right, was placed individually inside a chamber and sutured to the isometric force transducer. The force transducer, R = 360 Ω (Letica Scientific Instruments, Barcelona, Spain), transducer. The force transducer, R = 360 Ω (Letica Scientific Instruments, Barcelona, was adjusted to a mass of 0.2 g and allowed to stabilize for 15 min before starting the exper- iment. Data were registered using an analog-to-digital converter (ADC) (PowerLab 800, AD Instruments, Dunedin, New Zealand) with a sampling frequency of Fs = 1000 Hz. Electrical Stimulation Protocol Stimulation was performed using an analog stimulator (Grass S48, Astro-Med, Inc., West Warwick, RI, USA). Electrical stimulation of left atrial samples was applied using a current pulse of 40 mA for 0.4 ms at different stimulation frequencies of 1, 2, 3, 5 and 10 Hz. The frequency of electrical stimulation was related to the frequency with which the cardiac tissue contracted. Right atrial samples were not stimulated as the right atrium has a natural physiological pacemaker. Mathematics 2022, 10, 2786 5 of 21 2.3. Synthetic Dataset The experiment was based on the creation of a synthetic dataset of mechanical con- tractions of cardiac tissue. The synthetic dataset includes normal contraction and 7 types of different anomalies and includes a total of 100 signals for each category (i.e., n = 800 signals in total), every one with a signal length of 100,000 samples. For the creation of synthetic anomaly signals, different phenomena of anomalous contractions observed in the experi- mental recordings have been considered, including the presence of pulses with anomalous amplitude, duration, morphology or timing. The types of anomalies in the synthetic dataset are shown in Table 3. Normal signals were used to train the model and abnormal signals were used to validate both the anomaly classifier and the anomaly detector. In order to make the signals more realistic, noise was added to the signals following a normal dis- tribution with zero mean and a standard deviation σ = 0.005, robustly estimated from experimental signals. Table 3. Types of signals in the synthetic dataset. The x means not used parameter. Anomalies Repetitions Signal Type Signal Description Category Nanom Nrep 1. Normal signal Normal contraction Normal 6 x 2. Lower amplitude Phenomena of low amplitude contraction Anomalous 6 x Phenomena of large amplitude 3. Higher amplitude Anomalous 6 x contractions 4. Missing pulse Absence of contraction pulses Anomalous 1 x Periods of contraction with dynamics of 5. Slow pulse decay Anomalous 1 x slow decline Periods of contraction with dynamics of 6. Fast pulse decay Anomalous 1 x fast decline 7. Early or anticipated pulse Periods of anticipated contraction Anomalous 1 x 8. Pause/block Periods of absence of contraction pulses Anomalous 1 10 Different synthetic data generated at a stimulation frequency of 6 Hz, i.e., with the application of a stimulation pulse at approximately every 167 samples, are shown in Figure 2. 2.4. Architecture of Anomaly Detectors In the framework of mechanical signaling of cardiac tissue contraction, we are inter- ested in being able to localize the exact temporal position with a probability of having an anomaly. For this purpose, two detectors based on deep learning were implemented (Figure 3). As the first detector, a recurrent neural network (RNN) structure of LSTM type was used, whereas an autoencoder neural network was used as the second detector. The advantage of such anomaly detectors is that they can be trained with one-class data obtained from normal behavior without pathological conditions—that is, anomaly data are not necessary. Mathematics 2022, 10,2022, Mathematics x FOR10, PEER 2786 REVIEW 6 of6 21 of 22 * * ** ** Figure 2. Synthetic Figure signals 2. Synthetic representing signals representingeach anomalyatata stimulation each anomaly a stimulation frequency frequency of 6The of 6 Hz. Hz.asterisk The as- terisk(*)(*)indicates indicates a zoom a zoom intointo figure. figure. The triangle The triangle indicates indicates the anomaly the anomaly location.location. 2.4. Architecture of Anomaly Detectors In the framework of mechanical signaling of cardiac tissue contraction, we are inter- ested in being able to localize the exact temporal position with a probability of having an anomaly. For this purpose, two detectors based on deep learning were implemented (Fig- ure 3). As the first detector, a recurrent neural network (RNN) structure of LSTM type was used, whereas an autoencoder neural network was used as the second detector. The ad- vantage of such anomaly detectors is that they can be trained with one-class data obtained Mathematics 2022, 10, x FOR PEER REVIEW 7 of 22 Mathematics 2022, 10, 2786 7 of 21 from normal behavior without pathological conditions—that is, anomaly data are not nec- essary. A Signal D.1 x MAE Threshold - u |u | E Anomaly Subsequence LSTM y + E > T hr location B Signal D.2 x MAE Threshold - u |u | E Anomaly Subsequence Autoencoder y + E > T hr location Figure3. Figure 3. Anomaly Anomaly detectors detectorsbased basedon ondeep deeplearning. learning.(A)(A) Anomaly detector Anomaly based detector on an based on LSTM net- an LSTM work. (B) Anomaly detector based on an autoencoder. network. (B) Anomaly detector based on an autoencoder. Inthe In theuseuseofofautoencoders autoencodersand andLSTM LSTMnetworks networksas asanomaly anomalydetectors, detectors,when whennormal normal dataare data areintroduced, introduced,the thedetector detectorcan canreconstruct reconstructthetheinput, input,and andthen thenthetheerror errorbetween betweenthe the inputsignal input signalandandthetheoutput outputsignal signalisissmall. small.Conversely, Conversely,when whendatadatacontaining containinganomalies anomalies areentered are enteredinto intothe thenetwork, network,the thenetwork networkcannot cannotreconstruct reconstructthetheinput inputsignal signaland andthe theerror error becomeslarger. becomes larger.The TheMATLAB MATLABr2022a r2022asoftware software(Mathworks, (Mathworks,Boston,Boston,MA, MA,USA) USA)waswasused used to program the code for detecting anomalies, with which custom code to program the code for detecting anomalies, with which custom code was developed and was developed and specificlibraries specific librarieswere werealso alsoused, used,such suchas as Signal Signal Processing Processing Toolbox Toolboxand andMachine MachineLearning Learning Toolbox.AAcomputer Toolbox. computer with with an an 11th-generation 11th-generation IntelIntel CoreCore i9 processor i9 processor (CPU), (CPU), 32 GB32of GB RAM of RAM and anand Nvidiaan Nvidia GeForce GeForce GTX X GTX Titan Titan X graphic graphic card (GPU) card (GPU) was used was used in order in order to train to train the the deep deep learning learning models.models. The detailed The detailed steps steps necessary necessary to develop to develop the anomaly the anomaly detector detector based on an autoencoder based on an autoencoder(D.1) or an LSTM (D.1) neuralneural or an LSTM network (D.2) are network (D.2)described as follows are described (see as follows Supplementary Material S1 and (see Supplementary Material S1 and S2): S2): (a) (a)Organization Organizationof ofdata data First First of mechanicalcontraction of all, mechanical contractiondata datacorresponding corresponding to to normal normal andand anomalous anomalous cat- categories egories were were obtained. obtained. (b) (b)Normalization Normalizationof ofdata data Signals Signals were then normalized;centered were then normalized; centered and andstandardized standardized with thethe with z-score transforma- z-score transfor- tion, defined as follows: mation, defined as follows: x − µ( x ) xc = z-score( x ) = (1) σ( x) x − μ ( x) = z-score(data, where xc represents centered andxstandardized x ) = µ(x) is the mean of signal x and σ(x) (1) is c σ ( x) the standard deviation of signal x. (c) Construction of the subsequences matrix where xc represents Signals centered were arranged as aand standardized subsequences data,X, matrix, μ(x) an is thematrix LxK mean obtained of signal by x and σ(x) arrang- is the standard deviation of signal x. ing K signal windows of size L samples, as in Figure 4. The number of windows, K, was defined(c) as follows: Construction of the subsequences matrix K = S − matrix, Signals were arranged as a subsequences L X, an LxK matrix obtained by (2) ar- ranging K signal windows of size L samples, as in Figure 4. The number of windows, K, where K is the number of windows, S is the stride between windows and L is the window size. was Subsequences defined as follows: were organized in an LxK subsequences matrix, X, as follows: K =S− L (2) x c (1) x c (1 + S ) x c (1 + 2 · S ) · · · x c (1 + K · S )   where K is the number  x (of  2) windows, xc (2 + SS) is the xc (2stride + 2 · Sbetween ) · · · windows xc (2 + K ·and S) L is the window size.   X =  x c (3) x c (3 + S ) x c (3 + 2 · S ) . . . x c (3 + K · S )  (3)    . .. .. ..   .  . . . .   xc ( L) xc ( L + S) xc ( L + 2 · S) . . . xc ( L + K · S) where xc is the normalized signal with N samples; 1 × N. Mathematics 2022, Mathematics 10, x2786 2022, 10, FOR PEER REVIEW 8 8ofof 22 21 Figure 4. Creation of the subsequence matrix, X, from the contraction time series, x. Figure 4. Creation of the subsequence matrix, X, from the contraction time series, x. Subsequences (d) Training with were normalorganized data in an LxK subsequences matrix, X, as follows: Then, the autoencoder or LSTM network was trained with the sequences of the training  xc (1) data of the normal signal xc (1 + S) xc (1 + 2 ⋅ S)  xc (1 + K ⋅ S)  category. (e) Reconstruction  with test data  x(2) xc (2 + S) xc (2 + 2 ⋅ S)  xc (2 + K ⋅ S)  Then, validation or test data with the corresponding subsequences were taken, and the model trained X=  xcnormal with (3 + S)wasxused (3) xsignals (3 +to2 ⋅reconstruct S)  xthe + K ⋅ S)  (3signal. (3)  error  c c c (f) Reconstruction The reconstruction  error between    the input signal, x, and the reconstructed  signal, y  x ( (LSTM output), was calculated: L) x ( L + S) x ( L + 2 ⋅ S )  x ( L + K ⋅ S)  c c c c  Ei = |yi − xi | (4) where xc is the normalized signal with N samples; 1 × N. Then, the mean absolute error (MAE) was defined as follows: (d) Training with normal data Then, the autoencoder or LSTM network1wasNtrained with the sequences of the train- N i∑ ing data of the normal signal category. MAE = · Ei (5) =1 (e) Reconstruction with test data whereThen, xi isvalidation the input or test data signal and with yi is the the corresponding reconstructed signalsubsequences produced wereby taken, and the LSTM the model trained with network or the autoencoder. normal signals was used to reconstruct the signal. (f) (g)Reconstruction error Fitting the reconstruction error to a Gaussian distribution The The reconstruction error was reconstruction error between the input adjusted signal, x,normal to a univariate and thedistribution reconstructed signal,of by means y (LSTM output),likelihood the maximum was calculated: estimation (MLE): EiE∼=N y (µ̂,−σ̂x) (6) (4) i i i where µ̂ is the estimated mean and σ̂ is the estimated standard deviation. Then, the mean According absolute to [32], error (MAE) for simplicity, it is was defined as follows: recommendable to fit a Gaussian distribution, which is used to assess the likelihood of anomalous behavior. The maximum likelihood estimation (MLE) was calculated as follows [33]. 1 N (h) Measure of error distance MAE = ⋅ Ei N  i =1 (5) When error is found at the tail of a normal distribution, it is probable that an anomaly has occurred, so the rarity of the event can be measured according to how far the location of where xi is the error Eithe inputthe is from signal anddistribution, normal yi is the reconstructed signal N(µ̂, σ̂). For produced distance by theasLSTM calculation, net- proposed work or the autoencoder. by different authors [5,32], measurement of the Mahalanobis distance was used since it is a (g) Fitting measure the reconstruction of the distance between anerror errorto a Gaussian point, Ei , and distribution a statistical distribution. Therefore, The way, in some reconstruction error was the Mahalanobis adjusted distance to a univariate is being used as annormal anomaly distribution by means score. In this study, of thethe maximum univariate likelihood estimation Mahalanobis (MLE): distance, dM, was calculated as follows: 2 1 dM= dist( E, N N Ei(µ̂, σ̂ ))( μˆ=, σ(ˆ E ) − µ̂) · (7) (6) σ̂2 (i) Envelope extraction where μ̂ is the estimated mean and σˆ is the estimated standard deviation. Mathematics 2022, 10, 2786 9 of 21 The signal envelope, ζ n , was calculated by interpolating the Mahalanobis distance, dM, with a spline function, where the local maxima separated by at least np = 30 samples were selected. (j) Calculation of the threshold of anomalous samples Using the Kolmogorov–Smirnov one-sample normality test, it was confirmed that the validation error signals did not follow a normal distribution. Therefore, the threshold, Thr, for the detection of anomalous samples was defined as the upper quartile, Q3, of the envelope of the validation error plus five-fold the interquartile range (IQR): Thr = Q3(ζ n ) + 5 · IQR(ζ n ) (8) (k) Temporal localization of the anomalous samples An anomalous sample is considered when the value of the envelope of the Maha- lanobis distance, ζ n , is greater than the threshold, Thr: ζ n ≥ Thr (9) The entire methodological framework described from sections (a) to (k) is represented graphically in the Supplementary Materials. 2.4.1. Detector Based on an LSTM Network (D.1) The neural network used in the study was based on simple long short-term memory cells (vanilla LSTM), the output gate of which is described by the following recurrence functions [34]: ht = ot tanh(ct ) (10) ot = f t ct−1 + it cet where ht and ht −1 are the new state and the previous state of the hidden layer, respectively; ct is the current state of the cell and cet is the cell state candidate. The symbol in the previous equations refers to element-wise multiplication (Hadamard product). LSTM cells use additional gates for controlling which information of the hidden cell arrives to the input and to the following hidden state. This allows more efficient Mathematics 2022, 10, x FOR PEER REVIEW 10 of 22 network learning of the long-term relationship of data. The input gate, it , controls the update level of the cell state, ct . The forget gate controls the extent of forgetfulness of the cell state. Finally, the output gate controls the level of the cell state that has been added to to the hidden state, ht (see Figure 5). Moreover, the cell state candidate adds information the hidden state, ht (see Figure 5). Moreover, the cell state candidate adds information to to the cell state. Additionally, LSTM cells are able to learn long-term dependencies. the cell state. Additionally, LSTM cells are able to learn long-term dependencies. Figure5.5.Control Figure Controlgates gatesof ofaasimple simpleLSTM LSTMcell cell(vanilla (vanillaLSTM). LSTM). Table Table44summarizes summarizesthe thegovernance governanceequations equationsofofaacell cellininthe theLSTM LSTMnetwork network[34]. [34]. Table 4. Governance equations of an LSTM cell. Part Equation Hidden state ht = ot  tanh( c t ) Cell state (cell output) ct = ft  ct −1 + it  ct Mathematics 2022, 10, 2786 10 of 21 Table 4. Governance equations of an LSTM cell. Part Equation Hidden state ht = ot tanh(ct ) Cell state (cell output) ct = f t c t −1 + i t cet Cell state candidate cet = tanh(Uc · xt + Wc · ht−1 + bc ) Input gate it = σ(Ui · xt + Wi · ht−1 + bi ) Forget gate f t = σ (U f · x t + W f · h t −1 + b f ) Output gate ot = σ(Uo · xt + Wo · ht−1 + bo ) In the table, xt is the current input gate; Ui , Uf and Uo are the weights of the hidden layer for the input gate, forget gate and output gate, respectively; Wi , Wf and Wo are the weights of the hidden state for the input gate, forget gate and output gate, respectively; bi , bf and bo represent the bias for the input gate, forget gate and output gate, respectively; ht is the new state (hidden); ht −1 is the previous state (hidden); σ is the logistic sigmoid activation function and tanh is the hyperbolic tangent function. Sigma, σ, indicates that the gates’ activation function is defined as a logistic sigmoid function (logsig) as follows: 1 f ( xt ) = σ ( xt ) = (11) 1 + e− xt In a vanilla LSTM, the sigmoid function is used as the main activation function for the input gate, the output gate and the forget gate because it outputs a value between 0 and 1, and it can either let complete information flow or not throughout the gates. Additionally, to overcome the vanishing gradient problem caused by the derivative of the activation function, which is very close to zero in the gradient descent methods, we needed a function whose second derivative can be sustained for a long time before going to zero. Then, in the cell state candidate recurrence formula, the hyperbolic tangent function (tanh) ensures that the gate values are scaled between −1 and 1, thus increasing stability during gradient descent. Bidirectional LSTM (BiLSTM), which is simply the union of two conventional LSTM cells, allows cells to retrieve sequence information in both ways, either backward or Mathematics 2022, forward, to REVIEW 10, x FOR PEER continue the information transmission [35]. Figure 6 compares a conventional 11 of 22 LSTM cell with a bidirectional LSTM cell. Figure 6. Comparison of a unidirectional LSTM cell with a bidirectional cell. (A) Unidirectional Figure 6. Comparison of a unidirectional LSTM cell with a bidirectional cell. (A) Unidirectional LSTM cell (LSTM). (B) Bidirectional LSTM cell (BiLSTM). LSTM cell (LSTM). (B) Bidirectional LSTM cell (BiLSTM). For the detector based on a LSTM network, a neuronal network with 9 layers and 8 connections was developed (see Figure 7). In addition to LSTM cells, rectified linear units (RELUs), fully connected (FC) layers and one last regression layer (RL) were used. RELU activation functions are responsible for setting the negative values between layer and layer to zero; they are not important for Mathematics 2022, 10, 2786 11 of 21 Mathematics 2022, 10, x FOR PEER REVIEW 12 of 22 For the detector based on a LSTM network, a neuronal network with 9 layers and 8 connections was developed (see Figure 7).  Sequence input  BiLSTM 1  Bilateral Long Short-Term Memory RELU 1 yt  BiLSTM BiLSTM 2 ht σ LSTM  Layers ht RELU 2 LSTM  BiLSTM 3 xt  RELU 3  FC  RL 7. Architecture Figure 7. Architectureofoflong longshort-term short-termmemory memory(LSTM) network (LSTM) created network for anomaly created detection. for anomaly detection. TableIn5.addition to LSTM defining Hyperparameters cells, rectified linear units the architecture (RELUs), of the fully connected (FC) layers LSTM network. and one last regression layer (RL) were used. RELU activation functions are responsible Layer Name forType setting the Activations negative values between layerLearnables States for and layer to zero; they are not important learning as they do not have saturation regions and they avoid stagnation during some Sequence 1. Input phases 1 - - inputof the network training, unlike other activation functions. The hyperparameters that define the architecture of the LSTM Inputnetwork weights:used for× detecting 2080 1 anomalies in the signals of the mechanical contraction of cardiac tissue are summarized in Table Hidden states: 520 5. The definition of × 1 2. BiLSTM 1 BiLSTM 520 Recurrent weights: 2080 × 260 the best hyperparameters was achieved using an empirical study. Cell states: 520 × 1 Bias: 2080 × 1 For the training of the LSTM network, as these are real values, the loss function has 3. RELU 1 RELU been 520 squared defined as a mean - error (MSE) function, which calculates - the error between Input weights: 1280 × 520 network predictions and target values for the regression task: Hidden states: 320 × 1 4. BiLSTM 2 BiLSTM 320 Recurrent weights: 1280 × 160 1 N Cell states: 320 × 1 Bias: J (W1280 ) = × 1· ∑ (yt − yet )2 (12) 5. RELU 2 RELU 320 - N t =1 - Input weights: 2080 × 320 where yt represents ground truth data and yet is the output value of theHidden states: 520 × 1 neural network. 6. BiLSTM 3 BiLSTM 520 Recurrent weights: 2080 × 260 Cell states: 520 × 1 Bias: 2080 × 1 7. RELU 3 RELU 520 Weights: 1 × 520 8. Fully Connected FC 1 - Bias: 1 × 1 9. Regression RL 1 Mean squared error - Layer For the training of the LSTM network, as these are real values, the loss function has been defined as a mean squared error (MSE) function, which calculates the error between Mathematics 2022, 10, 2786 12 of 21 Table 5. Hyperparameters defining the architecture of the LSTM network. Layer Name Type Activations Learnables States 1. Input Sequence input 1 - - Input weights: 2080 × 1 Hidden states: 520 × 1 2. BiLSTM 1 BiLSTM 520 Recurrent weights: 2080 × 260 Cell states: 520 × 1 Bias: 2080 × 1 3. RELU 1 RELU 520 - - Input weights: 1280 × 520 Hidden states: 320 × 1 4. BiLSTM 2 BiLSTM 320 Recurrent weights: 1280 × 160 Cell states: 320 × 1 Bias: 1280 × 1 5. RELU 2 RELU 320 - - Input weights: 2080 × 320 Hidden states: 520 × 1 6. BiLSTM 3 BiLSTM 520 Recurrent weights: 2080 × 260 Cell states: 520 × 1 Bias: 2080 × 1 7. RELU 3 RELU 520 Weights: 1 × 520 8. Fully Connected FC 1 - Bias: 1 × 1 9. Regression Layer RL 1 Mean squared error - Numerical optimization that minimizes the loss of function by means of successive iterations, l, starting from a random initial value was used. In our case, for training the LSTM-type neural network, the adaptative moment estimation (Adam) algorithm was used; this is a stochastic gradient descent method based on adaptative estimation of the learning rates of moments of first, ml , and second order, vl , to further improve convergence time [36]. Computationally, it is a very efficient method as it uses little memory. Classical gradient descent methods use a single learning rate for all parameters, whereas the Adam method aims to improve network training by using adaptative training rates that can be adapted automatically to the loss function that is being optimized. When updating an iteration with the Adam method, a moving average is used such as the following: Update formula : Wl +1 = Wl − √α·m̂l v̂l +ε ml m̂l = (1− β 1 l ) ml = β 1 · ml −1 + (1 − β 1 ) · ∇ J (Wl ) (13) vl v̂l = (1− β 2 l ) vl = β 2 · vl −1 + (1 − β 2 ) · [∇ J (Wl )]2 where ml and vl are first- and second-order moments of the current state, respectively; Wl , Wl+1 and Wl −1 are the weight vectors of the current, previous and next states, respectively; ∇J is the gradient of the loss function and β1 and β2 are the exponential decay rates of the first and second moments, respectively. If gradients contain mostly noise, updating the weights with the moving average allows the moving average of the gradient to become smaller, and therefore the weight updates also become smaller. The moments indicated by the symbol (ˆ) are the moments with correction of the bias that appears at the beginning of training. The use of a regularizer allows for applying penalties on layer parameters during optimization. These penalties add to the loss function that optimizes the network. In the present study, a regularization penalty, L2 , was used. Additionally, adding a regularization term to the loss function weights, J(W), is one way to reduce overfitting [37,38]. Then, the loss function after L2 regularization was expressed as follows: JL2 (Wl ) = J (Wl ) + L2 · W 2 (14) where JL2 (W) is the loss function with L2 regularization, W is the weight vector of the current state and L2 is the penalty of L2 regularization. Mathematics 2022, 10, 2786 13 of 21 Mathematics 2022, 10, x FOR PEER REVIEW 14 of 22 All parameters used in the training options of the LSTM neural networks are sum- marized in Table 6. Parameters related to the Adam optimization algorithm are based on [36]. Table 6. Training parameters for the LSTM network. Table 6. Training parameters for the LSTM network. Description Parameter Name Value Optimization algorithm Description Solver Parameter name Name AdamValue First moment rate Optimization algorithm Gradient Solver name decay factor β Adam 1 = 0.9 First moment rate rate Second moment Squared Gradient gradient decay factor decay factor β1β=2 =0.9 0.999 Second Epsilonmoment rate Squared Epsilon gradient decay factor β2ε==0.999 10−8 Epsilon Maximum number of epochs Epsilon Max epochs ε =3010−8 Maximum number of epochs Max epochs 30 Mini-batch size Mini batch size 10 Mini-batch size Mini batch size 10 Shufflethe Shuffle thedata data Shuffle Shuffle Once Once Initiallearning Initial learningrate rate Initial Initial learnrate learnrate α0α=0 =0.001 0.001 L regularization L2 regularization L2 L2 regularization regularization L2L=2 =0.0001 0.0001 2.4.2.Autoencoder-Based 2.4.2. Autoencoder-BasedDetector Detector(D.2) (D.2) Autoencoders Autoencodersare areunsupervised unsupervisedlearning learningalgorithms algorithmswith withunlabeled unlabeleddata datathat thatare are based basedononneural networks neural networks and consist and of two consist parts: of two an encoder parts: an encoder andand a decoder. a decoder. The The function func- oftion the encoder is to compress of the encoder data. Fordata. is to compress givenFor input datainput given x = {xdata1 , x2 ,xx3=, .{x . .1,, xx2N, },x3the , …,encoder xN}, the compresses the N-dimensional encoder compresses inputh{h input in a code the N-dimensional in 1a, hcode 2 , h3 ,h{h . . .1,, hh2M , h},3,where …, hM}, Mwhere < N, named M < N, latent named representation, which includes latent representation, most of the which includes input most of theinformation. input information. The architecture of a The architec- basic ture autoencoder is shown inisFigure of a basic autoencoder shown8.in Figure 8. h = fW ( e ⋅ x + be) Encoder Decoder be Hidden bdT + f() + x f() x + h + T W e W d x = fW ( T d ⋅ h + bdT ) Figure8.8.Basic Figure Basicautoencoder autoencoderstructure. structure. Thenetwork The networkof of thethe autoencoder autoencoder can can be mathematically be mathematically represented represented by theby the corre- correspond- sponding encoding and decoding equations (see Figure ing encoding and decoding equations (see Figure 8): 8): h =h f=(W f (e W ⋅ xb+ ·x+ e) b ) e e (15) xe = f (Wd T · hT + bdT ) T (15) x = f (W ⋅ h + b ) d d where x and xe are the input and output signal vectors of the autoencoder N × 1, respectively; h where x and x is the hidden  areMthe state × input and W 1; We and output signal vectors of the autoencoder N × 1, respec- d are the weight vectors of the encoder and decoder tively; h is the hidden state M × 1; W M × N, respectively; and be and bd are the bias e and Wd areof the the weight encodervectors and theofdecoder the encoder M ×and 1, decoder respectively.M × N, respectively; and b e and bd are the bias of the encoder and the decoder M × 1,In respectively. order to avoid the vanishing gradient problem, the hyperbolic tangent function (tanh) In andorder logistic sigmoid to avoid thefunction vanishing(logsig) are commonly gradient used. problem, the Both are s-shaped hyperbolic curves; tangent function the only difference is that logsig lies between 0 and 1 whereas tanh lies (tanh) and logistic sigmoid function (logsig) are commonly used. Both are s-shaped between 1 and −1. Incurves; this study, the logistic sigmoid activation function (logsig) has been used for the only difference is that logsig lies between 0 and 1 whereas tanh lies between 1 both the and −1. In this study, the logistic sigmoid activation function (logsig) has been used for both the encoder and the decoder because it is commonly used for models where proba- bilities have to be predicted. Since the probability of anything exists only between the Mathematics 2022, 10, 2786 14 of 21 encoder and the decoder because it is commonly used for models where probabilities have to be predicted. Since the probability of anything exists only between the range of 0 and 1, the tanh function does not have a clear relationship to modeling probability. For the definition of the loss function of the autoencoder, since the input, x, has a real value, a mean squared error sparse (msesparse) function was used. This function calculates the error between the input and the output value of the autoencoder as well as adding two terms of regularization: 1 N N i∑ J (W ) = · ( xi − xei )2 + Ω L2 + ΩKL (16) =1 where N is the number of samples of the input signal vector; xi and xei are the values of the sample at the input and output of the autoencoder, respectively; ΩL2 is the regularization term L2 and ΩKL is the sparsity regularization term. In Equation (11), it can be observed that two regularization terms are used with the objective of avoiding overfitting. L2 regularization affects weights, whereas dispersion regularization, W, affects the activity of the hidden layers, h. To minimize loss function, the stochastic gradient descent (SGD) was used. The regularization term L2 or ridge regularization is expressed as the i-th sum of the squares of the weights, W, multiplied by the weight decay parameter, λ, which tends to decrease the magnitude of weights and helps to avoid overfitting [39,40]: s 1 L s l l +1 2 Ω L2 = λ · · ∑ ∑ ∑ (Wjil ) (17) 2 l =1 i =1 j =1 where λ is the weight decay parameter, L is the number of layers, sl is the number of units in the layer l and Wji is the weight between unit j of layer l and unit i of layer l + 1. It is interesting to be able to achieve a mean activation, ρ, of each hidden neuron, j, close to the dispersion parameter; that is, ρ̂ j = ρ. Typically, the dispersion parameter, ρ, takes a small value close to zero. In the present study, ρ = 0.1 was selected. In order to meet this restriction, activations of the hidden unit should be mostly close to 0. To achieve this, the dispersion regularization term, ΩKL , penalizes when ρ̂ j is significantly deviated from ρ on the basis of the Kullback–Leibler divergence (KL) concept [41]: M ΩKL = β · ∑ DKL (ρ k ρ̂ j ) (18) j =1 where M is the number of hidden neurons; β is the dispersion regularization parameter; ρ̂ is the vector of the mean activity of the hidden layer, h; and ρ is the sparsity proportion parameter. The Kullback–Leibler divergence (KL) is a measure between two statistical distributions— specifically the distance between two distributions of Bernoulli random variables, one with a mean ρ and the other with a mean ρ̂. The dispersion regularization parameter, β, regulates the weight of the dispersion term, ΩKL. The parameters used in the training options of the autoencoder are summarized in Table 7. The values of the regularization parameters and loss function selection are based on [39]. Mathematics 2022, 10, x FOR PEER REVIEW 16 of 22 Mathematics 2022, 10, 2786 15 of 21 Table 7. Training parameters for the autoencoder detector. Description Parameter Name Value Number of hidden Table 7. Training nodes for the autoencoder parameters Hiddendetector. size 200 L2 regularization L2 weight regularization λ = 0.0001 RegularizationDescription of dispersion Parameter Sparsity Name regularization Value β = 0.0001 Sparsity Numberproportion of hidden nodes Sparsity Hidden sizeproportion 200ρ = 0.1 L2 regularization Maximum number of epochs L2Max weight regularization epochs λ =200 0.0001 Regularization of dispersion Sparsity regularization β = 0.0001 Loss function Loss function msesparse Sparsity proportion Sparsity proportion ρ = 0.1 Activation Maximumfunction number of(encoder) epochs Encoder Max epochstransfer function 200logsig Activation function (decoder) Loss function Decoder Loss functiontransfer function logsig msesparse Activation function (encoder) Encoder transfer function logsig 3. Activation Results function (decoder) Decoder transfer function logsig 3.1. Synthetic Data 3. Results The intersection over union (IOU) and Sørensen–Dice index (SDI) criteria are com- 3.1. Synthetic monly used toData calculate the accuracy of an anomaly detector in synthetic data. The intersection over union (IOU) and Sørensen–Dice index (SDI) criteria are com- monly used to calculate the intersection ROI D ∩ accuracy of an anomaly detector ROIGTin synthetic data. IOU = = ⋅ ∈ [0,1] (19) union intersection ROI ROI D D∩∪ ROI ROI GT GT IOU = = · ∈ [0, 1] (19) union ROID ∪ ROIGT where ROID is the anomalous region detected by the anomaly detector and ROIGT is the where ROID is the anomalous region detected by the anomaly detector and ROIGT is the ground truth anomaly region. ground truth anomaly region. The The region region of of interest (ROI) corresponds interest (ROI) correspondsto tothe thecardinalities cardinalitiesofofthe thetwo two regions regions (i.e., (i.e., the the number of samples in each set). Figure 9 shows the practical calculation of the number of samples in each set). Figure 9 shows the practical calculation of the intersection inter- section over(IOU) over union (IOU)the unionwhen when theof region region of interest interest of the truth of the ground ground (ROItruth (ROIGT) and the GT ) and the region of region of interest detected by the anomaly detector (ROI interest detected by the anomaly detector (ROID ) are known. D ) are known. Figure Figure9.9.Method Methodofof obtaining obtainingthethe intersection over intersection union over (IOU) union in signals. (IOU) Colors in signals. boxesboxes Colors in the in the formula indicates the different regions of interest (ROI), according to the figure on the left. formula indicates the different regions of interest (ROI), according to the figure on the left. Alternatively, Alternatively,the theIOU IOUcan canbebewritten writtenusing usingthe thedefinition definitionofoftrue truepositive positive(TP), (TP),false false positive positive(FP) (FP)and andfalse falsenegative negative(FN), (FN),as asfollows: follows: intersection TP TP IOU = intersection = = IOU = ∈∈[0, [0,1] 1] (20) (20) union union TP++FP TP FP++FN FN This index is known by several other names, especially F1 score. In a similar way, the Sørensen–Dice index (SDI) or F1 score is defined as follows: 2 · TP SDI = ∈ [0, 1] (21) 2 · TP + FP + FN This index is known by several other names, especially F1 score. In a similar way, the Sørensen–Dice index (SDI) or F1 score is defined as follows: Mathematics 2022, 10, 2786 16 of 21 2 ⋅ TP SDI = ∈ [0,1] (21) 2 ⋅ TP + FP + FN ItItdiffers differs from the theIOU IOUindex indexwhich whichonly counts only true true counts positives once in positives bothinthe once nu-the both merator and the denominator. Both indices, IOU and Sørensen–Dice, range in numerator and the denominator. Both indices, IOU and Sørensen–Dice, range in the closed the closed intervalbetween interval between00 and and 1. 1. They They can can be beinterpreted interpretedasasa asimilarity similaritymeasure measureover sets. over sets. Figure 10 shows the detection accuracy using the anomaly detectors Figure 10 shows the detection accuracy using the anomaly detectors based based on onthethe autoencoderand autoencoder andLSTM LSTM for for one one particular particularcase caseofofaasynthetic syntheticsignal: signal:higher amplitude. higher amplitude. Figure 10. Detection of anomalies in the synthetic signal type: higher amplitude. (A) Autoencoder- Figure based 10. Detection detector. of anomalies (B) LSTM-based in the synthetic signal type: higher amplitude. (A) Autoencoder- detector. based detector. (B) LSTM-based detector. Table 8 shows the accuracy obtained by the classifiers based on the autoencoder and the Table 8 shows LSTM-type thenetwork; neural accuracyweobtained byuse decided to thethe classifiers based the SDI to evaluate on accuracy. the autoencoder Syn- and the LSTM-type neural network; we decided to use the SDI to evaluate thetic data which model seven types of anomalies representative of the contractilethesignals accuracy. Synthetic data which of heart tissue, model excluding seven normal typeswere signals, of anomalies representative used for evaluation of the (see Table 5). contractile signals of heart tissue, excluding normal signals, were used for evaluation (see Table 5). Table 8. Anomaly detectors’ accuracy. Table 8. Anomaly detectors’ accuracy. Sørensen–Dice Index (SDI) Signal Type Sørensen–Dice Index (SDI) Autoencoder LSTM Signal Type 2. Lower amplitude 0.55 Autoencoder 0.83 LSTM 3. Higher amplitude 2. Lower amplitude 0.55 0.71 0.99 0.83 4. 3. Missing pulse Higher amplitude 0.710.63 0.91 0.99 5. 4. Slow pulse Missing decay pulse 0.630.67 0.95 0.91 5. Fast 6. Slowpulse pulse decay decay 0.670.61 0.95 0.88 6. Early 7. Fast pulse decay or anticipated pulse 0.610.51 0.88 0.77 7. Early or anticipated pulse 0.51 0.77 8. Pause/block 0.61 0.9 8. Pause/block 0.61 0.9 The LSTM neural-network-based framework presented the best accuracy. It is also The noted LSTM that neural-network-based the early pulse signal was the framework most difficultpresented the best signal to detect accuracy. in both It is also approaches. noted that the early pulse signal was the most difficult signal to detect in both approaches. 3.2. Experimental Data This section describes the application of the anomaly detection framework based on deep learning to experimental data. In the previous section, it was validated with synthetic data that LSTM offers better accuracy in detecting anomalies (see Section 3.1). Therefore, the LSTM-type anomaly detector was tested on experimental signals of mechanical contrac- tion of left atrium cardiac tissue (AE) stimulated at successive frequencies of 1, 2, 4, 6, 8 3.2. Experimental Data This section describes the application of the anomaly detection framework based on deep learning to experimental data. In the previous section, it was validated with syn- Mathematics 2022, 10, 2786 thetic data that LSTM offers better accuracy in detecting anomalies (see Section 3.1). 17 of 21 Therefore, the LSTM-type anomaly detector was tested on experimental signals of me- chanical contraction of left atrium cardiac tissue (AE) stimulated at successive frequencies of 1, 2, 4, 6, 8 and 10 Hz and also on experimental data from the unstimulated right atrium, and 10 Hz and also on experimental data from the unstimulated right atrium, since it is since it is controlled controlled by pacemaker by pacemaker cells. Forcells. For stimulated stimulated left atrialeft atria recordings, recordings, normal ex- normal experimental perimental recordings were used for network training, while for right atria recordings were used for network training, while for right atria recordings, normal experi- recordings, normal mental experimental signals at signals stimulated stimulated a frequencyat aof frequency 5 Hz were of used, 5 Hz were as theused, as the contraction observed observed contraction in the heart tissue of the left atrium is approximately 5 Hz. Specifically, in the heart tissue of the left atrium is approximately 5 Hz. Specifically, 16 records were 16 records were selected from 12 mice where representative selected from 12 mice where representative anomalies were observed. anomalies were observed. Figure Figure 1111 shows showsthe experimental the experimental signals signalsofof the normal the normal category category used forfor used the training the training ofof the fault detector used in the experimental signals. the fault detector used in the experimental signals. Experimental data Time (samples) Figure Figure11.11. Normal Normalexperimental experimentalsignals with signals different with stimulation different frequencies stimulation (1 (1 frequencies to to 1010Hz). Hz).These These signals were used for training the anomaly detector that was applied to experimental signals. signals were used for training the anomaly detector that was applied to experimental signals. Figure Figure1212shows showsthe theLSTM-type LSTM-typeanomaly anomalydetector detectorapplied appliedto to experimental experimental mechani- mechanical calsignals signalsfrom fromthethe cardiac cardiac tissue tissue ofof thethe left left andand right right atria atria of of CD-1 CD-1 mice. mice. It visually It visually con- confirms firms thatthat the detector the detector is able is able to satisfactorily to satisfactorily identify identify the anomalies the anomalies present present in signal. in the the signal. Figure Figure12A 12Ashows showsan ananomaly anomalythat that appeared appeared in aa leftleftatrium atriumwhen whenstimulated stimulatedatat5 Hz. 5 Hz.Figure 12B Figure shows 12B shows anomalies anomalies detected detectedin in a right atrium, a right atrium, which whichremained remainedunstimulated unstimulated as it ashas a natural it has physiological a natural physiological pacemaker. pacemaker.These atriaatria These that did thatnot didfollow a stimulation not follow pattern a stimulation were more prone to the genesis of anomalies. Specifically, two cardiac arrest phenomena appeared, which were correctly identified by the anomaly detector. Mathematics 2022, 10, x FOR PEER REVIEW 19 Mathematics 2022, 10, 2786 pattern were more prone to the genesis of anomalies. Specifically, two18cardiac of 21 arrest nomena appeared, which were correctly identified by the anomaly detector. A * B * Figure 12. Figure 12. Anomalies Anomalies detected detected in experimental in experimental signals. (A) signals. (A) Left Left atrium. (B)atrium. (B) RightThe Right atrium. atrium. The terisk (*) indicates a zoom into figure. The triangle indicates the anomaly asterisk (*) indicates a zoom into figure. The triangle indicates the anomaly location. location. 4. Discussion4. Discussion We have reported a method We have reportedfor the detection a method for of theanomalies detection in of cardiac anomalies tissue based on in cardiac tissue base deep learning, in which deep twoin learning, anomaly which twodetectors anomalyhavedetectors been evaluated. have been The first detector evaluated. The first det used a recurrentused neural network a recurrent (RNN) neural of long(RNN) network short-term of long memory short-term(LSTM) type, (LSTM) memory while type, w the second used an autoencoder. The proposed methodology is “normal-only” the second used an autoencoder. The proposed methodology is “normal-only” bec because training of thetraining algorithm does of the not require algorithm doeslabeling of thelabeling not require whole dataset, but only of the whole of normal dataset, but only of no data. This is data. an approach that uses the power of deep learning representation This is an approach that uses the power of deep learning representation without wit usual a labelingusualconstraint, a labelingjust normal just constraint, data.normal No previous data. Nostudies previous have been have studies found in found i been the literature literature where thewhere methodology of detecting the methodology anomalies of detecting is applied anomalies is to mechanical applied to mechanical contractions intractions cardiac intissue. Nevertheless, we can compare it with the accuracy cardiac tissue. Nevertheless, we can compare it with the accuracy achieve achieved by deep-learning-based anomaly detectors deep-learning-based in cardiacinelectrical anomaly detectors signals, signals, cardiac electrical which iswhich in theis in the r range between between 0.75 and 0.1 (Zhao et al., 2021) [42]. In our study, despite dealinga with a s 0.75 and 0.1 (Zhao et al., 2021) [42]. In our study, despite dealing with signal with a high with signal-to-noise ratio andratio a high signal-to-noise a lot of biological and a lot of variability, biological when we used variability, the we used when LSTM-based anomaly detector, we achieved a maximum accuracy of 0.99 LSTM-based anomaly detector, we achieved a maximum accuracy of 0.99 and an accu and an accuracy of 0.88 in the worst of 0.88case. in the worst case. A set of synthetic A setdata was constructed, of synthetic data waswhich allowed constructed, for accurately which allowed for describing accuratelythedescribin different anomalies observed in the experimental recordings of the mouse different anomalies observed in the experimental recordings of the mouse cardiac cardiac atria. These artificial data were used to validate the correct performance of the anomaly detection These artificial data were used to validate the correct performance of the anomaly d methodology that was subsequently applied to real experimental data. tion methodology that was subsequently applied to real experimental data. The results of the anomaly detector framework presented in this study have been The results of the anomaly detector framework presented in this study have exported to an html file with interactive graphics where the anomalous regions are high- exported to an html file with interactive graphics where the anomalous regions are lighted. In addition, this file incorporates basic visualization tools that allow the researcher lighted. In addition, this file incorporates basic visualization tools that allow th or clinician to perform actions such as zooming in, zooming out, scaling and downloading searcher or clinician to perform actions such as zooming in, zooming out, scaling the results as an image file (png file). This whole set of interactive tools allows to delve downloading the results as an image file (png file). This whole set of interactive too deeper into the analysis or explanation of anomalies which are associated with the gen- lows to delve deeper into the analysis or explanation of anomalies which are assoc eration of cardiac arrhythmias and to easily gain actionable insights. Files with an html with the generation of cardiac arrhythmias and to easily gain actionable insights. extension can be used in any operating system (OS), do not require computer resources and with an html extension can be used in any operating system (OS), do not require comp can be used in offline mode (without Internet connection). All html interactive graphics files are freely accessibleand resources on can a Gitberepository used in offline page mode (Github)(without under Internet the MIT connection). Open software All html int tive graphics files license (See supplementary dataare freely accessible on a Git repository page (Github) under the section). From a biological point of view,(See Open software license thissupplementary methodology may databe section). of interest in the fields of heart physiology, pharmacology or pathology. Although the aim of this work was primarily methodological, it is important to note that the high frequency of heartbeats occurring physiologically in a mouse can make quantification difficult. The detection of abnormalities Mathematics 2022, 10, 2786 19 of 21 is a first step in the detection of different cardiac pathologies associated with changes in contraction amplitude (inotropic), frequency (chronotropic) or rhythm (arrhythmias). This type of experimental approach, together with traditional methods of calculating the parameters associated with contractions (amplitude, duration, frequency, etc.), may be useful in future studies where these alterations are studied. Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/math10152786/s1, S1: Autoencoder-based anomaly de- tector framework; S2: LSTM-based anomaly detector framework. Additionally, interactive graph- ics files are freely accessible on a Git repository page under the MIT Open software license (https://github.com/xavierMarimon/CardiacTissueAnomaly/). Author Contributions: Conceptualization, X.M., R.B. and M.J.; data curation, X.M., S.T., M.J. and A.O.; formal analysis, X.M. and R.B.; funding acquisition, R.B. and M.J.; investigation, X.M., S.T., R.B., M.J. and A.O.; methodology, X.M., S.T., R.B. and M.J.; project administration, R.B. and M.J.; resources, R.B. and M.J.; software, X.M. and A.O.; supervision, R.B. and M.J.; validation, X.M., S.T., R.B., M.J. and A.O.; visualization, X.M. and A.O.; writing—original draft, X.M. and R.B.; writing—review and editing, X.M., S.T., R.B., M.J. and A.O. All authors have read and agreed to the published version of the manuscript. Funding: This research was funded by the Spanish Ministry of Science and Innovation, grant number PID2020-116927RB-C22. Institutional Review Board Statement: The animal study protocol was approved by the Institutional Ethics Committee of Universitat de Autónoma de Barcelona. Informed Consent Statement: Not applicable. Data Availability Statement: Codes are available from the authors upon request. Acknowledgments: We gratefully acknowledge the support of NVIDIA Corporation with the dona- tion of the Titan X Pascal GPU used for this research and Marta Pulido for editing the manuscript and editorial assistance. Conflicts of Interest: The authors declare no conflict of interest. References 1. Lindemann, B.; Fesenmayr, F.; Jazdi, N.; Weyrich, M. Anomaly Detection in Discrete Manufacturing Using Self-Learning Approaches. Procedia CIRP 2019, 79, 313–318. [CrossRef] 2. Gupta, M.; Gao, J.; Aggarwal, C.C.; Han, J. Outlier Detection for Temporal Data: A Survey. IEEE Trans. Knowl. Data Eng. 2014, 26, 2250–2267. [CrossRef] 3. Chandola, V.; Cheboli, D.; Kumar, V. Detecting Anomalies in a Time Series Database; Technical Report; University of Minnesota: Minneapolis, MN, USA, 2009. 4. Hawkins, D.M. Identification of Outliers. Chapman and Hall: London—New York 1980, 188 S., £ 14, 50. Biom. J. 1987, 29, 198. [CrossRef] 5. Keogh, E.; Lin, J.; Fu, A. HOT SAX: Efficiently Finding the Most Unusual Time Series Subsequence. In Proceedings of the Fifth IEEE International Conference on Data Mining (ICDM’05), Houston, TX, USA, 27–30 November 2005; pp. 226–233. 6. Lis, A.; Dworakowski, Z.; Czubak, P. An Anomaly Detection Method for Rotating Machinery Monitoring Based on the Most Representative Data. J. Vibroeng. 2021, 23, 861–876. [CrossRef] 7. Yang, H.; Kan, C.; Liu, G.; Chen, Y. Spatiotemporal Differentiation of Myocardial Infarctions. IEEE Trans. Autom. Sci. Eng. 2013, 10, 938–947. [CrossRef] 8. Yang, H. Multiscale Recurrence Quantification Analysis of Spatial Cardiac Vectorcardiogram Signals. IEEE Trans. Biomed. Eng. 2011, 58, 339–347. [CrossRef] 9. Li, H.; Boulanger, P. A Survey of Heart Anomaly Detection Using Ambulatory Electrocardiogram (ECG). Sensors 2020, 20, 1461. [CrossRef] 10. Berndt, D.J.; Clifford, J. Using Dynamic Time Warping to Find Patterns in Time Series. In Proceedings of the KDD Workshop, New York, NY, USA, 26 April 1994. 11. Sivaraks, H.; Ratanamahatana, C.A. Robust and Accurate Anomaly Detection in ECG Artifacts Using Time Series Motif Discovery. Comput. Math. Methods Med. 2015, 2015, 1–20. [CrossRef] 12. Chuah, M.C.; Fu, F. ECG Anomaly Detection via Time Series Analysis. In Frontiers of High Performance Computing and Networking ISPA Workshops; Parimala, T., Ed.; Springer: Berlin, Germany, 2007; pp. 123–135. [CrossRef] Mathematics 2022, 10, 2786 20 of 21 13. Pham, N.D.; Le, Q.L.; Dang, T.K. HOT ASAX: A Novel Adaptive Symbolic Representation for Time Series Discords Discovery. In Lecture Notes in Computer Science; (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics); Springer: Berlin, Germany, 2010; Volume 5990, pp. 113–121. ISBN 3642121446. 14. Veeravalli, B.; Deepu, C.J.; Ngo, D. Real-Time, Personalized Anomaly Detection in Streaming Data for Wearable Healthcare Devices. In Handbook of Large-Scale Distributed Computing in Smart Healthcare; Khan, S.U., Zomaya, A.Y., Abbas, A., Eds.; Springer: Berlin, Germany, 2017; pp. 403–426. [CrossRef] 15. Christov, I.; Gómez-Herrero, G.; Krasteva, V.; Jekova, I.; Gotchev, A.; Egiazarian, K. Comparative Study of Morphological and Time-Frequency ECG Descriptors for Heartbeat Classification. Med. Eng. Phys. 2006, 28, 876–887. [CrossRef] 16. de Chazal, P.; Reilly, R.B. A Patient-Adapting Heartbeat Classifier Using ECG Morphology and Heartbeat Interval Features. IEEE Trans. Biomed. Eng. 2006, 53, 2535–2543. [CrossRef] 17. Llamedo, M.; Martínez, J.P. Heartbeat Classification Using Feature Selection Driven by Database Generalization Criteria. IEEE Trans. Biomed. Eng. 2011, 58, 616–625. [CrossRef] 18. Gretton, A.; Smola, A.; Huang, J.; Schmittfull, M.; Borgwardt, K.; Schölkopf, B. Covariate Shift by Kernel Mean Matching. In Dataset Shift in Machine Learning; The MIT Press: Cambridge, MA, USA, 2008; pp. 131–160. 19. Li, K.; Du, N.; Zhang, A. Detecting ECG Abnormalities via Transductive Transfer Learning. In Proceedings of the ACM Interna- tional Conference on Bioinformatics, Computational Biology and Biomedicine-BCB’12, Orlando, FL, USA, 8–10 October 2012; pp. 210–217. [CrossRef] 20. Ye, C.; Kumar, B.V.K.V.; Coimbra, M.T. Heartbeat Classification Using Morphological and Dynamic Features of ECG Signals. IEEE Trans. Biomed. Eng. 2012, 59, 2930–2941. [CrossRef] 21. Zhang, Z.; Dong, J.; Luo, X.; Choi, K.-S.; Wu, X. Heartbeat Classification Using Disease-Specific Feature Selection. Comput. Biol. Med. 2014, 46, 79–89. [CrossRef] [PubMed] 22. Rajesh, K.N.V.P.S.; Dhuli, R. Classification of ECG Heartbeats Using Nonlinear Decomposition Methods and Support Vector Machine. Comput. Biol. Med. 2017, 87, 271–284. [CrossRef] [PubMed] 23. Thomas, M.; Das, M.K.; Ari, S. Automatic ECG Arrhythmia Classification Using Dual Tree Complex Wavelet Based Features. AEU. Int. J. Electron. Commun. 2015, 69, 715–721. [CrossRef] 24. Wess, M.; Manoj, P.D.S.; Jantsch, A. Neural Network Based ECG Anomaly Detection on FPGA and Trade-off Analysis. In Proceedings of the 2017 IEEE International Symposium on Circuits and Systems (ISCAS), Baltimore, MD, USA, 28–31 May 2017; pp. 1–4. 25. Kiranyaz, S.; Ince, T.; Gabbouj, M. Real-Time Patient-Specific ECG Classification by 1-D Convolutional Neural Networks. IEEE Trans. Biomed. Eng. 2016, 63, 664–675. [CrossRef] [PubMed] 26. Kiranyaz, S.; Ince, T.; Gabbouj, M. Personalized Monitoring and Advance Warning System for Cardiac Arrhythmias. Sci. Rep. 2017, 7, 9270. [CrossRef] 27. Rajpurkar, P.; Hannun, A.Y.; Haghpanahi, M.; Tison, G.H.; Bourn, C.; Turakhia, M.P.; Ng, A.Y. Cardiologist-Level Arrhythmia Detection and Classification in Ambulatory Electrocardiograms Using a Deep Neural Network. Nat. Med. 2019, 25, 65–69. [CrossRef] 28. Acharya, U.R.; Fujita, H.; Lih, O.S.; Hagiwara, Y.; Tan, J.H.; Adam, M. Automated Detection of Arrhythmias Using Different Intervals of Tachycardia ECG Segments with Convolutional Neural Network. Inf. Sci. 2017, 405, 81–90. [CrossRef] 29. Zhou, H.; Kan, C. Tensor-Based ECG Anomaly Detection toward Cardiac Monitoring in the Internet of Health Things. Sensors 2021, 21, 4173. [CrossRef] 30. Li, D.; Zhang, J.; Zhang, Q.; Wei, X. Classification of ECG Signals Based on 1D Convolution Neural Network. In Proceedings of the 2017 IEEE 19th International Conference on e-Health Networking, Applications and Services (Healthcom), Dalian, China, 12–15 October 2017; pp. 1–6. 31. Bailly, R.; Malfante, M.; Allier, C.; Ghenim, L.; Mars, J. Deep Anomaly Detection Using Self-Supervised Learning: Application to Time Series of Cellular Data; IEEE Access: Porto, Portugal, 2021. 32. Malhotra, P.; Vig, L.; Shroff, G.M.; Agarwal, P. Long Short Term Memory Networks for Anomaly Detection in Time Series; ESANN: Bruges, Belgium, 2015. 33. McLachlan, G.; Peel, D. Finite Mixture Models; John Wiley & Sons, Inc.: Hoboken, NJ, USA, 2000; ISBN 9780471721185. 34. Hochreiter, S.; Schmidhuber, J. Long Short-Term Memory. Neural Comput. 1997, 9, 1735–1780. [CrossRef] [PubMed] 35. Wang, J.; Zhang, J.; Wang, X. Bilateral LSTM: A Two-Dimensional Long Short-Term Memory Model With Multiply Memory Units for Short-Term Cycle Time Forecasting in Re-Entrant Manufacturing Systems. IEEE Trans. Ind. Inform. 2018, 14, 748–758. [CrossRef] 36. Kingma, D.P.; Ba, J. Adam: A Method for Stochastic Optimization. In Proceedings of the 3rd International Conference for Learning Representations, San Diego, CA, USA, 22 December 2015. 37. Bishop, C.M. Pattern Recognition and Machine Learning; Springer: New York, NY, USA, 2006. 38. Murphy, K.P. Machine Learning: A Probabilistic Perspective; The MIT Press: Cambridge, MA, USA, 2012. 39. Olshausen, B.A.; Field, D.J. Sparse Coding with an Overcomplete Basis Set: A Strategy Employed by V1? Vis. Res. 1997, 37, 3311–3325. [CrossRef] 40. Nair, V.; Hinton, G.E. 3D Object Recognition with Deep Belief Nets. In Advances in Neural Information Processing Systems; Bengio, Y., Schuurmans, D., Lafferty, J., Williams, C., Culotta, A., Eds.; Curran Associates, Inc.: Red Hook, NY, USA, 2009; Volume 22. Mathematics 2022, 10, 2786 21 of 21 41. Coates, A.; Ng, A.Y. The Importance of Encoding versus Training with Sparse Coding and Vector Quantization. In Proceedings of the 28th International Conference on International Conference on Machine Learning, Washington, DC, USA, 28 June–2 July 2011; Omnipress: Madison, WI, USA, 2011; pp. 921–928. 42. Zhao, X.; Yan, H.; Hu, Z.; Du, D. Deep Spatio-Temporal Sparse Decomposition for Trend Prediction and Anomaly Detection in Cardiac Electrical Conduction. IISE Trans. Healthc. Syst. Eng. 2022, 12, 150–164. [CrossRef]