MODEL-BASED ECG ANALYSIS:TOWARDS PATIENT-SPECIFIC WEARABLE ECG MONITORING: MODEL-BASED ECG ANALYSIS:TOWARDS PATIENT-SPECIFIC WEARABLE ECG MONITORING

(1)

IT 20 020

Examensarbete 30 hp

Mars 2020

MODEL-BASED ECG ANALYSIS:

TOWARDS PATIENT-SPECIFIC

WEARABLE ECG MONITORING

Adnan Albaba

(2)

(3)

Teknisk- naturvetenskaplig fakultet UTH-enheten Besöksadress: Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0 Postadress: Box 536 751 21 Uppsala Telefon: 018 – 471 30 03 Telefax: 018 – 471 30 00 Hemsida: http://www.teknat.uu.se/student

Abstract

MODEL-BASED ECG ANALYSIS:

TOWARDS PATIENT-SPECIFIC

WEARABLE ECG MONITORING

Adnan Albaba

In this thesis, model-based analysis approach is considered as a possible solution towards a patient-specific point-of-care device for the purpose of electrocardiogram monitoring. Two novel methods are proposed, tested, and quantitatively evaluated. First, a method for estimating the instantaneous heart rate using the morphological features of one electrocardiogram beat at a time is proposed. This work is not aimed at introducing an alternative way for heart rate estimation, but rather illustrates the utility of model-based electrocardiogram analysis in online individualized monitoring of the heart function. The heart rate estimation problem is reduced to fitting one parameter, whose value is related to the nine parameters of a realistic nonlinear model of the electrocardiogram and estimated from data by nonlinear least-squares optimization. The method feasibility is evaluated on synthetic electrocardiogram signals as well as signals acquired from MIT-BIH databases at Physionet website. Moreover, the performance of the method was tested under realistic free-moving conditions using a wearable electrocardiogram and heart monitor with encouraging results. Second, a model-based method for patient-specific detection of deformed electrocardiogram beats is proposed. Five parameters of a patient-specific nonlinear electrocardiogram model are estimated from data by nonlinear least-squares optimization. The normal variability of the model parameters is captured by estimated probability density functions. A binary classifier, based on stochastic anomaly detection methods, along with a pre-tuned classification threshold, is employed for detecting the abnormal electrocardiogram beats. The utility of the proposed approach is tested by validating it on annotated arrhythmia data recorded under clinical conditions.

Handledare: Alexander Medvedev Ämnesgranskare: Hans Rosth Examinator: Phillipp Rummer IT 20 020

(4)

(5)

I

Acknowledgements

First and foremost, I must thank God for all the blessings in my life.

I am extremely grateful to my dear parents for always believing in my dreams and constantly supporting my ideas in every way possible.

I thank Prof. Alexander Medvedev, for accepting to supervise my thesis, for his guiding inputs and vital ideas, and, most of all, for teaching me how to do academic research. I thank Dr. Hans Rosth for his constructive feedback while reviewing my thesis. I thank the academic staff at the department of Information Technology, the department of Physics and Astronomy, and the department of Engineering Sciences at Uppsala University, for facilitating my journey during the last two years. Special thanks to Liselott Dominicus, simply for being the most helpful person at our department. Also, thanks to my colleagues at the Embedded Systems master program, and to my SensUs 2018 teammates and mentors.

In Sweden; Rami, without you, this dream would have never become true. Thank you for everything. Ammar, thanks for cracking the enigma of the Swedish bureaucracy! Nour, thanks for being a good subject for my experiments, yet a terrible TRIX partner. Ahmad, thanks for being a good TRIX partner. Deyaa and Sofiene, thanks for being such nice guys. Isabel, thanks for your friendliness and smiles.

In Kuwait; Hamza, your support will never be forgotten. Yassino, Sammano, and Almonther, thank you brothers for listening to my endless stories. Hafez, thanks for the IELTS. Naif, Bitar, Shurbaji, Morsy, Tahina, Okasha, Alhalaki, Saif, Arafat and Katot, thanks for allways being there for me. Bilal, thanks for being such a FIFA loser bro. In Syria; My family, thanks for always believing in me and praying for me.

In Jordan; My Professors and colleagues at JUST, my flatmates, and my dear friends, thanks for the great 5 years. Majd, thanks for keeping in touch. Shimaa, thanks for the unforgettable memories.

In Palestine; Farrah, thanks for the long night technology talks and jokes. In UAE; Bisher, thanks for keeping me an overseas company.

In Germany; Alhalabi’s, thank you guys for the Physiology talks.

In USA; Thaer, a conversation with you is like traveling to my near future, thank you. In Canada; Mouayad, thanks for the political and spiritual talks.

For those who I forgot to mention, last-minute submission guys! Hopefully, you will be acknowledged in my PhD dissertation.

(6)

(7)

Table of Contents II

4 Online Model-Based Beat-by-beat Heart Rate Estimation 16 4.1 Theory . . . 16 4.2 Methods . . . 16 4.2.1 Preprocessing . . . 17 4.2.2 Initial Settings . . . 18 4.2.3 Optimization . . . 18 4.2.4 Ground Truth . . . 19 4.3 Results . . . 19 4.3.1 ECGSYN Signals . . . 21 4.3.2 MIT-BIH Database . . . 21 4.3.3 MAX-ECG-MONITOR . . . 23

(8)

Table of Contents III 4.3.4 Quantitative Evaluation . . . 23 4.4 Discussion . . . 24 4.4.1 ECGSYN Signals . . . 24 4.4.2 MIT-BIH Database . . . 25 4.4.3 MAX-ECG-MONITOR . . . 27 4.4.4 Limitations . . . 28

5 Patient-Specific ECG Monitoring by Model-Based Stochastic Anomaly De-tection 30 5.1 Theory . . . 30 5.2 Methods . . . 31 5.2.1 Preprocessing . . . 31 5.2.2 Period Normalization . . . 32 5.2.3 Optimization . . . 33 5.2.4 Initial Settings . . . 34

5.2.5 Stochastic Anomaly Detection . . . 34

5.2.6 Decision Making . . . 35

5.2.7 Feature Selection and Classification-Threshold Tuning . . . 37

5.2.8 Ground Truth . . . 38 5.3 Results . . . 38 5.3.1 Quantitative Evaluation . . . 39 5.4 Discussion . . . 41 5.4.1 Methodology . . . 41 5.4.2 Results . . . 42 5.4.3 Limitations . . . 42 6 Conclusion 43 Literature 44

(9)

List of Tables IV

List of Tables

Table 4.1: Initial values for the parameters of the model . . . 18 Table 4.2: RMSE for ECGSYN signals with initial Ehr=15 . n(t)[mV] and Fs[Hz]

were varied. RMSE for heart rates within the range 15–120, the range 120–160 and total RMSE. . . 24 Table 4.3: Root Mean Squared Error for every MIT-BIH Database individual record

with initial E_hr= 15 (RMSE15) and 60 (RMSE60). . . 24 Table 5.1: Qualitative Evaluation of the Results for the Decisions in the

(10)

List of Figures V

List of Figures

Figure 2.1: Normal ECG beat. [25] . . . 6 Figure 2.2: (a) and (b) illustrate the incidence of a PVC and PAC episodes

re-spectively, where red (X) mark the abnormal beats and red filled (O) mark the normal beats; (c) and (d) show a segment of AFL and AFIB episodes, respectively; (e) and (f) illustrate the incidence of a VT and RBBB episodes respectively, where red (X) mark the abnormal beats and the red filled (O) marks the normal beats. . . 8 Figure 2.3: NUUBO’s wearable ECG platform. [40] . . . 10 Figure 3.1: Normal ECG beat generated using ECGSYN tool . . . 12 Figure 3.2: (a) depicts the waveform of x, y, and z. (b) depicts the phase plot of x

and y. . . 13 Figure 3.3: Trajectory generated by the ECG dynamical model in three-dimensional

state-space . . . 14 Figure 4.1: The effect of the sampling frequency Fson the performance of the HR

estimator. (a) Results for HR estimation using synthetic ECG signals sampled with Fs=512Hz, (b) Fs=360Hz and (c) Fs=60Hz. . . 20 Figure 4.2: The effect of adding uniformly distributed measurement noise n(t)to

the ECG beat data on the performance of the HR estimator, Fs=360Hz. (a) Results for HR estimation using synthesised ECG signals with n(t)

= 0 mV, (b) n(t)= 0.1 mV and (c) n(t)= 0.2 mV. . . 20 Figure 4.3: Results of the model-based HR estimation of the beats acquired from

MIT-BIH Normal Sinus Rhythm Database (a) and MIT-BIH Arrhyth-mia Database (b). Points marked with (o) represent the true HR values, while (x) and (+) represent the HR estimates with initial value of Ehr equal to 15 and 60, respectively. (c) Results of the model-based HR es-timation of the beats acquired from the MAX-ECG-MONITOR. Points marked with (o) represent the true HR values, while (x) represent the HR estimates with initial value of Ehrequal to 60. The tests were per-formed in three different physical statuses which are Sitting, Standing, and Running as illustrated by the vertical lines. (d) The HR estimates for a segment of the ECG test signal 106 acquired from the MIT-BIH Arrhythmia Database. The correlation coefficient between the HR estimates and true values is 0.6337. . . 22

(11)

List of Figures VI

Figure 4.4: The MAX-ECG-MONITOR deployed on the chest of the subject for wearable ECG acquisition. . . 23 Figure 4.5: The effects of local minima on the performance of the HR estimator.

(a) the estimated heart rates for 51 beats with true HR between 96 and 101 covered with the step 0.1, the initial value of Ehrwas set to 15, Fs= 360 Hz, and n(t)= 0 mV. Points marked with (o) represent the true HR values, while (x) represent the HR estimates. (b) shows the resultant squared norm of the residual for optimization procedures in (a). . . . 26 Figure 4.6: Tuning the initial value of Ehrimproves the performance of the HR

estimator by overcoming local minima issue illustrated in Fig. 4.5. (a) the estimated heart rates for 51 beats with true HR between 96 and 101 covered with the step 0.1, Fs = 360 Hz and n(t)= 0 mV. Points marked with (o) represent the true HR values, while (x) represent the HR estimates. (b) the resulting squared norm of the residual for optimization procedures in (a). . . 26 Figure 4.7: The difference in morphology between one ECG beat taken from the

MAX-ECG-MONITOR sitting test signal with true HR equals to 86 BPM (in Blue), and a synthetic ECG beat with true HR equals to 86 BPM generated by the ECGSYN model (in orange). It also shows the same ECG beat generated by ECGSYN model after being modified by tuning its parameters to fit with the MAX-ECG-MONITOR beat (dashed line in yellow). . . 27 Figure 4.8: (a) A one-minute segment of the ECG record 106 acquired from the

MIT-BIH Arrhythmia Database, where the amplitude of the signal changes during the process of the ECG acquisition. (b) HR estimates vs. true values plot for the test sample shown in (a). A significant change in the performance of the HR estimator after the point where the amplitude of the ECG signal changes which denoted by the black arrow is shown. (c) Amplitude scaling of the ECG segment normalizing the amplitude scale of the ECG segment shown in Fig. 4.8(a) cancels out the effect of the changing amplitude shown in (b). . . 29 Figure 5.1: A flow chart illustrating the main steps of the proposed period

nor-malization technique. . . 32 Figure 5.2: The results for the period normalization method; (a) depicts the

stretch-ing of a normal ECG beat with HR = 130 BPM, the resultant ECG beat is equivalent to a normal ECG beat with HR = 60 BPM. (b) depicts the shrinking of a normal ECG beat with HR = 60 BPM, the resultant ECG beat is equivalent to a normal ECG beat with HR = 40 BPM. Signals in orange represent the original ECG beats, while those in blue represent the scaled ones. . . 33

(12)

List of Figures VII

Figure 5.3: The ROC curves and their corresponding AUC for three different combinations of features of the presented classifier. . . 36 Figure 5.4: The distributions of the features aP, aR, aT, corresponding to P (X-axis),

R (Y-axis) and T-waves (Z-axis), for different patients in both normal and abnormal cases; (a), (b) and (c) depict the distribution of the es-timated parameters in normal case, AFIB and AFL respectively; (d) and (e) depict the distribution of the estimated parameters in case of normal and PAC beats, respectively; (f) and (g) depict the distribu-tion of the estimated parameters in case of normal and RBBB beats, respectively; (h) and (i) depict the distribution of the estimated values parameters in case of normal and PVC beats, respectively. Estimates for normal beats (training) are marked with blue (O), and estimations of normal and abnormal testing beats are marked with orange (X). . . 40 Figure 5.5: (a) and (b) depict The confusion matrices for the PR2T-based

classi-fier in both cases of tuning the classification threshold and the final evaluation of the presented algorithm, respectively. . . 41

(13)

List of Acronyms VIII

List of Acronyms

CVD Cardiovascular disease

ECG Electrocardiograph

EKG Electrocardiograph

IoT Internet of Things

HR Heart Rate

SA Sinoatrial

AV Atrioventricular

BPM Beat per Minute

BP Blood Pressure

PPG Photoplethysmography

SVD Singular Value Decomposition

PDF Probability Distribution Function

PVC Premature Ventricular Contraction

PAC Premature Atrial Beat

AFL Atrial Flutter

AFIB Atrial Fibrillation

VT Ventricular Tachycardia

RBBB Right Bundle Branch Block

RMSE Root Mean Square Error

RESNORM Resultant Squared Norm of the Residual

SNR Signal-to-Noise Ratio

OSA Orthogonal Series Approximation

AUC Area Under Curve

ROC Receiver Operating Characteristics

TP True Positive

FP False Positive

TN True Negative

FN False Negative RTP True Positive Rate RFP False Positive Rate

(14)

Introduction 1

1 Introduction

1.1 Overview

No one can deny the fact that cardiovascular disease (CVD) has recently become one of the biggest threats for mankind. On top of being the major cause of death globally, patients suffering from CVD are accounting for the highest percentage of hospitalized patients [1]. Locally speaking, CVD remains the main cause of death in Sweden [2]. CVD is the name for the group of disorders of heart and blood vessels, and include:

• Hypertension (high blood pressure (BP))

• Coronary heart disease (heart attack)

• Cerebrovascular disease (stroke)

• Peripheral vascular disease

• Heart failure

• Rheumatic heart disease

• Congenital heart disease

• Cardiomyopathies.

Electrocardiography (ECG) monitoring is one of the well-known techniques used by the cardiologists for collecting information about the structure and functionality of the cardiovascular system. Early recognition of cardiac-related issues such as angina, dyspnea, and syncope can be critical, and CVD patients need to visit their doctors for medical tests on a regular basis. ECG is frequently used in hospitals and clinics for its usefulness in diagnosing heart diseases. Monitoring of cardiac electrical activity for long time intervals can also be useful for detection of arrhythmias. In fact, “The diagnostic yield is increased by 15% to 39% by a 24-hour recording” [3].

Conventionally, interpreting the ECG signal is done by a specialist with the help of an ECG recording machine. This task requires intellect, training, and an organized approach. During the past 50 years, extensive research has been devoted to developing algorithms for computerized and automated ECG processing and analysis. As a result, computerized ECG monitoring is now a well-established approach, and around 100 million computerized ECG interpretations are being recorded every year in the United States alone [4].

(15)

Introduction 2

1.2 Problem Formulation and My Contributions

Computerized monitoring of ECG is now a well established method, thanks to the significant progress derived by many proposed algorithms over the years, see Section 1.3. Different ECG beat detection and classification techniques were presented in the literature to support wearable ECG monitors see Section 1.3. Such algorithms can either be used for processing the ECG signal in real-time or off-line. While the design criteria may differ between on-line and off-line ECG monitors, one common problem faced in their development is the inter-patient and intra-patient variation in the ECG morphology. The QT-interval has a strong inverse relation with the HR. Comparing QT-intervals between different controls with different HR would require the use of a QT-interval rate correction formula [5].

Researchers have paid less attention to the inter-patient and intra-patient variations in the QT-RR relationship. Such variation also limits the reliability as well as the robustness of any universal clinical decision support algorithm. It also contributes to the errors of any fixed HR correction formula [6].

Therefore, a beat detector and classifier producing good results for a certain training database can perform inconsistently, when dealing with a different patient cohort thus ef-fectively preventing reliable solutions for automated ECG processing. Individualization of the medical solution may be a possible solution to address this problem.

In this thesis, Model-Bsaed ECG analysis approach is considered for the previously mentioned problem. Two novel methods for estimating the HR and detecting the abnormalities in the ECG pattern are introduced, tested, and discussed.

The outcome of this work is summarized below:

• Online Model-Based Beat-by-beat Heart Rate Estimation: A novel method for estimating the instantaneous HR using the morphological features of one ECG beat at a time. This tool is not aimed at introducing an alternative way for estimating the HR, but rather illustrates the utility of model-based ECG analysis in online individualized monitoring of the heart function, and demonstrates the possibility of estimating the HR from the data of one single beat.

• Patient-Specific ECG Monitoring by Model-Based Stochastic Anomaly Detection: A novel model-based method for patient-specific detection of deformed ECG beats. Five parameters of a patient-specific nonlinear ECG model are estimated from data by nonlinear least-squares optimization. The normal variability of the model parameters is captured by estimated probability density functions. A binary classifier, based on stochastic anomaly detection methods, along with a pre-tuned classification threshold, is employed for detecting the abnormal ECG beats. Both tools are addressed separately in the following chapters, see Section 1.4.

(16)

Introduction 3

1.3 Related Work

Various methods and techniques were proposed for automatic ECG heartbeat classi-fication in literature. Quantitative ECG analysis was used by D. E. Krummen, et al. to demonstrate that new computerized algorithms can detect heart arrhythmia with relatively high accuracy [7]. Discrete Wavelet Transform was employed by L. Senhadji, et al. to act as a feature extraction approach and combined it with a linear discriminant classifier [8], while Shyu, et al. classified PVC beats using a combination of wavelet fea-ture extraction and Fuzzy Neural Network [9]. M. Lagerholm, et al. clustered ECG beats using Hermite functions along with Self-Organizing Maps [10]. Hidden Markov models were used by R.V. Andreao, et al. for ECG signal analysis [11]. M.I. Owis, et al. used the correlation dimension and largest Lyapunov exponent to detect ECG arrhythmia [12]. A number of additional ECG beat detection and classification techniques were presented in the literature to support wearable ECG monitors [13]–[16].

Whilst non of the previously mentioned works provides a truly patient-specific approach, Y.H. Hu, et al. brought up the topic of patient-adaptable ECG beat classifiers, by utilizing a mixture of experts approach [17].

As for the work been done on ECG mathematical models, A. Ruha, et al. preduced a synthetic ECG waveform generator to test their Real-Time Microprocessor QRS Detector system [18]. However, this generator was not meant to be used for generating highly realistic ECG waveform, Moreover, the synthetic ECG does not include P-wave or beat-to-beat variations.

On the other hand, ECG nonlinear dynamical model [19], introduced by P.E. McSharry, et al. and better known as ECGSYN, is based on time-varying differential equations and includes beat-to-beat variations in morphology as well as inter-beat timing variations, which makes it capable of generating highly realistic synthetic ECG waveform. G.D. Clifford, et al. used the ECGSYN model for filtering, compression and classification of the ECG [20].

It is worth mentioning that the ECGSYN model was utilized in this work, see Section 3.2.1 for more details.

The relation between the morphology of the ECG beat and the corresponding HR is taken into consideration especially in articles related to the topic of "ECG as a Biometric". This is because normalization is needed when developing an ECG based identification system [21], [22].

Stochastic anomaly detection was proposed in [23] and successfully applied to the analysis of eye-tracking data in [24].

To my knowledge, no previous work was done on the concept of model-based HR estimation using the information of one ECG beat, which supports the novelty of the proposed method in this thesis. Moreover, model-based stochastic anomaly detection for a patient-specific ECG monitor is also a new and novel approach to detect and identify abnormal ECG beats.

(17)

Introduction 4

1.4 Structure of The Thesis

The rest of this thesis is organized as follows; Chapter 2 presents a theoretical background about the ECG from physiological and pathological points of view. Chapter 3 gives a brief overview about the ECG modeling approaches and sets up the mathematical model of the ECG waveform, which is used in this work. Chapter 4 summarizes the theory and methods for the Online Model-Based Beat-by-beat Heart Rate Estimation method. It also presents the results of algorithm performance evaluation and discusses them. Chapter 5 summarizes the theory and methods for the Patient-Specific ECG Monitoring by Model-Based Stochastic Anomaly Detection method. It also presents the results of algorithm performance evaluation and discusses them. Conclusions are drawn in Chapter 6.

(18)

Background 5

2 Background

2.1 Electrocardiography (ECG)

An electrocardiogram - abbreviated as EKG or ECG - is a test that measures the electrical activity of the heartbeat, see Fig. 2.1. With each beat, an electrical impulse travels through the heart. This impulse causes the muscle of the heart to squeeze and pump blood from the heart.

The sinoatrial (SA) node, located in the right atrium, is the natural pacemaker of the heart and initiates the heart beat in normal cases. Electrical impulses originating from the SA node spread throughout both atria and stimulate them to depolarize. The atrial depolarization is represented by the P-wave. The PQ-segment represents the time it takes for the electrical impulses to travel from the SA node to the atrioventricular (AV) node located near the AV valve. The AV node serves as an electrical gateway to the ventricles and delays the electrical impulses, which cause the atria relax or repolarize. The atrial repolarization is masked by the QRS-complex. After that, the AV node passes the electrical impulses to the bundle of His that is then divided into right and left bundle branches. Those branches conduct the electrical impulses towards the apex of the heart. The electrical impulses are passed further onto the Purkinje fibers and spread throughout the ventricular myocardium, which causes the ventricles to depolarize. The ventricular depolarization is represented by the QRS-complex. The next step is the contraction of the ventricles, which causes a plateau in the myocardium action potential and is represented by the ST-segment. The ventricular repolarization then occurs and is represented by the T-wave.

Therefore, a normal ECG signal can be broken down into the following fiducial points:

• P-wave is the depolarization of the atria

• Q-wave is associated with the deflection immediately before ventricular depolar-ization

• R-wave is associated with the peak of the ventricular depolarization

• S-wave is associated with the deflection proceeding the ventricular depolarization

• T-wave is the repolarization of the ventricles

An ECG gives two major kinds of information: First, by measuring time intervals on the ECG, a doctor can determine how long the electrical wave takes to pass through the

(19)

Background 6

Figure 2.1: Normal ECG beat. [25]

heart. Finding out how long a wave takes to travel from one part of the heart to the next shows if the electrical activity is normal or slow, fast or irregular. Second, by measuring the amount of electrical activity passing through the heart muscle, a cardiologist may be able to find out if parts of the heart are too large or are overworked [26].

2.2 Heart Rate (HR)

The HR is the number of times in which the heart contracts or beats per minute [27]. Along with body temperature, blood pressure, and breathing rate, HR is one of the primary vital signs which indicate the status of the body0s life-sustaining functions. On top of clinical applications such as sleep test and stress test, HR is used in fitness and activity tracking solutions, as it varies according to the body0s physical needs for oxygen. HR can be measured manually by feeling specific points on the body, where the artery’s pulsation is transmitted. It can also be measured electrically by means of ECG and estimated from the electrocardiogram pattern. The most prominent feature in one ECG cycle is the R-wave, which is used as a marker for every beat and the time interval between two R-waves (RR-interval) can be used for estimating the HR.

From a physiological point of view, two main types of action potential are the pacemaker action potential which is generated spontaneously, i.e. the sinoatrial (SA) node located in the posterior wall of the right atrium, and the non-pacemaker action potential which is generated by depolarization current from adjacent cell (the Ventricle muscles). The

(20)

Background 7

HR is normally determined by the pacemaker action potential of the SA node at a rate which is determined by the spontaneous changes in conductances of potassium, sodium and calcium. If not affected by neurohumoral factors, this intrinsic automaticity exhibits a spontaneous firing rate of 100-115 beats per minute (BPM) and decreases with age [28]. Activating the vagus nerve decreases the HR and results in a normal resting HR between 60 and 80 BPM because of the significant vagal tone on the SA node at rest. On the contrary, an activation of sympathetic nerves leads to a HR increase above the intrinsic rate and explains the HR increment during physical exercise because of the reciprocal change in sympathetic and parasympathetic activity [28].

There are three main methods that are currently used for obtaining HR information: 1. ECG-based techniques capture the electrical impulses of the heart [29]. The

de-tection of QRS-complex is the most direct method for ECG-based heat rate mea-surement. Algorithms based on digital filters [30], derivative-based algorithms [31], wavelet-based algorithms [32], neural network-based approaches [33] and model-based methods [34] are used for QRS-complex detection.

2. Photoplethysmography (PPG) relies on the blood perfusion and the change of its volume. Optical HR monitoring is an example of PPG-based technologies implemented in wearable HR monitors [35].

3. Techniques based on arterial blood pressure that utilize the pulsatile waveform of BP [36].

Since HR measurements are subject to noise and artifact sources, several techniques have been developed for improving the process of HR estimation, e.g. by combining the signal quality indices with a Kalman filter [37].

2.3 Pathophysiology and ECG Patterns

In this Section, six types of arrhythmia (i.e. disorder in heart rhythm), which are investigated in this thesis, are briefly introduced:

Premature Ventricular Contraction (PVC)

A PVC occurs when the heart beat is initiated by Purkinje fibers instead of the SA node. ECG monitoring usually makes it possible to distinguish a PVC from a normal heart beat. According to [27], the specific effects caused by PVCs in the ECG (see Fig. 2.2(a)) are:

• Prolongation of the QRS-complex because the impulse is not conducted through the Purkinje system, but through the slow conducting muscles of the ventricles.

• High voltage level for the QRS-complex because the impulse travels in one direc-tion unlike in normal case where the depolarizadirec-tion waves of the two sides of the heart partially neutralize each other in the ECG.

(21)

Background 8

(a) (b)

(c) (d)

(e) (f)

Figure 2.2: (a) and (b) illustrate the incidence of a PVC and PAC episodes respectively, where red (X) mark the abnormal beats and red filled (O) mark the normal beats; (c) and (d) show a segment of AFL and AFIB episodes, respectively; (e) and (f) illustrate the incidence of a VT and RBBB episodes respectively, where red (X) mark the abnormal beats and the red filled (O) marks the normal beats.

(22)

Background 9

• The T-wave has an electrical potential polarity opposite to that of the QRS-complex because of the slow conduction of the impulse causing the muscle fibers to repolar-ize first.

Premature Atrial Beat (PAC)

PAC - and Aberrated Atrial Premature Contraction (APAC) - takes place when a region other than the atria depolarizes before the SA node and triggers a premature beat. The P-wave of the PAC occurs earlier in the heart cycle. Moreover, the PR-interval is shortened thus indicating that the ectopic origin of the beat is in the atria near the AV node. The interval between the PAC and the next contraction is slightly prolonged, which phenomenon is called a compensatory pause [38]. Fig. 2.2(b) illustrates a PAC episode.

Atrial Flutter (AFL)

A circus movement in the atria is usually the cause of the AFL condition, in which the electrical impulse travels as a single large wave in one direction around the atrial muscle. This leads to a rapid contraction rate of the atria at 200–350 BPM. In a typical ECG featuring AFL (see Fig. 2.2(c)), the P-waves are strong but a QRST-complex follows a P-wave only once for every two-three cycles of the atria [38].

Atrial Fibrillation (AFIB)

AFIB occurs when the atria beat rapidly and irregularly leading in an abnormal heart rhythm. An ECG during AFIB condition (see Fig. 2.2(d)) shows either no P-waves or only very fine wavy ones. On the other hand, the QRST-complexes are normal but irregular in timing. This irregularity in ventricular rhythm is due to the rapid and irregular arrival of impulses from the atrial muscle at the AV node [38].

Ventricular Tachycardia (VT)

VT is an arrhythmia which is caused by abnormal impulses in the ventricles. In VT cases, the heart rate is usually 100 BPM and out of sync with the atria. VT appears as a series of PVCs occurring one after another (see Fig. 2.2(e)) without any normal beats interspersed [38].

Right Bundle Branch Block (RBBB)

During an RBBB, the right ventricle is not activated by impulses which travel through the right bundle branch. However, the left ventricle is still activated by the left bundle branch. These impulses travel through the myocardium of the left ventricle to the right ventricle to depolarize them. Since conduction through the Bundle of His-Purkinje fibres is faster than conduction through the myocardium, the QRS-complex is prolonged [38].

(23)

Background 10

Figure 2.3: NUUBO’s wearable ECG platform. [40]

The QRS-complex usually shows extra deflection that reflects the rapid depolarisation of the left ventricle followed by the slower depolarisation of the right ventricle. Fig. 2.2(f) illustrates an RBBB episode.

2.4 Wearable ECG Technology

Wearable medical technology can be a practical solution for long term ECG monitoring, see Fig. 2.3 for an example of an industrial wearable solution for ECG monitoring. Wearable ECG monitoring techniques have been a hot area of research for the past decade and still. It is worth mentioning that wearable technology has begun dominating the market of consumable electronic gadgets along with internet of things (IoT). The number of connected wearable devices worldwide is expected to jump from an estimate of 325 million in 2016 to over 830 million in 2020 [39].

The increasing numbers of deployed wearable and IoT systems motivated researchers around the world to develop new methods which fit with the standards of wearable devices. However, many barriers are needed to be overcome in order to reach a true clinical adoption of the wearable medical monitors. Issues like motion artifacts, network connectivity, data processing, integration and clinical decision support are among the challenges faced when considering a wearable solution for medical application [41]. The main focus of this thesis is put on the aspect of clinical decision support for wearable ECG monitors, see Section 1.2.

(24)

Modeling of Electrocardiography 11

3 Modeling of Electrocardiography

3.1 Introduction

The field of physiological modelling attempts to understand living systems from a qualitative perspective. This is done by building mathematical models of how these living systems work.

Thanks to the significant advancement in the field of scientific computing over the past decades, building highly efficient and detailed models of organ systems has been made possible.

Clinical and biomedical research provide large amount of experimental data which can be interpreted using mathematical models.

Mathematical models can be used for reproducing pathological conditions by altering the parameters of these models.

Organ systems in the human body are considered to be highly complex in terms of structure as well as functionality. They usually interact with other organ systems and behave non-linearly.

Physiological models can be divided into two main categories:

• Large-scale/complex models

• Small/simple models

A major obstacle when it comes to physiological models -especially small models- is to decide the level of granularity. In order to estimate the parameters from clinical data, it is necessary to build models with appropriate order of parameters for that clinical data. Control theory and system dynamics areas propose methods for model order reduction. However, these methods may lack model structure preservation, which is essential in physiological models as each parameter has a specific interpretability. The process of translating physiological models into clinical applications takes into account the interpretability of each parameter of the model because clinicians are used to thinking in terms of certain variables such as BP and HR, and it is important to maintain this interpretability.

(25)

Figure 3.1: Normal ECG beat generated using ECGSYN tool

3.2 Models for ECG

Biomedical signals like ECG can be described by nonlinear dynamical systems0methods. Bio-electrical models of the heart are studied in the following three levels:

1. Single-Cell Model: in this level, the ionic current flow through the membranes of the myocardial cell is described [42].

2. Tissue Model: in this level, the tissue, or the cell network, model describes the ionic current flow between different myocardial cells. This process is done in temporal and spatial domains [42].

3. Whole-Heart Model: the whole-heart-model describes three main mechanisms: how the activation currents propagate in the heart, the electrical sources in the heart, and the extracellular electrical potentials on and within the surface of the body. Such model is capable of relating the ECG waveforms to the conduction ve-locity of heart tissue, the action potential and other electrophysiological properties of the cardiac system. This leads to a clinically comparable ECG signals [42]. It is possible to build a mathematical model of the electrical activity of the heart using a different approach. This can be done by modeling the morphology of the ECG signal rather than modeling the physiological process behind it. Such models can be used as a quantitative monitoring tool for the ECG, by studying how the parameters of these models change in time in the parameter space. Section 3.2.1 provides a detailed information on a nonlinear mathematical model used for ECG.

(26)

(a) (b)

Figure 3.2: (a) depicts the waveform of x, y, and z. (b) depicts the phase plot of x and y.

3.2.1 ECGSYN

"ECGSYN is a dynamical model for generating synthetic ECG signals with arbitrary morphologies where the user has the flexibility to choose the operating characteristics" [43].

Presented in [19], the ECGSYN model is based on three coupled ordinary differential equations ˙x=αx−ωy, (3.1) ˙y=αy+ωx, (3.2) ˙z= −

∑

i∈P,Q,R,S,T ai∆θiexp −∆θ 2 i 2b2 i ! − (z−z0), (3.3) where α=1− q (x2₊_y2₎_, _(3.4) ∆θi= (θ−θi)mod(2π), (3.5) −π≤θ=tan−1(y, x) ≤π, (3.6) z0(t) =A sin(2π f2t), (3.7)

ωrepresents the angular velocity of the trajectory. A is the amplitude of the baseline

wander. The variables x, y and z comprise the state vector of the ECG dynamical model that possesses an attracting limit cycle with one period of the solution corresponding to one heart beat, see Fig. 3.3.

(27)

Figure 3.3: Trajectory generated by the ECG dynamical model in three-dimensional state-space

The power spectrum S(f) of the RR-interval is given by Eq (3.8), which takes into account the effects of the respiratory sinus arrhythmia (RSA) as well as the Mayer waves.

S(f) = σ 2 1 2πc2₁exp (f− f1)2 2c2₁ ! (3.8) + σ 2 2 2πc2 2 exp (f− f2) 2 2c2 2 ! ,

where f1and f2are the means and c1and c2are the standard deviations of the distribu-tions.

A time series T(t), with S(f)power spectrum of an RR-interval, can be generated by taking the inverse Fourier transformation of a sequence of complex numbers, which havepS(f)as amplitudes and random phases distributed within 0 and 2π. Then, the time-dependent ω(t)will be given by Eq. (3.10).

T(t) = 60 Ψ +ζ∗ 60δ Ψ2 ∗ 1 $, (3.9) ω(t) = 2π T(t), (3.10)

whereΨ is the mean HR of the generated ECG, ζ is the inverse Fourier transformation of a sequence of complex numbers, which havepS(f)as amplitudes and random phases distributed within 0 and 2π, δ is the standard deviation of the HR, and $ is the standard deviation of ζ.

(28)

The ECG dynamical model reduces every electrocardiogram cycle into a set of fifteen parameters that describe the five fiducial points – P, Q, R, S and T, as shown in Fig. 3.1. Each point is characterized with three parameters ai, bi, θi, i=P, . . . , T, see Table 4.1.

3.3 Model-Based ECG Analysis

Unlike most of the existing techniques, model-based approaches take into account the nature of the underlying dynamics that generate the signal and/or the noise that affects it. As a result to that, once a model has been fitted to a segment of ECG, it can produce a filtered version of the waveform. Moreover, it can derive wave onsets and offsets using the estimated parameters of the model, thus, compressing the ECG and classifying beats. Another advantage for using model-based approaches is that the quality of the fit can be used as a confidence measure with respect to the filtering methods.

Existing techniques for filtering and segmenting ECGs are limited by the lack of an explicit patient-specific model to help isolate the required signal from contaminants. Only a vague knowledge of the frequency band of interest and almost no information concerning the morphology of an ECG are generally used [43].

(29)

Online Model-Based Beat-by-beat Heart Rate Estimation 16

4 Online Model-Based Beat-by-beat Heart

Rate Estimation

This Chapter is dedicated to present a novel method for estimating the instantaneous HR by using the morphological features of one ECG beat at a time. In this work, the utility of the model-based ECG analysis technique is illustrated in online personalized monitoring of the HR. The problem of estimating the HR is reduced to fitting one parameter. The value of this parameter is related to the values of nine parameters of an ECG realistic nonlinear model. The HR is estimated from information of one ECG beat using nonlinear least-squares optimization. The feasibility of the presented method is evaluated on synthetic ECG signals and signals acquired from MIT-BIH databases at Physionet website. In addition, the performance of the presented method is tested under realistic free-moving conditions using a wearable ECG and heart monitor.

4.1 Theory

It is believed that the morphological features of the ECG beats are influenced by the change in the HR. This is because of the relation between the HR and the intra-ventricular conduction [44], also see Section 1.2. By analyzing the ECG signals, which were acquired from healthy subjects at different values of HR, it is notable that the width of the QRS-complex decreases when the HR increases. These results are expected as the speed of the conduction across the ventricles increases together with an augmented HR as the sympathetic tone increases. Therefore, the time for the ventricular depolarization, represented by the QRS complex of the ECG, will become shorter [43].

4.2 Methods

The ECGSYN model, presented in Section 3.2.1, is considered in this Chapter.

Based on the fundamental relation between the width of the QT-interval and the change in the HR [44], also see Section 1.2, it is assumed that the coordinate z of each fiducial point (i.e. ai) is independent of the change in HR and can be assigned a fixed value in the parameter estimation procedure. Moreover, since the parameter set{θi}describes the displacement of each fiducial point from the R-peak, the parameter{θR}is set to zero. Therefore, the parameters are reduced to nine parameters that can be estimated from measured data using a nine-dimensional gradient descent in the parameter space.

(30)

The complexity along with the computational time of the optimization problem were further reduced by introducing a new parameter Ehr, representing the estimated HR, that is mathematically related to the nine parameters (i.e. bP, θP, bQ, θQ, bR, bS, θS, bT, and θT), see (4.3)–(4.5) b0= h bP bQ bR bS bT i , (4.1) θ0= h θP θQ θR θS θT i , (4.2) β=p60/E_hr, (4.3) bβ=β∗b0, (4.4) θβ=θ0 h β12 _β 1 β β12 i , (4.5)

Note the element-wise product in Eq. (4.5). The estimated HR (EHR) is used in Eq. (4.3)– (4.5) to updated the parameters b0and θ0to bβand θβ. The updated parameters are then

used to generate synthetic ECG cycle z(Ehr). Next, z(Ehr)is compared with the observed ECG cycle s(t)by computing the loss function using Eq. (4.6)–(4.7).

This parametrization, along with appropriate upper and low bounds, drastically reduces the computational overhead and is inspired by the initial settings recommended in [20]. The described ECG dynamical model is implemented in MATLAB by the function ecgsyn.m [45]. The MATLAB function ode45.m is used for numerically integrating the set of the ordinary differential equations in (3.1)–(3.3), producing a simulated ECG waveform in Fig. 3.1.

4.2.1 Preprocessing

The first step before estimating the model parameters is to segment the ECG data and isolate every cycle with the corresponding main features, i.e. P, Q, R, S, and T-waves. The number of samples in the isolated cycle does not matter as it can be easily compensated for in a later step by zero padding or truncation. Then, RR-intervals are calculated by subtracting the occurrence times of every two consecutive R-peaks. The next step is to perform beat isolation by considering the starting point for one beat cycle to be the number of samples equal to a half of the corresponding RR-interval before the R-peak and vise-versa.

All isolated beats are then be de-trended to set the baseline and comply with the assump-tion of zero z-offset. Lastly, to eliminate the parameter θR, the R-peak of the observed ECG cycle and the R-peak of the estimated model output are centered at the same point (say at the point zero).

(31)

Table 4.1: Initial values for the parameters of the model

P Q R S T

ai 1.2 -5.0 30.0 -7.5 0.75

bi 0.25 0.1 0.1 0.1 0.4

θi -1.2217 -0.2618 0 0.2618 1.7453

4.2.2 Initial Settings

Before starting the optimization process, the initial value for the HR needs to be specified. It is recommended to use an initial value that corresponds to the lower or higher bound which have been set for the algorithm (see Section 4.2.3 for more details). The sampling frequency is set according to the input signal while the internal sampling frequency is double of the ECG sampling frequency Fs. The means f1and f2are set to 0.1 and 0.25, respectively. The standard deviations c1and c2are both set to 0.01 in Eq. (3.8). The ratio

σ₁2/σ₂2(the LF/HF ratio) is set to 0.5 in Eq. (3.8). A is set to 0.15. Ψ is E_hr, and δ is set to 1

in Eq. (3.9).

The initial values for the fifteen model parameters are summarized in Table 4.1.

In this application, it is assumed that z-offset of the ECG signal is eliminated, which can be achieved by setting(z−z0) =0. Both x and y are set to one.

4.2.3 Optimization

Nonlinear least-squares method solving the optimization problem (4.6), (4.7) is utilized for model parameter estimation

min Ehr kf(Ehr)k22=min Ehr (f₁2(Ehr) + f22(Ehr) + · · · + fn2(Ehr)), (4.6) f(Ehr) =       s1(t1) −z1(Ehr) s2(t2) −z2(Ehr) .. . sn(tn) −zn(Ehr)       , (4.7)

where s(t)is the observed ECG beat, z(Ehr)is the ECGSYN model fit to the ECG signal, and n is the number of samples. The MATLAB function lsqnonlin.m is employed for solving this problem. Choosing the trust-region-reflective algorithm and specifying proper constraints for the optimizer (e.g. upper and lower bounds) appears to be critical for obtaining acceptable performance. It is noticed that setting the maximum and minimum change in variables for finite-difference gradients in a proper way resulted in a significant performance improvement of the optimizer. Adding a feedback from the previous iterations to influence the initial values of the parameters, see Section 4.3, also enhances the performance of the optimization process.

(32)

4.2.4 Ground Truth

Both MIT-BIH Arrhythmia Database [46] and the MIT-BIH Normal Sinus Rhythm Database [47] contain RR-intervals information as part of their annotation files. MIT-BIH Arrhythmia Database Directory provides HR ranges for every entire individual record but not on beat-by-beat basis.

To construct a set of true values of the heart rates, the RR-intervals information are used to calculate the HR for beat i according to

Hi=60FsR_int−1, (4.8) where Hiis the HR for beat i and Rintis the difference between the incidence time of the R-peak for beat i and the incidence time of the R-peak for beat i+1.

For the ECG signals acquired using the MAX-ECG-MONITOR, the true HR values are established using the RR-interval information. This was done by applying the R-peaks detector in [48]. Moreover, five HR measurements of the subject were taken before each test using a Samsung S8 plus mobile phone to measure the resting HR.

4.3 Results

In order to evaluate the performance of the presented method, the algorithm is first tested on realistic synthetic ECG signals at different HRs produced by the ECGSYN tool. The next step is to test the method on real-life ECG signals, and, for that end, the large online bank of physiological signal datasets, provided at Physionet website [47], is used. Finally, ECG signals acquired using the MAX-ECG-MONITOR were used to evaluate the performance of the presented method. Each ECG cycle was isolated to perform a beat-by-beat HR estimation test in order to capture the HR variability.

(33)

Figure 4.1: The effect of the sampling frequency Fson the performance of the HR estimator. (a) Results for HR estimation using synthetic ECG signals sampled with Fs=512Hz, (b) Fs=360Hz and (c) Fs=60Hz.

Figure 4.2: The effect of adding uniformly distributed measurement noise n(t)to the ECG beat data on the performance of the HR estimator, Fs=360Hz. (a) Results for HR estimation using synthesised ECG signals with n(t)= 0 mV, (b) n(t)= 0.1 mV and (c) n(t)= 0.2 mV.

(34)

4.3.1 ECGSYN Signals

ECGSYN tool gives the researcher the control over several parameters such as the mean HR and the standard deviation of HR. the method was tested using ECG signals with HR in the interval 15–160 BPM.

Every test suite consists of a series of 25 ECG beats (averaged to one beat in later step) per every specific heart rate from 15 to 160 BPM. For every test suite, one of the two parameters (the sampling frequency Fsor the additive uniformly distributed measurement noise n(t)) is varied, and the other one is fixed. The standard deviation of HR was set to zero, and the initial value of Ehrwas set to 15 for all tests.

Results are illustrated graphically in Fig. 4.1 and Fig. 4.2, with the numerical values summarized in Table 4.2, and discussed in Section 4.4.

4.3.2 MIT-BIH Database

The method is tested on real-life ECG signals by using the MIT-BIH Arrhythmia Database [46] that contains 48 recordings of two-channel ambulatory ECG obtained from 47 different subjects for around 30 minutes long and studied by the BIH Arrhythmia Laboratory. The MIT-BIH Normal Sinus Rhythm Database [47] was also used for testing the method. It includes 18 recordings of long-term ECG obtained from subjects who were found to have had no significant arrhythmia.

Beats with morphological defects, e.g. left and right bundle branch block beats, prema-ture beats, ectopic beats and premaprema-ture ventricular contractions, were excluded from the test suites. Moreover, beats with arrhythmia, such as atrial fibrillation, were also excluded from the test suites. On the other hand, beats with rate abnormalities, such as tachycardia or bradycardia, were included. Tests were performed on the modified Lead II ECG signals only, as the morphologies in the used ECG model are modeled after Lead II. It is worth mentioning that ECG Lead II is widely used in wearable ECG solutions as it gives a good view of the P-wave and it is most commonly used to record the rhythm strip [49].

Thanks to the provided annotation files, every beat was isolated, segmented, and pre-processed, as described in Section 4.2.1, before undergoing a beat-by-beat HR estimation. For every heart beat, the model-based HR estimation was performed twice, with different initial values for the parameter Ehr.

Overall, a total number of 1171 individual beats acquired from 17 different patients were processed. HR estimation results for signals acquired from MIT-BIH Normal Sinus Rhythm Database are depicted in Fig. 4.3(a), while Fig. 4.3(b) illustrates the results for signals acquired from MIT-BIH Arrhythmia Database. Results for both databases are illustrated numerically in Table 4.3 and discussed in Section 4.4.

(35)

(a) (b)

(c) (d)

Figure 4.3: Results of the model-based HR estimation of the beats acquired from MIT-BIH Normal Sinus Rhythm Database (a) and MIT-BIH Arrhythmia Database (b). Points marked with (o) represent the true HR values, while (x) and (+) represent the HR estimates with initial value of E_hrequal to 15 and 60, respectively. (c) Results of the model-based HR estimation of the beats acquired from the MAX-ECG-MONITOR. Points marked with (o) represent the true HR values, while (x) represent the HR estimates with initial value of Ehrequal to 60. The tests were performed in three different physical statuses which are Sitting, Standing, and Running as illustrated by the vertical lines. (d) The HR estimates for a segment of the ECG test signal 106 acquired from the MIT-BIH Arrhythmia Database. The correlation coefficient between the HR estimates and true values is 0.6337.

(36)

Figure 4.4: The MAX-ECG-MONITOR deployed on the chest of the subject for wearable ECG acquisition.

4.3.3 MAX-ECG-MONITOR

The MAX-ECG-MONITOR evaluation and development platform is a wearable solution for analyzing and accurately tracking ECG signals to provide valuable insight for clinical and fitness applications [50]. Fig. 4.4 illustrates the wearable setup of MAX-ECG-MONITOR.

Three tests were performed on one subject while sitting, standing-up, and running, by recording the ECG signals for one minute per test. A total number of 294 ECG beats were processed for HR estimation using the model-based method. The ECG signals were acquired with a sampling frequency Fs=128Hz. Fig. 4.3(c) illustrates the results of the three ECG tests.

4.3.4 Quantitative Evaluation

As for the approach of the performance evaluation, the performance indices in this work are not compatible to the positive predictability and sensitivity parameters [51]. Therefore, to quantify the performance of the presented method, the Root Mean Square Error (RMSE) was used:

V=r ∑

n

i=1(Ei−Oi)2

n (4.9)

as performance index. In (4.9), E refers to the estimates, O refers to the observations or the true values of the heart rates, and n refers to the total number of beats.

Table 4.2 and Table 4.3 show the RMSE values for ESGSYN signals as well as MIT-BIH database records, respectively. In the MAX-ECG-MONITOR tests, the RMSE values for the sitting, standing-up, and running tests are V=6.9174, V=3.7804, and V=9.5156, respectively.

(37)

Table 4.2: RMSE for ECGSYN signals with initial Ehr=15 . n(t)[mV] and Fs[Hz] were varied. RMSE for heart rates within the range 15–120, the range 120–160 and total RMSE.

ECGSYN Signals

Sampling rate Fs n(t) RMSE15−120 RMSE120−160 RMSET

512 0 1.2907 1.3733 1.3257

360 0 2.1831 2.0798 2.1506

360 0.1 2.9516 46.0554 24.8495

360 0.2 3.9610 54.6491 29.4973

60 0 2.3542 4.5238 3.0929

Table 4.3: Root Mean Squared Error for every MIT-BIH Database individual record with initial Ehr= 15 (RMSE15) and 60 (RMSE60).

Arrhythmia Database Normal Sinus Rhythm Database Record RMSE15 RMSE60 Record RMSE15 RMSE60

103 20.2833 25.2434 16272 12.0266 11.2904 105 28.3592 18.3208 16483 12.5322 7.2598 106 8.9791 13.9306 16539 11.4925 7.6650 113 21.3918 21.4623 16773 7.7009 8.7237 116 5.9964 6.2839 16786 3.7392 -119 14.9614 16.8384 16795 11.8204 -123 15.1424 20.5047 18184 13.1047 10.1677 209 8.8182 7.6326 - - -215 23.2076 11.3170 - - -220 33.3000 37.4901 - -

-4.4 Discussion

4.4.1 ECGSYN Signals

The test results of the HR estimation algorithm on the ECGSYN signals, illustrated graphically in Fig. 4.1 and Fig. 4.2 and numerically in Table 4.2, show an accurate and stable performance for the simulated HR within the range 15–160 BPM with V≤2.4 in all cases, when additive uniformly distributed measurement noise n(t)is zero.

Fig. 4.1, as well as Table 4.2, show that a decrease in the ECG sampling frequency leads to a slight drop in the performance of the HR estimator. Moreover, introducing additive uniformly distributed measurement noise n(t)into the ECG data results in a significant drop in the performance of the HR estimator, especially for the HR values higher than 120 BPM, as shown in Fig. 4.2 and Table 4.2.

(38)

It is worth mentioning that the optimization method used for the model-based HR estimation was found to be highly susceptible to local minima. Fig. 4.5 illustrates the effect of local minima on the performance of the model-based HR estimation, and how a slight change in the HR of the input ECG signal can result in a poor HR estimate. However, the problem of local minima can be overcome by carefully choosing the initial value of the parameter Ehr. To do that, the initial value of the parameter Ehris changed continuously to be equal to the arithmetic mean of the previous HR estimates for the test input ECG signal.

Fig. 4.6 confirms the significant improvement in the HR estimates as well as a huge reduction in the resultant squared norm of the residual (RESNORM) for the same test ECG input as in Fig. 4.5.

Continuously setting the initial value of Ehr to be equal to the arithmetic mean of the previous HR estimates might not produce the best estimates, as the effect of previous poor estimates might accumulate and affect the following estimation process over time.

4.4.2 MIT-BIH Database

The results of MIT-BIH Arrhythmia Database, illustrated graphically in Fig. 4.3(b) and numerically in Table 4.3, show inter-subject variation in performance while maintaining a strong correlation with the true HR values. Fig. 4.3(d) illustrates a part of the ECG test signal 106 acquired from the MIT-BIH Arrhythmia Database, where the HR of the patient increases gradually. It is notable that the results of the HR estimator, although not very accurate, correlate with the true values of HR with the correlation coefficient of 0.6337.

On the other hand, the results of MIT-BIH Normal Sinus Rhythm Database ECG signals, illustrated graphically in Fig. 4.3(a) and numerically in Table 4.3, show higher accuracy as well as robustness to inter-subject variability compared to the MIT-BIH Arrhythmia Database test signals.

One factor that affects the performance of the HR estimator is the amplitude scale of the ECG signal. Different signal acquisition units may produce results with different amplitude levels defined by the amplifier gain. While the gain is included in the description of each record within the MIT-BIH Databases, a deviation of the amplitude gain may occur within the same record as illustrated in Fig. 4.8(a). Fig. 4.8(b) shows how a variation in the amplitude gain affects the performance of the HR estimator. However, this effect can be neutralized by normalizing the amplitude scale of the ECG beat as demonstrated in Fig. 4.8(c).

(39)

Figure 4.5: The effects of local minima on the performance of the HR estimator. (a) the estimated heart rates for 51 beats with true HR between 96 and 101 covered with the step 0.1, the initial value of E_hrwas set to 15, Fs= 360 Hz, and n(t)= 0 mV. Points marked with (o) represent the true HR values, while (x) represent the HR estimates. (b) shows the resultant squared norm of the residual for

optimization procedures in (a).

Figure 4.6: Tuning the initial value of E_hr improves the performance of the HR estimator by overcoming local minima issue illustrated in Fig. 4.5. (a) the estimated heart rates for 51 beats with true HR between 96 and 101 covered with the step 0.1, Fs= 360 Hz and n(t)= 0 mV. Points marked with (o) represent the true HR values, while (x) represent the HR estimates. (b) the resulting squared norm of the residual for optimization procedures in (a).

(40)

Figure 4.7: The difference in morphology between one ECG beat taken from the MAX-ECG-MONITOR sitting test signal with true HR equals to 86 BPM (in Blue), and a synthetic ECG beat with true HR equals to 86 BPM generated by the ECGSYN model (in orange). It also shows the same ECG beat generated by ECGSYN model after being modified by tuning its parameters to fit with the MAX-ECG-MONITOR beat (dashed line in yellow).

4.4.3 MAX-ECG-MONITOR

The results of the HR estimation tests on data from the MAX-ECG-MONITOR, see Fig. 4.3(c) and numerical values in Section 4.3.4, show a good accuracy, especially in the stand-up position with an RMSE value V=3.7804. However, the level of accuracy drops in both sitting and running tests, scoring V=6.9174 and V=9.5156, respectively. While this drop can be attributed to the presence of motion artifacts in the case of the running test and the lack of digital denoising, it does not hold in the case of the sitting test.

A combination of the following two points could explain the accuracy drop:

• As illustrated in Fig. 4.4, the MAX-ECG-MONITOR uses a ground-free ECG elec-trodes configuration. This means that the produced ECG waveform is different from the conventional ECG Lead II waveform. The ground-free ECG setup results in an ECG waveform with suppressed P and T-waves, as illustrated in Fig. 4.7.

• The parameters ai, i=P, . . . , T related to the amplitude of the ECG fiducial points are assigned fixed values and not estimated in the model-based HR estimator. Consequently, the algorithm will attribute a change in the amplitude to a change in HR.

(41)

A possible solution would be to personalize the model metrics by assigning relative values to the parameters ai, i=P, . . . , T as illustrated in Fig. 4.7 by the ECG signal in dashed line. Unfortunately, this was proven inefficient for the following two reasons:

• As explained in Section 4.1, the QT-interval is related to the change in HR. However, the QRS complex tends to be least affected by the HR and therefore, shows small changes with the change in HR. This is due to the relatively short time that it takes for the action potential to travel through the Purkinje fibers between the endocardium and the epicardium, which is represented by the QRS complex. This was further investigated in [21].

• The sampling frequency of the MAX-ECG-MONITOR is Fs= 128 Hz, which is relatively low and narrows the bandwidth of the HR-related QRS changes in morphology. This leads to the optimization procedure mostly diverging to the lower or the upper bounds of the fitted parameters, which in this case is the HR. To address this issue, a model-based HR estimation method, tailored to the MAX-ECG-MONITOR by utilizing all the 15 parameters grouped in a suitable way, can be explored in a future work.

4.4.4 Limitations

Currently, the model-based HR estimation method, as presented, is limited by the following factors:

1. Lead configuration: the performance of the presented method is highly dependent on the ECG morphology represented by a specific lead configuration. Therefore, different lead configurations, as well as different lead placement relative to the heart, would lead to different HR estimates. The method was tested on modified Lead II ECG signals only. Moreover, the ECGSYN tool generates ECG signals similar to those acquired using Lead II configuration illustrated in Fig. 3.1. 2. Noise component: the signal-to-noise ratio (SNR) was found to be inversely

pro-portional to the RMSE, as shown in Table 4.2.

3. Abnormal morphology: since the model-based HR estimation method makes use of the slight variations in the morphological features of the tested ECG beat, it is basically assumed that these features are not affected by any source other than the change in HR. Therefore, pathological-related abnormalities are not taken into account so far.

4. Beat isolation: it is essential to isolate every ECG beat before performing the mea-surements. Isolation process can be performed using the RR-interval information or ECG segmentation methods [52].

(42)

(a)

(b) (c)

Figure 4.8: (a) A one-minute segment of the ECG record 106 acquired from the MIT-BIH Arrhythmia Database, where the amplitude of the signal changes during the process of the ECG acquisition. (b) HR estimates vs. true values plot for the test sample shown in (a). A significant change in the performance of the HR estimator after the point where the amplitude of the ECG signal changes which denoted by the black arrow is shown. (c) Amplitude scaling of the ECG segment normalizing the amplitude scale of the ECG segment shown in Fig. 4.8(a) cancels out the effect of the changing amplitude shown in (b).

(43)

Patient-Specific ECG Monitoring by Model-Based Stochastic Anomaly Detection 30

5 Patient-Specific ECG Monitoring by

Model-Based Stochastic Anomaly

Detection

This Chapter is dedicated to present a novel model-based stochastic anomaly detection approach to ECG monitoring. The utility of the presented method is demonstrated on patient data recorded in clinical setting.

A beat detector and classifier producing good results for a certain training database can perform inconsistently, when dealing with a different patient cohort thus effectively preventing reliable solutions for automated ECG processing. Individualized detection of ECG beat morphology deformation addresses this issue and can be implemented using a combination of model-based techniques and stochastic anomaly detection methods.

5.1 Theory

In order to individualize the ECG analysis process, patient-specific profiles are required. Given a mathematical model, these profiles comprise the estimated parameter distribu-tions corresponding to normal ECG beats. The calibration of these profiles is done by estimating the amplitude-related parameters of the ECG dynamical model, see Section 3, by fitting it to the ECG beats registered under normal conditions. A probability distri-bution function (PDF) is then estimated for the specific patient describing how those parameters vary.

Initial estimates of time and amplitude-related parameters as well as HR are required for generating synthetic ECG signals, see Section 3.2.1. To reduce the overhead of the optimization problem, only the amplitude-related parameters are fitted and tested for anomaly detection. This was made possible by assigning the time-related parameters fixed values. However, fixing the HR value renders the fitting process inefficient due to the relation between the HR and the width of the ECG cycle, see Section 4.1, hence the true z coordinates for each fiducial point will be shifted from the model-based ones. To overcome this issue, period normalization is employed, see Section 5.2.2.

Given a patient-specific profile, isolated ECG beats can be tested to decide whether they are normal or not by utilizing an anomaly detection algorithm. Anomaly detection is the process of identification of rare items (anomalous observations) in a data set, i.e. data points that are suspicious as being significantly different from the majority of the data [53], [54]. Typically, the detected anomaly will correspond to some kind of

(44)

Patient-Specific ECG Monitoring by Model-Based Stochastic Anomaly Detection 31

problem which in the present case is a pathological condition. In this work, a stochastic anomaly detection algorithm proposed in [23] and successfully applied to the analysis of eye-tracking data in [24] is utilized.

5.2 Methods

The ECGSYN model, presented in Section 3.2.1, is considered in this Chapter.

The ECG signal provides direct information on the conduction of the electrical impulses in the heart. In other words, the incidence of any type of abnormality in the ECG indicates either an artifact or a pathological episode, see Section 2.3. The timing and amplitude features of ECG beats are used as indicators for such abnormal activities. In this study, the coordinate z of each fiducial point (i.e. ai) is utilized to capture the morphological changes in the ECG beat. Therefore, the model individualization is reduced to estimating five parameters from the measured ECG by means of five-dimensional gradient descent in the parameter space. The complexity and, consequently, the computational time of the optimization problem are further reduced by assigning the time-related parameters of each fiducial point (i.e. biand θi) to fixed values, which opera-tion is performed through period normalizaopera-tion, see Secopera-tion 5.2.2. This parametrizaopera-tion, along with appropriate upper and lower bounds, drastically reduces the computational overhead.

As mentioned in the previous Chapter, the described ECG dynamical model is imple-mented in MATLAB by the function ecgsyn.m [45]. The MATLAB function ode45.m is used to numerically integrate the set of the ordinary differential equations in (3.1)–(3.3), producing a simulated ECG waveform in Fig. 3.1.

5.2.1 Preprocessing

The first step before estimating the model parameters is to segment the ECG data and isolate every cycle with the corresponding main features, i.e. P, Q, R, S, and T-waves. The number of samples in the isolated cycle is of minor importance as it can easily be compensated for in a later step by zero padding or truncation. Then, RR-intervals are calculated by subtracting the occurrence times of every two consecutive R-peaks. The next step is to perform beat isolation by considering the starting point for one beat cycle to be the number of samples equal to a half of the corresponding RR-interval before the R-peak and vise-versa. All isolated beats must then be de-trended to set the baseline and comply with the assumption of zero z-offset. Lastly, the R-peak of the observed ECG cycle and the R-peak of the estimated model output must be centered at the same point, e.g. at zero.

Now the isolated ECG beats undergo both amplitude and period normalization. The amplitude normalization step results in each isolated ECG beat having the peak-to-peak voltage Vppequal to 1.573 V, which is the Vppvalue for an ECG beat generated by the