Artificial Intelligence for Multiclass Rhythm Analysis for Out-of-Hospital Cardiac Arrest During Mechanical Cardiopulmonary Resuscitation

1. Introduction

Out-of-hospital cardiac arrest (OHCA) is one of the leading causes of death worldwide, with an annual incidence of 67 to 170 per 100 000 inhabitants in Europe and survival rates at hospital discharge of around 8% (0% to 18%) [1]. OHCA survival depends on several crucial factors, including bystander cardiopulmonary resuscitation (CPR) with emphasis on chest compresions (CCs), early defibrillation and the overall standard of medical care provided by the emergency medical services (EMS) [2].

Recognizing the patient’s cardiac rhythm throughout resuscitation is crucial for two key reasons: first, to guide therapy according to the treatment pathways defined by the international guidelines; and second, to retrospectively analyze the patient’s response to treatment. Regarding the first case, resuscitation guidelines emphasize the need of discriminating between shockable (Sh) rhythms, comprising ventricular fibrillation (VF) and pulseless ventricular tachycardia (VT), and non-shockable (NSh) rhythms, which include both organized rhythms (OR) and asystole (AS). The Sh/NSh discrimination is the most crucial decision during resuscitation, as defibrillation is the only treatment capable of restoring the normal function of the heart when a Sh rhythm is present [2,3]. A finer classification of the cardiac rhtyhm may also be needed, especially within the NSh rhyhm group, to determine other decisive therapies. For instance, the recommended treatment for AS consists of high-quality CPR, early administration of adrenaline and identification of the underlying cause of the arrest [4], whereas the presence of an OR could be an indicative of the return of spontaneous circulation, in which case the patient should be transported to hospital for post-resuscitation care and recovery [2]. As for the restrospective debriefing of resuscitation episodes, knowledge of the patient’s rhythm throughout the episode may offer valuable insights on the interaction between therapy and physiological response [5,6,7]; this may help identify optimal treatment strategies or clinical interventions that improve OHCA survival. One of the limitations for such retrospective studies is the lack of OHCA databases including cardiac rhythm annotations by expert clinicians, mostly due to the expensive and time-consuming manual labor required. Given all this, there is a clear need for the development of multiclass algorithms that automatically identify the patient’s cardiac rhythm, both in real time and retrospectively.

The state-of-the-art OHCA rhythm classification algorithms are mainly based on the analysis of the ECG, typically consisting of an ECG feature extraction stage followed by a machine learning (ML) classifier. They mostly consist of an ECG feature extraction stage followed by a ML classifier. ECG feature extraction has been approached in time [8,9], frequency [10,11], combined time-frequency [12] and complexity domains [13,14]. The ML approaches explored for the classification stage include K-nearest neighbors [15,16], support vector machines [17,18], artificial neural networks [19] and ensembles of decision trees [20]. Recently, OHCA rhythm classifcation has shifted towards deep learning (DL) techniques, such as convolutional neural networks (CNN) [21,22] or residual networks (ResNet) [23] that avoid the knowledge-based feature extraction process of traditional ML models. As determining the need for defibrillation is crucial in OHCA, the discrimination between Sh and NSh rhythms has been the most commonly addressed classification problem by the aforementioned algorithms. However, a less simplistic cardiac rhythm classification is needed to determine other decisive therapies during CPR [24,25,26]. To address this, Rad et al. [16] introduced the first multiclass OHCA rhythm classifier, where a set of features derived from the discrete wavelet analysis of the ECG were fed into different ML-based classifiers.

These ML and DL-based binary and multiclass algorithms have primarily focused on rhythm classification during interruptions in chest compressions (CCs), as the mechanical activity during CPR introduces ECG artifacts that hinder accurate rhythm detection. As a result, current commercial defibrillators require rescuers to pause CCs every 2 minutes for rhythm analysis. However, these interruptions reduce blood flow to vital organs, decreasing the chances of survival [27]. Over the past few decades, efforts have been made to develop more accurate rhythm analysis algorithms that could be applied during CCs, and thus help minimize CC interruptions [10,28]. These algorithms typically follow a similar structure to those used during non-CC intervals, but include a preliminar filtering step to remove CC-induced artifacts. Such methods have proven successful in classifying both Sh/NSh and multiclass rhythms during OHCA [29,30,31].

While earlier methods primarily focused on manually delivered CCs, the use of mechanical compression devices in OHCA assistance has notably increased. These devices provide CCs at a constant rate and depth, and although evidence supporting improved survival remains inconclusive [32,33,34], their growing adoption underscores the potential benefits they offer. Mechanical devices help ensure the quality of CCs in line with current resuscitation guidelines, even in situations where manual CPR might be compromised, such as during transport [35,36], in confined spaces, or during prolonged resuscitation efforts when rescuer fatigue may impact performance [37]. In addition, these devices alleviate the physical burden on healthcare providers, allowing them to focus on other aspects of patient care. While some proposals have addressed the Sh/NSh rhythm classification in the context of mechanical CCs [38,39,40], no solutions have yet been developed for multiclass rhythm classification.

This study introduces a robust ML/DL framework for a reliable multiclass classification of the cardiac rhythm to distinguish Sh, AS or OR rhythms during mechanical CPR. Both DL and ML architectures are preceded by an adaptive filtering stage to remove mechanical CC artifacts from the ECG. The ML framework includes a feature extraction stage and a random forest (RF) classifier. In the DL framework, the filtered ECG is directly analyzed using CNN- and ResNet-based classification models.

2. Materials

The data used in this study were collected from the Circulation Improving Resuscitation Care (CIRC) trial designed to compare automated load distributing band CPR (LDB-CPR) with high-quality manual CPR (M-CPR), in terms of survival [32,41]. Data were gathered between March 5, 2009 and January 11, 2011 in a randomized, unblinded, controlled group sequential trial of OHCA patients by three US (Fox Valley, Hilsborough, Houston) and two European (Vienna, Nijmegen) EMS. After EMS providers initiated manual CCs, patients were randomized to receive either LDB-CPR or M-CPR. The LDB device (AutoPulse, ZOLL Medical, Chelmsford, MA, USA) delivered CCs in a fixed position, with constant depth of 20% of the patient’s anterior posterior diameter of the chest and at a constant rate of 80 $\min$^-1 ($f_{\mathrm{LDB}}=1.33$ $\mathrm{Hz}$).

Anonymized data from Lifepak 12 and 15 monitor-defibrillators were exported to Matlab (MathWorks Inc., Natick, MA), using Physio-Control’s CODE-STAT data review software, and resampled to a sampling frequency of 250 Hz. The data included the ECG and thoracic impedance (TI) signals of each episode together with the CC instants detected by the CODE-STAT software. Figure 1 corresponds to a 70 $\mathrm{s}$ interval from an OHCA episode showing the ECG (corrupted by CCs) and TI signals in the panels a) and c), respectively. The blue circles on the TI signal indicate the CC instants detected by the CODE-STAT software. As can be seen, each fluctuation in the TI signal corresponds to a CC administered by the EMS. So, in Figure 1, two series of CCs can be clearly distinguished in the $0-15$ $\mathrm{s}$ and $47-70$ $\mathrm{s}$ time intervals corresponding to M-CPR and LDB-CPR, respectively. The middle interval, from 15 $\mathrm{s}$ to 47 $\mathrm{s}$, corresponds to a segment without CCs and therefore free of CC artifacts. Note the clear difference in the TI pattern during LDB-CPR and M-CPR, with much larger amplitude and a more regular pattern for LDB-CPR due to the constant depth and rate of the mechanical CCs. The panel b) of Figure 1 shows the instantaneous CC rate derived from the CC instants marked in panel c). It can be observed that the CC rate for manual CPR was variable and fluctuated around 140 $\min$^-1, but when the LDB device was applied the CC frequency stabilized in 80 $\min$^-1.

Episodes where LDB-CPR was administered were used to conduct this study and the application of the LDB device was identified when the CC rate stabilized at the device’s fixed rate of 80 $\min$^-1 for at least 20 $\mathrm{s}$ (notice the activation of the LDB device in panel b) of Figure 1). Then, 22 $\mathrm{s}$ signal segments were extracted, corresponding to a single cardiac rhythm, and comprising a 6 $\mathrm{s}$ CC-free interval and a 16 $\mathrm{s}$ corrupted by CC artifacts (refer to the highlighted segment in Figure 1). The intervals during CCs were used as inputs for the multiclass decision algorithms, while the artifact-free intervals were employed to annotate the real underlying rhythm of the patient. Rhythms were annotated as Sh, AS or O.R. The segments corresponding to the EMS of Hilsborough, Nijmegen and Vienna were annotated by consensus of three biomedical engineers and subsequently audited by a clinician specialized in the resuscitation field. These segments were used to evaluate the performance of the multiclass decision algorithms. The segments corresponding to the remaining two EMS were not audited by a clinician and were therefore used to train the algorithms.

The final database consisted of 15479 segments extracted from 2058 patients, whereof 9666 segments (1252 Sh, 3865 AS, 4549 OR) from 1178 patients were used to train/develop the multiclass decision algorithms and 5813 segments (1154 Sh, 1616 AS, 3043 OR) from 880 patients were used to test the performance.

3. Methods

This study proposes and evaluates different algorithms for the three-class (Sh, AS, OR) classification of ECG segments corrupted by LDB CC artifacts. All the algorithms were composed of two main stages: 1) an adaptive filter based on a Recursive Least Squares (RLS) to eliminate CC artifacts from the ECG and 2) a classification stage, for which RF, CNN and ResNet models were optimized to classify the rhythm in the filtered ECG as Sh, AS or OR. All the classification models were designed to analyze the filtered ECG in the interval from 2-14 $\mathrm{s}$ (see the highlighted interval in Figure 2) which may help avoid the filtering transients of the RLS filter. In what follows $t=n·T_s$, where $T_s=4\mathrm{m}\mathrm{s}$ is the sampling period ($f_s=250\mathrm{Hz}$), and n is the sample index.

3.1. CPR Artifact Suppressing Filter

During CPR, the corrupted ECG signal, $s_{\mathrm{cor}}\left(n\right)$, recorded by the defibrillator, can be expressed as

$$s_{\mathrm{cor}}\left(n\right)=s_{\mathrm{ecg}}\left(n\right)+s_{\mathrm{cc}}\left(n\right),$$

(1)

where $s_{\mathrm{ecg}}\left(n\right)$ is the patient’s uncorrupted ECG, reflecting the actual underlying heart rhythm and, $s_{\mathrm{cc}}\left(n\right)$ represents the artifact introduced by CCs.

An adaptive RLS filter, tailored for removing periodic interferences [38,42,43],was used to obtain an estimate of the CC artifact ${\widehat{s}}_{\mathrm{cc}}\left(n\right)$, which was then substracted from $s_{\mathrm{cor}}\left(n\right)$ to obtain the filtered ECG $s_{\mathrm{ecg}}\left(n\right)$, i.e an estimate of the true underlying heart rhythm. In this approach, the artifact is assumed to be quasi-periodic, and modeled as a truncated Fourier series of N terms:

$$s_{\mathrm{cc}}\left(n\right)=\sum _{k=1}^N{c}_k\left(n\right)cos(k{\omega }_{0}n+{\theta }_k\left(n\right))=\sum _{k=1}^N\left(a_k\left(n\right)cos\left(k{\omega }_{0}n\right)+b_k\left(n\right)sin\left(k{\omega }_{0}n\right)\right),$$

(2)

where ${\omega }_0$ is the fundamental discrete frequency of CCs which for a LDB device is constant at ${\omega }_0=2\pi f_{\mathrm{LDB}}{T}_s$ with $f_{\mathrm{LDB}}=1.33\mathrm{Hz}\equiv 80{\min }^{-1}$ and $T_s$ the sampling period. Based on this model, the estimated artifact can be expressed in vector format as

$${\widehat{s}}_{\mathrm{cc}}\left(n\right)={\mathbf{\Theta }}^{\mathbf{T}}\left(n-1\right)\mathbf{\Phi }\left(n\right),$$

(3)

where vectors $\mathbf{\Theta }\left(n\right)$ and $\mathbf{\Phi }\left(n\right)$, respectively, define the time-varying coefficients and the in-phase and quadrature reference signals of the Fourier series:

$$\begin{array}{cc}\hfill \mathbf{\Theta }\left(n\right)& ={[a_1\left(n\right)b_1\left(n\right)\cdots a_N\left(n\right)b_N\left(n\right)]}^T\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill \mathbf{\Phi }\left(n\right)& ={[cos\left({\omega }_{0}n\right)sin\left({\omega }_{0}n\right)\cdots cos\left(N{\omega }_{0}n\right)sin\left(N{\omega }_{0}n\right)]}^T\hfill \end{array}$$

(5)

The RLS filter estimates the $a_k\left(n\right)$ and $b_k\left(n\right)$ coefficients adaptively over time so that the error between the corrupted ECG, $s_{\mathrm{cor}}\left(n\right)$, and the estimated artifact, ${\widehat{s}}_{\mathrm{cc}}\left(n\right)$, is minimized at each iteration at the harmonics of the LDB device compression frequency, $f_{\mathrm{LDB}}$. Note that in this configuration (Figure 2), the error signal corresponds to the filtered ECG, ${\widehat{s}}_{\mathrm{ecg}}\left(n\right)$. The update equations of the filter are given by

$$\begin{array}{cc}\hfill {\widehat{s}}_{\mathrm{ecg}}\left(n\right)& =s_{\mathrm{cor}}\left(n\right)-{\widehat{s}}_{\mathrm{cc}}\left(n\right)\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill \mathbf{F}\left(n\right)& ={\displaystyle \frac{1}{\lambda }}\left[\mathbf{F}\left(n-1\right)-{\displaystyle \frac{\mathbf{F}\left(n-1\right)\mathbf{\Phi }\left(n\right){\mathbf{\Phi }}^T\left(n\right)\mathbf{F}\left(n-1\right)}{\lambda +{\mathbf{\Phi }}^T\left(n\right)\mathbf{F}\left(n-1\right)\mathbf{\Phi }\left(n\right)}}\right]\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \mathbf{\Theta }\left(n\right)& =\mathbf{\Theta }\left(n-1\right)+\mathbf{F}\left(n\right)\mathbf{\Phi }\left(n\right){\widehat{s}}_{\mathrm{ecg}}\left(n\right),\hfill \end{array}$$

(8)

where the gain matrix, $\mathbf{F}\left(n\right)$, and the coefficient vector were initialized to $\mathbf{F}\left(0\right)=0.03{\mathbf{I}}_{2N}$ and $\mathbf{\Theta }\left(0\right)={\mathbf{0}}^T$. As shown in the previous equations, there are two configurable parameters in the RLS filter: the number of harmonics, N, that model the artifact, and the forgetting factor, $\lambda$, which provides a trade-off between the filter’s adaptability and stability. These values were fixed to $N=35$ and $\lambda =0.989$, based on the optimal RLS configuration identified in [39], where the RLS filter was evaluated in terms of the performance of a Sh/NSh decision algorithm applied to the filtered signal. Panel a) of Figure 2 shows the adaptive RLS filtering schema, while panel b) displays the input and output signals, from top to bottom: the corrupted ECG, $s_{\mathrm{cor}}\left(n\right)$, the estimated CC artifact, ${\widehat{s}}_{\mathrm{cc}}\left(n\right)$, and the filtered ECG revealing the underlying rhythm of the patient, which in this case corresponds to an OR rhythm.

3.2. Optimization and Evaluation

The training set, composed of 9666 segments derived from 1178 patients, was used in a 10-fold cross validation (CV) approach for hyperparameter optimization and, in the case of RF models, for feature selection; details on the hyperparameters and feature selection are provided in the upcoming sections corresponding to each model. Data were divided patient-wise and ensuring that every fold retained a prevalence of each rhythm comparable to that of the entire dataset.

The remaining 5813 segments, corresponding to 880 patients, were used to evaluate the performance of the classifiers. For each class $i\in \left\{\mathrm{Sh},\mathrm{AS},\mathrm{OR}\right\}$, the sensitivity (${\mathrm{Se}}_i$), positive predictive value (${\mathrm{PPV}}_i$) and F₁-Score_i were computed, and the unweighted mean of all sensitivities (UMS), total accuracy (ACC) and unweighted mean of all F₁-Scores (UMFS) were used as summarizing metrics:

$$\begin{array}{cccccc}\hfill {\mathrm{Se}}_i& ={\displaystyle \frac{{\mathrm{TP}}_i}{{\mathrm{TP}}_i+{\mathrm{FN}}_i}},\hfill & \hfill {\mathrm{PPV}}_i& ={\displaystyle \frac{{\mathrm{TP}}_i}{{\mathrm{TP}}_i+{\mathrm{FP}}_i}},\hfill & \hfill {\mathrm{F}}_1\text{-}{\mathrm{Score}}_i& =2·{\displaystyle \frac{{\mathrm{PPV}}_i·{\mathrm{Se}}_i}{{\mathrm{PPV}}_i+{\mathrm{Se}}_i}}\hfill \\ \hfill \mathrm{UMS}& ={\displaystyle \frac{1}{3}}\sum _{i=1}^3{\mathrm{Se}}_i,\hfill & \hfill \mathrm{ACC}& ={\displaystyle \frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{TN}+\mathrm{FP}+\mathrm{FN}}},\hfill & \hfill \mathrm{UMFS}& ={\displaystyle \frac{1}{3}}\sum _{i=1}^3{\mathrm{F}}_1\text{-}{\mathrm{Score}}_i\hfill \end{array}$$

where ${\mathrm{TP}}_i$, ${\mathrm{TN}}_i$, ${\mathrm{FP}}_i$ and ${\mathrm{FN}}_i$ are the true positives, true negatives, false positives and false negatives for class i.

3.3. Classical Machine Learning Architecture

The classical ML algorithm is based on a state-of-the-art solution originally designed for multiclass OHCA rhythm classification during manual CPR [31]. In essence, the algorithm integrates a multiresolution ECG analysis approach, employing the Stationary Wavelet Transform (SWT) for feature extraction and a RF classifier for the subsequent classification. The SWT decomposes the 12 $\mathrm{s}$ window into 7 detail coefficients (d₁-d₇) sub-bands using a Daubechies 4 mother wavelet. A denoised version of the ECG is also reconstructed using the detail coefficients d₃ to d₇, corresponding to an analysis frequency band of 0.98–31.25 $\mathrm{Hz}$.

From the denoised ECG and the detail coefficients d₃-d₇, ninety-three features were extracted to characterize OHCA rhythm subtypes, representing over 25 years of research in the field. These features were divided into five analysis domains. Time domain features included characteristics like the mean and the standard deviation of the heart rate [25]; spectral features included classical measures like VF leakage [44] or the power proportion concentrated around the VF-fibrillation band [29]; complexity analysis covered entropy measures like sample or Shannon entropy [45]; statistical analysis measures involved amplitude distribution characteristics; finally, phase space features utilized time-delay embedding to extract dynamics in the ECG. The detailed description of each of the 93 features can be found in [31].

The training set was divided using a 10-fold CV approach to select the best subset of K features and to optimize RF hyperparameters. First, the optimal set of features was selected for each of the 10 training folds that constitute the 10-fold CV. Feature selection was based on a recursive feature elimination (RFE) approach, using the out-of-bag permutation importance as a ranking criterion [46,47]. Permutation importance is an inherent characteristic of the RF classifier that evaluates the significance of each feature by randomly shuffling its values in the training data of each tree in the forest and then measuring the resulting change in the out-of-bag error. In each iteration of the RFE algorithm, features were ranked, and the least important 3% were removed. This process was repeated until the optimal sets of K features, with $K=\left\{1,3,5,7,10·j\right\}$ and $j=1,...,9$, were selected for classification. Once the best K-feature subsets were selected in each of the 10 CV training folds, the RF classifier was optimized. Only one parameter of the RF classifier was considered for optimization: the minimum number of observations per leaf, $l_{\mathrm{size}}$, which controls the depth of the trees and was identified in [31] as essential for preventing overfitting. For every CV training fold and subset of K features, different $l_{\mathrm{size}}$ values were trained in the range $1\le l_{\mathrm{size}}\le 200$, and evaluated in the corresponding testing fold. The number of trees was set to $B=500$ and the number of predictors per split was fixed to the default value $\sqrt{K}$ for both feature selection and $l_{\mathrm{size}}$ optimization. $B=500$ is sufficient to stabilize accuracy without causing overfitting [48] and the default value of the number of predictors per split achieved by far the best performance in [31]. Finally, $l_{\mathrm{size}}$ was fixed to the default value $l_{\mathrm{size}}=1$ during the feature selection process.

3.4. Algorithm Based on CNNs

Figure 3 a) shows in blue the architecture of the multiclass OHCA rhythm classification algorithm based on CNNs. This architecture is based on the one proposed in [30] to discern Sh and NSh rhythms during manual CPR. The filtered ECG (1-D signal of $N=3000$ samples), is introduced to a CNN composed of B convolutional blocks (Figure 3 shows a 3-block as an example, B will be a trainable parameter) that extract the high level features of the signal, followed by two fully connected layers for the three-class classification. The b-th convolutional block is comprised of a one-dimensional convolutional layer (Conv1D) with $J_b$ filters of width $I_b$ , followed by a batch normalization (BN) layer, a rectified linear unit (ReLU), a max-pooling layer and a dropout layer.

The input to the first convolutional block is defined as ${\mathrm{s}}_0(n,1)={\widehat{s}}_{\mathrm{ecg}}\left(n\right)$. The expression ${\mathrm{s}}_{b-1}(n,m)$ will refer to the input of block b, or equivalently, the output of block $b-1$, where n and m represent the time and filter index, respectively. The output of the Conv1D at b-th convolutional block can be formulated as:

$$c_b(n,m)=b^m+\sum _{\rho =1}^{J_{b-1}}\sum _{i=1}^{I_b}{w}_{\rho ,i}^m{s}_{b-1}(n+i-1,\rho )$$

(9)

where the filter weights $w_{\rho ,i}^m$, and the biases for channel shifting $b^m$ are the learnable parameters adjusted during training.

BN layers modify the output of the preceding layer to prevent complex weight interactions from altering the data distribution. This accelerates training by allowing the use of larger learning rates and improves generalization while reducing overfitting [49]. For every training mini-batch $\mathcal{B}$, a BN layer calculates the channel-wise means ${\mu }_{\mathcal{B},m}$ and variances ${\sigma }_{\mathcal{B},m}^2$ and then normalizes each channel following

$${\widehat{c}}_b(n,m)={\displaystyle \frac{c_b(n,m)-{\mu }_{\mathcal{B},m}}{\sqrt{{\sigma }_{\mathcal{B},m}^2+ϵ}}},$$

(10)

where $ϵ$ is a small value included to ensure numerical stability. The normalized channels are then adjusted through scaling and shifting to optimize the final ReLu layer. As a result, the outputs, $z_b(n,m)$, can be expressed as

$$z_b(n,m)={\gamma }_m·{\widehat{c}}_b(n,m)+{\beta }_m,$$

(11)

where ${\gamma }_m$ and ${\beta }_m$ are trainable parameters. On inference, a moving average of the mini-batch means ${\mu }_{\mathcal{B},m}$ and variances ${\sigma }_{\mathcal{B},m}^2$ observed during training is typically applied in (9).

Max-pooling layers downsample input data by taking the maximum value from each block of K elements along the time dimension n, so that the output for block b can be represented as:

$$s_b(n,m)=\max {\left\{z_b(k,m)\right\}}_{k=(n-1)·K+1,...,n·K}$$

(12)

Finally, the ReLU layers add nonlinearity to the network through the activation function $f\left(x\right)=\max \left\{0,x\right\}$, enabling the model to learn intricate nonlinear mappings.

Zero-padding was applied before the convolution operations, so the only reduction of dimensionality was due to the max-pooling layers ($K=3$). The dropout layer at the end of each block serves as a regularization mechanism, operating exclusively during training to prevent overfitting. This layer temporarily disables a randomly chosen fraction of the network’s adjustable parameters. The output of the last convolutional block was flattened and fed to a dense network composed of two fully connected layers with 10 and 3 neurons, respectively. Finally, a softmax layer transformed the output of the final 3 neurons into values ranging from $\left[0,1\right]$, representing the probabilities that a given segment corresponds to a Sh ($p_{\mathrm{Sh}}$), AS ($p_{\mathrm{AS}}$) or OR ($p_{\mathrm{OR}}$) rhythm.

The CNN parameters adjusted during the training phase were the following: the number of convolutional blocks $B=\left\{3,4,5,6\right\}$, the width of the filters $I=\left\{2,4,8,16,32,64,92\right\}$ (the same filter width is considered along the B blocks) and the number of filters which varies from block to block $\mathbf{L}=\left(J_1,J_2,...,J_B\right)$. Six filter configurations, with increasing number of filters (from sparse to dense), were studied: ${\mathbf{L}}_{\mathbf{1}}=(1,2,4,8,16,32)$, ${\mathbf{L}}_{\mathbf{2}}=(2,4,8,16,32,64)$, ${\mathbf{L}}_{\mathbf{3}}=(4,8,16,32,64,128)$, ${\mathbf{L}}_{\mathbf{4}}=(6,12,24,48,96,192)$, ${\mathbf{L}}_{\mathbf{5}}=(8,16,32,64,128,256)$, ${\mathbf{L}}_{\mathbf{6}}=(10,20,40,80,160,320)$. The values in parentheses correspond to the number of filters $J_b$ for blocks $b=1,...,6$. For architecures with $B<6$ blocks, central values (with upwards bias) were selected, so for 3, 4 and 5 blocks the ${\mathbf{L}}_{\mathbf{2}}$ configuration would be for instance: (8,16,32), (4,8,16,32) and (4,8,16,32,64), respectively.

3.5. Algorithm Based on ResNets

The third architecture explored was a ResNet, which mitigates the issue of performance deterioration as layers are added, enabling deeper networks [50]. The main components of a ResNet are residual blocks, which comprise a main path, including convolutional, batch normalization and other typical CNN layers, and a shortcut path, which directly connects the input and the main path output. Be x the input to a residual block, and $H\left(x\right)$ the desired data transformation; instead of learning this transformation directly, residual blocks focus on learning the difference between the input and the output, called the residual $F\left(x\right)=H\left(x\right)-x$. This is achieved by simple addition of the main path and shortcut path outputs, and makes it easier for the network to learn by focusing on refining the input rather than completely transforming it.

Figure 3 b) shows the layout of the ResNet architecture [50] which, similarly to in Jaureguibeitia et al. [23], intended to replicate that of the CNN, deepening the network while maintaining a coherent structure. Similarly to the CNN, the network was composed of $B=3,4,5,6$ blocks, each consisting of two residual blocks following the main path pre-activation configuration (conv-BN-ReLu-conv-BN) proposed by Han et al. [51]; the first block of the network was an exception to this rule, and consisted of a single, much simpler conv-BN-ReLU configuration with no shortcut path. Pooling layers were replaced by strided convolutions, which skip every other step in the filtering process. When adjustments to length and depth were needed, the shortcut path of the first residual block included a strided convolution to create a linear projection of the input. As in the CNN, each convolutional layer in the network used an identical filter size selected from $I=2,4,8,16,32,64,92$. Similarly, the possible configurations of number of filters were selected from ${\mathbf{L}}_{\mathbf{1}}-\mathbf{L}\mathbf{6}$, with the number of filters per block applied to each convolutional layer within that block. Finally, the hidden fully connected layer was replaced by a global average pooling layer, which outputs the mean value of each input channel [52]. The final values of B, I, and configuration $\mathbf{L}$, were all optimized during the training phase.

Both in the CNN and ResNet architectures, the weights and biases of every layer were optimized as to minimize categorical crossentropy, using stochastic gradient descent with a momentum of 0.8. The initial learning rate was fixed at 0.02 and it was reduced by a factor of 0.8 at every epoch. The training process was conducted for 20 epochs, with a batch size of 256 samples [53].

4. Results

Figures 4 and 5 show the results obtained by the classical ML algorithm on the training data. The left panel of Figure 4 shows the mean performance metrics (UMS, UMFS, and average between UMS and UMFS) as a function of the number of features K; these were calculated as the average, across all $l_{\mathrm{size}}\left(i\right)$ values considered, of the CV performance metrics obtained for K-feature and $l_{\mathrm{size}}=l_{\mathrm{size}}\left(i\right)$ RF models. The best compromise between model simplicity and performance was obtained for $K=20$ as metrics barely increase for a greater value of K. The right panel of Figure 4 shows the CV performance metrics as a function of $l_{\mathrm{size}}$; RF models of $K=20$ features were considered, based on the previous results. In terms of average between UMS and UMFS, the optimal range for $l_{\mathrm{size}}$ was 1 $<=l_{\mathrm{size}}<=9$, with a significant decline in performance observed for larger values of $l_{\mathrm{size}}$. The value $l_{\mathrm{size}}=1$ was chosen as optimal, as it produced the most balanced UMS and UMFS results. Figure 5 shows the 30 ECG features with higher probability of selection for the $K=20$ and $l_{\mathrm{size}}=1$ CV RF models. These probabilities were estimated by counting the number of times the features were selected in the 10 iterations of the feature selection algorithm in the 10-fold CV loop. The most important 9 parameters, i.e. those that were selected in the 100% of iterations, were highly heterogeneous, deriving from all detail coefficients as well as from the denoised ECG, and corresponding to complexity (SampEn), time (bCP), statistical (IQR, StdAbs1, Hmb, Hcmp) and phase space domains (SkewPSD). The acronyms used in Figure 5 are the same as those found in [31], where a detailed description of each parameter is provided. Given these results, a single optimal RF configuration was defined, using $l_{\mathrm{size}}=1$ and the most important $K=20$ features as per Figure 5. On a 10-fold CV loop over training data, this configuration obtained identical UMS and UMFS scores of 81.3%. As a single model trained on the complete training data set and evaluated on the test set, it obtained a UMS and UMFS of 85.3% ad 85.1%, respectively (all the details can be cound in Table 1).

The impact of altering the main parameters of the CNN architecture is depicted in Figure 6. The top row shows the CV performance metrics (UMS, UMFS, and average between UMS and UMFS) for a varying filter size I and a fixed filter configuration ${\mathbf{L}}_3$. The best performance in terms of average between UMFS and UMS (third column on both rows) was obtained for a filter width of $I=64,32,32$ and 16 when $B=3,4,5$ and $B=6$ were used. The bottom row studies the effect of changing the filter configuration when the filter size is fixed at/to these optimal values. The ${\mathbf{L}}_3$ configuration for $B=5$ blocks and the ${\mathbf{L}}_5$ configuration with $B=6$ blocks achieved the best and very similar performances, with the 6-block one being slightly better. Therefore, the optimal model for the CNN architecture was composed of 6 convolutional blocks with 8, 16, 32, 64, 128 and 256 filters, all of them of width 16. In the training set, this model achieved a UMFS and UMS of 83.0% and 84.8%, respectively. In the testing set, the model obtained a UMFS and UMS of 86.1% and 87.5%, respectively (all the details can be found in Table 1).

Figure 7 analyzes the effect of changing the parameters of the ResNet architecture. Similar to CNNs, the filter width was optimized first with the filter configuration fixed at ${\mathbf{L}}_{\mathbf{3}}$ (first row) for fixed optimal filter sizes and then followed by the selection of filter configuration (second row). The best classification results were obtained for 4 blocks. Adding a fifth block increased the complexity of the network (number of trainable parameters) and slightly decreased performance. Using only 3 blocks resulted in a large decrease in performance, or an overly simplistic model. As the third column of the first row shows, the best performance in terms of average between UMFS and UMS was obtained for $I=32$ when $B\le 5$ and $I=16$ when $B=6$. With the filter width fixed, the second row of Figure 7 shows that the best performance was achieved for 4 blocks and ${\mathbf{L}}_{\mathbf{4}}$ filter configuration. Thus, the optimal configuration of the ResNet architecture consisted of 4 blocks containing 6, 12, 24 and 48 filters respectively, all of them of width 32. That is, the ResNet input block used 6 filters and the remaining 3 blocks used 12, 24 and 48 filters in each of the residual blocks. In the training set, this configuration obtained a UMFS and UMS of 86.9% and 86.6%, respectively. In the testing set, the ResNet model obtained a value of 88.3% for both the UMFS and the UMS (all the details can be found in Table 1).

Table I shows the results of the OHCA rhythm classification algorithms based on classical ML, CNN and ResNet models when the aforementioned optimal configurations were applied in to the testing set. As the summarizing metrics demonstrate, the DL-based algorithms performed better than the traditional ML ones. The CNN model outperformed the RF model by 1 percentage point in UMFS and 2 percentage points in UMS and ACC. The ResNet was the best-performing model, outperforming the CNN model by approximately 2 percentage points in UMFS and ACC and one percentage point in UMS.

Given the importance of the Sh/NSh discrimination during resuscitation therapy and as additional experiment, the best performing method (ResNet) was adapted for the binary Sh/NSh classification task. For this experiment, AS and OR rhythms were grouped into the NSh category, using 9666 segments (1252 Sh, 8414 NSh) to train and 5813 segments (1154 Sh, 4659 NSh) to test the model. The optimal ResNet architecture for this problem was selected in an analogous manner to that for the 3-class problem. Performance was evaluated in terms of sensitivity (SE, the proportion of correctly classified Sh rhythms) and specificity (SP, the proportion of correctly classified NSh rhythms), in line with the minimum performance requirements for Sh/Nsh discrimination recommended by the AHA; the balanced accuracy (the mean of SE and SP) was chosen as summarizing metric. Figure 8 shows the performance metrics obtained across the 10-fold CV loop in the training set. As in previous figures, the top row analyzes the impact of altering the width of the filters, I, while the bottom one explores the impact of altering the configuration of the filters, $\mathbf{L}$, for a fixed filter size. As shown in the third panel of the bottom row, the ${\mathbf{L}}_6$ configuration with 3 and 4 blocks and $I=32$ achieved the maximum BAC, with the 4-block one being slightly higher. This architecture obtained a SE/SP/BAC scores of 85.5%/98.6%/92.0% and 90.6%/98.5%/94.6% in the training and testing sets, respectively.

5. Discussion

The adoption of mechanical CPR devices has significantly increased in recent years, primarily through two main technologies: LDB and piston-driven mechanical CC devices. Mechanical CPR guarantees high-quality CCs when manual CPR is i) subjected to fatigue, ii) practically challenging or iii) cannot be delivered safely.

This is, to the best of our knowledge, the first study that covers multiclass OHCA rhythm classification during mechanical CPR. The proposed architectures integrated the collective knowledge derived from active research in the rhythm classification field, considering an adaptive filter for CPR artifact suppression, followed by a multiclass rhythm classifier for which the performance of a traditional ML algorithm and two DL architectures based on CNNs and ResNets were compared.

DL-based algorithms outperformed the classical ML algorithm by at least 2 percentage points in UMS and ACC. Considering that the classical ML algorithm relies on over 20 years of expert knowledge in ECG feature engineering for OHCA rhythm classification, these results highlight the power of DL algorithms to learn discriminative features by leveraging all the hidden information in the ECG. This simplifies the feature extraction process, saving time and, more importantly, improving the quality of the features extracted.

The algorithm based on ResNets offered the best performance, achieving a UMFS, UMS and ACC of 88.3%, 88.3% and 88.2% for the 3-class classification task, respectively. This performance is similar to that obtained in [31], the only study in the literature that analyzes multiclass OHCA rhythm classification during manual CPR. The characteristics of manual compressions are rescuer dependant, which means the variability of the resulting artifacts anticipates a more complex filtering challenge. However, manual artifacts showed significantly smaller artifact amplitudes and less harmonic components (smaller bandwidths) compared to LDB artifacts, and this balance resulted in similar levels of accuracy [40,54,55]. For the Sh/NSh problem, the BAC was 94.6%, with a SE of 90.6% and a SP of 98.5%. This is a very important problem since it addresses shock advice decisions during CPR. Shock advice algorithms for defibrillators are normally tested on artifact-free data. In that scenario, the AHA requires a minimum SE and SP of 90% and 95%, respectively[56]. Our solution met those requirements.

The database used in this study was fully derived from OHCA cases. It is unclear whether the proposed algorithms would perform differently for in-hospital data; however, given that in-hospital resuscitation does not entail differences in LDB-CPR, it is tempting to argue that the analysis and results we have presented for OHCA patients also will be valid for in hospital cardiac arrest and CPR.

In conclusion, this study presents the first method for an automatic and clinically safe multiclass OHCA rhythm analysis during LBD-CPR. The proposed solution together with already available solutions for piston-driven [38] and manual CCs [31] would cover rhythm analysis in every CPR scenario. This may open the possibility of a reliable multiclass OHCA rhythm analysis during CPR, contributing to guide therapy while reducing no-flow intervals, and thereby improving survival in OHCA. In addition, this method could be used to automatically annotate large OHCA datasets in order to retrospectively analyze the effects of therapy on OHCA outcome, contributing to identify those treatments that increase OHCA survival.