# N-ACT: An Interpretable Deep Learning Model for Automatic Cell Type
and Salient Gene Identification

A. Ali Heydari  Oscar A. Davalos  Katrina K. Hoyer  Suzanne S. Sindi

###### Abstract

Single-cell RNA sequencing (scRNAseq) is rapidly advancing our understanding
of cellular composition within complex tissues and organisms. A major
limitation in most scRNAseq analysis pipelines is the reliance on manual
annotations to determine cell identities, which are time consuming,
subjective, and require expertise. Given the surge in cell sequencing,
supervised methods–especially deep learning models–have been developed for
automatic cell type identification (ACTI), which achieve high accuracy and
scalability. However, all existing deep learning frameworks for ACTI lack
interpretability and are used as “black-box” models. We present _N-ACT_
(Neural-Attention for Cell Type identification): the first-of-its-kind
interpretable deep neural network for ACTI utilizing neural attention to
detect salient genes for use in cell-types identification. We compare N-ACT to
conventional annotation methods on two previously manually annotated data
sets, demonstrating that N-ACT accurately identifies marker genes and cell
types in an unsupervised manner, while performing comparably on multiple data
sets to current state-of-the-art model in traditional supervised ACTI.

Interpretable Deep Learning, Deep Learning, single-cell RNA sequencing,
Attention, Machine Learning, Computational Biology, ICML

AReferences (Appendix)

##  1 Introduction

Single-cell RNA-sequencing (scRNAseq) technologies allow for measuring
transcriptome-wide gene expression at the single-cell level. In contrast to
bulk-RNA sequencing, scRNAseq can elucidate dynamic expression patterns
between different cellular populations, providing a tremendous advantage when
studying organisms as well as delineating intra-population heterogeneities
(Erfanian et al., 2022).

Accurate identification of cell types in scRNAseq studies remains a
challenging and time-consuming task (Pasquini et al., 2021). Cell type
annotation is often performed manually by experts – a long, laborious, and
subjective process (Pasquini et al., 2021; Clarke et al., 2021). To mitigate
these challenges, researchers have developed automatic cell type
identification (ACTI) pipelines, including deep learning (DL) models such as
ACTINN (Ma & Pellegrini, 2019), which learn from datasets that have _already
been annotated from the same (or similar) populations_ (see (Pasquini et al.,
2021) for a review on ACTI). However, the supervision required for these
approaches limits their utility for studies which lack prior knowledge of
tissue- or sample-specific cell types. Such models can not provide the much-
needed biological interpretability on the algorithm’s decision-making needed
to translate scRNAseq findings to inform experimental design.

In this work, we introduce _Neural-Attention for Cell Type identification
(N-ACT)_ : An interpretable _unsupervised_ DL model that employs _attention_
to detect salient genes, and applies this information to identify specific
cell-types. Our results on multiple datasets show that N-ACT accurately
predicts preliminary annotations with no prior knowledge about the system,
providing a valuable complementary framework to experimental studies and
computational pipelines.

##  2 Methods and Approach

Figure 1: N-ACT, an Interpretable Model for Unsupervised/Supervised Cell-Type
Identification. N-ACT consists of flexible stages and modules that can be
modified for different objectives. We focus on unsupervised ACTI, which
consists of the following stages: _(Stage I)_ Generating labels through graph-
based clustering followed by identifying five thousand highly-variable genes
for efficient training (standard practice in scRNAseq pipelines), _(Stage II)_
Training the N-ACT DL core to predict generated (or true) labels, with the
goal of identifying salient genes to predict the correct cell types (labels),
and _(Stage III)_ Interpreting model predictions by extracting attention
values, thus constructing a ranked list of “attentive genes” used to compare
to existing literature (referred to as “Querying Salient Genes” in the figure)
and thus predicting cell types. We describe each component in more detail in
the Appendix.

N-ACT’s framework consists of three stages (shown in Fig. 1): (I) Assigning
labels (if no labels are available) and identifying highly-variable genes
(HVGs), (II) a deep neural network (NN), which is the core of N-ACT, and (III)
interpretation. In this section, we present stage II, i.e. N-ACT’s DL core,
and later describe stages (I) and (III) in Section 4. N-ACT’s DL core has
three modules: (1) Additive attention module, responsible for learning the
importance of each gene; (2) Multi-headed projection module tasked with
learning a set of non-linear operators mapping attention outputs to the last
layer, i.e. (3) Task module, which is designed to fulfill a specific
downstream task.

###  2.1 Additive Attention

Attention (or neural-attention) is a recent DL mechanism that has transformed
computer vision and natural language processing research (refer to (Chaudhari
et al., 2021) for a review on attention in DL). Attention networks aim to
mimic the way humans understand “context” in sentences or details in images by
focusing on a subset of significant features for a given objective. The use of
attention-based NN for scRNAseq analysis is still in its infancy, with only a
few successful applications to date. To identify salient genes, we use an
additive attention module (Bahdanau et al., 2015) in a feed-forward NN
(similar to (Raffel & Ellis, 2015)) aiming to learn the optimal weighting
(_importance_) of all genes for each cell, given a downstream task.

The first step in the DL core (attention) is used to calculate a gene-score
matrix (weighted version of scRNAseq count matrix), representing expression
data in later layers. These importance scores enable gene prominence
quantification for the downstream task, allowing for interpretation of model
decision making. Given a gene expression matrix $X\in\mathbb{R}^{C\times N}$,
where $C$ and $N$ denote the number of cells and genes, respectively, we
define the gene-score matrix $\Gamma$ and the attention weights $A$ as shown
in Eq. (1):

$\Gamma=A\odot X,\text{where
}A_{i,j}=\frac{e^{L_{i,j}}}{\sum_{j=1}^{N}e^{L_{i,j}}},$ (1)

with $L=NN(X)$ denoting a linear NN111$NN:\mathbb{R}^{C\times
N}\to\mathbb{R}^{C\times N}$ is a linear operator of the form $NN(X)=XW+B$,
with input $X$ and biases $B\in\mathbb{R}^{C\times N}$ and weights
$W\in\mathbb{R}^{N\times N}$.. The learned operator $A$ is leveraged after
training to identify important (or “attentive”) genes for interpretability.

###  2.2 Multi-Headed Projections

Our _Projection_ mechanisms are intermediate layers between the attention
layer and the downstream task module. The goal in using projection modules is
to strike a balance between model capacity and efficiency: Too much capacity
(e.g. too many non-linear layers) could lead to significant over-fitting,
while insufficient capacity prevents the model from learning the correct
representations. We design the projection blocks to be _multi-headed_
(consisting of $h\in\mathbb{N}$ separate linear operators), a concept shown by
(Vaswani et al., 2017) as effective in learning different representations.
Such design allows efficient consideration of different gene subsets and
improves model performance without the need for numerous non-linear layers.
N-ACT consists of $k$ projection blocks, each consisting of $h$ heads. Outputs
from each head is then concatenated and inputted to a point-wise feed-forward
network (equivalent to $1\times 1$ convolution layer) with a Rectified Linear
Unit (ReLU) (Nair & Hinton, 2010) that adds non-linearity. Through careful
ablation studies, we found that $k=2$ projection blocks with $h=10$ heads
provides the appropriate balance of accuracy and efficiency. Details on
projection blocks and relevant ablation studies are provided in Appendix F.

###  2.3 Task Module and Architecture Choices

The last stage of our DL core is the _task module_ , which can be adjusted
based on desired objectives. Given our ACTI goal, we chose a non-linear
mapping between the projection block’s output and the labels (either provided
[supervised] or generated in the earlier stages [unsupervised]).

Our task module connects the projections to the number of labels, followed by
Leaky ReLU activation (Xu et al., 2015) (depicted in Appendix F). N-ACT
minimizes a standard cross entropy loss (Appendix C) using the Adam gradient-
based optimizer (Kingma & Ba, 2014) at a learning rate $\alpha=10^{-4}$ for 50
epochs. Additional information on N-ACT’s training scheme is provided in
Appendix F.

##  3 Datasets Studied

We tested the model’s ability to learn cell-types on four datasets (two for
supervised and two for unsupervised ACTI). Results are presented on different
datasets for each learning setting to showcase N-ACT’s versatility and utility
for different systems and species (similar results were achieved on all
datasets in all tasks). Data were minimally pre-processed (only quality-
controlled) and divided roughly 85%:15% for training and testing using
“balanced split” (described in (Heydari et al., 2022)). A brief description of
each dataset is provided below, with more details on pre-processing and the
data presented in Appendix B.

Datasets for supervised training: (1) Mouse HDF (PubMed ID: 34548614) consists
of scRNAseq of murine aortic cells of mice on a normal diet versus mice on a
high-fat diet, resulting in 24K cells and 10 annotated populations after
processing. (2) Immune CSF (PubMed ID: 33382973) profiles single-cells from
cerebrospinal fluid (CSF) of immune CSF, viral encephalitis, and non-
inflammatory and autoimmune neurological disease. Cells were isolated from 31
patients, resulting in 80K cells (70K cells after processing) and 15 annotated
populations.

Datasets for unsupervised training: (1) COVID PBMC (PubMed ID: 33357411)
profiles the transcriptional immune dysfunction triggered in moderate and
severe COVID-19 patients using scRNAseq. Peripheral blood mononuclear cells
(PBMC) were isolated from 20 patients and were sequenced resulting in 69K
cells (64K cells after processing) with 9 cell populations. (2) Immune cSCC
(Pub Med ID: 32579974) consists of scRNAseq from healthy skin and cutaneous
squamous cell carcinoma (cSCC) tumors. 10 patients with cSCC tumors had
healthy skin and tumor cells sequenced, resulting in 48K cells (47K cells
after processing) and 14 cell populations.

##  4 Results

In this section, we provide the results of using N-ACT for datasets described
in Section 3. Standard evaluation tools were used to measure model performance
with each metric detailed in Appendix C.

###  4.1 Supervised ACTI Performance

We benchmark N-ACT against the current state-of-the-art supervised model,
ACTINN, though the goal of this work is to generate unsupervised annotations.
To show the importance of our novel architecture in effective attention
utilization, we added the same feed-forward attention module to ACTINN
(denoted by “ACTINN+ATTN”), which significantly hindered the model’s
performance (Table 1). Given the comparable performance of our model to the
state-of-the-art DL algorithm, our results show that N-ACT can effectively
learn supervised ACTI. Moreover, the poor performance of ACTINN + ATNN
highlights the importance of appropriate architectures needed for effective
use of attention for interpretability. We next consider unsupervised cell-type
annotation on two previously annotated datasets.

Table 1: Benchmarking N-ACT on Supervised ACTI. Although the main goal of N-ACT is interpretable ACTI in an unsupervised manner, our model can be used in a supervised setting as well, while still providing biological interpretability. W-F1: Weighted F1 score, NW-F1: Non-Weighted F1 score, Interp: Interpretability  Model | W-F1 | NW-F1 | Interp?
---|---|---|---
| Mouse HDF |  |
ACTINN | 0.9703 | 0.9677 | No
ACTINN+Attn | 0.8759 | 0.7438 | Yes
N-ACT (Ours) | 0.9681 | 0.9712 | Yes
| Immune CSF |  |
ACTINN | 0.9357 | 0.8898 | No
ACTINN+Attn | 0.7528 | 0.2413 | Yes
N-ACT (Ours) | 0.9285 | 0.8963 | Yes

###  4.2 Unsupervised ACTI Performance

Figure 2: Evaluating N-ACT’s Unsupervised ACTI. [(A), (C)]: To evaluate the
interpretability and the performance of the attentive genes, we calculated the
Pearson correlations of mean attention values in each population (per gene)
for all cell types in (A) Immune cSCC and (C) COVID PBMC. [(B), (D)]: To
assess N-ACT’s predictive capabilities, we queried the TF-IDF-ranked attentive
genes from CellMeSH and recorded the placement when the predicted cell type
matched the actual annotation (called a “hit”). To interpret our model’s
mistakes, we also retrieved the first prediction from CellMeSH, which we
include as “Top Retrieval”.

Unsupervised label generation: Given that the cell types are not known a
priori, labels must be generated before training the DL core. To do so, we
choose to perform unsupervised clustering using the Leiden algorithm (Traag et
al., 2019), a standard scRNAseq clustering technique in many pipelines
(Heydari & Sindi, 2022), allowing label generation without supervision. All
results shown in this section follow this approach. However, given the
requisite comparison against the true annotations for this work, we ensure the
number of clusters generated by our model is equal to the number of annotated
populations from the original publication, without enforcing any additional
constraints.

Finding salient genes: To identify _attentive_ (salient) genes in each
population, we leverage the attention scores calculated by our model: First,
we calculate the mean attention score per cell type for each gene. To analyze
the performance of our model in assigning importance to various genes, we
investigate the correlation between the mean attention scores for all genes in
each cell type (Fig. 2(A),(C)). We find a low correlation for those
populations that are not closely related (e.g. dendritic and endothelial cells
in Fig. 2(A)), and a high correlation between those that are related (e.g.
CD4+ and CD8+ T cells in Fig. 2(C)). These results indicate that N-ACT has
accurately learned to assign importance to the same gene sets across similar
populations, as desired.

Next, we use mean attention scores per cluster and selected the top 100 genes
with the highest average values, constructing an object that includes the top
100 genes for each cell per cluster. Using this object, we calculate term
frequency (TF)-inverse document frequency (IDF) for each gene in the object.
TF-IDF is a standard natural language processing technique that weights the
importance of a word in a document (see Appendix C for more information). In
this formulation, each row of the object (containing top genes for all cells
in a cluster) is a document used to calculate TF and IDF. The TF-IDF values
for each gene were then multiplied by the original attention values, providing
the final saliency scores. TF-IDF normalization down-weights common
housekeeping genes that frequently appear in each cell population but are less
useful in identifying cell-types. Lastly, we re-rank genes based on these
weighted scores and select the top 25 as the attentive genes, which we use for
cell-type identification.

Identifying cell types: Once attentive genes are identified, various
techniques for finding their corresponding cell-types can be applied. To
showcase the accuracy and automation capabilities of our model, we queried the
attentive genes using CellMeSH (Mao et al., 2021), a probabilistic cell type
querying tool that uses a database built from indexed literature to map marker
genes to probable cell types (CellMeSH details provided in Appendix D). It is
important to note that the bias in each database can affect results, and that
other specialized databases or methods can further improve the identification
process (see Appendix G). Fig. 2(B) and 2(D) present the accuracy of cell-type
prediction using N-ACT-identified salient genes for Immune cSCC and for COVID
PBMC, respectively. These results demonstrate that N-ACT accurately identifies
attentive genes that are known markers for the underlying populations without
any prior knowledge of the system or species.

##  5 Conclusions and Discussion

In this work, we presented N-ACT, the first-of-its-kind interpretable DL model
for ACTI. We show that N-ACT effectively identifies cell-types, in a
supervised and, more importantly, unsupervised manner. N-ACT is a first
attempt at providing interpretability in this context, and we believe our
improvements and developments reduce subjectivity while significantly
minimizing the time needed for annotating scRNAseq datasets. Optimizing the
scRNAseq annotation process will accelerate translational and basic research
by enabling scientists to focus on the underlying biological questions. Our
results demonstrate that N-ACT accurately identifies salient genes that are
known markers for the underlying populations, without prior knowledge of the
system or species. Moreover, the interpretability of our framework is useful
for predicting the correct cell-types, and for better understanding the data
when there is ambiguity (e.g. Appendix G) or when the model makes mistakes. As
such, N-ACT provides a powerful tool for facilitating discovery, even if its
top prediction is ultimately incorrect. Using our model, cells can be assigned
cell-types by new users without prior expertise in the given system, within
minutes. From these conclusions, we hypothesize that attention can be further
utilized to identify unique relationships between different genes and cells,
which would not otherwise be apparent. Despite successful application of DL in
scRNAseq space, most DL models are not interpretable. Biologically-
interpretable DL models, such as N-ACT, can provide crucial information on the
algorithm’s decision making, while assisting scientists in understanding
underlying complex biological networks.

## Code and Data Accessibility

All source code and reproducibility/tutorial notebooks, alongside download
links to trained models and datasets, are available at
https://github.com/SindiLab/NACT.

## Acknowledgements

The authors were supported from the National Institutes of Health
(R15-HL146779 & R01-GM126548), National Science Foundation (NSF) (DMS-1840265)
and University of California (UC) Office of the President and UC Merced
COVID-19 Seed Grant. Computational resources were supported by the NSF, Grant
No. ACI-2019144 and in part through NSF awards CNS-1730158, ACI-1540112 and
ACI-1541349. Due to space constraints, complementary citations are included in
the Appendix.

## References

  * Bahdanau et al. (2015) Bahdanau, D., Cho, K., and Bengio, Y.  Neural machine translation by jointly learning to align and translate.  In Bengio, Y. and LeCun, Y. (eds.), _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_ , 2015.  URL http://arxiv.org/abs/1409.0473.
  * Chaudhari et al. (2021) Chaudhari, S., Mithal, V., Polatkan, G., and Ramanath, R.  An attentive survey of attention models.  _ACM Trans. Intell. Syst. Technol._ , 12(5), oct 2021.  ISSN 2157-6904.  doi: 10.1145/3465055.  URL https://doi.org/10.1145/3465055.
  * Clarke et al. (2021) Clarke, Z. A., Andrews, T. S., Atif, J., Pouyabahar, D., Innes, B. T., MacParland, S. A., and Bader, G. D.  Tutorial: guidelines for annotating single-cell transcriptomic maps using automated and manual methods.  _Nature Protocols_ , 16(6):2749–2764, 2021.  doi: 10.1038/s41596-021-00534-0.  URL https://doi.org/10.1038/s41596-021-00534-0.
  * Erfanian et al. (2022) Erfanian, N., Heydari, A. A., Iañez, P., Derakhshani, A., Ghasemigol, M., Farahpour, M., Nasseri, S., Safarpour, H., and Sahebkar, A.  Deep learning applications in single-cell omics data analysis.  _bioRxiv_ , 2022.  doi: 10.1101/2021.11.26.470166.  URL https://doi.org/10.1101/2021.11.26.470166.
  * Heydari & Sindi (2022) Heydari, A. A. and Sindi, S. S.  Deep learning in spatial transcriptomics: Learning from the next next-generation sequencing.  _bioRxiv_ , 2022.  doi: 10.1101/2022.02.28.482392.  URL https://doi.org/10.1101/2021.11.26.470166.
  * Heydari et al. (2022) Heydari, A. A., Davalos, O. A., Zhao, L., Hoyer, K. K., and Sindi, S. S.  ACTIVA: realistic single-cell RNA-seq generation with automatic cell-type identification using introspective variational autoencoders.  _Bioinformatics_ , 38(8):2194–2201, 02 2022.  ISSN 1367-4803.  doi: 10.1093/bioinformatics/btac095.  URL https://doi.org/10.1093/bioinformatics/btac095.
  * Kingma & Ba (2014) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization, 2014.  URL https://arxiv.org/abs/1412.6980.
  * Ma & Pellegrini (2019) Ma, F. and Pellegrini, M.  ACTINN: automated identification of cell types in single cell RNA sequencing.  _Bioinformatics_ , 36(2):533–538, 07 2019.  ISSN 1367-4803.  doi: 10.1093/bioinformatics/btz592.  URL https://doi.org/10.1093/bioinformatics/btz592.
  * Mao et al. (2021) Mao, S., Zhang, Y., Seelig, G., and Kannan, S.  CellMeSH: probabilistic cell-type identification using indexed literature.  _Bioinformatics_ , 38(5):1393–1402, 12 2021.  ISSN 1367-4803.  doi: 10.1093/bioinformatics/btab834.  URL https://doi.org/10.1093/bioinformatics/btab834.
  * Nair & Hinton (2010) Nair, V. and Hinton, G. E.  Rectified linear units improve restricted boltzmann machines.  In _Proceedings of the 27th International Conference on International Conference on Machine Learning_ , ICML’10, pp. 807–814, Madison, WI, USA, 2010. Omnipress.  ISBN 9781605589077.
  * Pasquini et al. (2021) Pasquini, G., Rojo Arias, J. E., Schäfer, P., and Busskamp, V.  Automated methods for cell type annotation on scrna-seq data.  _Computational and Structural Biotechnology Journal_ , 19:961–969, 2021.  ISSN 2001-0370.  doi: https://doi.org/10.1016/j.csbj.2021.01.015.  URL https://www.sciencedirect.com/science/article/pii/S2001037021000192.
  * Raffel & Ellis (2015) Raffel, C. and Ellis, D. P. W.  Feed-forward networks with attention can solve some long-term memory problems, 2015.  URL https://arxiv.org/abs/1512.08756.
  * Traag et al. (2019) Traag, V. A., Waltman, L., and van Eck, N. J.  From louvain to leiden: guaranteeing well-connected communities.  _Scientific Reports_ , 9(1):5233, 2019.  doi: 10.1038/s41598-019-41695-z.  URL https://doi.org/10.1038/s41598-019-41695-z.
  * Vaswani et al. (2017) Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, L. u., and Polosukhin, I.  Attention is all you need.  In Guyon, I., Luxburg, U. V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems_ , volume 30. Curran Associates, Inc., 2017.  URL https://proceedings.neurips.cc/paper/2017/file/3f5ee243547dee91fbd053c1c4a845aa-Paper.pdf.
  * Xu et al. (2015) Xu, B., Wang, N., Chen, T., and Li, M.  Empirical evaluation of rectified activations in convolutional network.  _arXiv preprint arXiv:1505.00853_ , 2015.

Appendix and Supplementary Material

##  Appendix Appendix A Computational Environment

Model development and testing was performed in Python (v 3.9.7) and data
processing was performed in R (v. 4.1.2) (detailed in section B). Our models
were developed and tested on A100 Nvidia GPUs, and data pre-processing on a
14-inch Macbook Pro with Apple M1 Max and 64 GB RAM. Data I/O was done in
Scanpy (v. 1.7.0) scanpy. The DL core of our model was developed in PyTorch
(v. 1.9.1); however, developing/testing N-ACT on A100 GPUs required
installation of a specific version of PyTorch (torch==1.9.0+cu111), which is
provided in N-ACT’s package repository. A complete list of requirements of
Python packages is also available in N-ACT’s GitHub repository:
https://github.com/SindiLab/NACT.

##  Appendix Appendix B  Datasets Studied and Pre-Processing Workflow

###  B. 1 Pre-Processing

All data (count matrices and manual annotations) are publicly available from
NCBI gene expression omnibus (GEO) and the Broad Institute Single Cell Portal
(SCP), with links provided below. Datasets were processed using the Seurat
package (v. 4.1.0) in R Seurat. Manual annotations were merged with count
matrices using a variety of tidyverse (v. 1.3.1) functions, and subsequently
added to Seurat object as metadata. Data filtering consisted of the standard
practice of removing cells with fewer than 200 expressed genes and removing
genes present in fewer than 3 cells. Next, we retained cells with less than
10% mitochondrial reads to mitigate cellular debris. Lastly, cell types
containing less than 100 cells were removed and excluded from the dataset.

After filtering, we identified the top 5,000 highly variable genes (HVG’s) for
each dataset (HVG procedure is detailed in Appendix C. 3). To minimize the
biological and technical effects in each dataset based on patient and/or
biological conditions such as normal versus disease state, we utilized Harmony
(v 0.1.0) to perform integration Harmony. Resulting Uniform Manifold
Approximation and Projection (UMAP) UMAP plots for each dataset with original
annotations are provided in Fig. A.3. SeuratDisk (v. 0.0.0.9019) was used to
convert the Seurat data object into an AnnData object compatible with Scanpy.
To perform clustering and generate cell labels (in the unsupervised case), we
used Scanpy’s pipeline for clustering (consisting of dimensionality reduction
using principal component analysis (PCA), followed by Leiden clustering). As
mentioned in the main manuscript, we found Leiden resolutions that led to the
same number of clusters as the annotated populations in order to compare our
predictions to the ground truth labels.

###  B. 2 Datasets

To evaluate N-ACT capabilities, we chose four large datasets relevant to
immune researchers in light of the current pandemic. Three datasets are
comprised of human cells, and one dataset consists of mouse cells. The three
human datasets were chosen due to the distinct disease conditions being
evaluated in the original studies. These included two COVID viral infection
studies and human cSCC cancer study. We included a mouse dataset to
demonstrate our model can be effectively used on non-human and non-immune
datasets. All datasets were generated using the 10x Genomics platform (see
10XplatformReview, DLinSpatialTranscriptomicsReview-App for a review of next-
generation scRNAseq). A summary for each dataset is provided in subsequent
subsections (see Fig. A.3).

Figure A.3: A UMAP-reduced plot of cell populations present in each dataset
investigated in this work.

####  B. 2.1 Mouse-HDF

Mouse-HDF MouseHDF [SCP1361] consists of scRNAseq of aortic cells in mice
given a normal diet and mice given a high-fat diet (HFD), resulting in 24K
cells (24K cells after processing). The authors identified 27 clusters, for 10
different cell populations.

Download Link:
https://singlecell.broadinstitute.org/single_cell/study/SCP1361/single-cell-
transcriptome-analysis-reveals-cellular-heterogeneity-in-the-ascending-aorta-
of-normal-and-high-fat-diet-mice

####  B. 2.2 Immune-CSF

Immune-CSF ImmuneCSF [GSE163005] profiles single cells in cerebrospinal fluid
(CSF) of Neuro-COVID, non-inflammatory, autoimmune neurological diseases, and
viral encephalitis patients. Cells isolated from 31 patients: 8 COVID
patients, 9 non-inflammatory, 9 autoimmune, 5 viral encephalitis resulting in
a total 80K cells. After processing the dataset contains 70K cells with 15
populations.

Download Link: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE163005

####  B. 2.3 COVID PBMC

COVID-PBMC COVIDPBMC [GSE163005] consists of the evaluation of the
transcriptional immune dysfunction triggered during moderate and severe
COVID-19 patients using scRNA-seq. Peripheral blood mononuclear cells (PBMC)
were isolated from 20 patients and were sequenced. Patients included in the
study ranged from healthy ($n=3$), moderate COVID ($n=5$), acute respiratory
distress syndrome (ARDS-Severe, $n=6$), and recovering ARDS-Recovering
($n=6$), resulting in 69K cells. After pre-processing the dataset, we retained
64K cells. The authors identified 9 populations in the dataset.

Download Link: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE154567

####  B. 2.4 Immune-cSCC

Immune-cSCC ImmunecSCC [GSE144236] evaluates scRNAseq of normal skin and
patient cutaneous squamous cell carcinoma (cSCC) tumors. Patients (n=10) with
the cSCC tumors had normal skin and tumor cells sequenced, resulting in 48K
cells. We retained 47K cells after pre-processing. The authors identified 7
major cell populations. The myeloid cell population (CD14+Hi) is composed of
various subpopulations.

Download Link: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE144236

##  Appendix Appendix C : Definition of Standard Mathematical Expressions and
Evaluation Metrics

###  C. 1 Cross Entropy Loss

The downstream task was to predict the correct labels (original labels or
generated labels), with the learning objective being a standard cross-entropy
loss:

$H(\mathbf{y},\hat{\mathbf{y}})=\sum_{j=1}^{k}y_{j}\log(\hat{y}_{j})+(1-y_{j})\log(1-\hat{y}_{j}),$
(A.2)

where $\mathbf{y}$ and $\hat{\mathbf{y}}$ denote the original/generated labels
and model predicted labels, respectively. In this formulation, $\mathbf{y}$ is
a one-hot encoded vector of the cell types; that is

$y_{j}=\begin{cases}1\hskip 7.22743pt\text{if cell type is }j\\\ 0\hskip
7.22743pt\text{otherwise.}\end{cases}$

###  C. 2 Formulation of TF-IDF for Cell Type Querying

In this work, we leverage attention scores generated for each gene in all
clusters. First, we calculate the average attention score for each gene per
cluster. Next, take the top 100 genes with the highest attention scores per
cluster and create an object containing only gene names. Using the top 100
genes per cluster object, we calculate the term frequency(TF)-inverse document
frequency(IDF) for each gene in the matrix. TF-IDF is a standard tool in
natural language processing (NLP) often used to weight the importance of words
appearing in documents, as shown in (A.3).

$\displaystyle TF(g,P)=\frac{f_{P}(g)}{|P|}$ (A.3a)

$\displaystyle IDF(g,D)=-\log p(g|D)=\log\left(\frac{|D|}{|\\{P\in D,g\in
P\\}|}\right)$ (A.3b)

$\displaystyle TF\text{-}IDF(g,P,D)=TF(g,P)\cdot IDF(g,D)$ (A.3c)

where $f_{P}$ is the raw count of term (gene) $g$ in the document (population)
$P$, with a corpus of documents $D$. We consider each row in the matrix to be
a document in which we calculate gene frequencies, and then calculate the IDF,
which down-weights more common genes in the matrix. We multiply both TF and
IDF for each gene to obtain the TF-IDF score. The generated TF-IDF score down-
weights common genes in the top 100 by multiplying the TF-IDF scores by the
average attention scores. Once the attention scores have been weighted, the
top 25 genes per cluster are submitted to CellMeSH (CITE: PMID: 34893819),
which generates a prediction cell type.

###  C. 3 Selecting Highly Variable Genes

To select Highly Variable Genes (HVGs), we utilized Seurat’s
FindVariableFeatures() function, described in HVG_V3. In this step, the aim is
to calculate a measure of single-cell dispersion while including the mean
expression. Step one learns a mean-variance relationship from the data, which
is then used to calculate an expected standard deviation for each gene. To do
so, the mean $\mu_{g}$ and variance $\sigma_{g}$ of each gene is computed from
raw data. Next, a curve fitting using a $2^{\text{nd}}$ degree polynomial is
applied to learn $f(\mu_{g})=\hat{\sigma}_{g}$, providing a regularized
estimator of variance given the mean of a gene $g$. Using this, we perform
feature transformations without removing higher-than-expected variations. That
is, given the raw counts $X_{ig}$ for gene $g$ in cell $i$, the mean raw value
for gene $g$, $\mu_{g}$, and $\sigma_{g}$, the expected standard deviation of
gene $g$ from the global fit, we have the transformed count value of gene $g$
as:

$z_{ig}=\frac{X_{ig}-\mu_{g}}{\sigma_{g}},$ (A.4)

for all gene $g$ across all cell $i$. Finally, variance across the new
standard feature values is calculated and used to rank the genes. Although the
convention is to select 2,000 HVGs, we chose to select 5,000 HVGs to increase
the complexity and truly test model interpretability.

###  C. 4 Hit@k

Cell type retrieval position is measured using N-ACT attentive genes for each
cell population using _Hit@k_. Hit@k is a metric measuring retrieval of a
target value among the top $k$ retrievals. In our model, a “hit” is the
published cell type annotation (ground truth label) existing in the database.
Rankings for $k=\\{1,3,5,10\\}$ are reported in Fig. 2, and in Fig. A.4.

###  C. 5 F1 Score

F1 score is a standard metric for evaluating a classifier. F1 score is the
harmonic mean of precision and recall, which is shown in Eq. (A.5),

$F1=2\left(\frac{\text{precision}\cdot\text{recall}}{\text{precision}+\text{recall}}\right).$
(A.5)

For more information on weighted and non-weighted F1 score, see
https://scikit-
learn.org/stable/modules/generated/sklearn.metrics.f1_score.html.

##  Appendix Appendix D CellMeSH Cell Type Predictions

CellMeSH CellMesh is a probabilistic method that generates a cell type
prediction based on a gene list. CellMeSH aims to alleviate two common issues
with querying gene lists from the literature: 1) publication bias, which
relates to specific genes/cell types that are studied more often than other
types, and therefore have more literature associated with them; and 2) Noise
present in gene/cell-type mapping and the corresponding database. CellMesh
uses a database built from the National Library of Medicine (NLM) MEDLINE
indexed records called Medical Subject Headings. In building the CellMeSH
database, Mao et al. use TF-IDF to reduce publication bias thus addressing
publication bias.

To address the second issue, CellMeSH utilizes a probabilistic querying
method. Given gene query list $Q$, CellMeSH assumes a probabilistic model for
a query genes $g$ being obtained from a cell-type $C$, shown in Eq. (A.6):

$P(g|C)=\begin{cases}\alpha\cdot w_{C}(g)\hskip 14.45377pt&g\in Q\cap C\\\
(1-\alpha)\cdot\frac{1}{N_{g}-K_{C}}\hskip 14.45377pt&g\in
Q\cap\overline{C}\end{cases}$ (A.6)

with $w_{C}(g)$ being the adjusted weight of gene $g$ in cell type $C$ (using
TF-IDF), $N_{g},K_{C}$ denoting the total number of genes and total number
genes with non-zero weight in $C$. $\alpha$ is a parameter which aims to help
with noise in the database. For each candidate cell type $C$, CellMeSH
calculates a log-likelihood score:

$L(Q|C)=\log P(Q|C)=\sum_{g}\log P(g|C),$ (A.7)

which are then used to rank $C$ based on their values. Lastly, CellMeSH uses a
maximum likelihood-based estimation to generate predictions. In this
framework, the top cell type, $C^{*}$, can be found as:

$C^{*}=\text{argmax}_{C}\log L(Q|C).$ (A.8)

##  Appendix Appendix E : Main Manuscript Results (Larger Figures)

Figure A.4: (For the ease of the reader) An enlarged version of Hit@K results,
initially presented in Fig. 2 of the main manuscript. Figure A.5: (For the
ease of the reader) An enlarged version of correlation results presented in
Fig. 2 of the main manuscript.

##  Appendix Appendix F N-ACT Architecture Details

Figure A.6: A Detailed View of N-ACT and the Multi-Headed Projection Module.
We dive deeper into the projection block of N-ACT in Appendix F. An
implementation view of the DL architecture is presented in Fig. A.7.

###  F. 1 Projection Block

The idea behind the N-ACT projection block is to learn various representations
for different gene subsets in each cell. Projection block design was inspired
by the multi-head attention architecture presented in AttentionIsAllYouNeed-
App; however, the projections are not multi-head attention mechanisms. One way
to think about the multi-head projection block is to view it as a set of $h$
linear projections, with each $l_{h}:\mathbb{R}^{B\times
N}\to\mathbb{R}^{B\times d}$ ($B$ is the number of samples and $N$ is the
number of genes) done sequentially and independent of one-another (e.g. in a
for loop). However, as noted by AttentionIsAllYouNeed-App, these projections
can be done more efficiently through creating a tensor
$L\in\mathbb{R}^{B\times h\times d}$, which acts the same way as the
collection of the individual linear operators. In this formulation, $0\equiv
N(\text{mod }h)$, and the last projection component, the “concatenation” layer
(depicted in Fig. A.6) allows us to reshape the model output which can be
passed along to remaining layers.

To increase model capacity and allow N-ACT to learn complex mappings, outputs
are non-linearly “activated” through a Point-Wise Feed Forward Neural Network
(which can be thought of as 1$\times$ 1 convolution). However, there is the
possibility of projection blocks learning a non-linear mapping unrelated to
interpretability, resulting in loss of interpretability. Therefore, we
hypothesized that adding residual connection between output of the attention
module and input of each subsequent layer would improve performance and
interpretability. Our ablation studies show that adding a residual connection
improves model performance (Table A.3). Skip connections are added to
projection block outputs and normalized using Layer-Norm LayerNorm) before
being used as input for the next layer.

Figure A.7: An Implementation View of the N-ACT DL Core. In this illustration,
$C$ corresponds to the number of cells inputted, $N$ denotes the number of
features ($N=5000$ in our results), $h$ is the number of projection heads, $k$
is the number of projection blocks(in our model, $h=10$ and $k=2$) and $T$
denotes the number of classes (distinct labels) present in the data.

###  F. 2 Multi-Head Ablation Study

We investigated the effect of head number in the multi-head projection block
and found that 10 heads provided the best results (Table A.2). We also tested
our model accuracy when using one, two and three projection blocks and found
that two projection blocks provided the best balance of accuracy and
efficiency. To test our hypothesis regarding residual connections, we tested
N-ACT (10 heads, 2 projection blocks) with and without skip connections. As
shown in Table A.3, model accuracy drops without the residual connections.
Residual connections are identity mappings added to the output of each layer,
with this quantity normalized and used as inputs for subsequent layers.

Table A.2: Ablation Study on the Number of Projection Heads. We fixed all other hyperparameters and studied the effect of head number on accuracy and interpretability. Training times are the average of 5 training runs on an A100 GPU. Accuracy of each model remained the same across different training settings, since all random parameters were initialized with the same random seed. Number of Heads | W-F1 | NW-F1 | Hit@5 | Avg. Training Time
---|---|---|---|---
1 | 0.9278 | 0.8968 | 0.85 | $9.11\pm 0.14\text{ (Min)}$
5 | 0.9317 | 0.9077 | 0.85 | $9.38\pm 0.31\text{ (Min)}$
8 | 0.9307 | 0.9156 | 1.00 | $9.31\pm 0.22\text{ (Min)}$
10 | 0.9322 | 0.9173 | 1.00 | $9.56\pm 0.27\text{ (Min)}$
20 | 0.9324 | 0.9156 | 1.00 | $9.87\pm 0.24\text{ (Min)}$
Table A.3: Effectiveness of Residual Connections on N-ACT. Residual Connection? | W-F1 | NW-F1
---|---|---
Yes | 0.9322 | 0.9173
No | 0.9243 | 0.8725

###  F. 3 Training Scheme

As mentioned in the main manuscript, we train N-ACT for 50 epochs using Adam
Adam-App optimizer, with a fixed learning rate. We tested training N-ACT for
more epochs and found that the accuracy of predictions increases; however, a
learning rate scheduling is beneficial in avoiding overfitting (when training
over 100 epochs). When training for longer epochs, we employed an exponential
learning rate decay with $\gamma=0.95$ and a decay schedule of every 10
epochs. We chose 50 epochs to balance accuracy and training time efficiency.

###  F. 4 Supervised ACTI Comparison

Table A.4: Benchmarking N-ACT on Supervised ACTI. Although the main goal of N-ACT is interpretable ACTI in an unsupervised manner, our model can be used in a supervised setting as well, while still providing biological interpretability. The reported training times are the median of 5 runs for each model. The same random seed was used to initialize all random parameters to ensure reproducibility of results across different runs. Training/Testing Support: the number of samples used in training/evaluation, W-F1: Weighted F1 score, NW-F1: Non-Weighted F1 score Model | W-F1 | NW-F1 | Training Support | Testing Support | Med. Training Time
---|---|---|---|---|---
|  |  | Mouse HDF |  |
ACTINN | 0.9703 | 0.9677 | 21,003 | 2,998 | 1.5 Minutes
N-ACT (Ours) | 0.9681 | 0.9712 | 21,003 | 2,998 | 5.4 Minutes
|  |  | Immune CSF |  |
ACTINN | 0.9357 | 0.8898 | 66,549 | 8,991 | 3.8 Minutes
N-ACT (Ours) | 0.9285 | 0.8963 | 66,549 | 8,991 | 11.2 Minutes
|  |  | COVID PBMC |  |
ACTINN | 0.9323 | 0.9148 | 55,005 | 9,864 | 3.1 Minutes
N-ACT (Ours) | 0.9322 | 0.9173 | 55,005 | 9,864 | 9.4 Minutes
|  |  | Immune cSCC |  |
ACTINN | 0.9646 | 0.9315 | 40,027 | 6,994 | 2.4 Minutes
N-ACT (Ours) | 0.9684 | 0.9455 | 40,027 | 6,994 | 7.7 Minutes

##  Appendix Appendix G  Utility of N-ACT for Disambiguation of Broad
Annotations

Figure A.8: N-ACT Used to Disambiguate Broad Annotations. Here, we use N-ACT-
identified salient genes to query a specialized database (Azimuth) to
disambiguate original broad manual annotations.

Given that many manual annotations are typically performed with only a few
genes (often two or three BestPracticesTheis), it is possible to have
populations that are broad and ambiguous. This was the case with two
populations of COVID PBMC data COVIDPBMC, namely the original annotations
“Proliferating Lymphocytes” and “ Unidentified Lymphocytes” (Fig. A.8). Here
we utilize N-ACT to disambiguate these broad annotations without re-clustering
or further complex analyses.

In the results shown in the main manuscript, we utilized CellMesh CellMesh.
However, using specialized databases could provide refined predictions. To
demonstrate this, we investigated the two broad annotations in the COVID PBMC
data and queried attentive genes from the Azimuth Cell Type 2021 database
tailored for immune cells Azimuth [using Enrichr R Package (v 3.0) Enrichr1,
Enrichr2, Enrichr3], to perform an enrichment analysis. We query our top 25
genes against the Azimuth Cell Types 2021 reference for the enrichment
analysis and select the top 25 attentive genes (as described in 4.2 and
Appendix C) for “proliferating lymphocytes” and “unidentified lymphocytes.”

Our enrichment analysis for “proliferating lymphocytes” resulted in multiple
significant (significant adjusted $p-$values) populations, with the most
probable types being CD8+ proliferating T cells (shown in A.8). Intuitively,
these results are expected given the quantitative and qualitative similarity
of this population with the CD8+ T cell population. It is important to note
that the difference in gene overlap percentages of top three predictions,
namely CD8 proliferating T, CD4 proliferating T, and proliferating natural
killer populations is very small, which may explain why the original
annotations were left broadly as proliferating lymphocytes. Enrichment
analysis for “unidentified lymphocytes” yielded natural killer cells as the
most probable cell type, with other viable populations also being
statistically significant. Similar to the previous population, our results are
possibly intuitive given the closeness of the proliferating lymphocytes to
CD8+ T cells and natural killer cells. Additionally, there is known overlap in
the gene sets between of CD8+ T cells and natural killer cells. Lastly, we
note that the second most probable population based on attentive genes are
CD8+ effector memory T, suggesting that the “unidentified lymphocyte”
population could likely be refined into two or more populations. These results
further signify the utility of N-ACT for unsupervised annotation, and the
applicability of our framework in tandem with other annotation forms to
provide interpretability and validation.

StyleFiles/icml2022 A


# Validating an SVM-based neonatal seizure detection algorithm for
generalizability, non-inferiority and clinical efficacy

Karoliina T. Tapani<EMAIL_ADDRESS>Päivi Nevalainen  Sampsa
Vanhatalo  Nathan J. Stevenson  BABA Center, Helsinki University Hospital,
University of Helsinki, Finland  Department of Neuroscience and Biomedical
Engineering, Aalto University, Espoo, Finland1 Brain Modelling Group, QIMR
Berghofer Medical Research Institute, Australia  Department of Medical
Physics, Kymenlaakso Central Hospital, Kotka, Finland  Epilepsia Helsinki and
Department of Clinical Neurophysiology, HUS diagnostic center, Helsinki
University Hospital and University of Helsinki, Finland

###### Abstract

Neonatal seizure detection algorithms (SDA) are approaching the benchmark of
human expert annotation. Measures of algorithm generalizability and non-
inferiority as well as measures of clinical efficacy are needed to assess the
full scope of neonatal SDA performance. We validated our neonatal SDA on an
independent data set of 28 neonates. Generalizability was tested by comparing
the performance of the original training set (cross-validation) to its
performance on the validation set. Non-inferiority was tested by assessing
inter-observer agreement between combinations of SDA and two human expert
annotations. Clinical efficacy was tested by comparing how the SDA and human
experts quantified seizure burden and identified clinically significant
periods of seizure activity in the EEG. Algorithm performance was consistent
between training and validation sets with no significant worsening in AUC
($p>0.05$, $n=28$). SDA output was inferior to the annotation of the human
expert, however, re-training with an increased diversity of data resulted in
non-inferior performance ($\Delta\kappa$=0.077, 95% CI: -0.002-0.232, $n=18$).
The SDA assessment of seizure burden had an accuracy ranging from 89-93%, and
87% for identifying periods of clinical interest. The proposed SDA is
approaching human equivalence and provides a clinically relevant
interpretation of the EEG.

###### keywords:

neonatal EEG , EEG monitoring , neonatal intensive care unit , seizure ,
support vector machine

††journal: Computers in Biology and Medicinemyfootnotemyfootnotefootnotetext:
Department of Neuroscience and Biomedical Engineering Aalto University School
of Science P.O. Box 12200 FI-00076 AALTO, Finland

##  1 Introduction

Clinical suspicion of neonatal seizures or an identified risk of seizures due
to e.g. hypoxic–ischaemic encephalopathy (HIE), meningitis or stroke are key
indications for electroencephalography (EEG) monitoring in the neonatal
intensive care units [1]. To guide interventions and inform prognosis,
reliable and accurate detection of seizures is essential [2, 3]. Recent
guidelines insist that neonatal seizures must have an EEG correlate as
clinical recognition of seizures is confounded by a range of factors [4].
Neonatal EEG monitoring generates a considerable amount of data and its
interpretation is a challenge that has led to the development of several
automated seizure detection algorithms (SDA) [5, 6, 7, 8, 9, 10, 11, 12, 13,
14, 15, 16, 17, 18].

While the performance of many of these algorithms has thoroughly been assessed
within the training data, few of them have been validated on an independent
data set [19, 8, 20, 21, 22, 17]. This step is essential as it ensures that
SDAs built on relatively small data sets of neonatal EEG generalize to the
larger population of neonatal EEG recordings. While generalizability is
important, it alone does not suffice to define SDA performance. An SDA must
also be accurate and clinically useful. The complete assessment of SDAs,
therefore, must address four important questions: 1) does the SDA performance
generalize?; 2) does the SDA generate an annotation that is non-inferior to
the human expert, taking into account ambiguity in the annotations of human
experts?; 3) does the annotation of the SDA accurately reflect the seizure
burden experienced by a neonate over a clinically relevant period of time
(e.g. birth to 72-96 h post-natal age)?; 4) will the implementation of the SDA
improve clinical practice?

The SDA developed by Temko et al. (2011) has the most comprehensive assessment
to date [23]. Promising early performance was validated on an independent data
set and the annotation of the SDA was shown to provide an accurate reflection
of seizure burden [22]. Non-inferiority of the SDA annotation to the
annotation of the human expert was, however, not comprehensively answered and
the performance of the SDA when clinically trialled was underwhelming [24].
Recent advances in SDAs and subsequent improvements in accuracy, hold out hope
that modern SDAs will soon achieve a level of performance that will,
ultimately, improve clinical practice.

In this paper we 1) assess the generalizability of our previously developed
SDA on an unseen validation data set; 2) determine if the SDA generates an
annotation that is non-inferior to the human expert; and 3) determine if the
annotation of the SDA accurately estimates the seizure burden on a neonate.
While we do not have access to the implementation of the SDA of Temko et al.
(2011), we will use our implementation of their feature set, trained on our
data, as a benchmark.

##  2 Materials and Methods

This study utilizes two cohorts of EEG recordings. The first data set
(training data set) is used to develop and train the neonatal SDA. The second
data set (validation data set) is used to validate the SDA. To test the
ability of the SDA to generalize to unseen data, several performance measures
evaluating differences between the SDA and human annotation of the EEG, were
compared between the training and validation data sets. The smaller the
difference in these performance measures the more generalized the SDA. The
non-inferiority of the SDA annotation compared to the human annotation of the
validation set was examined using measures of inter-observer agreement (IOA).
Clinically relevant measures of seizure burden were compared between the SDA
and human expert annotations to assess the clinical efficacy of the algorithm.
Finally, the validation data set was split in two, with the first third
($n=10$) added to the training set and the latter two thirds ($n=18$) saved as
a final test set. This process helps determining if more diverse training data
allows improving SDA generalizability or non-inferiority.

###  2.1 EEG acquisition and annotation

Training data set: This data set consists of EEG recordings from a cohort of
79 neonates with a variety of aetiologies. These neonates were recorded with a
short duration paradigm and a whole scalp, 10-20 EEG system between 2010 and
2014. The EEG was recorded using a Nicolet One EEG monitor (Nicolet One, Natus
Medical Incorporated, USA) with 19 scalp electrodes placed on positions Fp1,
Fp2, F3, F4, F7, F8, T3, T3, C3, C4, T5, T6, P3, P4, O1, O2, Fz, Cz, Pz with
sampling rate of 256Hz. The median duration of each recording was 74 min (IQR:
64 to 96 min). A bipolar derivation (double banana) was annotated by three
clinical experts and used to train the SDA. The data set contains 460 seizures
from 39 neonates (in consensus). More details can be found in [25].

Validation data set: This data set consists of EEG recordings from a cohort of
28 neonates with stroke or HIE (Table 1). This validation set is a subset of
the cohort first studied in [26].

Table 1:  Table of patient characteristics for the validation data set ($n=28$). Data shown as $n$ (%), or mean (standard deviation). Patient characteristics ($n=28$) |
---|---
Gestation age (weeks) | $39.2(\pm 2)$
Birth weight (g) | $3300(\pm 550)$
pH | $7.2(\pm 0.11)$
Etiologies |
Stroke | 25 (89%)
Hypoxic-ischemic encephalopathy | 3 (11%)

These neonates were recorded with long duration brain monitoring between 2011
and 2016. The signal was acquired at 250Hz using a Nicolet One EEG monitor
(Nicolet One, Natus Medical Incorporated, USA) and four needle electrodes
placed on positions F3, F4, P3, P4. The total length of all recordings was 86
days 16 hours (median 77h 50min, IQR: 48h 20min - 94h 35min). A three channel,
bipolar derivation (F3-P3, F4-P4, P3-P4) was annotated by two clinical experts
(PN, SV: both employed at Helsinki University Hospital) and processed by the
SDA. The start and stop times of each seizure were noted. For annotations, the
internationally approved neonatal seizure definition by Clancy was used [27].
Clinical experts were blinded to each others annotations.

The training and validation sets were both collected in the neonatal intensive
care unit (NICU) of the Children’s Hospital, Helsinki University Hospital,
Finland. The EEG was recorded as part of the standard care in our NICU. The
use of the training recordings were approved by the Ethics Committee of the
Helsinki University Children’s Hospital, Finland. The Institutional Research
Review Board at Helsinki Children’s Hospital approved the extraction and
collation of the validation data set from hospital records and waived consent
due to the study’s retrospective and observational nature.

###  2.2 SDA

This SDA is based on a ’bag of features’ calculated on 16s epochs of EEG. The
combination of features to form a decision statistic was performed with a
support vector machine (SVM) [15]. We used 22 summary EEG measures/features to
represent a 16 s epoch of EEG. One additional feature to the previously
published algorithm was maximum amplitude, which we used to discard epochs
with high amplitude. To improve the robustness of the SDA, we incorporated an
outlier detection system into our SDA. This was due to the relatively small
size of our database and the expected difference in EEG quality between high
density, short duration training data and low density, long duration
validation data. Therefore we used a modified version of our previous system
that was first proposed in [15, 28]. Post-processing of the initial SDA output
included eliminating outliers (we also eliminated EEG epochs detected as ’bad
electrode’ by the EEG monitor), applying a temporal moving average to each
channel, taking the maximum value across channels and thresholding to form a
binary (seizure/no seizure) decision. The binary decision was then extended in
time using a collar. The parameters for outlier determination and seizure
detection were optimized on training data by maximizing the agreement between
the SDA and human annotation.

####  2.2.1 Outlier detection

Outlier detection identifies periods of the EEG recording that fall outside
its experience; in our case the span of the feature set calculated on the
training data. These methods are important to incorporate into clinical
decision algorithms as the full variety of EEG patterns may never be captured
in training data sets and when the detection of the rarer, clinically
important class (seizure) is critical. The implementation of outlier
detection, therefore, ensures that unusual EEG patterns (typically associated
with artefact or rare neurological conditions) which fall outside the span of
training data are dealt with differently. For our SDA, we adhere to the
principle of primum non nocere and annotate ’no seizure’ when an outlier is
detected.

The effect of outlier detection on seizure annotation is shown in Figure 1.
The example is based on a 2 feature SVM. In 1A, the original SVM shows the
separation between the two classes, ’seizure’ and ’no seizure’. In 1B the
black outlier boundary is included and is based on distances of the 3 nearest
neighbours within the training data. Any data falling outside the decision
boundary is labelled as ’no seizure’ resulting in a smaller decision space for
’seizure’.

Figure 1: The effect of outlier detection on SDA annotation. The blue dots
denote most of the training samples in the database. A) The original SVM based
SDA. The thick red line is the SVM decision boundary for determining ‘seizure’
or ‘no seizure’ and the yellow circles represent incoming, unseen data
(potential false detections). B) SVM based SDA with outlier detection. The
black boundary is based on kNN distances of the training data; the updated
decision boundary reduces the size of the feature space assigned to ’seizure’.

As not all features in our original feature set are best suited for outlier
detection, we use both general feature space distances, $k$-nearest neighbors
(kNN) search, and feature specific (EEG amplitude) distances. If an incoming
data point is either greater than a threshold Euclidean distance from the
k-nearest neighbours in the training data set or greater than a threshold EEG
amplitude, then an outlier is detected and no seizure annotated.

We used the knnsearch function in Matlab ([29]), which outputs ascending
Euclidean distances of an incoming data point to an array of training data
points. In the case of our SVM-based detector, this data point is represented
by a vector of features computed on a single epoch of EEG. The number of
distances measured, that is the number of neighbors, is defined by parameter
$k$.

We complemented our thresholding of feature space distances with a stricter,
separate threshold on EEG amplitude. The number of neighbours, distance
threshold for the kNN search, EEG amplitude threshold, decision collar and
decision threshold were, therefore, optimized on the training data set.

###  2.3 Analysis of algorithm performance

The most commonly used measure in SDA performance assessment is the area under
the curve of the receiver operator characteristics (AUC), a second by second
comparison of the algorithm annotation to the human annotation. It combines
two important measures, true positive rate (or sensitivity) and false
detection rate (or specificity) on multiple detection thresholds. These
measures can be defined with the numbers of true positives/seizure (TP), true
negatives/non-seizure (TN), false positives (FP) and false negatives (FN).
Sensitivity is defined by TP/(TP+FN) and specificity by TN/(TN+FP). In
addition to these temporal measures we assessed seizure detection rate (SDR),
that is the number of correctly detected seizure events, and false detections
per hour (FD/h), which is the number of falsely detected seizure events
divided by the duration of the recording in hours. These measures can be
supplemented with measures of inter-observer agreement (IOA) such as the Kappa
family of statistics; we use Cohen’s kappa as we only compare between two
annotations at a time [30].

We calculated performance measures on individual neonates and also on a
single, concatenated, consensus of human expert annotation for comparison. As
each data set contained neonates with no seizure, a concatenated annotation
avoids complications of per neonate assessments (sensitivity/AUC and SDR
cannot be calculated on neonates who do not have seizure). It also provides a
more useful assessment of event based measures as several recordings in the
training data set had limited events on relatively short recordings resulting
in poor quantization per neonate.

####  2.3.1 Implementation of Outlier Detection

The cross-validation results of the original SDA and SDA with outlier
detection were compared using statistical testing. Differences between AUC,
sensitivity and specificity calculated on a per neonate basis were tested
using Wilcoxon signed rank test (for paired data) due to non-normal
distribution of these values across training and validation cohorts. AUC and
sensitivity can only be calculated on neonates with seizure, while specificity
was calculated on all neonates in the training and validation cohorts. All
tests were two-sided and the level of significance was 0.05.

Analysis of concatenated recordings was also performed using bootstrapped
statistics. In this case, a concatenated annotation from each data set was
generated using a random sampling of neonate annotations (sampling with
replacement). The difference between performance measures were calculated on
each re-sampled annotation. The difference in performance measures was
calculated for 1000 re-samplings and used to estimate a 95% confidence
interval of differences. If the difference spanned zero, there was no
significant difference in performance between data sets. This statistical
method was used to compare performance where applicable (we cannot compare
bootstrapped distributions between different data sets).

As outlier parameters were optimized using IOA between the algorithm and human
experts, we used Cohen’s kappa as the primary measure of performance
differences between the SDA, with and without, outlier detection.

####  2.3.2 Generalizability of the SDA

We measured generalization by comparing perfmance measures between cross-
validation assessments of the training dataset and the validation dataset. If
there was no statistically significant difference in these measures between
the training and validation sets then we claim that the SDA generalizes to a
larger population of neonates. The primary measure of generalizability was
statistical testing of AUC.

Differences between AUC, sensitivity and specificity calculated on a per
neonate basis were tested similarly to the previous section, with the
exception that a Mann-Whitney U tests (for unpaired data) was used. If the
difference in AUCs is not significant, the SDA generalizes on unseen data.

####  2.3.3 Non-inferiority to human annotation

The limitation of AUC as a measure of performance is that it does not
inherently determine if the SDA annotation of the EEG is ’good enough’. This
is particularly valid in instances where there is ambiguity in the ’gold
standard’ used to calculate the AUC such as the manual annotation of the EEG
by human experts. We, therefore, determined the non-inferiority of the SDA to
the annotations of the human expert using bootstrapped measures of IOA. Each
expert and the SDA were compared to each other using Cohen’s kappa statistic.
The difference between the kappa of SDA/human expert pairs and human/human
pairs was estimated using a bootstrap (1000 iterations). If the difference
spanned zero then the SDA/human kappa was assumed to be equivalent/non-
inferior to the human/human kappa. This analysis uses the complete annotation
of each human expert and avoids the requirement of consensus annotations.

We also investigated SDA errors at the annotation and feature level. At the
annotation level, we compared general summary measures of the annotation such
as the total seizure burden, seizure number and seizure duration. At the
feature level, we determined if some features were systematically different
between the training and validation data sets.

####  2.3.4 Clinical efficacy

While non-inferiority provides an effective means to compare the annotation of
an SDA to the annotation of human experts, it does not evaluate the practical
clinical impact of any discrepancies. To this end, we determined the clinical
efficacy of SDA by comparing clinically relevant interpretations of SDA and
human annotations of the EEG. The first set of measures compares a continuous
assessment of seizure burden. This measure is considered clinically important
as it provides a clear quantitative target for treatment.

While some clinicians are aggressively treating any amount of seizures with a
variety of anti-epileptic medications [1, 31, 32], there is also substantial
variance among clinicians with respect to the level of seizure burden that is
considered to prompt treatment. We compared hourly assessments of seizure
burden using the correlation coefficient and a diagnostic binary decision of
periods of clinical targets and periods of no clinical target in an EEG
recording. We define a clinical target as a two-hour period with two 30s
seizures or 3 minutes of accumulated seizure [4]. We then compare the number
of clinical targets defined by the human annotation to the number defined by
SDA annotation using sensitivity and specificity.

We also evaluate the prognostic utility of the SDA. Recent work has shown that
high total seizure burden ($>$ 45 minutes) and high maximum hourly seizure
burden ($>$ 13 minutes per hour) in a long duration recording of EEG are
associated with poor neurodevelopmental outcome in neonates with HIE [33].
Note, that a high hourly seizure burden is typically used to define status
epilepticus ($>$30 minutes per hour). We use 13 minutes per hour as an
alternate definition of high seizure burden due to its association with poor
developmental outcome [33]. We compare the SDA annotation to the human
consensus annotation for the prediction of high or low seizure burden in each
neonate using measures of sensitivity and specificity. We also show how
differences in human interpretations affect these measures by averaging
sensitivity and specificity between experts, using expert 1 as the ’gold
standard’ and using expert 2 as the ’gold standard’.

###  2.4 Re-training

Finally, we investigated the incorporation of additional data from the
validation set into the training set on SDA performance. We extracted seizure
and non-seizure epochs from the first 10 neonates in the validation data set.
The epochs for each class were randomly sampled from consensus segments, so
that the size of this new training set was comparable to the size of the
original training set. The non-seizure epochs were randomly sampled from the
second half of the recordings only, as the majority of false detections
occurred in this time period (see Results). To increase diversity in the
training data set, not the total number of samples, every second sample from
the new and from the original training data was used to re-train the SVM.
Generalizability and non-inferiority measures were calculated on the re-
trained SVM to determine the effectiveness of added diversity in the training
data.

##  3 Results

###  3.1 Implementation of outlier detection and generalization of the SDA

The incorporation of an outlier detection significantly improved the agreement
between the SDA output and the consensus annotation of the human expert (Table
2). Slight, but significant reductions in per infant AUC were offset by
significant improvements in specificity and false detections per hour.

The SDA had consistent performance between training and validation sets with
our primary measure for generalizability, the AUC, suggesting that the
algorithm generalizes to unseen data (Table 2). That is, there was no
significant reduction in performance between the cross-validated training data
set and the validation data. The decrease in sensitivity was significant,
while there was no statistical difference for specificity.

The proposed SDA outperformed an SDA using the Temko feature set (Table 2).

Table 2:  The performance of a neonatal SDA with outlier detection. The neonatal SDA with outlier detection was applied to the original dataset (10-fold cross-validation) and compared to the original SDA without outlier detection. The neonatal SDA was also applied to an independent validation set and compared to an SDA based on the Temko feature set (Temko FS). Median (IQR) are presented for AUC, sensitivity (Sens) and specificity (Spec). The value with 95%CIs in brackets are presented for cAUC, cSDR, cFD/h and cKappa (estimated on concatenated recordings). a and b denote a significant increase and decrease, respectively, between the original SDA and the SDA with outlier detection. c denotes a significant increase in SDA performance between the training and validation data sets. d denotes significant decrease in performance between proposed SDA and SDA based on Temko features. Primary measures of performance are shown in bold. $n_{t}$ is the total number of neonates in each dataset and $n_{s}$ is the number of neonates with consensus seizure in each dataset. | Original Training | Original Training | Validation | Validation
---|---|---|---|---
| (10-fold CV) | (10-fold; outlier detection) | (outlier detection) | (Temko FS)
| $n_{t}=79,n_{s}=39$ | $n_{t}=79,n_{s}=39$ | $n_{t}=28,n_{s}=24$ | $n_{t}=28,n_{s}=24$
AUC | 0.992 (0.938-0.998) | 0.986 (0.908-0.995)b | 0.960 (0.896-0.980) | 0.918 (0.823-0.961)d
Sens | 0.823 (0.341-0.975) | 0.784 (0.363-0.958) | 0.605 (0.218-0.753)c | 0.385 (0.061-0.674)d
Spec | 0.993 (0.979-1) | 0.996 (0.982-1)a | 0.998 (0.997-0.999) | 0.998 (0.994-0.999)d
cAUC | 0.956 (0.931-0.979) | 0.955 (0.928-0.977) | 0.967 (0.953-0.974) | 0.949 (0.931-0.961)d
cSDR | 0.750 (0.582-0.851) | 0.685 (0.531-0.791) | 0.761 (0.631-0.841) | 0.657 (0.522-0.750)d
cFD/h | 2.290 (1.312-3.897) | 1.967 (0.964-3.478)b | 0.357 (0.167-0.690) | 0.302 (0.135-0.513)
cKappa | 0.668 (0.464-0.749) | 0.672 (0.474-0.764)a | 0.729 (0.656-0.789) | 0.623 (0.533-0.694)d

###  3.2 Non-inferiority to human annotation

The SDA provides an annotation of the EEG that is slightly inferior to human
experts (Table 3; $\Delta\kappa$ fails to span 0). Example EEG segments that
were correctly detected as seizure, missed or falsely labelled as seizure are
shown in Figure 3.

Table 3:  Inter-observer agreement between human and algorithm based annotations of seizure/non-seizure in the validation set of EEG. E1 is human expert 1, E2 is human expert 2, Consensus is the consensus annotation of human experts, $\Delta\kappa$ is difference in IOA between a single human expert and the SDA and the IOA between human experts (results are presented as median, 95% confidence interval), FS is feature set. Inter-observer Agreement | E1 | E2 | Consensus
---|---|---|---
SDA (proposed FS) | 0.591 | 0.703 | 0.729
SDA (Temko FS) | 0.511 | 0.600 | 0.623
E1 | - | 0.792 | 0.805
E2 | - | - | 0.982
SDA $\Delta\kappa$ (proposed FS) | 0.197 (0.149-0.269) | 0.085 (0.011-0.176) | -
SDA $\Delta\kappa$ (Temko FS) | 0.277 (0.220-0.364) | 0.187 (0.113-0.289) | -

At the feature level, the major differences between training and validation
data sets were between the smoothed nonlinear energy operator feature within
EEG epochs without seizure (training data: median [IQR] 39.2 [15.2-74.8],
validation data: median [IQR] 69.9 [43.5-111.4]; values are not normalized).
Increases in the nonlinear energy imply increases in high frequency energy in
seizure signals.

At the annotation level, the SDA has a noticeable reduction in seizure
detection rate and, in general, it generates an excessive number of seizures
compared to the consensus seizure annotation. This is in the context of an SDA
annotation that is, nevertheless, more conservative and diffuse with its
annotation; the SDA has lower estimates of seizure burden measures (maximum
and total, Table 4). The SDA also tends to generate more false detections
during the second half of recordings compared to the first. A representative
example is shown in Figure 2. In general, the median false detection rate was
0.13 per hour (IQR: 0.03-0.17) during the first half of the recording and 0.20
per hour (IQR: 0.06-0.42) during the second half of the recording. The second
half of recordings is more likely to contain EEG patterns associated with
recovery and, therefore, more normal EEG such as sleep states. These patterns
may not be prevalent in the training data.

Table 4:  Seizure summary statistics of seizures annotated by SDAs and the consensus of human experts. Results are presented as median (IQR) summarized across neonates with seizures. Seizure duration is presented as median (IQR) over all annotated seizures. 1 24 neonates had seizures annotated by both human consensus and the SDA (proposed FS) and 2 23 by consensus and SDA (Temko FS). Seizure statistics | Consensus | SDA (proposed FS) | SDA (Temko FS)
---|---|---|---
Seizure patients | 24 | 281 | 272
Seizure number | 33 (11-68) | 29 (17-119) | 14 (8-120)
Seizures (total) | 1362 | 2313 | 1749
Seizure duration (s) | 96 (55-167) | 44 (32-84) | 64 (52-104)
Burden (max; min/h) | 13 (3-29) | 8 (3-20) | 9 (2-21)
Burden (total; min) | 42 (11-163) | 28 (10-138) | 24 (8-201)
Figure 2: An example aEEG and seizure burden time course. Note the appearance
of sleep-wake-cycling after 50 h. This will introduce normal EEG patterns such
as mixed frequency activity, high voltage synchronous activity and trace
alternant, that will not be prevalent in the original training data as those
recordings were typically taken earlier in the seizure burden time course.

(a) True detection on patient 2. The seizure is on all channels, but most
prominent on channel F4-P4.

(b) False detection on patient 34. Channel F3-P3 raises the false alarm.

(c) Missed seizure on patient 1. Seizure is on all channels, but predominantly
on channel P3-P4.

Figure 3: Examples of a true detection, a false detection and a missed
seizure.

###  3.3 Clinical efficacy

To assess the clinical efficacy of the SDA, we used more clinically focused
metrics (Table 5(c)). The key component of this assessment is evaluation of
the short-term seizure burden per hour (example in Figure 4). The correlation
between human consensus and SDA for short-term seizure burden per hour on
concatenated data was 0.89 (95% CI: 0.80-0.96). The correlation between human
experts assessment of short-term seizure burden was 0.95 (95% CI: 0.93-0.98).

Table 5: Agreement in clinically relevant measures between the SDA and human
consensus annotation of the EEG on all neonates. POI - period of clinical
interest (hourly epochs of EEG that exceed a treatment threshold for
seizures), nPOI - not a POI, H - human annotation (consensus), S - algorithm
annotation, HB - high seizure burden, LB - low seizure burden.

(a) Period of Clinical Interest (POI) | POI (S) | nPOI (S)
---|---|---
POI (H) | 151 | 30
nPOI (H) | 103 | 716

(b) Total Seizure Burden | HB (S) | LB (S)
---|---|---
HB (H) | 11 | 2
LB (H) | 0 | 15

(c) Maximum Hourly Seizure Burden | HB (S) | LB (S)
---|---|---
HB (H) | 11 | 3
LB (H) | 0 | 14

Figure 4: Hourly (short-term) seizure burden evolution over time on human
consensus and SDA annotation. Each marker represents an hour segment and the
marker thickness the hourly seizure burden. Time series were cut after no more
events occurred for visualization purposes. Red dots stand for no seizure.

The hourly seizure burden is commonly used as a target to initiate treatment
in clinical studies. The SDA identifies epochs that exceed a target threshold
with an accuracy of 87% (Table 5(a)). The sensitivity of the SDA is 83% and
the specificity is 87% which are in a similar range to inter-observer
agreement between human experts (sensitivity 89% and specificity 97%) for this
measure.

The sensitivity and specificity of the SDA for detecting high total seizure
burden ($>45$ minutes; compared to the consensus annotation) are 85% and 100%,
respectively (Table 5(b)). The sensitivity and specificity of human experts
are 85% and 79%, respectively. Finally, for detecting high maximum hourly
seizure burden ($>13$ minutes per hour; compared to consensus), the
sensitivity and specificity for the SDA are 79% and 100% (Table 5(c)), while
the sensitivity and specificity of human experts are 94% and 92%,
respectively. Importantly, the SDA does not detect high total or hourly
seizure burden in a low seizure burden neonate.

###  3.4 Re-training

Incorporating seizure and non-seizure data from the first 10 recordings in the
validation set into the training data set improved AUC and Kappa on the
remaining 18 neonates in the validation set (Table 6). Improved sensitivity
was offset by a slight reduction in specificity. These improvements were
sufficient to generate an annotation that was non-inferior to the human
annotation ($\Delta\kappa$ = 0.077, 95% CI: -0.002-0.232, $n=18$). In
contrast, the original SDA, when applied to the truncated validation data set
retained its inferior annotation of seizures ($\Delta\kappa$ = 0.084, 95%CI:
0.003-0.226, $n=18$). Re-training did not significantly change primary
performance measures when evaluated within a 10-fold cross-validation on the
original cohort of 79 neonates (Table 6). Specificity was slightly, but
significantly, reduced.

Re-training an SDA based on the Temko feature set showed slight, but not
significant, improvements in performance compared to a version trained on the
original cohort of 79 neonates. Compared to the proposed, retrained SDA, SDA
retrained with Temko features had inferior performance with all measures
except specificity and false detections per hour (Table 6). The annotation of
the EEG by the retrained SDA with the Temko feature set remained inferior to
the human expert annotation ($\Delta\kappa$ = 0.145, 95%CI of 0.083-0.287,
$n=18$).

Table 6:  The performance of the re-trained neonatal SDA. The re-trained neonatal SDA was applied to the original dataset (10-fold cross-validation) and the truncated validation set ($n=18$). These values were compared to the original SDA applied to the truncated validation set and a re-trained SDA based on the Temko feature set (Temko FS). Median (IQR) results are presented for AUC, sensitivity (Sens) and specificity (Spec). The value with 95%CIs in brackets are presented for cAUC, cSDR, cFD/h and cKappa (estimated on concatenated recordings). a and b denote a significant increase and decrease, respectively, between Validation and Validation (retrained), c and d denote a significant increase and decrease, respectively, between Validation (retrained) and Validation (retrained Temko FS), and f denotes a significant decrease between Training (10-fold CV) and Training (outlier detection; 10-fold CV) in Table 2. $n_{t}$ is the number of infants in the cohort and $n_{s}$ is the number of infants with consensus seizure in the cohort. | Original Training | Validation | Validation | Validation
---|---|---|---|---
| (10-fold CV) |  | (retrained) | (retrained Temko FS)
| $n_{t}=79,n_{s}=39$ | $n_{t}=18,n_{s}=14$ | $n_{t}=18,n_{s}=14$ | $n_{t}=18,n_{s}=14$
AUC | 0.979 (0.930-0.996) | 0.965 (0.880-0.980) | 0.976 (0.932-0.983)a | 0.953 (0.867-0.974)d
Sens | 0.816 (0.433-0.947) | 0.519 (0.156-0.759) | 0.648 (0.311-0.837)a | 0.395 (0.057-0.645)d
Spec | 0.991 (0.973-1)f | 0.998 (0.997-0.999) | 0.997 (0.996-0.998)b | 0.998 (0.997-0.999)c
cAUC | 0.964 (0.944-0.980) | 0.964 (0.932-0.978) | 0.974 (0.953-0.984)a | 0.956 (0.912-0.971)d
cSDR | 0.667 (0.515-0.778) | 0.679 (0.447-0.868) | 0.729 (0.447-0.876) | 0.600 (0.307-0.763)
cFD/h | 1.943 (0.918-3.331) | 0.414 (0.145-0.976) | 0.332 (0.147-0.610) | 0.191 (0.092-0.303)
cKappa | 0.680 (0.507-0.769) | 0.722 (0.571-0.795) | 0.735 (0.560-0.808) | 0.662 (0.484-0.726)

##  4 Discussion

In this study, we have evaluated the performance of our neonatal SDA on an
independent validation data set consisting of limited channel, long duration
EEG recordings. The performance of the SDA was not significantly different
between cross-validation on the training data set and the validation data set,
using the primary generalizability measure, the AUC. While the overall
performance of the SDA remained consistent, it was, nevertheless, inferior to
the human expert annotation of the validation set. The inferiority of the SDA
annotation was mitigated by retraining on a larger data set containing data
portions of the validation data set.

There is a limited number of validation studies performed on published
neonatal SDAs [19, 8, 20, 21, 22]. These studies provide seizure detection
rates of 53-66% and false detections of 0.04-2.3 per hour; values that are
compatible with our SDA performance 76% at 0.35 false detections per hour. The
limitations of these studies are that 1) they do not test the generalizability
of their SDAs and 2) they do not evaluate the non-inferiority of the SDA
annotation compared to that of the human expert taking into account the
inherent ambiguity in the annotation of multiple human experts. The algorithm
with the highest congruence between training and validation sets is that of
Temko et al˙ with a seizure detection rate of 71% at 0.25 false detection per
hour on the training set and a seizure detection rate of 75% at 0.36 false
detection per hour on a validation set [22]. This algorithm was further shown
to potentially aid clinical seizure recognition with earlier detection of
seizure onset and improved temporal localization of clinical target periods
[22]. It has recently been clinically trialled and shown to improve the
recognition of hourly seizures, but not the identification on neonates with
seizure, which may help optimize the timing of treatment [24]. We have shown
that our feature set is capable of outperforming the feature set of Temko et
al˙ that underlies this algorithm, when trained on identical data.

Initially, the SDA provided an inferior annotation of seizures compared to
human experts on the validation data set. The most obvious failing of the SDA
was an excessive annotation of seizures. The SDA annotated twice as many
seizures as the human experts although the majority of these seizures were
multiple detections within a single human annotated seizure. Nevertheless, a
considerable proportion were periods annotated by human experts as ’no
seizure’. Reassuringly, visual inspection of false detections showed EEG with
repetitive characteristics, activity that is commonly associate with seizure
(see Figure 3).

There are two key differences between the training and validation data sets
that may explain increased errors in the validation set 1) reduced electrode
density and increased recording duration of EEG in the validation set, and 2)
different brain injury etiologies in the two data sets. Firstly, EEG
monitoring with a reduced electrode set is commonly used in the NICU as it
offers less disruption for neonates, is compatible with visualizations such as
the aEEG and reduces the demand of cot-side resources. While these recordings
capture the majority of seizures and have electrodes located on regions that
are more robust to artefact, we have shown that low density recordings result
in a different annotation to standard EEG recordings and have a reduced
prevalence of seizure [34, 35, 36, 25]. This means that the manifestation of
neonatal seizure may be more diffuse as the distance between seizure loci and
recording electrodes will be higher. Long duration EEG also contains a wider
variety of normal EEG patterns that may not be seen in short duration, higher
density recordings performed early to confirm the presence of seizures. This
explains why the ’no seizure’ feature values within the validation data set
have a larger spread compared to the training data. Secondly, the validation
set contains seizures mostly from a single aetiology (stroke) whereas the
training set contained seizures from several aetiologies (HIE, stroke,
meningitis) [25]. Seizures in neonates with stroke may have different
morphology to seizures from other aetiologies and stroke was poorly
represented in the training data set. This also means that recordings in the
validation set will not contain prolonged periods of isoelectric EEG and
burst-suppression that are apparent in the EEG recordings of neonates with
conditions such as HIE. Despite these potential limitations, the SDA had
consistent, general performance between the data sets and non-inferiority to
human annotation could be achieved by incorporating a subset of the validation
data set into the training set.

Nevertheless, the SDA annotation provides an accurate ($\sim 90\%$)
representation of key interpretations of the annotations such as identifying
periods of EEG that meet treatment thresholds and definitions of high seizure
burden that are associated with poor prognosis. These values may exceed
current clinical practice when on-demand EEG annotation is not available,
however, the risk of false identifications may outweigh increased sensitivity.
A risk that is minimized for a non-inferior annotation.

Bag of features approaches to neonatal seizure detection are being superseded
by deep neural networks [17]. We and others have shown that deep convolutional
neural networks provide similar performance to the proposed ’bag of features’
approach when applied to our training set [28, 17] with deep convolutional
networks providing more consistent performance over several data sets [17].
These methods will, nevertheless, face similar challenges such as lack of
diversity within training data sets, subjectivity in the gold standard of
human expert annotation (the training target) and performance assessment.

A potentially controversial component of our study is the use of 10s as a
minimum duration of seizure. The original choice of a 10s minimum seizure
duration was always arbitrary and acknowledges that as evolution is an
important component of determining the presence of seizure on an EEG
recording, there must be some sort of minimum duration to judge such an
evolution [27]. This minimum seizure duration has been abolished in the latest
recommendations on the Task Force for Neonatal Seizure [4] but its abolition
is disputed [37]. With the lack of biophysical evidence for a minimum seizure
duration, it is best to turn to other forms of evidence such as studies of
inter-observer agreement. These works provide the only data-driven evidence
base for a minimum seizure duration and suggest that 30s may be a feasible
option [36]. This, more lenient, definition of minimum seizure duration has
the potential to improve the accuracy of SDAs. Furthermore, it is a simple way
to achieve SDA annotations that generalize to the larger population of
neonates and are non-inferior to the visual interpretation of the human
expert.

##  5 Conclusion

We validated a neonatal SDA on an independent data set of neonatal EEG. We
show consistent performance between training and validation data sets, but an
SDA annotation that is inferior to the human expert. This annotation is,
nevertheless, sufficient to provide highly accurate estimation of seizure
burden and identification of clinical targets. Increasing the diversity of the
training set resulted in a non-inferior SDA annotation for the left-out part
of the validation data set. This raises key open questions for neonatal SDAs:
1) how much data is required and from what population of neonates should it be
sampled to train neonatal SDAs that are non-inferior to the human expert, 2)
is non-inferior annotation compared to human experts the ultimate target and
3) if and when validated non-inferiority is achieved, will it improve the
clinical recognition of seizures and improve health outcomes for critically
ill neonates.

##  6 Declaration of competing interest and conflict of interest

The authors have no competing interests to declare and no conflicts of
interest.

##  7 Acknowledgements

Funding for this research was provided by the Finnish Cultural Foundation
(00181077), Finnish Pediatric Foundation, the Finnish Academy (313242, 288220,
321235), Aivosäätiö, Neuroscience Center at University of Helsinki, Helsinki
University Hospital, Helsinki University Hospital Research Funds, HUS and
University of Helsinki Researcher position, Sigrid Juselius Foundation and
European Union’s Horizon 2020 Research and Innovation Programme (H2020-MCSA-
ITN-813483 and H2020-MCSA-IF-656131). The study sponsors had no involvement in
the study design, nor in the collection, analysis or interpretation of data.
Additionally we acknowledge the computational resources provided by the Aalto
Science-IT project.

##  8 Author contribution statement

KTT contributed to study design, developed analysis tools, performed data
processing and scripting for results. PN contributed to data collection,
collation and annotation of the EEG recordings. SV contributed to study
design, data collection, collation and annotation of the EEG recordings. NJS
contributed to study design and algorithm development. KTT, NJS, PN, and SV
contributed to the writing and editing of the paper manuscript.

## References

  * [1] J. Rennie, G. Boylan, Treatment of neonatal seizures, Archives of Disease in Childhood-Fetal and Neonatal Edition 92 (2) (2007) F148–F150.
  * [2] S. Björkman, S. Miller, S. Rose, C. Burke, P. Colditz, Seizures are associated with brain injury severity in a neonatal model of hypoxia–ischemia, Neuroscience 166 (1) (2010) 157–167.
  * [3] S. Miller, J. Weiss, A. Barnwell, D. Ferriero, B. Latal-Hajnal, A. Ferrer-Rogers, N. Newton, J. Partridge, D. Glidden, D. Vigneron, et al., Seizure-associated brain injury in term newborns with perinatal asphyxia, Neurology 58 (4) (2002) 542–548.
  * [4] R. M. Pressler, M. R. Cilio, E. M. Mizrahi, S. L. Moshé, M. L. Nunes, P. Plouin, S. Vanhatalo, E. Yozawitz, L. S. de Vries, K. Puthenveettil Vinayan, et al., The ILAE classification of seizures and the epilepsies: Modification for seizures in the neonate. position paper by the ILAE task force on neonatal seizures, Epilepsia 62 (3) (2021) 615–628.
  * [5] A. Liu, J. Hahn, G. Heldt, R. Coen, Detection of neonatal seizures through computerized EEG analysis, Electroen Clin Neuro 82 (1) (1992) 30–37.
  * [6] J. Gotman, D. Flanagan, J. Zhang, B. Rosenblatt, Automatic seizure detection in the newborn: Methods and initial evaluation, Electroen Clin Neuro 103 (3) (1997) 356–362.
  * [7] P. Celka, P. Colditz, A computer-aided detection of EEG seizures in infants: A singular-spectrum approach and performance comparison, IEEE T Bio-Med Eng 49 (5) (2002) 455–462.
  * [8] L. Smit, R. Vermeulen, W. Fetter, R. Strijers, C. Stam, Neonatal seizure monitoring using non-linear EEG analysis, Neuropediatrics 35 (06) (2004) 329–335.
  * [9] M. A. Navakatikyan, P. B. Colditz, C. J. Burke, T. E. Inder, J. Richmond, C. E. Williams, Seizure detection algorithm for neonates based on wave-sequence analysis, Clinical Neurophysiology 117 (6) (2006) 1190–1203.
  * [10] W. Deburchgraeve, P. Cherian, M. De Vos, R. Swarte, J. Blok, G. Visser, P. Govaert, S. Van Huffel, Automated neonatal seizure detection mimicking a human observer reading EEG, Clinical Neurophysiology 119 (11) (2008) 2447–2454.
  * [11] J. Mitra, J. R. Glover, P. Y. Ktonas, A. T. Kumar, A. Mukherjee, N. B. Karayiannis, J. D. Frost Jr, R. A. Hrachovy, E. M. Mizrahi, A multi-stage system for the automated detection of epileptic seizures in neonatal EEG, J Clin Neurophysiol 26 (4) (2009) 218–226.
  * [12] A. Temko, E. Thomas, W. Marnane, G. Lightbody, G. Boylan, EEG-based neonatal seizure detection with support vector machines, Clinical Neurophysiologyphysiol 122 (3) (2011) 464–473.
  * [13] A. Ansari, P. Cherian, A. Dereymaeker, V. Matic, K. Jansen, L. De Wispelaere, C. Dielman, J. Vervisch, R. Swarte, P. Govaert, et al., Improved multi-stage neonatal seizure detection using a heuristic classifier and a data-driven post-processor, Clinical Neurophysiology 127 (9) (2016) 3014–3024.
  * [14] A. O’Shea, G. Lightbody, G. Boylan, A. Temko, Neonatal seizure detection using convolutional neural networks, in: 2017 IEEE 27th International Workshop on Machine Learning for Signal Processing (MLSP), IEEE, 2017, pp. 1–6.
  * [15] K. T. Tapani, S. Vanhatalo, N. J. Stevenson, Time-varying EEG correlations improve automated neonatal seizure detection, International Journal of Neural Systems 29 (04) (2019) 1850030.
  * [16] A. H. Ansari, P. J. Cherian, A. Caicedo, G. Naulaers, M. De Vos, S. Van Huffel, Neonatal seizure detection using deep convolutional neural networks, International Journal of Neural Systems 29 (04) (2019) 1850011.
  * [17] A. O’Shea, G. Lightbody, G. Boylan, A. Temko, Neonatal seizure detection from raw multi-channel EEG using a fully convolutional architecture, Neural Networks 123 (2020) 12–25.
  * [18] D. Y. Isaev, D. Tchapyjnikov, C. M. Cotten, D. Tanaka, N. Martinez, M. Bertran, G. Sapiro, D. Carlson, Attention-based network for weak labels in neonatal seizure detection, Proceedings of machine learning research 126 (2020) 479.
  * [19] J. Gotman, D. Flanagan, B. Rosenblatt, A. Bye, E. Mizrahi, Evaluation of an automatic seizure detection method for the newborn EEG, Electroencephalography and Clinical Neurophysiology 103 (3) (1997) 363–369.
  * [20] R. Lawrence, A. Mathur, S. N. T. Tich, J. Zempel, T. Inder, A pilot study of continuous limited-channel aEEG in term infants with encephalopathy, The Journal of Pediatrics 154 (6) (2009) 835–841.
  * [21] P. Cherian, W. Deburchgraeve, R. Swarte, M. De Vos, P. Govaert, S. Van Huffel, G. H. Visser, Validation of a new automated neonatal seizure detection system: a clinician’s perspective, Clinical Neurophysiology 122 (8) (2011) 1490–1499.
  * [22] S. R. Mathieson, N. J. Stevenson, E. Low, W. P. Marnane, J. M. Rennie, A. Temko, G. Lightbody, G. B. Boylan, Validation of an automated seizure detection algorithm for term neonates, Clinical Neurophysiology 127 (1) (2016) 156–168.
  * [23] A. Temko, E. Thomas, W. Marnane, G. Lightbody, G. Boylan, Performance assessment for EEG-based neonatal seizure detectors, Clinical Neurophysiology 122 (3) (2011) 474–482.
  * [24] A. M. Pavel, J. M. Rennie, L. S. de Vries, M. Blennow, A. Foran, D. K. Shah, R. M. Pressler, O. Kapellou, E. M. Dempsey, S. R. Mathieson, et al., A machine-learning algorithm for neonatal seizure recognition: a multicentre, randomised, controlled trial, The Lancet Child & Adolescent Health 4 (10) (2020) 740–749.
  * [25] N. Stevenson, K. Tapani, L. Lauronen, S. Vanhatalo, A dataset of neonatal EEG recordings with seizure annotations, Scientific data 6 (1) (2019) 1–8.
  * [26] P. Nevalainen, M. Metsäranta, S. Toiviainen-Salo, T. Lönnqvist, S. Vanhatalo, L. Lauronen, Bedside neurophysiological tests can identify neonates with stroke leading to cerebral palsy, Clinical Neurophysiology 130 (5) (2019) 759–766.
  * [27] R. R. Clancy, A. Legido, The exact ictal and interictal duration of electroencephalographic neonatal seizures, Epilepsia 28 (5) (1987) 537–541.
  * [28] N. Stevenson, K. Tapani, S. Vanhatalo, Hybrid neonatal EEG seizure detection algorithms achieve the benchmark of visual interpretation of the human expert, in: 2019 41st Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), IEEE, 2019, pp. 5991–5994.
  * [29] J. H. Friedman, J. L. Bentley, R. A. Finkel, An algorithm for finding best matches in logarithmic expected time, ACM Transactions on Mathematical Software (TOMS) 3 (3) (1977) 209–226.
  * [30] J. Cohen, A coefficient of agreement for nominal scales, Educational and psychological measurement 20 (1) (1960) 37–46.
  * [31] L. G. van Rooij, L. Hellström-Westas, L. S. de Vries, Treatment of neonatal seizures, in: Seminars in Fetal and Neonatal Medicine, Vol. 18, Elsevier, 2013, pp. 209–215.
  * [32] R. Wickström, B. Hallberg, M. Bartocci, Differing attitudes toward phenobarbital use in the neonatal period among neonatologists and child neurologists in Sweden, European Journal of Paediatric Neurology 17 (1) (2013) 55–63.
  * [33] L. Kharoshankaya, N. J. Stevenson, V. Livingstone, D. M. Murray, B. P. Murphy, C. E. Ahearne, G. B. Boylan, Seizure burden and neurodevelopmental outcome in neonates with hypoxic–ischemic encephalopathy, Developmental Medicine & Child Neurology 58 (12) (2016) 1242–1248.
  * [34] N. J. Stevenson, L. Lauronen, S. Vanhatalo, The effect of reducing EEG electrode number on the visual interpretation of the human expert for neonatal seizure detection, Clinical Neurophysiology 129 (1) (2018) 265–270.
  * [35] L. Webb, M. Kauppila, J. A. Roberts, S. Vanhatalo, N. Stevenson, Automated detection of artefacts in neonatal EEG with residual neural networks, Computer Methods and Programs in Biomedicine (2021) 106194.
  * [36] N. J. Stevenson, R. R. Clancy, S. Vanhatalo, I. Rosén, J. M. Rennie, G. B. Boylan, Interobserver agreement for neonatal seizure detection using multichannel EEG, Ann Clin Transl Neurol 2 (11) (2015) 1002–1011.
  * [37] R. A. Shellhaas, Neonatal seizures reach the mainstream: The ILAE classification of seizures in the neonate, Epilepsia 62 (3) (2021) 629–631.