# Compelling new electrocardiographic markers for automatic diagnosis

Cristina Rueda<EMAIL_ADDRESS>Itziar Fernández  Yolanda Larriba
Alejandro Rodríguez-Collado  Christian Canedo  Department of Statistics and
Operations Research, Universidad de Valladolid, Paseo de Belén 7, 47011,
Valladolid, Spain

###### Abstract

Background and Objective: The automatic diagnosis of heart diseases from the
electrocardiogram (ECG) signal is crucial in clinical decision-making.
However, the use of computer-based decision rules in clinical practice is
still deficient, mainly due to their complexity and a lack of medical
interpretation. The objetive of this research is to address these issues by
providing valuable diagnostic rules that can be easily implemented in clinical
practice. In this research, efficient diagnostic rules friendly in clinical
practice are provided. Methods: In this paper, interesting parameters obtained
from the ECG signals analysis are presented and two simple rules for automatic
diagnosis of Bundle Branch Blocks are defined using new markers derived from
the so-called FMMecg delineator. The main advantages of these markers are the
good statistical properties and their clear interpretation in clinically
meaningful terms. Results: High sensitivity and specificity values have been
obtained using the proposed rules with data from more than 35000 patients from
well known benchmarking databases. In particular, to identify Complete Left
Bundle Branch Blocks and differentiate this condition from subjects without
heart diseases, sensitivity and specificity values ranging from 93% to 99% and
from 96% to 99%, respectively. The new markers and the automatic diagnosis are
easily available at https://fmmmodel.shinyapps.io/fmmEcg/, an app specifically
developed for any given ECG signal. Conclusions: The proposal is different
from others in the literature and it is compelling for three main reasons. On
the one hand, the markers have a concise electrophysiological interpretation.
On the other hand, the diagnosis rules have a very high accuracy. Finally, the
markers can be provided by any device that registers the ECG signal and the
automatic diagnosis is made straightforwardly, in contrast to the black-box
and deep learning algorithms.

###### keywords:

FMM model , ECG waves , Diagnostic rule , Bundle Branch Block

##  1 Introduction

The relevance of the automatic diagnosis of diseases is beyond the debate. In
the case of cardiovascular diseases, it can be a matter of life and death.
Many researchers have developed automatic rules for diagnoses, but their use
in clinical practice is still deficient. One reason is that many cases are
based on black-box algorithms that physicians do not trust because they lack a
medical interpretation. Furthermore, the validation process is often not
clean, as it depends heavily on the selection of signals and patients
(training and test) and on the preprocessing stage. In addition to these
limitations of most studies, there are those derived for the lack of consensus
on the definition of electrocardiographic features and criteria to diagnose
certain diseases using such features.

The proposal in this paper is to prove that we can get around all these
problems by the definition of new markers using the FMMecg delineator [1]. The
FMMecg delineator describes the fragment of an ECG signal corresponding to a
heartbeat as the combination of five waves corresponding to the fundamental
waves in a heartbeat: P, Q, R, S, T. The details of the model specifications,
such as the restrictions and the estimation algorithm are given in [1]. In
addition, that paper reveals the methodology’s potential to describe a variety
of non-pathological, noisy, and pathological ECG patterns. The FMMecg wave
decomposition for a typical non-pathological pattern, from a left lead signal,
is shown in Figure 1.

Figure 1: The five waves: P, Q, R, S, T identify by the FMMecg approach for a
typical non-pathological pattern. Data from patient $37$ from the Georgia
database.

In particular, using the FMMecg model parameters, two powerful markers, omeR
and omeS, are presented here and are used as rules for the diagnosis of the
Bundle Branch Block (BBB), which is a defect in the electrical conduction
system of the heart. In such cases, ventricular enlargement or hypertrophy
occurs, and the QRS complex will widen, deform and prolong [2].

On the one hand, a left BBB (LBBB) happens when the activation of the left
ventricle of the heart is delayed, which causes the left ventricle to contract
later than the right ventricle. In particular, LBBB is associated to a higher
risk of suffering different cardiovascular diseases, including a high
mortality ascribed to acute myocardial infarction. Some authors also noted the
crucial need for a prompt identification of LBBB in determining the best
health strategy and future for patients. Specifically relevant is the
detection of patients with complete LBBB (CLBBB) who may respond positively to
cardiac resynchronization therapy. [3] is a good review of the state of the
art in the LBBB disease.

On the other hand, the right BBB (RBBB) is a blockage of electrical impulses
to the heart’s right ventricle and it is one of the most common
electrocardiographic abnormalities that is often detected in asymptomatic
patients [4]. Although the RBBB has been associated with fewer complications
for cardiovascular disease in comparison to LBBB [5], the complete RBBB
(CRBBB) has an impact on patients with other diseases, including acute
myocardial infarction, and the appearance of CRBBB in patients with other
cardiovascular diagnoses worsen their prognosis [6], so it is advisable to
develop tools for early detection.

The electrographic diagnosis criteria for CLBBB and CRBBB depends on the
location of the lead analyzed. In this paper, we initially consider only left-
sided leads, specifically: I, II and V5. The typical ECG signals, in the
presence of complete BBB (CBBB), are shown in Figure 2 for lead I. In
agreement with the figure, the widely mentioned CLBBB diagnostic criteria for
such leads can be summarized in four points as follows [7, 8, 9, 3]:

  1. 1.

Prolonged QRS complex duration. Even longer for male than female.

  2. 2.

Notched and slurred R-waves.

  3. 3.

Prolonged R wave peak time.

  4. 4.

T wave inversions and ST-segment depression.

It worth mentioning the Strauss criteria [9] as it is widely used and closely
related to the criteria recommended by different institutions as the American
Heart Association Electrocardiography and Arrhythmias Committee, the Council
on Clinical Cardiology/American College of Cardiology Foundation/Heart Rhythm
Society [8]. The Strauss criteria is based on the first two points in the
enumeration above.

Figure 2: Typical ECG morphology in complete BBB. ECG signals are represented
as grey points and FMMecg predictions as a red line. Left panel (A): data from
patient 2087 (CPCS database) with CLBBB diagnosis. Right panel (B): data from
patient 10131 (PTB-XL database) with CRBBB diagnosis.

Concerning CRBBB, although there is no consensus in its diagnosis criteria
[6], most researchers assume the following points as typical CRBBB signs for
left-sided leads:

  1. 1.

Large QRS complex duration.

  2. 2.

Wide, negative, and slurred S wave.

  3. 3.

Large S wave duration.

The ECG findings described above for CLBBB or CRBBB diagnoses are apparently
simple but nothing is further from the truth. First, there is no consensus on
the definition of the QRS complex duration or the R wave peak time. The T wave
inversion or the ST-segment depression are also vague descriptions [10]. There
is no standard definition of the QRS notch and slur patterns [11], and no
simple criteria exist to determine their number and location. Moreover, many
of the well-known ECG features, such as the QRS complex duration, strongly
depend on the heart rate [12]. In particular, this association confounds the
relationship between QRS durations and sex, which affects the definition of
the criteria for diagnosis. In addition, a considerable QRS duration is a
typical characteristic in both CLBBB and CRBBB which makes the difference a
problematic task.

Finally, significant systematic differences among engineering solutions of
manufacturers of automated ECG have been detected due to preprocessing and
denoising stages and the algorithms used. Specifically, such points as the end
of the R wave and the end of the QRS complex have no precise definition.
Consequently, different results might be expected of the QRS duration and
other features of the same underlying ECG waveform by different algorithms
[10].

In the last few years, together with a considerable increase in electronic
health databases, such as Physionet’s Physiobank [13], machine learning
algorithms, including deep learning, have become popular in the automatic
diagnosis of cardiac abnormalities, see [14, 15, 16] among others. In
particular, few studies have focused exclusively on the automated
identification of BBB, we have found [17] and [18]. In those cases, the
classifiers achieved sensitivity and specificity values above 95% for RBBB and
98% for LBBB. However, these algorithms have been developed and tested using
ECG recordings from only $47$ subjects from the MIT-BIH Arrhythmia database
[19]. Furthermore, there has been recently a proliferation of methods for the
automatic diagnosis of CLBBB patients known to have a good clinical response
to cardiac resynchronization therapy. For example, [20] proposed a method
based on wavelet analysis, obtaining high and moderate values of sensitivity
and specificity (92.9% and 65.1%, respectively). A more accurate diagnosis was
achieved in [21], which combines a random forest classifier and a neural
network with sensitivity and specificity values of 91.7% and 88.7%,
respectively. Both algorithms were trained and validated using data on 600
patients from the Multicenter Automatic Defibrillator Implantation
Trial—Cardiac Resynchronization Therapy (MADIT-CRT) database [22]. It is
important to note that these classifiers use the same database for training,
test and validation and that the number of patients is limited, which
increases the risk of overfitting and failure when used in any other database.
Another well-known drawback of machine learning is the lack of any
physiological interpretation of the rules. These limitations have prevented
its use in clinical practice.

The new markers, $omeR$ and $omeS$, presented here have a consistent
mathematical and morphological meaning, with a definition independently of the
electrocardiograph or expert, and this consequently makes them universal and
quite useful in diagnosis. In particular, omeR captures the width of the R
wave, is sex-independent, is closely related to the duration of QRS complex
and the R peak time, and emerges as a measure of the degree of LBBB severity.
Moreover, omeS captures the relevance of negative S waves. A complementary
rule combining omeR and omeS is useful to diagnose CBBB.

The potential of new markers and the corresponding diagnostic rules are
validated from such benchmarking databases as PTB-XL [23], Georgia (see the
2020 PhysioNet/Computing in Cardiology Challenge [24] for details) and CPCS
[25]. The distribution of omeR by diagnosis and database, as well as by sex
and age, is described. In particular, the focus is centered on the PTB-XL
results, as it is one of the few databases where BBB are well classified as
complete/incomplete, as well as left and/or right. Finally, to facilitate the
use of the new tools for the reader, an app has been developed that has the
registered fragment as input and the new indexes as output. It is available on
https://fmmmodel.shinyapps.io/fmmEcg/.

##  2 Methods

omeR and omeS are defined as functions of the FMMecg basic parameters. The
FMMecg model is a particular element of the family of Frequency Modulated
Möbius (FMM) models, which was developed to analyze oscillatory signals and is
briefly presented below.

###  2.1 The FMM approach

Oscillatory signals are defined in the time domain and, without loss of
generality, it is assumed that the time points are in $[0,2\pi)$. In any other
case, transform the time points $t^{\prime}\in[t_{0},T+t_{0}]$ by
$t=\frac{(t^{\prime}-t_{0})2\pi}{T}$.

A mathematical term describing an FMM wave is defined as follows:

$W(t,A,\alpha,\beta,\omega)=A\cos(\beta+2\arctan(\omega\tan(\frac{t-\alpha}{2}))),$

where, $A\in\Re^{+}$ is a scale parameter measuring the wave’s amplitude,
$\alpha\in[0,2\pi)$ is a location parameter, while $\beta\in[0,2\pi)$ and
$\omega\in(0,1]$ are parameters describing the shape. $\beta$ is a skewness
parameter that also indicates upward and/or downward peak direction, and
$\omega$ measures the width.

The FMMm model is defined as an additive signal plus error model, where the
signal is defined as a sum of waves as follows:

$X(t_{i})=M+\sum_{J=1}^{m}W_{J}(t)+e(t_{i});$

where $W_{J}(t)=W(t,A_{J},\alpha_{J},\beta_{J},\omega_{J});J=1,\dots,m$,
$M\in\Re$ is an intercept, and the error term, which accounts for the noise,
is assumed Gaussian with a common variance.

In addition, restrictions among the parameters are incorporated depending on
the application. The estimators of the parameters are obtained using maximum
likelihood. Refer to [26] for the details.

Furthermore, other important parameters of practical use are the peak and
trough times for each wave $K$, denoted by $tmax_{K}$ and $tmin_{K}$,
respectively, and their predicted values denoted by $fmax_{K}$ and $fmin_{K}$,
respectively, as follows:

$\displaystyle tmax_{K}$
$\displaystyle=\alpha_{K}+2\arctan(\frac{1}{\omega_{K}}\tan(\frac{-\beta_{K}}{2}))$
$\displaystyle tmin_{K}$
$\displaystyle=\alpha_{K}+2\arctan(\frac{1}{\omega_{K}}\tan(\frac{\pi-\beta_{K}}{2}))$
$\displaystyle fmax_{K}$ $\displaystyle=M+\sum_{J=1}^{m}W_{J}(tmax_{K})$
$\displaystyle fmin_{K}$ $\displaystyle=M+\sum_{J=1}^{m}W_{J}(tmin_{K})$

and the distances between waves by

$d_{JK}=1-\cos(\alpha_{J}-\alpha_{K});J,K=1,...,m$

In the case in which models describing different leads are considered
simultaneously, a superscript indicating the lead will be incorporated to the
notation.

###  2.2 Morphological interpretation of the FMMecg parameters

A typical non-pathological ECG signal, such as that shown in Figure 1, has
positive P, R and T waves and non positive Q and S waves, which is well
described by the beta FMM parameters: $\beta_{J}\in
arc(2\pi/3,4\pi/3),J=P,R,T$ and $\beta_{J}\notin arc(2\pi/3,4\pi/3),J=Q,S$.
Moreover, the omega and amplitude FMM parameters capture other aspects of the
shape of a typical pattern, such as, a prominent and sharp R wave ($A_{R}$
high and $\omega_{R}$ low); a prominent and flat T wave ($A_{T}$ high and
$\omega_{T}$ moderate); a less prominent and not sharp P wave ($A_{P}$ low and
$\omega_{P}$ moderate) and sharp and not prominent Q and S waves
($A_{Q},A_{S},\omega_{Q},\omega_{S}$ low values).

Specific values for FMM basic parameters are given for different patterns,
including that one, in the supplementary material of [1].

Typical CLBBB and CRBBB patterns can also be described in terms of FMM
parameters. In particular, Table 1 includes the description of the criteria
for BBB diagnosis enumerated in the introduction. Examples that illustrate the
correspondence in Table 1 are included in the A. Besides, Figure 3 illustrates
the configuration of waves and parameters that corresponds to the CLBBB and
CRBBB signals presented in Figure 2.

Brief description | FMM parameter representation
---|---
Broad R wave or prolonged R-peak time | high values of $\omega_{R}$
Prolonged QRS duration | high values of $d_{QS}$
M-shaped, qR-wave |  $fmax_{Q}$ or $fmin_{Q}$ close to $fmax_{R}$
M-shaped, Rs-wave |  $fmax_{S}$ or $fmin_{S}$ close to $fmax_{R}$
Negative S-wave |  $\beta_{S}\in arc(5\pi/3,\pi/3)$.
Broad S-wave | high values of $\omega_{S}$
Large S-wave | high values of $A_{S}$
ST segment depression | high values of $tmin_{S}$.
T-wave inversion |  $\beta_{T}\in arc(5\pi/3,\pi/3)$.
Table 1: Features for BBB diagnostic criteria in leads I, II and V5.

Alternative parameter configurations to those in Table 1 capture also
morphologies with notched R waves. Anyway, notched R waves are often observed
as a consequence of $\omega_{R}$ being high. This property, together with the
potential of $\omega_{R}$ to measure the R peak time and the QRS complex
duration, prompts this latter parameter to be an excellent predictor of CLBBB
by itself, since many experts associate this disease with the simultaneous
combination of this disease occurrence of notches and a prolonged QRS complex.

Figure 3: The five waves of FMMecg fitted to data of complete BBB patients.
Left panel (A): data from patient 2087 with CLBBB diagnosis in CPCS database.
P wave (blue) $\omega_{P}=0.110$ and $\beta_{P}=4.197$; Q wave (green)
$\omega_{Q}=0.070$ and $\beta_{Q}=5.794$; R wave (red) $\omega_{R}=0.077$ and
$\beta_{R}=2.930$; S wave (orange) $\omega_{S}=0.045$ and $\beta_{S}=3.816$; T
wave (purple) $\omega_{T}=0.092$ and $\beta_{T}=1.291$. Right panel (B): data
from patient 10131 with CRBBB diagnosis in PTB-XL database. P wave (blue)
$\omega_{P}=0.116$ and $\beta_{P}=4.654$; Q wave (green) $\omega_{Q}=0.033$
and $\beta_{Q}=0.892$; R wave (red) $\omega_{R}=0.045$ and $\beta_{R}=3.352$;
S wave (orange) $\omega_{S}=0.118$ and $\beta_{S}=0.042$; T wave (purple)
$\omega_{T}=0.198$ and $\beta_{T}=3.675$.

###  2.3 CLBBB and CBBB diagnostic rules

In this paper, ECG fragments from leads I, II, and V5 are considered as input.
A preprocessing stage is done where basic checks are performed in order to
discard high noise signals, correct the trend and perform a simple
normalization of the data. Also, a Pan Tompkins algorithm [27] is used to
locate the QRS complexes used to divide the signal into beats. Moreover, the
imputation of missing values has been done. Details of this preprocessing
stage are given in the B.

Now, using an updated algorithm similar to that described in [1], the signal
is analyzed using the FMMecg. A representative heartbeat is obtained for each
patient and lead, using the medians of the parameters. For such a
representative heartbeat, let us define:

CLBBB RULE: omeR $>0.06$ $\displaystyle
omeR=\max\limits_{L=\\{I,II,V5\\}}\\{\omega_{R}^{L}\\}$

The index developed resembles the Strauss criteria [9] as it accounts for the
QRS duration and R peak time. In fact, omeR is the width of the main positive
wave in the QRS. More sophisticated rules can be defined using parameters that
describe other parts of the morphology, such as the ST-segment depression.
However, other diseases, such as RBBB, also has these characteristic patterns,
and the specificity decreases while the complexity of the index definition
increases. omeR itself is a measure of the likelihood for a given patient to
have the LBBB disorder.

The typical pattern for a CRBBB is a pronounced negative S wave, which is
described with a high value of $\omega_{S}$ and $\beta_{S}$ close to $2\pi$.
However, in our databases, more than 5% of CRBBB patients verify $omeR>0.06$
and would therefore be diagnosed as CLBBB. Then, although more sophisticated
rules that increase the accuracy in the sacrifice of simplicity could be
defined for specific CRBBB diagnoses, we have opted for a simple rule for a
combined CBBB diagnosis that complements the rule defined above. The combined
rule is given as follows:

CBBB RULE: omeR $>0.025$ and omeS $>0.05$ $\displaystyle
omeR=\max\limits_{L=\\{I,II,V5\\}}\\{\omega_{R}^{L}\\}$ $\displaystyle
omeS=\max\limits_{L=\\{I,II,V5\\}}\\{\omega_{S}^{L}\cdot I_{[\beta_{S}^{L}\in
arc(5\pi/3,\pi/3)]}\\}$

where $I_{[\cdot]}$ is an indicator taking a value of $1$ if the argument is
true and $0$ otherwise.

It is interesting to note that the use of median parameters and the maximum of
three values prevent the occasional influence of incorrect parameter
estimators that are derived from individual heartbeats and a single lead.
These errors may happen due to an incorrect definition of the fragment limits
corresponding to a specific heartbeat, signal noise, and a lesser extent, a
misidentification of the R wave.

##  3 Results

A total of $106233$ ECG fragments corresponding to $35411$ patients over
eighteen and with diverse diagnoses, have been analyzed. The data are from the
benchmarking databases PTB-XL (patients with likelihood diagnosis $\geq 80$),
Georgia, and CPCS. These databases have been used widely in the literature, to
be more precise, in the challenge of Computing in Cardiology 2020 [24]. The
duration of the ECG signal is quite variable across databases; only the
fragment corresponding to the first 14 seconds has been considered for each
patient. Only 348 patients (less than 1%), were discarded due to noise
artifacts in the signal. The distribution across databases is as follows: 159
(0.86%) from PTB-XL, 90 (0.88%) from Georgia and 99 (1.48%) from CPCS.

Two numerical analyses are presented in this section. The first is dedicated
to omeR. The distribution of this fundamental marker across diagnosis, sex,
and age is explored. The relation of omeR with the RR interval, the time
elapsed between two successive R waves, is also examined. The second study is
the validation of the rules defined above.

A diagnostic label has been assigned to each patient, using SNOMED CT ontology
(bioportal.bioontology.org/ontologies/SNOMEDCT), as is done for the challenge
of Computing in Cardiology 2020 [24]. Thus, patients with a BBB diagnostic are
labeled with CLBBB, ILBBB, CRBBB, or IRBBB according to the side (Left/Right)
and the degree of the defect (Complete/Incomplete). Note that in the CPCS
database, patients with a code LBBB are labeled as CLBBB, as the code for
ILBBB is not considered. Furthermore, the label CD is assigned to those with a
block diagnostic different to BBB; the NORM label is assigned to patients
diagnosed as non-pathological and OTHER to the rest.

All the patients are classified in one and only one category, even those with
several diagnoses (except for a patient with CLBBB and RBBB labelled as LBBB).

###  3.1 omeR distribution across diagnostic, age and sex

In order to consider omeR as a gold standard to detect anomalous ECG is
necessary to know reference values and normal ranges across diagnostics, sex
and age.

Median values and the percentage interval for omeR are given for patients in
the PTB-XL, Georgia and CPCS databases, across diagnostics in Table 2.

Database | Diagnosis | N | $P_{5}$ | $P_{50}$ | $P_{95}$
---|---|---|---|---|---
PTB-XL | CLBBB | 503 | 0.080 | 0.148 | 0.243
| BBB | 2168 | 0.028 | 0.047 | 0.190
| CD | 2529 | 0.029 | 0.048 | 0.109
| OTHER | 5275 | 0.026 | 0.043 | 0.088
| NORM | 8310 | 0.026 | 0.036 | 0.052
Georgia | CLBBB | 208 | 0.063 | 0.137 | 0.246
| BBB | 1170 | 0.027 | 0.052 | 0.180
| CD | 1024 | 0.026 | 0.045 | 0.099
| OTHER | 6253 | 0.025 | 0.044 | 0.084
| NORM | 1723 | 0.027 | 0.038 | 0.058
CPCS | CLBBB | 226 | 0.057 | 0.133 | 0.236
| BBB | 2044 | 0.028 | 0.043 | 0.144
| CD | 688 | 0.026 | 0.038 | 0.066
| OTHER | 3067 | 0.029 | 0.045 | 0.089
| NORM | 812 | 0.029 | 0.040 | 0.055
Table 2: Percentile range and median for omeR across diagnostic in the PTB-XL,
Georgia and CPCS databases.

Moreover, we have claimed in the introduction that an important property of
omeR is that it is independent of sex. Table 3 shows percentile ranges and
median values for the RR interval and omeR across databases, age and sex. Only
NORM patients are considered, as spurious associations could be found if other
diagnoses were included because they are related to omeR alterations.

The numbers in Table 3 show that the distribution of omeR does not depend on
sex or age, while the distribution of the RR interval does. Increasing values
of the RR interval medians are obtained with increasing age, and for males
against females, as is well documented in the literature [28, 29]. Moreover,
it is also shown that the omeR distributions are much more homogeneous across
databases than those of the RR interval.

|  | RR interval | omeR
---|---|---|---
|  | N | $P_{5}$ | $P_{50}$ | $P_{95}$ | $P_{5}$ | $P_{50}$ | $P_{95}$
PTB-XL | Male | 3774 | 648 | 899 | 1183 | 0.026 | 0.036 | 0.052
| Female | 4536 | 620 | 856 | 1125 | 0.026 | 0.036 | 0.051
| $<25$ | 461 | 636 | 864 | 1172 | 0.025 | 0.035 | 0.047
| $[25,50]$ | 3104 | 634 | 872 | 1146 | 0.026 | 0.036 | 0.051
| $[51,70]$ | 3495 | 645 | 890 | 1157 | 0.026 | 0.035 | 0.051
| $>70$ | 1250 | 592 | 842 | 1132 | 0.027 | 0.037 | 0.054
Georgia | Male | 792 | 618 | 819 | 988 | 0.028 | 0.038 | 0.060
| Female | 931 | 628 | 809 | 983 | 0.027 | 0.038 | 0.056
| $<25$ | 50 | 634 | 771 | 946 | 0.028 | 0.040 | 0.061
| $[25,50]$ | 605 | 621 | 800 | 984 | 0.028 | 0.039 | 0.057
| $[51,70]$ | 791 | 620 | 822 | 988 | 0.027 | 0.037 | 0.058
| $>70$ | 277 | 636 | 834 | 986 | 0.027 | 0.037 | 0.058
CPCS | Male | 322 | 640 | 831 | 1008 | 0.029 | 0.040 | 0.054
| Female | 490 | 629 | 809 | 970 | 0.030 | 0.040 | 0.055
| $<25$ | 52 | 657 | 800 | 944 | 0.035 | 0.042 | 0.054
| $[25,50]$ | 452 | 632 | 806 | 969 | 0.029 | 0.041 | 0.055
| $[51,70]$ | 244 | 654 | 846 | 1015 | 0.029 | 0.039 | 0.054
| $>70$ | 64 | 618 | 857 | 1027 | 0.028 | 0.038 | 0.053
Table 3: RR interval and omeR percentile ranges and median for NORM patients
across gender and age groups in the PTB-XL, Georgia, and CPCS databases.

###  3.2 Rule performance in the PTB-XL, Georgia and CPCS databases

Standard measures of test validation, such as sensitivity (SE) or specificity
(SP), have been calculated using both rules for different types of patients in
Table 4 and 5, respectively.

The target of the rules is patients with CLBBB or CRBBB. However, the
differentiation between complete or incomplete blocks is not homogeneous
across the databases, partially explaining the differences in sensitivity
between them. Moreover, the distribution of patients across other diagnoses
also explains the differences in specificity for the category ALL.

The PTB-XL(1) database is defined as the sub-base containing PTB-XL patients
with a single diagnosis. The PTB-XL(1) is a good reference, as patients have a
unique diagnostic and it discriminates between complete or incomplete blocks.
For this database, the sensitivity of CLBBB (CBBB) rules is 99% (98%), the
specificity is 99% (90%) for NORM patients, and 96% (87%) for patients without
an LBBB (CBBB) diagnostic. Sensitivity and specificity values are very similar
for PTB-XL.

Moreover, in the Georgia and CPCS databases, the slightly lower sensitivity
values obtained for the CLBBB rule, 96% and 93%, respectively, compared with
those of PTB-XL, may reflect the fact that the degree of injury is not well
differentiated in these databases. Specifically, some patients labeled as LBBB
could be ILBBB. Sensitivity values for the CBBB rule are similar to that of
PTB-XL.

On the other hand, the specificity is very high, especially for NORM patients
and the CLBBB rule, where the values range from 96% in Georgia to 99% in PTB-
XL.

| N | SE(CLBBB) | SE(ILBBB) | SE(LBBB) | SP(NORM) | SP(ALL)
---|---|---|---|---|---|---
PTB-XL(1) | 12393 | 99% | 67% | 99% | 99% | 96%
PTB-XL | 18282 | 98% | 88% | 97% | 99% | 88%
Georgia | 10170 | 96% | 68% | 88% | 96% | 84%
CPCS | 6611 | 93% | — | 93% | 98% | 86%
Table 4: CLBBB rule: Sensitivity (SE) and specificity (SP) across subgroups of patients by diagnostic. The ALL category includes all the patients except those with LBBB labels. | N | SE(CBBB) | SE(IBBB) | SE(BBB) | SP(NORM) | SP(ALL)
---|---|---|---|---|---|---
PTB-XL(1) | 12393 | 98% | 27% | 67% | 90% | 87%
PTB-XL | 18282 | 96% | 34% | 74% | 90% | 79%
Georgia | 10170 | 96% | 53% | 72% | 87% | 72%
CPCS | 6611 | 95% | — | 74% | 92% | 77%
Table 5: CBBB rule: Sensitivity (SE) and specificity (SP) across subgroups of
patients by diagnostic. The ALL category includes all the patients except
those with BBB labels.

###  3.3 FMMecg app

The app is freely available on https://fmmmodel.shinyapps.io/fmmEcg/. The
instructions for use are very simple. The app requires the ECG fragment from
either I, II, or V5 leads, including one or multiple heartbeats, recorded, as
input. The output includes: plots of the observed and predicted values, a plot
of the FMM median waves, tables with the median values of the estimated FMM
parameters, the $omeR$ and $omeS$ markers, and the BBB diagnosis. Note that
the values of $omeR$ and $omeS$ calculated in the app are obtained from a
single lead. When that information on multiple leads is available, the values
of the markers can be obtained for each lead independently with the app, and
then the maximum of the marker values is calculated by hand.

In the C provides further details about the app, including multiple-use
examples.

##  4 Discussion

In this paper, new ECG markers are proposed, defined using the FMMecg
delineator. In the main place, omeR, which is related to the QRS duration, a
widely used index that is relevant to know. In particular, the literature has
shown that a large QRS duration is a characteristic associated to different
diseases such as ischemic cardiomyopathy [30], myocardial infarction [31] or
sudden death [32].

Compared with the latter index, omeR has a more precise meaning, is measured
in a normalized scale, is independent of sex, and does not depend on the
measuring device or the researcher, as it is a parameter of a statistical
model that is estimated by maximum likelihood. Normal ranges for omeR have
been provided, as well as ranges across diagnostics, sex, and age. Relevant
differences have been observed only across diagnostics. Then, the association
of the QRS duration with age and sex, which is documented in the literature
[28, 29], can be explained by the association of these three variables with
the RR interval length.

Furthermore, the results in this and previous works [33, 1, 26], have shown
that the FMMecg parameters are estimated with high accuracy, and they are
consistent and robust. In particular, omeR because it is based on the main
median parameter of the most prominent wave and is less affected by noise.
Besides, the FMMecg parameters have much potential for detecting causes and
diagnoses derived from a prolonged QRS as they describe the morphology of the
QRS complex in an exact way, with four parameters describing the shape and
length of each of the waves Q, R, or S. In particular, while a wide R wave is
associated with left blocks, a prominent and negative S wave describes the
typical morphology of a right block. In both cases, the QRS duration is large,
and omeR and omeS differentiate these two diseases automatically from other
patterns.

Despite their simplicity, the ability of the two classifiers proposed in this
paper to correctly identify a BBB is highly satisfactory. The results are
particularly good for the CLBBB rule, which confirms the potential of omeR in
diagnosis. Specifically, the sensitivity and specificity values obtained for
the CLBBB classification are similar or even higher than those obtained from
approaches proposed by other authors, such as [19, 17, 18, 20, 21].
Nevertheless, comparisons with other rules are not fair as none of them has
such universal character or they are not possible to be applied automatically.

It is noteworthy that the rules can be adapted to the information available.
On the one hand, they can be used with data from only one or two leads, and on
the other hand, when the omeR and omeR baseline values are known in advance
for a given patient, the thresholds in the rules can be redefined to be
percentages of these values.

It is also worth pointing out that the databases used to validate the rules
are not explicitly designed for the problem at hand. They are databases
collected and labelled for very different purposes and consequently exhibit an
entirely different distribution of patients across diagnostics, which may
explain the difference in accuracy values across databases. However, they do
coincide in showing much better results in identifying complete left blocks
compared to the complete right or the incomplete ones. Nevertheless, omeR and
omeS could provide an opportunity to study the evolution of the degree of BBB
and even prevent more severe blocks in the future. In fact, according to some
authors, ILBBB (IRBBB) might be a precursor of CLBBB (CRBBB) [8, 34].

The usefulness of omeR, omeS and other markers derived from the FMMecg basic
parameters is promising and vast and goes much further than the diagnosis of
BBB and anything we can say in this paper. In particular, these markers are
much more easily registered than serial ECG, so current values are compared
with those obtained at other moments in time that can run from one hour to
many years before. On the other hand, the automatic registration of the
markers, not necessarily supervised by an expert, facilitates many tasks. For
instance, monitoring omeR after a trascatheter aortic valve replacement could
be essential to identifying an LBBB that persists more than 72 hours and
assist in prophylactic decision making [3]. Alternatively, it could also
identify the presence of a new BBB in patients with an acute myocardial
infarction, in which pacemaker insertion may not be not beneficial [35].

Finally, the derivation of consistent markers for the automatic interpretation
of the ECG would reduce the observed differences in the interpretation of the
electrocardiographic abnormalities among health professionals, the first cause
being a manual reading of the ECGs. Unifying the criteria would help improve
the competence of non-cardiologist physicians and achieve a better and
automatic diagnosis.

Several lines of work for future research emerge from this paper. On the one
hand, the definition of rules for the diagnosis of other diseases.
Specifically, the automatic diagnosis of myocardial infarction, using FMMecg
parameters, will be our immediate challenge. On the other hand, the
simultaneous estimation of signals from different leads using multivariate
models would reduce the already rare cases of the incorrect identification of
waves or anomalous parameter estimators.

## Acknowledgements

The authors gratefully acknowledge the financial support received by the
Spanish Ministerio de Ciencia e Innovación [PID2019-106363RB-I00] to CR, IF,
YL as well the Call for predoctoral contracts of the UVa 2020 to AR-C.

## Author contributions

CR: Conceived aims, theoretical proposal, conceptual design, wrote the
manuscript. CR,IF: revised specialized bibliograhy and propose LBBB criteria.
CR, YL design the preprocesing stage. YL, AR-C, CC: developed the
computational code and the app. All the authors approved the manuscript.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

##  Appendix A FMMecg Examples for BBB Features

This appendix shows examples that illustrate how FMMecg parameters describe
BBB features as those given in Table 1 of the main document. Row 1 in Figure 1
corresponds to a typical (non-pathological) pattern simulated using the median
FMMecg parameters obtained from the analysis of NORM patients in PTB-XL.

The rows 2-5 show patterns using alternative parameter configurations, which
are derived from that of row 1. The inner tables in each of the plots in
Figure 1 show the specific FMMecg parameter configuration that differs from
the typical pattern.

Specifically, $\omega_{R}$ is related to the width of the R wave, as the plots
in row 2 show. Besides, row 3 illustrates that a higher distance between
$\alpha_{Q}$ and $\alpha_{S},$ is related with a prolonged QRS duration, and
also that changes in $A_{Q}$ and $\beta_{Q}$ (resp. $A_{S}$ and $\beta_{S}$)
yield a qR (resp. Rs) notched pattern. Parameter configuration related to S
wave is illustrated in row 4, together with a specific configuration for ST
segment depression. Finally, T wave inversions are observed in row 5.

Figure 1: Simulated examples to illustrate BBB diagnostic features showed in
Table 1 of the main text

##  Appendix B Data preprocessing

ECG data preprocessing is required to remove high noise signals, or other
motion artifacts in the recording, to normalize the signal and to correct the
trend on ECG heartbeats. Conventionally, ECG data come from multiple leads and
their preprocessing is based on the following stages: ECG denoising, QRS
detection, ECG segmentation, scaling and trend correction [36, 37]. Most of
the preprocessing procedures proposed in the literature cover these stages. A
brief discussion of that is given below.

Most of noise suppression techniques are often based on band-pass filterings,
such as the Butterworth filter [38, 39]. They highly rely on a fixed cut-off
frequency and cannot track the changing characteristics of the time-varying
ECG signal. As a result, the ST segment or the QRS complex may be notably
distorted [40].

Regarding QRS detection, several methods have prolifered in the last decade
[41, 42]. However, Pan-Tompkins algorithm [27], reported in 1985, is
considered as the benchmark in QRS detection. Once these complexes are set,
ECG heartbeats must be delimeted, existing a wide range of strategies to cover
this goal. Among others, there are those based on P onset and T offset
locations [43], although these marks are not precisely defined in practice.
While in others [1, 44], the cutting is based on the length of the RR
intervals, i.e. the time elapsed between two successive QRS detections.

ECG segmentation often displays a high degree of uncertainty, due to the noise
in the device, missing sections of the ECG signals, or because the signal
contains a transient rise of signal amplitude [45]. In such cases, Pan
Tompkings usually fails resulting in missed peaks or false detections [46,
47]. Therefore, an ECG segmentation revision is desirable. Most of these
procedures are based on the RR interval properties such as length or density
[48, 44].

The ECG signal normalization is commonly based on linearly scaling methods
[49, 50]. While polynomial regression is widely extended for ECG detrending
[38, 36]. Finally, ECG fragment distortions, caused by sudden body movements,
or other noise artifacts, are usually removed if the heartbeat’s amplitude is
relatively high with regard to that observed in the signal [51, 52].

A relatively simple data preprocessing algorithm is defined here. First, noise
suppression is dismissed in this work since FMMecg is robust against noise
artifacts [1]. For QRS detection, Pan Tompkins algorithm is applied on raw ECG
data, from one or three leads. ECG heartbeats are delimited based on the RR
intervals. ECG signals are normalized using a linearly scaling and heartbeat
detrending is based on linear regression. The details of the preprocessing are
described below. An overview of this algorithm is illustrated in Figure 1.

Three-lead algorithm: The raw ECG signals from the leads I, II and V5 are the
inputs. The five steps of the algorithm are:

  1. 1.

QRS detection and ECG segmentation: Raw ECG signals from leads I, II and V5
are independently analyzed with Pan Tompkins. Then, ECG hearbeats are
delimited from the QRS annotations ($t^{QRS}$) and RR interval lengths as
follows: [$t^{QRS}-40\%RR,$ $t^{QRS}+60\%RR$], see [1] for details. Lead II is
used as reference [53, 54].

  2. 2.

ECG segmentation revision: Median RR length is used to identify and remove too
short/long heartbeats and those with distant R peaks or misplaced $t^{QRS}$
locations, following the ideas given in [55, 48, 1].

  3. 3.

Scaling and detrending: ECG recording is scaled into [-1,1] using min-max
normalization as in [50]. Next, as noted in [56, 46, 50], those heartbeats
with significant trends or remarked differences at the heartbeat’s boundaries,
are detrended by using linear regression [57, 58, 59].

  4. 4.

Remove ECG distortions: Heartbeats whose amplitudes are considerably larger
with regard to the QRS amplitudes in the signal are removed as done in [51,
52].

  5. 5.

QRS annotations checking: The QRS annotations from two leads are considered to
be in agreement if their distance is less than that corresponding to 4%
frequency (Hz). If no two leads coincide in at least three QRS annotations,
the patient is discarded. Otherwise, in the case the reference lead does not
match with any of the others in at least three QRS annotations, the reference
derivation is changed to that with lower RR length variation coefficient (vc)
among I and V5. In this latter case the algorithm goes back to Step2.
Otherwise, the preprocessing has finished.

For each lead, the outputs are ECG fragments with valid or invalid (removed)
heartbeats. Any patient with less than three valid heartbeats in a given lead
is discarded.

Removed heartbeats are imputed after the FMMecg analysis. The imputed values
are the median FMMecg parameter values obtained from the valid heartbeats in
the corresponding fragment.

Single-lead algorithm: The raw ECG signals one leads are the inputs. The
algorithm given above must be adapted as follows. Step 1 is conducted just
focusing on the input lead. Then, Steps 2-4 are conducted similarly. Step 5
reduced to discard patients with less than three valid heartbeats.

Figure 1: Outline of the three-lead preprocessing algorithm

##  Appendix C FMMecg Analysis app

###  C.1 Implementation

The application has been developed in the programming language R, using the
package shiny [60]. Shiny applications are mainly composed of two modules: the
User Interface (UI), designed with HTML widgets and CSS elements, and the
Server, responsible for the computational tasks implemented in R.

The UI module makes use of several R packages to create a complete and
flexible interface. In particular, the dashboard appearance is given by
shinydashboard [61] and flexdashboard [62] packages. The dashboard has been
implemented with the packages shinyWidgets [63], shinyjs [64], and ggplot2
[65], which provide a collection of control elements, JavaScript interactions,
and a visualization engine, respectively.

The data analysis using the FMM approach is implemented using the FMM package
[66, 67]. Detailed information on the model specifications and the estimation
algorithm are given in [1].

###  C.2 Application structure

An overview of the app workflow is shown in Figure 1. Flexibility is granted
to the app by its different input elements, which are the following:

  * 1.

Main input. An uploaded tabular data file, which can either have a single beat
or multiple recorded in one of the I, II, or V5 leads. It must contain a
single column with a column header name and be one of the following
extensions: .csv, .xls, or .xlsx.

  * 2.

Data options. User may indicate whether the data is composed by a single beat
or by multiple beats. The user must insert the observation number of the QRS
annotation, in the case of single beat analysis, and the sampling rate in Hz,
for multiple beat analysis.

  * 3.

Algorithm options: number of backfittings. More backfitting iterations imply a
better prediction at the cost of a higher computational time.

Once the previous options have been correctly specified, the FMMecg analysis
can be started by pressing the prediction button. The following outputs are
given:

  * 1.

The predicted FMM signal plotted along with the input data.

  * 2.

Accuracy measures: $var_{J};J=P,Q,R,S,T$ are measures accounting for the
percentage of the variability explained for each wave relative to the previous
fitted waves. Let J be the k-th wave fitted, then:

$var_{J}=R^{2}_{1,...,k}-R^{2}_{1,...,k-1},$

where $R^{2}_{1,...,k}$, is the proportion of variance explained by the FMM
model defined by the first k waves, out to the total variance.

  * 3.

Individual waves plot. Plot of estimated $W_{J}(t),\,J\in\\{P,Q,R,S,T\\}$,
$t\in(0,2\pi]$.

  * 4.

Median wave parameters. Median estimated values of
$A_{J},\alpha_{J},\beta_{J},\omega_{J}$, $J\in\\{P,Q,R,$ $S,T\\}$.

  * 5.

omeR and omeS estimated markers and diagnostic button. The diagnostic is made
using the omeR and omeS rules to detect CLBBB and CBBB pathologies,
respectively.

Predicted signal and individual waves plots are shown simultaneously when data
corresponds to a single beat. In the case of analysing multiple beat signal,
the app provides a button to shift between the two plots.

For demonstration purposes, example data is also provided. In particular, a
single healthy beat, a single beat with CLBBB and a multiple beat signal with
CBBB.

Figure 1: FMMecg Analysis app workflow chart

###  C.3 app examples

Three examples from PTB-XL database have been included in the app. The first
example (healthy) is the ECG signal recorded from patient 14, lead II, the
second example (CLBBB) is the ECG signal recorded from patient 796, lead I,
and the third example (CBBB) is a multiple beat signal recorded from patient
10131, lead V5. Upon initialization, the app analyzes the first example. The
app automatically provides the data options if example data is analyzed, which
for this case are the ‘single beat’ option and 68 for the QRS annotation.
Lastly, the user can choose the number of backfittings of the FMM algorithm,
which can be initiated by pressing the prediction button. The app shows the
predicted signal and individual waves plots, the high prediction accuracy, the
estimated parameter values for each wave, and the omeR and omeS markers. By
clicking the diagnostic button, the pathology rules are evaluated. In the
first example the patient is diagnosed as healthy.

Similarly, in the second example, the data inputs are analyzed with the
options ‘single beat’ and 70 for the QRS annotation. While for the third
example, the options ‘multiple beats’ and $200$ Hz for the sampling rate must
be provided. Also in these two cases, pathologies are diagnosed, in
particular, CLBBB in the second example and CBBB in the third.

Users can upload their own ECG fragment to be analyzed through the FMM
approach with the corresponding options/buttons on the left side of the
screen. After uploading the data, whether the ECG fragment contains a single
or multiple beats must be indicated and also the QRS annotation, in the former
case, or the sampling rate, in the latter case. Finally, the user should
choose the number of backfittings with the correponding button.

## References

  * Rueda et al. [2021] C. Rueda, Y. Larriba, A. Lamela,  The hidden waves in the ECG uncovered revealing a sound automated interpretation method,  Scientific Reports 11 (2021) 1–11.
  * McAnulty and Rahimtoola [1984] J. H. McAnulty, S. H. Rahimtoola,  Bundle branch block,  Progress in Cardiovascular Diseases 26 (1984) 333–354.
  * Tan et al. [2020] N. Y. Tan, C. M. Witt, J. K. Oh, Y.-M. Cha,  Left bundle branch block: current and future perspectives,  Circulation, arrhythmia and electrophysiology 13 (2020) e008239.
  * De Bacquer et al. [2000] D. De Bacquer, G. De Backer, M. Kornitzer,  Prevalences of ECG findings in large population based samples of men and women,  Heart 84 (2000) 625–633.
  * Breithardt and Breithardt [2000] G. Breithardt, O.-A. Breithardt,  Left bundle branch block, an old-new entity,  Journal of Cardiovascular Translational Research 5 (2000) 107–116.
  * Alventosa-Zaidin et al. [2019] M. Alventosa-Zaidin, G. Pera, C. R. Saumell, N. M. Miralles, M. Z. Sanchez, T. G. Garcia, L. G. Font, M. B. Camps, J. Francisco-Pascual, J. B. Terradellas,  Diagnosis of right bundle branch block: a concordance study,  BMC Family Practice 20 (2019) 58.
  * Willems et al. [1985] J. L. Willems, E. O. Robles de Medina, R. Bernard, P. Coumel, C. Fisch, D. Krikler, N. A. Mazur, F. L. Meijler, L. Mogensen, P. Moret, Z. Pisa, P. M. Rautaharju, B. Surawicz, Y. Watanabe, H. J. Wellens,  Criteria for intraventricular conduction disturbances and pre-excitation,  Journal of the American College of Cardiology 5 (1985) 1261–1275.
  * Surawicz et al. [2009] B. Surawicz, R. Childers, B. J. Deal, L. S. Gettes,  AHA/ACCF/HRS recommendations for the standardization and interpretation of the electrocardiogram: part III: intraventricular conduction disturbances a scientific statement from the American Heart Association Electrocardiography and Arrhythmias Committee, Council on Clinical Cardiology; the American College of Cardiology Foundation; and the Heart Rhythm Society endorsed by the International Society for Computerized Electrocardiology,  Journal of the American College of Cardiology 53 (2009) 976–981.
  * Strauss et al. [2011] D. G. Strauss, R. H. Selvester, G. S. Wagner,  Defining left bundle branch block in the era of cardiac resynchronization therapy,  The American Journal of Cardiology 107 (2011) 927—934.
  * Kligfield et al. [2018] P. Kligfield, F. Badilini, I. Denjoy, S. Babaeizadeh, E. Clark, J. De Bie, B. Devine, F. Extramiana, G. Generali, R. Gregg, E. Helfenbein, J. Kors, R. Leber, P. Macfarlane, P. Maison-Blanche, I. Rowlandson, R. Schmid, M. Vaglio, G. van Herpen, J. Xue, B. Young, C. L. Green,  Comparison of automated interval measurements by widely used algorithms in digital electrocardiographs,  American Heart Journal 200 (2018) 1–10.
  * Xia et al. [2016] X. Xia, B. Wieslander, D. G. Strauss, G. S. Wagner, W. Zareba, A. J. Moss, J.-P. Couderc,  Automatic QRS Selvester scoring system in patients with left bundle branch block,  EP Europace 18 (2016) 308–314.
  * Mason et al. [2016] J. W. Mason, F. Badilini, M. Vaglio, R. L. Lux, B. Aysin, T. E. Moon, B. Heinz, I. Strachan,  A fundamental relationship between intraventricular conduction and heart rate,  Journal of Electrocardiology 49 (2016) 362–370.
  * Goldberger et al. [2000] A. L. Goldberger, L. A. Amaral, L. Glass, J. M. Hausdorff, P. C. Ivanov, R. G. Mark, J. E. Mietus, G. B. Moody, C.-K. Peng, H. E. Stanley,  PhysioBank, PhysioToolkit, and PhysioNet: components of a new research resource for complex physiologic signals,  Circulation 101 (2000) e215–e220.
  * Lyon et al. [2018] A. Lyon, A. Mincholé, J. P. Martínez, P. Laguna, B. Rodriguez,  Computational techniques for ECG analysis and interpretation in light of their contribution to medical advances,  Journal of The Royal Society Interface 15 (2018).
  * Hannun et al. [2019] A. Y. Hannun, P. Rajpurkar, M. Haghpanahi, G. H. Tison, C. Bourn, M. P. Turakhia, A. Y. Ng,  Cardiologist-level arrhythmia detection and classification in ambulatory electrocardiograms using a deep neural network,  Nature Medicine 25 (2019) 65–69.
  * Ribeiro et al. [2020] A. H. Ribeiro, M. H. Ribeiro, G. M. Paixão, D. M. Oliveira, P. R. Gomes, J. A. Canazart, M. P. Ferreira, C. R. Andersson, P. W. Macfarlane, W. Meira Jr, T. B. Schön, A. L. Ribeiro,  Automatic diagnosis of the 12-lead ECG using a deep neural network,  Nature Communications 11 (2020).
  * Hao et al. [2018] C. Hao, S. Sivanesan, M. Majmudar, K. S. Rajput,  Combinational feature based random forest classification for enhanced bundle branch block beat detection,  in: 2018 IEEE EMBS International Conference on Biomedical Health Informatics (BHI), pp. 319–322.
  * Allami et al. [2016] R. Allami, A. Stranieri, V. Balasubramanian, H. F. Jelinek,  A genetic algorithm-neural network wrapper approach for bundle branch block detection,  in: 2016 Computing in Cardiology Conference (CinC), pp. 461–464.
  * Moody and Mark [2001] G. B. Moody, R. G. Mark,  The impact of the MIT-BIH arrhythmia database,  IEEE Engineering in Medicine and Biology Magazine 20 (2001) 45–50.
  * Martín-Yebra and Martínez [2019] A. Martín-Yebra, J. P. Martínez,  Automatic diagnosis of strict left bundle branch block using a wavelet-based approach,  PLoS One 14 (2019) e0212971.
  * Yang et al. [2020] T. Yang, R. E. Gregg, S. Babaeizadeh,  Detection of strict left bundle branch block by neural network and a method to test detection consistency,  Physiological Measurement 41 (2020) 025005.
  * Moss et al. [2009] A. J. Moss, W. J. Hall, D. S. Cannom, H. Klein, M. W. Brown, J. P. Daubert, N. M. Estes III, E. Foster, H. Greenberg, S. L. Higgins, M. A. Pfeffer, S. D. Solomon, D. Wilber, W. Zareba,  Cardiac-resynchronization therapy for the prevention of heart-failure events,  New England Journal of Medicine 361 (2009) 1329–1338.
  * Wagner et al. [2020] P. Wagner, N. Strodthoff, R.-D. Bousseljot, D. Kreiseler, F. I. Lunze, W. Samek, T. Schaeffter,  PTB-XL, a large publicly available electrocardiography dataset,  Scientific Data 7 (2020) 1–15.
  * Alday et al. [2021] E. A. P. Alday, A. Gu, A. J. Shah, C. Robichaux, A.-K. I. Wong, C. Liu, F. Liu, A. B. Rad, A. Elola, S. Seyedi, Q. Li, A. Sharma, G. D. Clifford, M. A. Reyna,  Classification of 12-lead ECGs: the PhysioNet/Computing in Cardiology Challenge 2020,  Physiological Measurement 41 (2021) 124003.
  * Liu et al. [2018] F. Liu, C. Liu, L. Zhao, X. Zhang, X. Wu, X. Xu, Y. Liu, C. Ma, S. Wei, Z. He, J. Li, Y. Kwee, E. Ng,  An open access database for evaluating the algorithms of electrocardiogram rhythm and morphology abnormality detection,  Journal of Medical Imaging and Health Informatics 8 (2018) 1368–1373.
  * Rueda et al. [2021] C. Rueda, A. Rodríguez-Collado, Y. Larriba,  A novel wave decomposition for oscillatory signals,  IEEE Transactions on Signal Processing 69 (2021) 960–972.
  * Pan and Tompkins [1985] J. Pan, W. J. Tompkins,  A real-time QRS detection algorithm,  IEEE Transactions on Biomedical Engineering BME-32 (1985) 230–236.
  * Dhingra et al. [2006] R. Dhingra, M. J. Pencina, T. J. Wang, B.-H. Nam, E. J. Benjamin, D. Levy, M. G. Larson, W. B. Kannel, R. B. D’Agostino Sr, R. S. Vasan,  Electrocardiographic QRS duration and the risk of congestive heart failure: the Framingham Heart Study,  Hypertension 47 (2006) 861–867.
  * Badheka et al. [2013] A. O. Badheka, V. Singh, N. J. Patel, A. Deshmukh, N. Shah, A. Chothani, K. Mehta, P. Grover, G. T. Savani, S. Gupta, A. Rathod, G. R. Marzouk, R. D. Mitrani, M. Moscucci, M. G. Cohen,  QRS duration on electrocardiography and cardiovascular mortality (from the National Health and nutrition examination survey—III),  The American Journal of Cardiology 112 (2013) 671–677.
  * Iuliano et al. [2002] S. Iuliano, S. G. Fisher, P. E. Karasik, R. D. Fletcher, S. N. Singh, Department of Veterans Affairs Survival Trial of Antiarrhythmic Therapy in Congestive Heart Failure,  QRS duration and mortality in patients with congestive heart failure,  American Heart Journal 143 (2002) 1085–1091.
  * Fosbøl et al. [2007] E. L. Fosbøl, M. Seibæk, B. Brendorp, D. V. Møller, M. Ersbøll, C. Torp-Pedersen, L. Køber, Danish Investigations and Arrhythmia ON Dofetilide (DIAMOND) study group,  Differential prognostic importance of QRS duration in heart failure and acute myocardial infarction associated with left ventricular dysfunction,  European Journal of Heart Failure 9 (2007) 814–819.
  * Teodorescu et al. [2011] C. Teodorescu, K. Reinier, A. Uy-Evanado, J. Navarro, R. Mariani, K. Gunson, J. Jui, S. S. Chugh,  Prolonged QRS duration on the resting ECG is associated with sudden death risk in coronary disease, independent of prolonged ventricular repolarization,  Heart Rhythm 8 (2011) 1562–1567.
  * Rueda et al. [2019] C. Rueda, Y. Larriba, S. D. Peddada,  Frequency modulated möbius model accurately predicts rhythmic signals in biological and physical sciences,  Scientific Reports 9 (2019) 1–10.
  * Senesael et al. [2020] E. Senesael, S. Calle, V. Kamoen, R. Stroobandt, M. De Buyzere, F. Timmermans, J. De Pooter,  Progression of incomplete toward complete left bundle branch block: A clinical and electrocardiographic analysis,  Annals of Noninvasive Electrocardiology 25 (2020) e12732.
  * Epstein et al. [2008] A. E. Epstein, J. P. DiMarco, K. A. Ellenbogen, N. M. Estes, R. A. Freedman, L. S. Gettes, A. M. Gillinov, G. Gregoratos, S. C. Hammill, D. L. Hayes, M. A. Hlatky, L. K. Newby, R. L. Page, M. H. Schoenfeld, M. J. Silka, L. W. Stevenson, M. O. Sweeney, S. C. Smith Jr, A. K. Jacobs, C. D. Adams, J. L. Anderson, C. E. Buller, M. A. Creager, S. M. Ettinger, D. P. Faxon, J. L. Halperin, L. F. Hiratzka, S. A. Hunt, H. M. Krumholz, F. G. Kushner, B. W. Lytle, R. A. Nishimura, J. P. Ornato, R. L. Page, B. Riegel, L. G. Tarkington, C. W. Yancy, American College of Cardiology/American Heart Association Task Force on Practice Guidelines,  ACC/AHA/HRS 2008 guidelines for device-based therapy of cardiac rhythm abnormalities: a report of the American College of Cardiology/American Heart Association Task Force on Practice Guidelines (Writing Committee to Revise the ACC/AHA/NASPE 2002 Guideline Update for Implantation of Cardiac Pacemakers and Antiarrhythmia Devices) developed in collaboration with the American Association for Thoracic Surgery and Society of Thoracic Surgeons,  Circulation 117 (2008) 861–867.
  * Aspuru et al. [2019] J. Aspuru, A. Ochoa-Brust, R. A. Félix, W. Mata-López, L. J. Mena, R. Ostos, R. Martínez-Peláez,  Segmentation of the ECG signal by means of a linear regression algorithm,  Sensors 19 (2019) 775.
  * Tseng et al. [2016] K.-K. Tseng, D. Lee, C. Chen,  ECG identification system using neural network with global and local features,  in: International Association for Development of the Information Society, ERIC, 2016.
  * Friganovic et al. [2018] K. Friganovic, D. Kukolja, A. Jovic, M. Cifrek, G. Krstacic,  Optimizing the detection of characteristic waves in ECG based on processing methods combinations,  IEEE Access 6 (2018) 50609–50626.
  * Pandey and Tiwari [2014] H. Pandey, R. Tiwari,  An innovative design approach of Butterworth filter for noise reduction in ECG signal processing based applications,  Progress in Science in Engineering Research Journal (2014).
  * Sun et al. [2002] Y. Sun, K. L. Chan, S. M. Krishnan,  ECG signal conditioning by morphological filtering,  Computers in Biology and Medicine 32 (2002) 465–479.
  * Chen and Chuang [2017] C.-L. Chen, C.-T. Chuang,  A QRS detection and R point recognition method for wearable single-lead ECG devices,  Sensors 17 (2017) 1969.
  * Elgendi et al. [2014] M. Elgendi, B. Eskofier, S. Dokos, D. Abbott,  Revisiting QRS detection methodologies for portable, wearable, battery-operated, and wireless ECG systems,  PloS one 9 (2014) e84018.
  * Sahoo et al. [2018] P. K. Sahoo, H. K. Thakkar, W.-Y. Lin, P.-C. Chang, M.-Y. Lee,  On the design of an efficient cardiac health monitoring system through combined analysis of ECG and SCG signals,  Sensors 18 (2018) 379.
  * Tseng et al. [2016] K.-K. Tseng, F. Zeng, W. Ip, C.-H. Wu,  ECG sensor verification system with mean-interval algorithm for handling sport issue,  Journal of Sensors 2016 (2016) 1–12.
  * Nurmaini et al. [2021] S. Nurmaini, A. E. Tondas, A. Darmawahyuni, M. N. Rachmatullah, J. Effendi, F. Firdaus, B. Tutuko,  Electrocardiogram signal classification for automated delineation using bidirectional long short-term memory,  Informatics in Medicine Unlocked 22 (2021) 100507.
  * Arefin et al. [2015] M. R. Arefin, K. Tavakolian, R. Fazel-Rezai,  QRS complex detection in ECG signal for wearable devices,  in: 2015 37th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), pp. 5940–5943.
  * Mondelo et al. [2017] V. Mondelo, M. J. Lado, A. J. Méndez, X. A. Vila, L. Rodríguez-Liñares,  Combining 12-lead ECG information for a beat detection algorithm,  Theoretical and Applied Informatics 3 (2017) 5–9.
  * Sadhukhan and Mitra [2012] D. Sadhukhan, M. Mitra,  R-peak detection algorithm for ECG using double difference and RR interval processing,  Procedia Technology 4 (2012) 873–877.
  * Yang et al. [2013] M. Yang, B. Liu, M. Zhao, F. Li, G. Wang, F. Zhou,  Normalizing electrocardiograms of both healthy persons and cardiovascular disease patients for biometric authentication,  PloS one 8 (2013) e71523.
  * Liang et al. [2020] Y. Liang, S. Yin, Q. Tang, Z. Zheng, M. Elgendi, Z. Chen,  Deep learning algorithm classifies heartbeat events based on electrocardiogram signals,  Frontiers in Physiology 11 (2020) 569050.
  * Chu and Delp [1990] C. H. Chu, E. J. Delp,  Nonlinear methods in electrocardiogram signal processing,  Journal of Electrocardiology 23 (1990) 192–197.
  * Sörnmo and Laguna [2006] L. Sörnmo, P. Laguna, Electrocardiogram (ECG) signal processing, John Wiley & Sons, Inc. Hoboken, NJ, USA.
  * Meek and Morris [2002] S. Meek, F. Morris,  Introduction. I—Leads, rate, rhythm, and cardiac axis,  BMJ 324 (2002) 415–418.
  * Mondéjar-Guerra et al. [2019] V. Mondéjar-Guerra, J. Novo, J. Rouco, M. G. Penedo, M. Ortega,  Heartbeat classification fusing temporal and morphological information of ECGs via ensemble of classifiers,  Biomedical Signal Processing and Control 47 (2019) 41–48.
  * Rezk et al. [2011] S. Rezk, C. Join, S. El Asmi,  An algebraic derivative-based method for R wave detection,  in: 2011 19th European Signal Processing Conference, IEEE, Barcelone, Spain.
  * Sun et al. [2005] Y. Sun, K. L. Chan, S. M. Krishnan,  Characteristic wave detection in ECG signal using morphological transform,  BMC Cardiovascular Disorders 5 (2005) 28.
  * Mitov [1998] I. Mitov,  A method for assessment and processing of biomedical signals containing trend and periodic components,  Medical Engineering & Physics 20 (1998) 660–668.
  * Sharma and Sharma [2017] T. Sharma, K. K. Sharma,  QRS complex detection in ECG signals using locally adaptive weighted total variation denoising,  Computers in Biology and Medicine 87 (2017) 187–199.
  * Kontaxis et al. [2021] S. Kontaxis, E. Gil, V. Marozas, J. Lázaro, E. García, M. Posadas-de Miguel, S. Siddi, M. L. Bernal, J. Aguiló, J. M. Haro, C. de la Cámara, P. Laguna, R. Bailón,  Photoplethysmographic waveform analysis for autonomic reactivity assessment in depression,  IEEE Transactions on Biomedical Engineering 68 (2021) 1273–1281.
  * Chang et al. [2021] W. Chang, J. Cheng, J. Allaire, C. Sievert, B. Schloerke, Y. Xie, J. Allen, J. McPherson, A. Dipert, B. Borges, shiny: Web Application Framework for R, 2021\. R package version 1.6.0.
  * Chang and Borges Ribeiro [2018] W. Chang, B. Borges Ribeiro, shinydashboard: Create Dashboards with ’Shiny’, 2018\. R package version 0.7.1.
  * Iannone et al. [2020] R. Iannone, J. Allaire, B. Borges, flexdashboard: R Markdown Format for Flexible Dashboards, 2020. R package version 0.5.2.
  * Perrier et al. [2021] V. Perrier, F. Meyer, D. Granjon, shinyWidgets: Custom Inputs Widgets for Shiny, 2021. R package version 0.6.0.
  * Attali [2020] D. Attali, shinyjs: Easily Improve the User Experience of Your Shiny Apps in Seconds, 2020. R package version 2.0.0.
  * Wickham [2016] H. Wickham, ggplot2: Elegant Graphics for Data Analysis, Springer-Verlag New York, 2016\.
  * Fernández et al. [2021a] I. Fernández, A. Rodríguez-Collado, Y. Larriba, A. Lamela, C. Canedo, C. Rueda, FMM: An R package for modeling rhythmic patterns in oscillatory systems, 2021a. URL: https://arxiv.org/abs/2105.10168.
  * Fernández et al. [2021b] I. Fernández, A. Rodríguez-Collado, Y. Larriba, A. Lamela, C. Canedo, C. Rueda, FMM: Rhythmic Patterns Modeling by FMM Models, 2021b. R package version 0.3.0.


# Towards a Standardized Reinforcement Learning Framework for AAM Contingency
Management

Luis E. Alvarez  Surveillance Systems
MIT Lincoln Laboratory
Lexington, MA USA
<EMAIL_ADDRESS>Marc Brittain  Surveillance Systems
MIT Lincoln Laboratory
Lexington, MA USA
<EMAIL_ADDRESS>Kara Breeden  Surveillance Systems
MIT Lincoln Laboratory
Lexington, MA USA
<EMAIL_ADDRESS>
###### Abstract

Advanced Air Mobility (AAM) is the next generation of air transportation that
includes new entrants such as electric vertical takeoff and landing (eVTOL)
aircraft, increasingly autonomous flight operations, and small UAS package
delivery. With these new vehicles and operational concepts comes a desire to
increase densities far beyond what occurs today in and around urban areas, to
utilize new battery technology, and to move toward more autonomously-piloted
aircraft. To achieve these goals, it becomes essential to introduce new safety
management system capabilities that can rapidly assess risk as it evolves
across a span of complex hazards and, if necessary, mitigate risk by executing
appropriate contingencies via supervised or automated decision-making during
flights. Recently, reinforcement learning has shown promise for real-time
decision making across a wide variety of applications including contingency
management. In this work, we formulate the contingency management problem as a
Markov Decision Process (MDP) and integrate the contingency management MDP
into the AAM-Gym simulation framework. This enables rapid prototyping of
reinforcement learning algorithms and evaluation of existing systems, thus
providing a community benchmark for future algorithm development. We report
baseline statistical information for the environment and provide example
performance metrics.

††DISTRIBUTION STATEMENT A. Approved for public release. Distribution is
unlimited. This material is based upon work supported by the National
Aeronautics and Space Administration under Air Force Contract No.
FA8702-15-D-0001. Any opinions, findings, conclusions or recommendations
expressed in this material are those of the author(s) and do not necessarily
reflect the views of the National Aeronautics and Space Administration. © 2023
Massachusetts Institute of Technology. Delivered to the U.S. Government with
Unlimited Rights, as defined in DFARS Part 252.227-7013 or 7014 (Feb 2014).
Notwithstanding any copyright notice, U.S. Government rights in this work are
defined by DFARS 252.227-7013 or DFARS 252.227-7014 as detailed above. Use of
this work other than as specifically authorized by the U.S. Government may
violate any copyrights that exist in this work.

##  I Introduction

Advanced Air Mobility (AAM) is a rapidly evolving field of aviation that has
the potential to revolutionize air transportation. AAM encompasses
transportation of people and cargo (e.g., passenger shuttling and package
delivery) between local, regional, and urban locations with a diverse set of
vehicle types which includes small uncrewed aerial systems (sUAS) and larger
occupied vehicles such as autonomous conventional takeoff and landing (CTOL)
and vertical takeoff and landing (VTOL) aircraft. Due in part to the
envisioned high density of operations, AAM requires increasingly autonomous
operations to enable operational efficiency and reduce risk in the airspace.

Phases of AAM operations are described by the FAA and NASA as evolving in UAM
Maturity Levels (UML) that range from 1 to 6 with increasing levels of
expected densities and complexities of roles and responsibilities [1, 2, 3].
For example, UML 1 envisions a low-density human piloted airspace environment,
whereas UML 4+ envisions 100s of simultaneous operations over a local region
with a mix of human and autonomously piloted aircraft. UML 4+ expects high
density, dynamic, and autonomous operations enabled by artificial intelligence
(AI) based decision-making algorithms. Therefore, it is important to
understand the trade-offs of various AI algorithms by comparing their
performance across a range of safety-critical situations.

Contingency management (CM) is a capability within a larger umbrella referred
to by NASA and the National Academies as In-Time Aviation Safety Management
Systems (IASMS) [4, 5]. The scope of CM includes (a) planning a set
contingencies prior and/or during flight, (b) ability to decide which
contingency to execute based on available information (e.g., aircraft or
system states, hazard level or likelihood) and (c) execution of contingency.
In particular, IASMS CM capability has focused on path planning techniques
that have an online system monitor various hazards, namely collisions
(airborne or ground), loss of power, loss of control, and loss of flight
capability [4].

To reduce the probability of entering an unsafe state, the majority of these
online systems use path planning heuristic algorithms. Predominant algorithms
explored are Monte Carlo approaches of path planning where the algorithm finds
an optimal path within a time horizon or explores a set of alternate
trajectories given the aircraft’s dynamics constraints [6, 7, 8]. Example path
planning approaches utilize variations of Rapidly-Exploring Random Tree (RRT)
and A* to generate an optimized path based on distance or other environmental
metrics (e.g., surrounding air traffic). Parameter training for these
approaches leverage techniques from machine learning such as bagging methods
(e.g., bootstrap sampling) with uncertainties represented by random
distributions or pre-trained Bayesian belief networks [9, 10, 11]. A drawback
of these techniques is the requirement for subject matter experts to create
and bound the heuristics used by the model.

Reinforcement Learning (RL), a subset of machine learning, alleviates the
drawbacks of the above approaches by learning a policy (akin to heuristics)
through repeated interaction with an environment to achieve a specified
objective. Recently, RL has shown promise in many challenging games such as
Go, Atari, Warcraft, and, most recently, Starcraft II with beyond human-level
performance [12, 13, 14]. In addition, RL has been applied to many real-world
air transportation problems such as conflict detection and resolution [15, 16,
17], autonomous separation assurance [18, 19, 20], and collision avoidance
[21]. Furthermore, the Airborne Collision Avoidance System X (ACAS X) family
of systems proved RL agents can improve safety and reduce nuisance alerts in
real-world operation compared to the heuristic-based systems, in particular
the Traffic Alert and Collision Avoidance System (TCAS) that is currently
mandated on all commercial aircraft [22, 23, 24, 25]. These notable
advancements demonstrate the capability of computational learning algorithms
to potentially augment and facilitate human tasks in real-world environments.

In this work, a framework is proposed to enable the development and evaluation
of RL-based contingency management designs, particularly those that constrain
contingency actions to a limited and deterministic set of highly assured
maneuvers. This is a principle of the run-time assurance concept used for
safety-critical systems. A set of metrics are introduced to evaluate the
nominal (unequipped) performance of the framework, providing a baseline for RL
research.

##  II Background

RL is one type of sequential decision making where the objective is to learn a
policy in a given environment. RL requires the environment to be formulated as
a Markov Decision Process (MDP), a mathematical framework for modeling
decision making processes with stochastic transitions. An MDP can be defined
by the five-tuple $(S,A,R,T,\gamma)$, where $S$ is the state space, $A$ is the
action set, $R$ is the reward function $R(s,a)$ providing the reward for
taking action $a$ from state $s$, and $T$ is the probability transition
function, $T(s,a,s^{\prime})$, that produces the probability of transition
from state $s$ to the next state $s^{\prime}$ when taking action $a$. In RL,
the transition matrix $T$ is often unknown. The discount factor $\gamma$
determines how far in the future to look for rewards. As $\gamma\rightarrow$
0, immediate rewards are emphasized, whereas, when $\gamma\rightarrow$ 1,
future rewards are prioritized.

From this formulation, the RL agent is able to derive an optimal policy in the
environment by maximizing a cumulative reward function. Let $\pi$ represent
some policy and $\tau$ represent the total time for a given environment, then
the optimal policy can be defined as follows:

$\pi^{*}=\operatorname*{arg\,max}_{\pi}E[\sum_{t=0}^{\tau}(r(s_{t},a_{t})|\pi)].$
(1)

By designing the reward function to reflect the objective in the environment,
the optimal solution can be obtained by maximizing the total reward.

In environments with a low dimensional state-action space, or environments
with discrete state-actions values, the policy $\pi$ can be represented as a
look-up table and solved using dynamic programming approaches such as
Q-learning [26]. However, in many real-world environments the state-action
space can not be represented by discrete values, or the dimensionality of the
state-action space is too large to construct a look-up table. Deep
reinforcement learning (DRL) alleviates the aforementioned issues by
introducing a neural network to represent the policy $\pi$. Various approaches
exist for training DRL agents including value-based methods [13, 27, 28], and
policy-based optimization [29, 30, 31].

Large scale simulation capabilities are essential for RL algorithms as they
explore the transitions between states as a consequence of their actions
within the modeled environment. Recently, MIT Lincoln Laboratory introduced
the AI testbed for advanced air mobility, or AAM-Gym to provide a
representative fast-time simulation environment for research and development
of decision making tools for AAM concepts of operation [32]. By leveraging
AAM-Gym, new RL use-cases, such as contingency management, can be developed to
explore the trade-offs of AI/RL-based decision-making functions with existing
approaches. The proposed RL contingency management framework utilizes AAM-Gym
[32] with BlueSky [33] as the back-end simulation capability to explore the
impact of RL algorithm choice, hazards (e.g., wind fields, rate of loss of
navigation, battery discharge prediction failures), and aircraft parameters on
the overall training and effectiveness of agents. By utilizing AAM-Gym, the
state space of each agent can make use of all the aircraft and environment
information if desired, as seen in Tables I and II.

In order to provide representative real world scenarios to AAM-Gym, we
leverage the Lincoln Laboratory UAMToolKit framework; an event-driven
simulation environment that can generate spatio-temporal traffic operations
with a diverse set of flight routes [34]. When modeling Urban Air Mobility
(UAM) operations, time variant cycles of demand across the New York City
metropolitan area are provided by estimating the demand as a percentage of the
expected replacement by AAM of for-hire-vehicle (FHV) operations (e.g. Uber,
Lyft, taxi, limousine services). The event driven model outputs a complete
history of aircraft movements and metadata, which includes departure times,
cruising speeds, arrival times, flight route (i.e., flight level and
waypoints), origin-destination pair, and aircraft type. UAMToolKit can be
queried by AAM-Gym allowing for vast amounts of unique and diverse traffic
patterns for the RL agents to train and evaluate against.

TABLE I: Aircraft Information Available for Agent Training. Parameter | Range |
---|---|---
Heading | 0 - 360° |
Altitude | 0 - 5,000 ft |
Aircraft horizontal speed | 0 - 120 knots |
Horizontal acceleration | 0.1 - 0.5$g$ |
Flight path distance to destination | - |
Energy level | 0 - 350 KWh |
Flight route waypoints | 0 - complete route |
TABLE II: Environment Information Available for Agent Training. Parameter | Notes
---|---
Wind | North and east magnitudes
Vertiport locations | Latitude and longitude
Approximate population density | From NYC FHV database
Probability of loss of navigation ($P_{nav}$) | Guassian noise per aircraft
Aircraft in vicinity | All fields available per aircraft

##  III Approach

In this section, the expansions of AAM-Gym for the RL contingency management
use-case are described as well as the MDP formulation.

###  III-A Traffic Generation

When modeling UAM operations, UAMToolKit requires information about the fleet
size, aircraft capacity of each vertiport, vehicle passenger capacity,
cruising speeds, the corridor network for operation, distributions of turn-
around time, and wind data. The focus of this research is on contingency
management during flight and not scheduling impacts due to turn-around and
wind data. Therefore, operations are modeled with zero wind field and a one
minute of turn-around time, the minimum simulation step. Operations occur
between eight altitude lanes evenly spaced from 1,000–5,000 ft utilizing the
existing helicopter network as depicted in [34]. Passenger capacity is
maintained at four passengers, and 50–500 aircraft are attempting to meet 5%
of the FHV demand, with landing pad spaces in each vertiport maintained so all
aircraft have a landing pad at the beginning of the simulation and an
additional landing pad for emergency. For example, in the case of the 500
fleet size scenario, each of the 29 vertiports selected in [34] would have 18
landing pads.

###  III-B Algorithms

AAM-Gym allows researchers to develop, train, and evaluate RL algorithms that
adhere to the OpenAI Gym [35] protocol, allowing for rapid prototyping of new
and existing algorithms. AAM-Gym has support for both general RL based
frameworks such as Ray’s RL-Lib [36] and air transportation based RL
algorithms such as the Deep Distributed Multi-Agent Variable with Attention
(D2MAV-A) [15].

###  III-C Hazard Modeling

The hazards initially modeled by the RL framework include: uncertainty in
battery consumption, max charge capacity due to number of charge cycles, wind
information, and navigation capability. Other hazards will be modeled in the
future.

As part of the RL contingency management use-case, two energy consumption rate
models were developed. The energy consumption rate is modeled within the RL
framework as either a linear consumption model or a mixture of a BADA 3 and
Bada H [37, 38] based energy consumption model during cruising phases with a
linear model at speeds below 20 knots.

The linear energy consumption model is represented by (2) where $\alpha$ is
the constant energy consumption rate per time step, $\Delta{t}$ is the time
step between simulation iteration, $C$ is the total number of charge cycles an
aircraft has experienced, and $\beta$ is the number of charge cycles when
uncertainty in consumption begins to occur. The uncertainty in consumption
rate after $\beta$ charge cycles is modeled as a truncated Gaussian
distribution, $N$, of range $[0,1]$ for values above $\phi$.

$\Delta{E}=\begin{cases}\alpha\cdot\Delta{t}&\text{if $C\leq\beta$}\\\
\alpha\cdot\Delta{t}+\text{$N$}&\text{if $C>\beta$}\\\ \end{cases}.$ (2)

When each aircraft is created, it is randomly assign a value for the number of
charge cycles, $C$, by a uniform distribution of a given range. The model then
uses the defined minimum energy level ($E_{min})$, maximum energy level
($E_{max}$), minimum charge cycles ($C_{min}$), and maximum charge cycles
($C_{max}$) to create a linear fit and assign the maximum energy level
possible given the number of cycles the aircraft has experienced (example
shown in Fig. 1).

Figure 1: Example of randomly sampled energy assigned to aircraft at creation
(red circle) from a linear function of charge cycles and energy level when
$E_{max}=250$ KWh, $C_{min}=0$ cycles, $E_{min}=100$ KWh and $C_{max}=10,000$
cycles.

The wind fields can be provided as constant wind fields, input data from
Weather Research Forecast (WRF) model [39], or as multivariate distributions.
When providing discrete locations for wind field information, the BlueSky
simulation environment uses multi-linear interpolation to provide the wind
vector at specified locations.

The probability of navigation loss, $P_{nav}$, can also be modeled as a
multivariate distribution providing the probability of losing navigation at
specified latitude and longitudes or as a discrete event occurring with a set
probability.

###  III-D Reward Model and Terminal States

As previously presented, the RL framework uses an MDP to construct the
contingency management problem where the researchers define the state space,
actions, transitions, and rewards for each agent. In this framework the
transition model is unknown as its dictated by the simulation environment and
the hazard parameter chosen. Therefore, it is not explicitly defined by the
researchers. The action space of the MDP evaluated is shown in Table III. The
available actions are defined as follows:

  1. 1.

Heading change - The aircraft turns left or right by a pre-defined magnitude
(e.g., $\pm$5$\degree$, 0$\degree$) and continues on this heading until a new
action is chosen.

  2. 2.

Land now - The aircraft chooses to execute a controlled vertical decent from
present position.

  3. 3.

No alert - The aircraft takes no action and continues on its current heading
or active flight plan.

  4. 4.

Use assigned flight route - The aircraft continues on its planned flight route
or returns to the nearest waypoint, if prior actions caused a deviation, and
continues along flight route.

TABLE III: Actions available for agents Description | Value
---|---
Heading change | -5°,0,-5°
Land now | -
No alert | -
Use assigned flight route | -

The reward function utilized by the MDP is defined by

$R(s_{t},h_{t},a_{t})=R(s_{t})+R(h_{t})+R(a_{t})-\Omega,$ (3)

where $R(s_{t})$, $R(h_{t})$, and $R(a_{t})$ are defined by (4), (5), (7) and
$\Omega$ is a step penalty. $R(s_{t})$ is a state dependant reward function
based on the terminal states of the aircraft as defined by (4). The RL
framework allows for three types of terminal states to be defined as total
loss of energy, loss of navigation, and when an aircraft reaches an altitude
of zero feet, each respectively adding a penalty $\delta_{energy}$,
$\delta_{navigation}$, and $\delta_{range\\_to\\_destination}$.

$R(s_{t})=\delta_{energy}+\delta_{navigation}+\delta_{range\\_to\\_destination}$
(4)

A multivariate Gaussian distribution is used to define the rewards of landing
at destination and alternative vertiports. It is defined as

$R(h_{t})=\sum_{i=0}^{n}\delta_{vertiport}\cdot
e^{\frac{(X-\mu_{xi})^{2}+(Y-\mu_{yi})^{2}}{-2\sigma^{2}}}$ (5)

where $X$ and $Y$ are the latitude and longitude of the aircraft under
control. The latitude and longitude of each vertiport are defined by $\mu_{x}$
and $\mu_{y}$ respectively. The standard deviation $\sigma$ is set to 0.0005°.
The variable $\delta_{i}$ is defined by

$\delta_{vertiport}=\\\ \begin{cases}\delta_{vertiport\\_destination}&\text{if
vertiport = destination}\\\ \delta_{vertiport\\_other}&\text{if vertiport
$\neq$ destination}\\\ \end{cases}.$ (6)

Last, $R(a_{t})$ is defined by

$R(a_{t})=\\\ \begin{cases}\delta_{land}+\delta_{action\\_penalty}&\text{if
$a_{t}=$ ``Land now''}\\\ \delta_{action\\_penalty}&\text{if $a_{t}\neq$ ``No
alert''}\end{cases},$ (7)

where $\delta_{action\\_penalty}$ is a penalty applied when any action is
chosen other than “No alert” and $\delta_{land}$ is the penalty for choosing
to “Land now”.

###  III-E Metrics

The RL framework extends the capabilities of AAM-Gym and allows users to
define performance metrics to help assess the CM algorithm within a Hydra
configuration file [40] leveraging the aircraft and environment information,
as defined in Table I and II. In addition, the following information was
recorded: spatial distance from corridor, number of aircraft that reached
their destination, energy level of each aircraft, total reward experienced by
each aircraft, average wind in scenario, time spent in flight route, and
actions taken.

##  IV Numerical Experiment and Results

Preliminary experiments conducted to date are described in this section with
the parameters explored and their resulting impact on the operations.
Additional experiments are planned.

###  IV-A Parameter Sweep

The Hydra configuration interface of AAM-Gym allows for multiple experiments
to be launched simultaneously to explore the impact of environment decisions
and rewards. This experiment explores the impact of modeling an increased
consumption rate on the battery model, differences in maximum charge capacity,
and activating probability of loss of navigation, $P_{nav}$, per aircraft per
simulation step. The parameters that were explored and the range of values
chosen are shown in Table IV. For this paper we explore the impact of each of
these parameters on the likelihood that an operation will reach its
destination or other assigned vertiports as a function of training steps.

TABLE IV: Environment Parameters Explored and Values Parameter | Values Explored
---|---
Maximum initial energy level ($E_{max}$) | 50, 150, 250 KWh
Minimum initial energy level ($E_{min}$) | 100 KWh
Maximum battery cycles ($C_{max}$) | $10\,000$
Minimum battery cycles ($C_{min}$) | 0
Energy consumption uncertainty threshold ($\beta$) | $3000$
Gaussian truncation parameter ($\phi$) | 0.5
Probability of loss of navigation capability ($P_{nav}$) | 0.0, $1\text{\times}{10}^{-5}$

###  IV-B Results

When no RL-based CM agent is in control of the aircraft (i.e., unequipped) the
likelihood of reaching their destination is strongly correlated to the maximum
initial energy level, $E_{max}$. This trend is evident by the higher fraction
of aircraft reaching the destination, $\bar{P}_{dest}$, as seen in Fig. 2. Due
to all AAM-Gym parallel workers initializing their respective simulations at
the same time with a low number of aircraft, we see a spike in fraction of
aircraft reaching their destination in Fig. 2. As the AAM-Gym parallel workers
continue the simulation, the number of aircraft stabilizes, thus creating a
steady number of aircraft reaching their destination as a function of
$E_{max}$ and $P_{nav}$.

This is expected as operations with lower energy levels are at risk of quickly
depleting their energy reserves during flight without reaching any vertiport
if an agent does not dictate an action to undertake. The unequipped results
also show the number of aircraft reaching their final destination decreases
when the probability of navigation loss, $P_{nav}$, is above zero, as seen in
Table V. A change in $P_{nav}$ from 0 to $1\text{\times}{10}^{-5}$ results in
an decrease of aircraft arrivals by 0.12% to 0.25% as it increases the
likelihood of entering a terminal state.

Figure 2: Rolling mean of the fraction of aircraft reaching their destination as a function of simulation steps, probability of loss of navigation ($P_{nav}$), and maximum initial energy level ($E_{max}$). TABLE V: Simulation Parameters and Metrics of Interest $P_{nav}$ | $E_{max}~{}\text{(KWh)}$ | $\max\bar{P}_{dest}$
---|---|---
$1\text{\times}{10}^{-5}$ | 50.0 | 0.4639
0.0 | 50.0 | 0.4654
$1\text{\times}{10}^{-5}$ | 150.0 | 0.8719
0.0 | 150.0 | 0.8745
$1\text{\times}{10}^{-5}$ | 250.0 | 0.9414
0.0 | 250.0 | 0.9439

##  V Conclusion

In this work, a RL framework was developed to enable standardized research and
development in the area of contingency management for AAM. It is shown how the
problem is formulated as a Markov Decision Process and through using AAM-Gym
and UAMToolkit, a representative simulation environment is established for
training and validation of RL agents. Preliminary results are presented to
demonstrate the impact of various simulation parameters that induce a
contingency event.

In future work, further analysis of the reward model, MDP environment
structure, and algorithm design is required to assess the viability of RL
algorithms in contingency management. While baseline results were presented in
this study, future work will explore various state-of-the-art RL approaches to
understand the viability of learning-based algorithms. After a broad analysis
of the performance of these algorithms, a selected set will be tuned through
the use of surrogate optimization. In addition, the simulation environment
will be further refined with the introduction of crowd source population
databases for causality penalties, and use of eVTOL vehicle performance models
with improved battery models.

##  VI Acknowledgment

The authors would like to thank Steve Young of NASA for his guidance and
support of this effort. The work was funded by NASA’s System-Wide Safety
Project via an Interagency Agreement with the US Air Force (IA 80LARC23TA002,
Project 10384).

## References

  * [1] Federal Aviation Administration, “Urban Air Mobility Concept of Operations version 1.0,” 2020. [Online]. Available: https://nari.arc.nasa.gov/sites/default/files/attachments/UAM_ConOps_v1.0.pdf
  * [2] National Aeronautics and Space Administration, “National Aeronautics and Space Administration (NASA) UAM Vision Concept of Operations (ConOps) UAM Maturity Level (UML)4,” 2020. [Online]. Available: https://ntrs.nasa.gov/api/citations/20205011091/downloads/UAMVisionConceptofOperationsUML-4v1.0.pdf
  * [3] Federal Aviation Administration, “Urban Air Mobility Concept of Operations version 2.0,” 2022. [Online]. Available: https://www.faa.gov/sites/faa.gov/files/UrbanAirMobilityUAMConceptofOperations2.0_0.pdf
  * [4] K. K. Ellis, P. Krois, J. Koelling, L. J. Prinzel, M. Davies, and R. Mah, “A Concept of Operations (ConOps) of an in-time aviation safety management system (IASMS) for Advanced Air Mobility (AAM),” in _AIAA Scitech 2021 Forum_ , 2021, p. 1978.
  * [5] National Academies of Sciences, Engineering, and Medicine, _Advancing Aerial Mobility: A National Blueprint_. Washington, DC: The National Academies Press, 2020. [Online]. Available: https://nap.nationalacademies.org/catalog/25646/advancing-aerial-mobility-a-national-blueprint
  * [6] J. Kim, P. Sharma, E. Atkins, N. Neogi, E. Dill, and S. Young, “Assured contingency landing management for advanced air mobility,” in _2021 IEEE/AIAA 40th Digital Avionics Systems Conference (DASC)_ , 2021, pp. 1–12.
  * [7] H. Haghighi, D. Delahaye, and D. Asadi, “Performance-based emergency landing trajectory planning applying meta-heuristic and dubins paths,” _Applied Soft Computing_ , vol. 117, p. 108453, 2022. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1568494622000242
  * [8] P. Váňa, J. Sláma, J. Faigl, and P. Pačes, “Any-time trajectory planning for safe emergency landing,” in _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2018, pp. 5691–5696.
  * [9] E. Ancel, S. D. Young, C. C. Quach, R. F. Haq, K. Darafsheh, K. M. Smalling, S. L. Vazquez, E. T. Dill, R. C. Condotta, B. E. Ethridge, L. R. Teska, and T. A. Johnson, _Design and Testing of an Approach to Automated In-Flight Safety Risk Management for sUAS Operations_. [Online]. Available: https://arc.aiaa.org/doi/abs/10.2514/6.2022-3459
  * [10] S. Choudhury, S. Scherer, and S. Singh, “Rrt*-ar: Sampling-based alternate routes planning with applications to autonomous emergency landing of a helicopter,” in _2013 IEEE International Conference on Robotics and Automation_ , 2013, pp. 3947–3952.
  * [11] M. Warren, L. Mejias, J. Kok, X. Yang, F. Gonzalez, and B. Upcroft, “An automated emergency landing system for fixed-wing aircraft: Planning and control,” _Journal of Field Robotics_ , vol. 32, no. 8, pp. 1114–1140, 2015. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1002/rob.21641
  * [12] D. Silver, A. Huang, C. J. Maddison, A. Guez, L. Sifre, G. Van Den Driessche, J. Schrittwieser, I. Antonoglou, V. Panneershelvam, M. Lanctot _et al._ , “Mastering the game of Go with deep neural networks and tree search,” _nature_ , vol. 529, no. 7587, pp. 484–489, 2016.
  * [13] V. Mnih, K. Kavukcuoglu, D. Silver, A. A. Rusu, J. Veness, M. G. Bellemare, A. Graves, M. Riedmiller, A. K. Fidjeland, G. Ostrovski _et al._ , “Human-level control through deep reinforcement learning,” _nature_ , vol. 518, no. 7540, pp. 529–533, 2015.
  * [14] O. Vinyals, I. Babuschkin, W. M. Czarnecki, M. Mathieu, A. Dudzik, J. Chung, D. H. Choi, R. Powell, T. Ewalds, P. Georgiev _et al._ , “Grandmaster level in StarCraft II using multi-agent reinforcement learning,” _Nature_ , vol. 575, no. 7782, pp. 350–354, 2019.
  * [15] M. W. Brittain and P. Wei, _One to Any: Distributed Conflict Resolution with Deep Multi-Agent Reinforcement Learning and Long Short-Term Memory_. [Online]. Available: https://arc.aiaa.org/doi/abs/10.2514/6.2021-1952
  * [16] D.-T. Pham, P. N. Tran, S. Alam, V. Duong, and D. Delahaye, “Deep reinforcement learning based path stretch vector resolution in dense traffic with uncertainties,” _Transportation Research Part C: Emerging Technologies_ , vol. 135, p. 103463, 2022. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0968090X21004514
  * [17] P. N. Tran, D.-T. Pham, S. K. Goh, S. Alam, and V. Duong, “An Interactive Conflict Solver for Learning Air Traffic Conflict Resolutions,” _Journal of Aerospace Information Systems_ , vol. 17, no. 6, pp. 271–277, 2020. [Online]. Available: https://doi.org/10.2514/1.I010807
  * [18] M. W. Brittain, X. Yang, and P. Wei, “Autonomous Separation Assurance with Deep Multi-Agent Reinforcement Learning,” _Journal of Aerospace Information Systems_ , vol. 18, no. 12, pp. 890–905, 2021. [Online]. Available: https://doi.org/10.2514/1.I010973
  * [19] W. Guo, M. Brittain, and P. Wei, “Safety Enhancement for Deep Reinforcement Learning in Autonomous Separation Assurance,” in _2021 IEEE International Intelligent Transportation Systems Conference (ITSC)_ , 2021, pp. 348–354.
  * [20] M. Brittain and P. Wei, “scalable autonomous separation assurance with heterogeneous multi-agent reinforcement learning,” _IEEE Transactions on Automation Science and Engineering_.
  * [21] S. Li, M. Egorov, and M. Kochenderfer, “Optimizing collision avoidance in dense airspace using deep reinforcement learning,” _arXiv preprint arXiv:1912.10146_ , 2019.
  * [22] M. J. Kochenderfer and J. P. Chryssanthacopoulos, “Collision avoidance using partially controlled Markov decision processes,” in _Agents and Artificial Intelligence_ , J. Filipe and A. Fred, Eds. Springer, 2013, vol. 271, pp. 86–100.
  * [23] M. P. Owen, A. Panken, R. Moss, L. Alvarez, and C. Leeper, “ACAS Xu: Integrated collision avoidance and detect and avoid capability for UAS,” in _Digital Avionics Systems Conference (DASC)_. IEEE, 2019, pp. 1–10.
  * [24] L. E. Alvarez, I. Jessen, M. P. Owen, J. Silbermann, and P. Wood, “ACAS sXu: Robust decentralized detect and avoid for small unmanned aircraft systems,” in _Digital Avionics Systems Conference (DASC)_. IEEE, 2019, pp. 1–9.
  * [25] S. M. Katz, L. E. Alvarez, M. Owen, S. Wu, M. W. Brittain, A. Das, and M. J. Kochenderfer, _Collision Risk and Operational Impact of Speed Change Advisories as Aircraft Collision Avoidance Maneuvers_. [Online]. Available: https://arc.aiaa.org/doi/abs/10.2514/6.2022-3824
  * [26] C. J. Watkins and P. Dayan, “Q-learning,” _Machine learning_ , vol. 8, pp. 279–292, 1992.
  * [27] W. Dabney, M. Rowland, M. Bellemare, and R. Munos, “Distributional reinforcement learning with quantile regression,” _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 32, no. 1, Apr. 2018. [Online]. Available: https://ojs.aaai.org/index.php/AAAI/article/view/11791
  * [28] M. Andrychowicz, F. Wolski, A. Ray, J. Schneider, R. Fong, P. Welinder, B. McGrew, J. Tobin, O. Pieter Abbeel, and W. Zaremba, “Hindsight experience replay,” in _Advances in Neural Information Processing Systems_ , I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, Eds., vol. 30. Curran Associates, Inc., 2017. [Online]. Available: https://proceedings.neurips.cc/paper_files/paper/2017/file/453fadbd8a1a3af50a9df4df899537b5-Paper.pdf
  * [29] J. Schulman, S. Levine, P. Abbeel, M. Jordan, and P. Moritz, “Trust region policy optimization,” in _Proceedings of the 32nd International Conference on Machine Learning_ , ser. Proceedings of Machine Learning Research, F. Bach and D. Blei, Eds., vol. 37. Lille, France: PMLR, 07–09 Jul 2015, pp. 1889–1897. [Online]. Available: https://proceedings.mlr.press/v37/schulman15.html
  * [30] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov, “Proximal policy optimization algorithms,” _arXiv preprint arXiv:1707.06347_ , 2017.
  * [31] V. Mnih, A. P. Badia, M. Mirza, A. Graves, T. Lillicrap, T. Harley, D. Silver, and K. Kavukcuoglu, “Asynchronous methods for deep reinforcement learning,” in _Proceedings of The 33rd International Conference on Machine Learning_ , ser. Proceedings of Machine Learning Research, M. F. Balcan and K. Q. Weinberger, Eds., vol. 48. New York, New York, USA: PMLR, 20–22 Jun 2016, pp. 1928–1937. [Online]. Available: https://proceedings.mlr.press/v48/mniha16.html
  * [32] M. Brittain, L. E. Alvarez, K. Breeden, and I. Jessen, “Aam-gym: Artificial intelligence testbed for advanced air mobility,” in _2022 IEEE/AIAA 41st Digital Avionics Systems Conference (DASC)_ , 2022, pp. 1–10.
  * [33] J. Hoekstra and J. Ellerbroek, “BlueSky ATC Simulator Project: an Open Data and Open Source Approach,” in _Proceedings of the seventh International Conference for Research on Air Transport (ICRAT)_ , 06 2016.
  * [34] L. E. Alvarez, J. C. Jones, A. Bryan, and A. J. Weinert, _Demand and Capacity Modeling for Advanced Air Mobility_. [Online]. Available: https://arc.aiaa.org/doi/abs/10.2514/6.2021-2381
  * [35] G. Brockman, V. Cheung, L. Pettersson, J. Schneider, J. Schulman, J. Tang, and W. Zaremba, “OpenAI Gym,” 2016. [Online]. Available: https://arxiv.org/abs/1606.01540
  * [36] E. Liang, R. Liaw, R. Nishihara, P. Moritz, R. Fox, K. Goldberg, J. Gonzalez, M. Jordan, and I. Stoica, “RLlib: Abstractions for distributed reinforcement learning,” in _International Conference on Machine Learning_. PMLR, 2018, pp. 3053–3062.
  * [37] A. Nuic, D. Poles, and V. Mouillet, “Bada: An advanced aircraft performance model for present and future atm systems,” _International journal of adaptive control and signal processing_ , vol. 24, no. 10, pp. 850–866, 2010.
  * [38] V. Mouillet and D. Phu, “Bada family h-a simple helicopter performance model for atm applications,” in _2018 IEEE/AIAA 37th Digital Avionics Systems Conference (DASC)_. IEEE, 2018, pp. 1–8.
  * [39] J. C. Jones, T. Bonin, and E. Mitchell, _Evaluating Wind Hazards for Advanced Air Mobility Operations_. [Online]. Available: https://arc.aiaa.org/doi/abs/10.2514/6.2023-4104
  * [40] O. Yadan, “Hydra - a framework for elegantly configuring complex applications,” Github, 2019. [Online]. Available: https://github.com/facebookresearch/hydra