\section{Introduction}
\label{sec:introduction}

\IEEEPARstart{E}{pilepsy} is a chronic neurological disorder that affects about 1\% population in the world \cite{hosseini2017optimized,detti2018patient}. Such abnormal brain activity may cause seizures, unusual behavior, and sensations, and sometimes loss of awareness.
Thus, an accurate and timely diagnosis for epileptic seizure is crucial to reduce both life and financial cost.
Non-invasive Electroencephalogram (EEG) spontaneous is the most common method for seizure diagnosis \cite{zhang2019survey} but faces several challenges.

\subsection{Challenges}
\label{sub:challenges}

First, most of the existing epileptic seizure detection methods focus on patient-dependent scenario \cite{detti2018patient,altaf201516} but rarely consider the patient-independent situation\footnote{In some literature, the patient dependent/independent are called patient specific/non-specific.}.
The former refers to detect the epileptic seizure of a patient by learning from his own historical records, while the latter learn from the records from other patients.
Patient-dependent methods can achieve high accuracy because the training and testing samples are gathered from the same source and have similar distribution. Patient-dependent algorithms have been widely investigated in the last decade due to the good performance for recorded patients, however, that becomes deficient for new patients since the data distributions from different patients are much more diverse.
In contrast, patient-independent methods advance in alerting potential patients but are easily corrupted by inter-patient noises (e.g., gender, age, and epileptic type).
A majority of existing studies fail to eliminate such noises, and, to the best of our knowledge, few studies try to eliminate the inter-patient corruption by modeling the common features from the training samples.
Thus, in this study, we propose a robust method that can learn the common seizure pattern while also mitigate the influence of inter-patient factors.

Second, the non-invasive EEG signals are with low signal-to-noise ratio and are fragile to subjective factors (e.g., emotion and fatigue) and environmental factors (e.g., noise) \cite{li2015feature,goh2018spatio}.
To discover a latent and informative representation of the signals, most traditional seizure detection methods do feature engineering on raw EEG signals, which is time-consuming and highly dependent on expertise \cite{boashash2016automatic}.
Recently, deep learning has shown great success in lots of research topics (e.g., computer vision, natural language processing, and brain-computer interface)
because of the excellent automatic feature learning ability \cite{thodoroff2016learning}.
Therefore, deep learning algorithms have great potential in reducing the negative effects of manual feature engineering by automatically capturing the key difference between epileptic and non-epileptic patterns.

Another issue of current epileptic seizure diagnosis approaches is about the explainability.
Generally, we are interested in not only the seizure states of a patient but also the fine-grained observations.
For example, the epileptic seizure may occur in one specific brain region (e.g., temporal lobe) that is known as focal seizure or several brain regions (e.g., frontal, temporal, and occipital lobes) that is known as generalized epilepsy \cite{arunkumar2018entropy}.
 The current seizure tests can provide the overall detection results while are unable to show the specific active brain regions.
To this end, for a specific patient, we attempt to design a powerful model, involving attention mechanism \cite{yuan2019fusionatt}, which can not only recognize the seizure state but also discover which EEG channels are more important to the diagnosis decision.

\subsection{Motivation}
\label{sub:motivation}
The key target of this work is to develop an automatic representation learning algorithm that is dependent on neurological status (e.g., seizure or normal) while is independent of person states (e.g., age and gender). The representation is learned from non-invasive EEG signals and further be used to diagnose the epileptic seizure.
As mentioned above, a main challenge in patient-independent seizure feature learning is that the learning process can be easily affected by patient-related information (e.g., age, gender, and emotion). In order to eliminate such influences, we propose an adversarial method \cite{zhang2019adversarial} to decompose the EEG signals into a seizure-related component and a patient-related component.

In short, the intuitive idea of this work is separating the input raw EEG signals (denoted by $\mathcal{E}$) into two parts:
a seizure-related component, denoted by $\mathcal{S}$, that contains informative descriptions about seizure state and is insensitive to the patient identity; and a patient-related component, denoted by $\mathcal{P}$, that includes patient identity information and is insensitive to the seizure state.
The following two constraints should be satisfied during the process.
First, the sum of seizure-related component $\mathcal{S}$ and patient-related component $\mathcal{P}$ should be equal to the original EEG signal $\mathcal{E}$.
And second, the decomposed $\mathcal{S}$ and $\mathcal{P}$ are supposed to contain pure and informative seizure and patient features, respectively.

\subsection{Contributions}
\label{sub:contributions}

We propose an adversarial representation learning framework to construct a robust patient-independent detection algorithm that exterminates the corruptions by inter-patient noises. The framework harnesses both deep generative model and convolutional discriminative model to capture informative seizure and patient representations.

Moreover, we improve the explainability by incorporating attention mechanism. The model can not only enhance the contribution of important EEG channels but also visualize the determine brain regions.
The main contributions of this work are highlighted as follows:

\begin{itemize}
    \item We propose a robust framework targeting patient-independent epileptic seizure detection. The proposed approach is efficient in capturing the seizure-specific representations directly from the raw EEG signals. The resuable Python-based code can be find here\footnote{https://github.com/xiangzhang1015/adversarial\_seizure\_detection}.
    \item We propose a unified adversarial learning framework to extract the informative seizure-specific representation while eliminating the inter-patient noises. Moreover, the attention mechanism enables our approach to automatically explore the importance of each EEG channel, consequently, to bring fine-grained analysis for seizure diagnosis.
    \item We implement extensive experimental comparisons in order to evaluate the proposed approach over a benchmark seizure diagnosis dataset. The results demonstrate that our model outperforms competitive state-of-the-art methods effectively and efficiently.
\end{itemize}

\section{Related Work}
\label{sec:related_work}
To date, numerous analysis methods are proposed in order to automatically diagnose epileptic seizure. Non-invasive EEG signals provide important information about epileptogenic networks that must be analyzed and understood before the initiation of therapeutic procedures \cite{faust2015wavelet}. In order to discover the latent informative features, Bhattacharyya et al. \cite{bhattacharyya2017multivariate} adopted an empirical wavelet transform method to explore the multivariate signals in the time-frequency domain, in which the instantaneous amplitudes and frequencies were jointly considered.
Moreover, Fan et al. \cite{fan2018detecting} investigated the spatial-temporal synchronization pattern in epileptic brainwaves through spectral graphical representation and developed an efficient multivariate approach for detecting seizure in real-time. Detti et al. \cite{detti2018patient} attempted to find the distinguishable non-invasive EEG patterns in order to predict the epileptic seizure onset a short-term (e.g., a few minutes) ahead.

Apart from traditional statistic models, deep learning-based models have gained increasing attention. Schirrmeister et al. \cite{schirrmeister2017deep} adopted a deep learning-based architecture, ConvNets, to diagnose pathological patients from normal individuals based on the decoding of EEG samples. The experiments showed that the deep neural networks outperformed the conventional linear analysis classifiers.
Lin et al. \cite{lin2016classification} proposed a deep framework for the automatic detection of epileptic EEG by combining a stacked sparse autoencoder and a softmax classifier, which firstly learned the sparse and high-level representations from the preprocessed data and then sent these representations into a softmax classifier for automatic diagnosis.

\textbf{Summary.} However, all the mentioned seizure diagnosis methods (either traditional or deep learning models) are focusing on the patient-dependent scenario which means that the training EEG samples and testing samples are from the same or same group patients. Although very limited studies mentioned patient-independent validation like \cite{orosco2016patient}, they attempt to extract the general representations through enhancing training dataset instead of designing a robust patient-independent feature learning algorithm.
In this way, the diagnosis model can handle intra-patient factors but fail to eliminate inter-patient noise. The clinical therapy requires more complex situations (patient-independent) where the testing patient is unseen in the training stage. It is much difficult to perform accurate seizure diagnosis considering personal factors such as chronological age, gender, characteristics, and health state. In summary, proposing a robust method to deal with the patient-independent challenge is necessary.

\section{Methodology}
\label{sec:methodology}

\begin{figure*}[!t]
\centering
  \includegraphics[width=\linewidth]{method_overview.pdf}
  \caption{Workflow of the proposed patient-independent epileptic seizure diagnosis method. The collected raw EEG signals are first embedded into a latent space and be decomposed into seizure and patient representation. Afterward, the implicate representations are used to reconstruct the seizure and patient components, respectively.  At the meantime, the representations are sent to the seizure diagnosis and patient detection classifiers in order to diagnose the user's seizure stage and recognize the user's identity. The key branch is the seizure diagnosis while the patient detection branch is adopted to separate the patient-related information from the input EEG signals. The dashed lines represent the data flow of EEG decomposition.
  }
  \label{fig:method_overview}

\end{figure*}

\subsection{Overview}
\label{sub:overview}

The proposed patient-independent seizure representation learning method contains three sections (Figure~\ref{fig:method_overview}): the EEG decomposition, the seizure diagnosis, and the patient detection. These three sections are jointly trained in an adversarial way. To satisfy the first constraint, the EEG signals are mapped into a latent space and decompose into the seizure and patient representation corresponding to the seizure state and personal information, respectively. Then, the two latent representations are used to reconstruct the original EEG signal. If the reconstructed signal equal to the input signal, we regard the learned seizure and patient representation as effective features of the decomposed EEG data.

In addition, to satisfy the second constraint that guarantees the purification of the decomposed representations, we design a seizure-based diagnosis section and a patient classifier section following the latent representations.
The seizure diagnosis section, which is following the seizure representation, aims to recognize the epileptic seizure from a normal state. And a satisfactory seizure representation is learned until the diagnosis achieves competitive performance.
Similarly, the patient classifier that follows the patient representation targets recognizing the user's identity, and a well-performed diagnosis assures the pureness of the patient representation.

The EEG decomposition section focuses more on the reconstruction performance of the decomposed components than the pureness of the components; conversely, the seizure diagnosis section and patient classifier consider the pureness of the two components other than the effectiveness of the representations. In a whole, the two aspects are optimized toward opposite directions, which form an adversarial training situation.
It should be noted that the patient classifier is only adopted to improve the adversarial training and help to have a better signal decomposition performance.

Furthermore, we combine an attention mechanism with our seizure diagnosis section in order to pay more attention to the important EEG electrodes. The effectiveness of attention mechanism has been demonstrated by a number of research topics like natural language process \cite{zhou2016attention}. The modified attention-based seizure diagnosis is supposed to discover which electrodes make more contributions to the

seizure.

\subsection{EEG Decomposition}
\label{sub:signal_decomposition}
Assume the gathered EEG signals are denoted by $\mathcal{E} = \{\bm{E_i} \in \mathbb{R}^{M \times N} \}$, in which each EEG sample $\bm{E_i}$ has $M$ time series and $N$ channels. We suppose the EEG signals can be decomposed into the seizure-related component $\mathcal{S} = \{\bm{S_i} \in \mathbb{R}^{M \times N} \}$ and the patient-related component $\mathcal{P} = \{\bm{P_i} \in \mathbb{R}^{M \times N} \}$, described as\footnote{For simplify, we ignore the sub script.}
\begin{equation}
    \bm{E} = \bm{S} + \bm{P}
\end{equation}

Nevertheless, the $\mathcal{S}$ and $\mathcal{P}$ cannot be directly calculated through traditional methods.
To overstep this challenge, we proposed a novel decomposition model based on deep neural networks. In particular, the $\bm{E}$ is transformed into a latent space through convolutional operation followed by a max-pooling layer:
\begin{equation}
    \bm{\bar{S}} = ReLU(\bm{w_s} \circledast \bm{E} + \bm{b_s})
\end{equation}
\begin{equation}
    \bm{\bar{P}} = ReLU(\bm{w_p} \circledast \bm{E} + \bm{b_p})
\end{equation}
where $\bm{\bar{S}}, \bm{\bar{P}} \in \mathbb{R}^{J \times K \times H}$ denote the latent seizure and patient representation, in which $J$ and $K$ denote the rows and columns of the learned representation, respectively, along with $H$ denotes the number of convolutional filters. The $\circledast$ denotes convolution operation and $\bm{w_s}, \bm{b_s}, \bm{w_p}, \bm{b_p}$ represent the corresponding weights and biases. In each convolutional layer, the activation function is ReLU and the padding method is 'SAME'.

Then, in order to ensure the representations to contain enough discriminative information, we attempt to reconstruct the EEG signals through deconvolutional operations:
\begin{equation}
    \bm{S} = ReLU(\bm{w'_s} \circledast \bm{\bar{S}} + \bm{b'_s})
\end{equation}
\begin{equation}
    \bm{P} = ReLU(\bm{w'_p} \circledast \bm{\bar{P}} + \bm{b'_p})
\end{equation}
where $\bm{S}$ and $\bm{P}$ represent the decomposed seizure- and patient-related components, which have the same shape with the input EEG signal $\bm{E}$.
The reconstructed signal $\bm{E'}\in \mathbb{R}^{M \times N}$ can be calculated by combining the two decomposed signals
\begin{equation}\label{eq:w1}
    \bm{E'} = w_1\bm{S} + w_2\bm{P}; w_1+w_2 =1
\end{equation}
where $w_1$ and $w_2$ are weights.

To guarantee the decomposition performance and separate the useful seizure-related information while eliminating the intra-patient corruption, we force the reconstructed signal $\bm{E'}$ to approximate the original signal $\bm{E}$ by minimizing the Mean Square Error (MSE, a.k.a. $\ell_2$) loss between them. The MSE distance has been widely applied to design the reconstruction loss function due to its simplicity and efficiency \cite{amaral2013using}. The MSE loss indeed outperforms the cross-entropy and integral probability metrics (e.g., Wasserstein loss \cite{arjovsky2017wasserstein}) in our preliminary experiments. The MSE loss function can be calculated as
\begin{equation}
    \mathcal{L}_{D} = \left \|  \bm{E} - \bm{E'}\right \|_2
\end{equation}

The learned latent seizure $\bm{\bar{S}}$ and patient representation $\bm{\bar{P}}$, compared to the decomposed components $\bm{S}$ and $\bm{P}$, have lower dimension and represent high-level features.
Thus, the learned representations are used for the corresponding seizure diagnosis and patient detection in the next process.

\subsection{Attention-based Seizure Diagnosis}
\label{sub:task_specific_classifier}

The weights of different EEG electrodes are varying in seizure diagnosis \cite{temko2011instantaneous}; thus, we involve the attention mechanism to learn the channel importance and pay different attention to various signal channels. The attention mechanism, as mentioned above, allows modeling of dependencies among input sequences \cite{vaswani2017attention} and has shown success in some research topics \cite{zhou2016attention}. Meanwhile, the excellent latent feature learning ability of Convolutional Neural Networks (CNNs) has been widely used in research areas such as computer vision \cite{zhang2017learning} and natural language process \cite{zhang2016mgnc}.
Thus, we propose an attention-based CNN algorithm to automatically extract the distinctive information from the received seizure representation $\bm{\bar{S}}$ which has shape $[J, K, H]$.

As shown in Figure~\ref{fig:method_overview}, the attention-based seizure diagnosis contains four convolutional stages: the flatten layer, two fully-connected layers (FC) and an output layer. Each convolutional stage contains a convolutional layer and a max-pooling layer. We assume the input layer receives the learned seizure feature $\bm{\bar{S}}$ and sent to the convolutional layers:
\begin{equation}
    \bm{C^i} = ReLU(\bm{w^i_s} \circledast \bm{L^{i-1}} + \bm{b^i_s}), i \in \{1, 2, 3\}
\end{equation}
where $\bm{C^i}$ and $\bm{L^{i}}$ denote the $i$-th convolutional and pooling layer, respectively.
If $i=1$, $\bm{L^0}$ equals to the input $\bm{\bar{S}}$; otherwise, the pooling layer can be measured by:
\begin{equation}
    \bm{L^i_u} = \underset{u \in \mathcal{U}}{max}\{\bm{C^i_u}\}    , i \in \{1, 2, 3\}
\end{equation}
where $\mathcal{U}$ represents the max-pooling perception field and $u$ denotes the $u$-th element.

The forth pooling layer $\bm{L^4}$ is flattened to 1-D vector $\bm{F^0}$ and then fed into the FC layers:
\begin{equation}
    \bm{F^i} = sigmoid(\bm{w^i_f} \bm{F^{i-1}} + \bm{b^i_f}), f \in \{ 1, 2, 3\}
\end{equation}
where $sigmoid$ denotes the activation function and $\bm{F^i}$ denotes the $i$-th FC layer.

The predicted seizure state $\bm{\hat{y_s}}$ is related to the last FC layer $\bm{F^2}$ and the learned attention weights $\bm{a_{tt}}$:
\begin{equation}
    \bm{\hat{y_s}} = \bm{F^2} \cdot \bm{a_{tt}}
\end{equation}
where $\cdot$ denotes dot production and $\bm{a_{tt}}$ is directly learn from the EEG sample:
\begin{equation}
\label{eq:att}
    \bm{a_{tt}} = sigmoid(\bm{w_a} \bm{E} + \bm{b_a})
\end{equation}
where $\bm{w_a}$ and $\bm{b_a}$ denote the corresponding parameters.

The cross-entropy loss function of seizure diagnosis $\mathcal{L}_s$ is defined as:
\begin{equation}
    \mathcal{L}_s = -(\bm{y_s}log(\bm{p(\bm{\hat{y_s}})})+ (1-\bm{y_s})log(1-\bm{p(\bm{\hat{y_s}})}))
\end{equation}
where $\bm{y_s}$ and $\bm{p(\bm{\hat{y_s}})}$ denote the ground truth of seizure state and the predicted probability of the patient in seizure state, respectively.

\subsection{Patient detection}
\label{sub:subject_detection}
The architecture of the patient detection is almost identical to the attention-based seizure diagnosis except two main differences: 1) the patient detection component does not need the attention module; 2) the patient detection performs multi-class instead of binary classification.
Similar to $\mathcal{L}_s$, we can measure the multi-class cross-entropy loss function $\mathcal{L}_p$ of the patient identity detection classifier as:
\begin{equation}

    \mathcal{L}_p = - \sum_{c=1}^{C}  (\bm{y_{p}}log(\bm{\hat{y_{p}}}))
\end{equation}
where $\bm{y_p}$, $\bm{\hat{y_{p}}}, $and $C$ denote the ground truth of patient identity, the predicted identity, and the overall number of training patients, respectively.

\subsection{Training Details}
\label{sub:loss_function}
In the proposed approach, we have three loss functions which are the reconstruction loss $\mathcal{L}_{D}$ in signal decomposition and the classification loss $\mathcal{L}_s$ and $\mathcal{L}_p$ in the attention-based seizure diagnosis and patient detection, respectively.

We propose an adversarial training strategy to jointly train all the loss functions:
\begin{equation}
    \mathcal{L} = \mathcal{L}_{D} + \mathcal{L}_s + \mathcal{L}_p +\ell_2
\end{equation}
which
assures all losses work on the gradient at the same time and converge to a trade-off position of that balances the reconstruction performance and the purification of the decomposed signals. The $\ell_2$ represents the $\ell_2$ norm with coefficient of $0.0001$ to prevent overfitting.
Then, since the seizure diagnosis is the most crucial component in this approach, we train the $\mathcal{L}_s$ one more time in order to raise its priority. To sum up, in each training epoch, the $\mathcal{L}$ and $\mathcal{L}_s$ are optimized once in turn. We adopt the Adam optimizer with learning rate $e^{-4}$ for both loss functions. The algorithm is trained for $250$ epochs. A dropout layer with $0.8$ keep rate is added to the flatten layer in both classifiers in order to prevent overfitting. The overall flowchart is presented in algorithm~\ref{alg:method}.

\begin{algorithm}[!t]
\begin{small}
\caption{The Proposed Approach}
\label{alg:method}
\begin{algorithmic}[1]
\renewcommand{\algorithmicrequire}{\textbf{Input:}}
\renewcommand{\algorithmicensure}{\textbf{Output:}}
 \REQUIRE Raw EEG data $\bm{E}$
 \ENSURE  Seizure State: $\bm{\hat{y_s}}$
 \FOR{$\bm{E} \in \mathcal{E}$}
    \STATE \#EEG Decomposition:
    \STATE $\bm{\bar{S}} = ReLU(\bm{w_s} \circledast \bm{E} + \bm{b_s})$
    \STATE $\bm{\bar{P}} = ReLU(\bm{w_p} \circledast \bm{E} + \bm{b_p})$
    \STATE $\bm{S} = ReLU(\bm{w'_s} \circledast \bm{\bar{S}} + \bm{b'_s})$
    \STATE $\bm{P} = ReLU(\bm{w'_p} \circledast \bm{\bar{P}} + \bm{b'_p})$
    \STATE $\bm{E'} = (\bm{S} + \bm{P})/2$
    \STATE $\mathcal{L}_{D} = \left \|  \bm{E} - \bm{E'}\right \|_2$
    \STATE \#Attention-based Seizure Diagnosis:
        \STATE $\bm{F^2} \leftarrow \bm{\bar{S}}$
        \STATE $\bm{a_{tt}} = sigmoid(\bm{w_a} \bm{E} + \bm{b_a})$
        \STATE $\bm{\hat{y_s}} = \bm{F^2} \cdot \bm{a_{tt}}$
        \STATE $\mathcal{L}_s \leftarrow \bm{y_s}, \bm{\hat{y_s}}$
    \STATE \#Patient Detection:
        \STATE $\bm{\hat{y_p}} \leftarrow \bm{\bar{P}}$
        \STATE $\mathcal{L}_p \leftarrow \bm{y_p}, \bm{\hat{y_p}}$
\ENDFOR
\STATE Minimize $\mathcal{L} = \mathcal{L}_{D} + \mathcal{L}_s + \mathcal{L}_p +\ell_2$
\RETURN $\bm{\hat{y_s}}$
\end{algorithmic}
\end{small}
\end{algorithm}

\section{Experiments}
\label{sec:experiments}
In this section, we evaluate the proposed approach over a benchmark dataset of an epileptic seizure. First of all, we present the experimental setting and hyper-parameters in the training model. Then we report the performance comparison results among our method with several competitive state-of-the-art baselines. The comparison is assessed by metrics (such as accuracy, sensitivity, specificity, ROC curves, and AUC scores). Moreover, we provide the latency analysis to demonstrate the efficiency of our method. At last, we attempt to explore the importance of each brain region in seizure diagnosis and present our empirical hypothesis.

\subsection{Dataset}
\label{sub:dataset}
The proposed approach is evaluated over a benchmark dataset TUH corpus \cite{obeid2016temple} for patient-independent epileptic seizure diagnosis. The TUH is a neurological seizure dataset of clinical EEG recordings associated with $22$ channels from a $10/20$ configuration. Since the epileptic seizure may only occur on several channels, we define the seizure onset when more than $12$ channels having a seizure. In order to have enough seizure samples, we only select the subject who contains more than $250$ seconds of seizure state. There are 14 subjects meet the criteria.
For each subject, we use $500$ seconds (half normal and half seizure) of EEG signals with the sampling rate of $250$ Hz. Each EEG sample has a window size of $250$ (i.e., lasting for 1 second) with $50\%$ overlapping. There are 13,986 samples in total. The epileptic seizure state is labelled as 1 while the normal state is labelled as 0.
In the evaluation, we adopt Leave-One-Out strategy, i.e., 14-fold cross-validation, which is one subject work as the testing set while all the other subjects are used for training. For example, the accuracy of subject 0 is calculated when subjects 1 $\sim$ 13 are used for training while subject 0 is used for testing.

\begin{table*}[t]
\centering
\caption{The accuracy comparison among our methods with the baselines and state-of-the-art studies. We adopt Leave-One-Out strategy for patient-independent testing. }
\label{tab:comparison}
\resizebox{\textwidth}{!}{
\begin{tabular}{llllllllllllllll}
\hline
\multirow{2}{*}{\textbf{Methods}} & \multicolumn{15}{c}{\textbf{Subject ID}} \\ \cline{2-16}
 & \textbf{0} & \textbf{1} & \textbf{2} & \textbf{3} & \textbf{4} & \textbf{5} & \textbf{6} & \textbf{7} & \textbf{8} & \textbf{9} & \textbf{10} & \textbf{11} & \textbf{12} & \textbf{13} & \textbf{Average} \\ \hline
\textbf{SVM} & 0.786 & 0.358 & 0.842 & 0.684 & 0.67 & 0.61 & 0.824 & 0.868 & 0.252 & 0.558 & 0.898 & 0.4 & 0.614 & 0.528 & 0.643 \\
\textbf{RF} & 0.398 & 0.681 & 0.92 & 0.754 & 0.314 & 0.634 & 0.696 & 0.568 & 0.298 & 0.602 & 0.906 & 0.612 & 0.658 & 0.603 & 0.619 \\
\textbf{KNN} & 0.786 & 0.652 & 0.878 & 0.674 & 0.682 & 0.608 & 0.808 & 0.618 & 0.562 & 0.612 & 0.882 & \textbf{0.652} & 0.672 & 0.647 & 0.699 \\
\textbf{Schirrmeister \cite{schirrmeister2017deep}} & 0.793 & 0.743 & 0.965 & 0.758 & 0.789 & 0.665 & 0.813 & 0.871 & 0.619 & 0.634 & 0.919 & 0.571 & 0.744 & 0.711 & 0.760 \\
\textbf{Ansari \cite{ansari2019neonatal}} & 0.739 & 0.757 & 0.953 & 0.749 & 0.746 & 0.677 & 0.807 & 0.729 & 0.579 & 0.611 & 0.902 & 0.559 & 0.75 & 0.729 & 0.735 \\
\textbf{Lin \cite{lin2016classification}} & 0.848 & 0.601 & 0.832 & 0.654 & 0.832 & 0.603 & 0.802 & 0.868 & 0.668 & 0.611 & 0.876 & 0.603 & \textbf{0.784} & 0.714 & 0.736 \\
\textbf{Kiral \cite{kiral2018epileptic}} & 0.805 & 0.669 & 0.855 & 0.709 & 0.772 & 0.619 & 0.823 & 0.836 & \textbf{0.746} & 0.598 & 0.835 & 0.556 & 0.745 & 0.726 & 0.736 \\
\textbf{Ours} & \textbf{0.841} & \textbf{0.826} & \textbf{0.978} & \textbf{0.774} & \textbf{0.842} & \textbf{0.733} & \textbf{0.911} & \textbf{0.914} & 0.697 & \textbf{0.652} & \textbf{0.923} & 0.604 & 0.772 & \textbf{0.787} & \textbf{0.805} \\ \hline
\end{tabular}
}
\end{table*}

\subsection{Hyper-parameter Setting}
\label{sub:hyper_parameter_setting}
Next, we report the hyper-parameter settings in detail. In the signal decomposition, the input EEG sample has shape $[M=250, N=22]$, the convolutional layer has 4 filters with size $[3, 3]$, and $[2, 2]$ strides. The followed max-pooling has $[2,1]$ window with $[2,1]$ strides. All the padding methods in this work are 'SAME'. The deconvolutional layers have 1 filter and the other setting are the same.
The seizure- and patient- branches have identical hyper-parameters.

Based on our empirical hyper-parameter tuning, the attention-based seizure diagnosis and patient detection have same hyper-parameter settings: the four convolutional layers have 16, 32, 64, 128 filters with $[3,3], [3,3], [2,2], [2,2]$ sizes (all have $[1,1]$ strides), respectively. All the kernel sizes and strides of the first three max-pooling layers are $[2,2]$ and $[2,1]$ for the last max-pooling layer. The two FC layers have $300$ and $22$ hidden neurons, respectively. The attention layer has $22$ hidden neurons corresponding to $22$ input channels. The detailed information can be found in Figure~\ref{fig:method_overview}. We also investigated the weights in Equation~\ref{eq:w1}, which show that the performance is not sensitive to the averaging weights. The experimental performance ranges in $[0.784, 0.806]$ while $w_1$ changes from 0.1 to 0.9 (while $w_2$ varies from 0.9 to 0.1).  One possible reason is that other parameters within the de/convolutional operation will automatically adjust based on the change of $w_1$, which alleviated the influence brought by $w_1$. For simplicity, we set $w_1=w_2=0.5$ in this work.

\begin{figure}[t]
    \centering
    \begin{subfigure}[t]{0.23\textwidth}
        \centering
        \includegraphics[width=\textwidth]{TUH_cm.pdf}
        \caption{Confusion matrix}
        \label{fig:cm}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[t]{0.25\textwidth}
        \centering
        \includegraphics[width=\textwidth]{roc.pdf}
        \caption{ROC curves with AUC}
        \label{fig:roc}
    \end{subfigure}%
    \caption{Confusion matrix and ROC curves}
    \label{fig:cm_roc}

\end{figure}

\subsection{Overall Comparison}
\label{sub:overall_comparison}
We compare our approach with three widely-used baselines and a set of competitive state-of-the-art methods. The key parameters of the baselines are listed here: Linear SVM ($C = 1$), Random Forest (RF, $n = 50$), and K-nearest Neighours ($k = 3$)\footnote{We directly present the optimal hyper-parameter setting and omitted the tuning procedure to save space.}.
Moreover, the comparable state-of-the-art studies for patient-independent seizure diagnosis are listed as follows. In this work, we mainly compare with the deep learning-based advance methods since it has been demonstrated that deep neural networks generally perform better than traditional models (Section~\ref{sec:related_work}).
\begin{itemize}
    \item Schirrmeister et al. \cite{schirrmeister2017deep} apply convolutional neural networks to distinguish seizure segments by decoding task-related information from EEG signals.
    \item Ansari et al. \cite{ansari2019neonatal} adopt CNN to extract the latent features which are fed into a Random Forest (RF) classifier for the final epileptic seizure detection in neonatal babies.
    \item Lin et al. \cite{lin2016classification} propose a sparse deep autoencoder with three hidden layers to extract the representative features from epileptic EEG signals.
    \item Kiral et al. \cite{kiral2018epileptic} design a deep neural network for seizure diagnosis and further develop a prediction system on a wearable device.
\end{itemize}
Since the state-of-the-art algorithms are designed and evaluated in the context of patient-dependent situations, to make them comparable, we reproduce all the state-of-the-art models in python code in our work to perform a patient-independent task. Take \cite{schirrmeister2017deep} as an example, we built a CNN model with 7 convolutional layers (the architecture and hyper-parameters are according to Table 2 in \cite{schirrmeister2017deep}) to perform patient-independent epileptic seizure detection.
All the baselines and our method are evaluated over the same dataset with the Leave-One-Out strategy. The overall performance of all the compared methods is reported in Table~\ref{tab:comparison}. From this table, we can clearly observe that our approach outperforms all the baselines and state-of-the-art models by achieving the average accuracy of 80.5\%, illustrating the effectiveness of the proposed method in patient-independent seizure representation learning. The advantage is demonstrated under all of the 14 subjects.

Moreover, to have a closer observation of the diagnosis results, we report the confusion matrix and receiver operating characteristic (ROC) curves with area under the curve (AUC) scores. The results of one of the best-performed subjects are shown in Figure~\ref{fig:cm_roc}. From the confusion matrix, we can observe that our approach achieves a high sensitivity of 97.4\% and a slightly low but still competitive specificity of 88.1\%.
Moreover, the ROC curves and AUC score are presented in Figure~\ref{fig:roc}.

Furthermore, we present the converging trend with the increase of training iterations. Figure~\ref{fig:loss_acc} tracks the evolution of the overall loss function and the testing accuracy, which shows that our model can converge smoothly. In Figure~\ref{fig:loss}, we dig deeper to investigate the trend of each loss function component. It can be observed that the seizure detection loss $\mathcal{L}_s$ is the most important loss which takes about half proportion. This also demonstrates that our strategy of training $\mathcal{L}_s$ one more time is appropriate.

\begin{comment}
\begin{figure}[t]
    \centering
    \begin{subfigure}[t]{0.24\textwidth}
        \centering
        \includegraphics[width=\textwidth]{training_time.pdf}
        \caption{Training time}
        \label{fig:training}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[t]{0.24\textwidth}
        \centering
        \includegraphics[width=\textwidth]{testing_time.pdf}
        \caption{Testing latency}
        \label{fig:tesing}
    \end{subfigure}%
    \caption{Latency analysis}
    \label{fig:latency}
    \vspace{-3mm}
\end{figure}

\subsection{Latency Analysis}
\label{sub:latency_analysis}
In addition to the diagnosis performance, latency is another important metric for an EEG-based epileptic seizure diagnosis system.
In this section, we present the latency comparison among the proposed method and the state-of-the-art baselines. All the experiments are run in an NVIDIA TITAN X GPU platform with 10 Gbps memory speed and 3584 cores.

Figure~\ref{fig:latency} illustrates the training time and testing latency achieved by the proposed model in comparison with the baselines. It can be observed that our method takes $1417.4 $seconds for training but only require $0.06$ seconds for testing. In which, the EEG decomposition component costs $0.01$ seconds while the attention-based seizure diagnosis takes around $0.05$ seconds.
In summary, our method has competitive latency even though takes the longest training time compared with other methods. Since the input EEG signals are acquired from 1 second, the diagnosis latency of 0.06 seconds is acceptable in a potential online epileptic seizure diagnosis system.
\end{comment}

\begin{figure}[t]
    \centering
    \begin{subfigure}[t]{0.24\textwidth}
        \centering
        \includegraphics[width=\textwidth]{cost_acc.pdf}
        \caption{Overall Loss and Accuracy}
        \label{fig:loss_acc}
    \end{subfigure}%
        \hfill
    \begin{subfigure}[t]{0.24\textwidth}
        \centering
        \includegraphics[width=\textwidth]{training_cost.pdf}
        \caption{Loss components}
        \label{fig:loss}
    \end{subfigure}
    \caption{Convergence trend with the increase in the number of training iterations}
    \label{fig:}

\end{figure}

\subsection{Latency Analysis}
\label{sub:latency_analysis}
In addition to the diagnosis performance, latency is another important metric for an EEG-based epileptic seizure diagnosis system.
In this section, we present the latency comparison among the proposed method and the state-of-the-art baselines. All the experiments are run in an NVIDIA TITAN X GPU platform with 10 Gbps memory speed and 3584 cores.

We evaluate the training time and testing latency of our method. The results show that we only requires $0.06$s ($0.01$s for EEG decomposition while $0.05s$ for attention-based seizure diagnosis) for testing although $1417.4$s for training.
In summary, the diagnosis latency of our method is acceptable in a potential online epileptic seizure diagnosis system.

\begin{figure*}[t]
    \centering
    \begin{subfigure}[t]{0.24\textwidth}
        \centering
        \includegraphics[width=\textwidth]{figure/brain_lobes.png}
        \caption{Brain lobes distribution}
        \label{fig:lobes}
    \end{subfigure}%

    \begin{subfigure}[t]{0.24\textwidth}
        \centering
        \includegraphics[width=0.8\textwidth]{att_0_backup.pdf}
        \caption{Attention topography (S1)}
        \label{fig:attention_1}
    \end{subfigure}
    \begin{subfigure}[t]{0.24\textwidth}
    \centering
    \includegraphics[width=0.8\textwidth]{att_3_backup.pdf}
    \caption{Attention topography (S8)}
    \label{fig:attention_2}
    \end{subfigure}
    \begin{subfigure}[t]{0.24\textwidth}
    \centering
    \includegraphics[width=0.8\textwidth]{att_1_backup.pdf}
    \caption{Attention topography (S12)}
    \label{fig:attention_3}
    \end{subfigure}
    \caption{Demonstration of the learned channel attention topography of three randomly selected subjects (S1, S8, S12)
    . It can be inferred from the topography that the epileptic seizure occurs from the temporal lobe around the T5 channel and then spread to all the brain surfaces except the earlobes.
    }
    \label{fig:topography}

\end{figure*}

\subsection{Explainability Analysis}
\label{sub:attention}
Explainability is important to a deep learning framework in order to have a better understanding of how the system works. Our method is enabled to capture the importance of each brain region in the diagnosis of epileptic seizure. The proposed method is designed to automatically pay different attention to various EEG channels. For a well-trained and effective model, the learned attention weights can be regarded as the importance of the channels. We conduct leave-one-out cross-validation experiments to investigate the leaned attention topography of each subject (we only present the topographies of three randomly selected subjects due to page limitation). In this section we report which region in the brain hemisphere contributes more to the epileptic seizure.

Figure~\ref{fig:lobes} shows the diagrammatic representation of brain lobes\footnote{https://socratic.org/questions/what-lobe-of-the-brain-is-responsible-for-vision} including frontal, temporal, parietal, and occipital lobes while Figure~\ref{fig:attention_1} to Figure~\ref{fig:attention_3} illustrate the topographies of the learned channel attention from three different subjects based on Equation~\ref{eq:att}.
The red color represents higher weight while blue color indicates lower weight.
From the attention topographies, we can observe that the left temporal lobe (around $T5$ channel) always have the highest weights across all the subjects, which indicates the temporal lobe has a close relationship with epileptic seizure.

 This conclusion is consistent with some biological inference that a region of particular importance in adults with epilepsy is mesial of the temporal lobe \cite{lerche2013ion}. This region is comprised of such structures as the hippocampus and amygdala, which control emotion, and the uncus, which is responsible for processing smells. This is why seizures generally bring strong emotion effects (e.g., fear) and sensation of an acrid odor (e.g., burning rubber)\footnote{https://www.massgeneral.org/childhood-epilepsy/overview/brain.aspx}.
Second, the brain regions near the left earlobe (around $A1$ channel) and the left part of the occipital lobe (near $O1$) seem have a rather lower weights indicating the less significance in seizure diagnosis.
The patterns keeps consistent across different patients, which illustrates that the proposed attention model is reasonable and effective.

The above-mentioned hypotheses are empirical inference based on the designed explainable model and we hope these results could bring some inspiration to researchers in biological, brain science, and medical areas.

\section{Discussion and Future Work}
\label{sec:discussion_and_future_work}

In this work, we propose a novel deep learning framework aiming at patient-independent epileptic seizure diagnosis, which is demonstrated has competitive performance and high interpretability.

Our method can protect the patients' privacy. By extracting the seizure representation in the EEG decomposition part, the learned $\bar{\bm{S}}$ contains rare person-related information. Thus, it's difficult to infer the patient's identity from the $\bar{\bm{S}}$. In the real-world deployment, the patients' privacy is protected.

Next, we discuss the opening challenges of this work and the potential future work.
First, our model achieved the average patient-independent diagnosis accuracy of around 80\%, which outperforms the state-of-the-art baselines but still not enough in clinical deployment.
One feasible solution is adopt ensembling strategies such as voting, bootstrap, aggregating. Take voting strategy as an example, we select multiple EEG segments and predict the seizure state independently.
Then, the final decision is the votes from all the prediction results. Such strategy may increase the latency but improves the prediction accuracy. In the perspective of time efficiency, if we have five EEG segments where each segment costs one second, the data collection latency will be 5 seconds, which is acceptable in a real-world scenario. In the perspective of prediction performance, suppose the original probability of correct classification $p$ is $0.8$; then, after voting from all of the EEG segments, e.g. 5 segments, the probability becomes  $1-[C^4_5p(1-p)^4-C^5_5(1-p)^5] = 0.9933$, which is satisfactory in most of the application scenarios. Thus, we can choose an appropriate number of segments to assure the model achieve better detection accuracy while is time-efficient in real-world cases.

Moreover, the proposed approach has a relatively high complexity containing three computational components. As a result, the training time of our model is much longer than the state-of-the-art models. One of our future scopes is to develop a lightweight framework to save computational resources.

By the way, the proposed method contains a large set of hyper-parameters (such as the convolutional filter size) which may lead to heavy tuning workload. One possible method is to adopt some hyper-parameter optimization algorithm (e.g., orthogonal array tuning method \cite{zhang2019deep}) to efficiently learn the optimal hyper-parameter settings.

In addition, the seizure samples in the training stage are noisy because they may have different seizure EEG channels. For instance, some samples may have 12 seizure channels while some have 21 seizure samples, however, all of them are labelled as 'seizure'. This phenomenon corrupts the training data quality, disturbs the data distribution, and has a negative impact on the interpretation. Nevertheless, the dataset does not contain enough samples in which all the channels are seizure. Thus, to keep the balance between the number of seizure and normal samples, we have to set a lower threshold as `12 seizure channels' to form more seizure samples (Section~\ref{sub:dataset}).

\section{Conclusion}
\label{sec:conclusion}
In this paper, we propose a novel and generic deep learning framework aiming at patient-independent epileptic seizure diagnosis. The proposed approach refines the seizure-specific representation by eliminating the inter-subject noise through adversarial training. Moreover, we involve the attention mechanism to learn the contribution of each EEG channel in the epileptic seizure detection, which empowers our method with great explainability. We conduct extensive experiments over a benchmark dataset. The experimental results show that the proposed approach not only outperforms a batch of state-of-the-art methods by a large margin but also indicates low testing latency and highly explainability.

\balance
\bibliographystyle{IEEEtran}


\section{Introduction} \label{S:1}

Many high- and mid-Z nuclei exhibit neutron-induced resonances in the epithermal range (from 1 eV to 100 eV).  Neutron Resonance Transmission Analysis (NRTA) uses the isotopic uniqueness of these resonances in a transmission measurement mode to determine isotopic and elemental compositions of target materials.  Depending on the isotope, most of these interactions are the result of either $(n,\gamma)$, $(n,n')$, or $(n,$fission$)$ processes.  The transmitted spectrum exhibits attenuation dips at energies that uniquely map to the resonance energies of the nuclei in the target material.
In these measurements, a pulsed neutron beam (e.g. via spallation by a proton beam) is used in combination with a detector sensitive to epithermal neutrons, such as a lithium doped scintillator or a boron doped microchannel plate (MCP)~\cite{ref:tremsin}.  The time between the neutron generation pulse and the time at which the neutron is detected provides the TOF of individual neutrons, enabling a reconstruction of their energy and thus achieving spectroscopic radiography.  In combination with known cross section data, the transmission spectrum can be used to precisely and simultaneously determine the identity and the areal density of the multiple isotopes that make up the target object.  A pixelated detector can further provide two-dimensional coordinate information~\cite{ref:losko,ref:tremsin,vogel2017neutron}, which can then be used to determine the two-dimensional areal density of the isotopes, thus producing isotopic imaging of the target.

Other neutron transmission techniques exist as well.  In particular, work by~\citeauthor{ref:blackburn2007fast}\cite{ref:blackburn2007fast} and ~\citeauthor{perticone2019fast}~\cite{perticone2019fast} has shown that $\sim$MeV fast neutrons can be used to map the composition of a number of low-Z elements, such as carbon, oxygen, and nitrogen.   Most mid-Z to high-Z elements however do not exhibit experimentally distinguishable resonances at $\sim$MeV energies, however.  Other radiographic methods are available to determine the hydrogenous content of an object~\cite{rahon2016spectroscopic,rahon2020hydrogenous,sowerby}.  Additionally, other applications of epithermal neutrons exist.  In particular, Boron Neutron Capture Therapy uses intense fluxes of epithermal neutrons in oncological applications, as described in detail by \citeauthor{doi:harling}\cite{doi:harling} and \citeauthor{forton2009overview}\cite{forton2009overview}.

In recent years NRTA has been used in many areas of nuclear engineering~\cite{ref:chichester2012assessing,ref:chichester2012JNMM,bourke2016non} and archaeology~\cite{schillebeeckx2015neutron,andreani2009novel}.
Epithermal NRTA has also been proposed as a platform for warhead verification in arms control treaties~\cite{ref:hecla2018nuclear,ref:engel2019}. Thus the isotopic specificity of NRTA gives it a very broad range of applicability. For historical and practical reasons, NRTA experiments have primarily taken place at facilities which were designed and built for nuclear data measurements.  Most of these facilities possess kilometer-long accelerators, significantly limiting the application of an otherwise powerful technique. The utility of NRTA can thus be significantly extended by the design of compact, possibly mobile, and relatively low cost platforms. Pulsed deuterium-tritium (DT) neutron sources, augmented by carefully optimized moderation, shielding, and collimation, can produce a compact and precise beam, suitable for mobile applications.  When combined with a shielded neutron detector, such a beam can constitute a novel instrument for isotope-specific radiography of objects composed of mid- and high-Z elements. In this Monte Carlo-based study we focus on the optimization of compact DT-based platforms, showing the feasibility and practicality of such platforms when it comes to materials analysis.  This technique can in the future benefit from ongoing research in high intensity, fast-pulsed DT and deuterium-deuterium (DD) sources~\cite{podpaly2018environment,ref:starfire}.  In this work, we perform Monte Carlo based feasibility and optimization studies of such an instrument, presenting an optimal configuration for high-Z isotopic imaging and demonstrating the feasibility of such a platform for mobile nuclear security applications. While recent work was performed on the development of compact thermal and epithermal sources~\cite{hasemi2017optimization, paradela2016neutron}, no design studies have been performed on the development of a compact and complete NRTA system.

\section{Monte Carlo Simulations} \label{S:2}
Monte Carlo simulations were performed to model two effects: the impact of moderation process on the accuracy of the TOF energy reconstruction; the spectral sensitivity to changes in isotopic composition of the target.  The purpose of the simulations was to independently model shielding configuration and moderator thickness and to optimize these criteria.  Optimization of the moderator improved the neutron epithermalization efficiency and optimization of the shielding geometry improved the TOF energy reconstruction and reduced the background.

\subsection{Simulations}
Simulations of viable shielding geometries were performed with a Geant4 based Monte Carlo application~\cite{allison2016_geant4,ref:grasshopper}.  Two series of simulations were performed to optimize both the moderator thickness and shielding geometries for the compact source.  All simulations used a 14.1 MeV isotropic neutron source to simulate neutrons produced by a compact DT source.  Neutron generation was modeled as an instantaneous pulse.   At later stages in the optimization analysis the time distribution of the output was smeared with a square function of a width 5~$\mu$s, which is the shortest possible pulse width of the commercially available Thermo Electron P325 neutron generator~\cite{ref:p325}.  In order to increase computational efficiency, a dedicated and separate simulation was performed to determine the optimal moderator thickness.
The shielding optimization was performed by simulating a realistic experimental geometry, composed of: a concrete room filled with ambient air; the DT source shielding, collimation, and moderation; a tungsten-silver target as a test object; and an idealized detector with 100\% detection efficiency, surrounded by its own shielding to reduce room return.  Since the optimal moderator thickness was determined in the separate epithermalization simulations (see Section~\ref{sec:moderator_optimization}),
these new simulations only iterated over the shielding and collimation geometries.  Moderator thicknesses were held constant at 10 cm for all geometries.  To suppress room return and collimator in-scattering, shielding geometry was modeled as a sphere of  75 cm radius consisting of of 5\% borated polyethylene.  This allowed to attenuate the off-axis neutron flux in order to reduce thermal and epithermal backgrounds.  To minimize the impact of in-scattering, the collimator walls and detector shielding was modeled to consist of boron carbide.  A schematic of the optimized shielding geometry is shown in Figure \ref{fig:portgeo}.  The target was chosen to consist of a mix of tungsten, molybdenum, and silver. These elements have very clear resonances in the $\sim$eV range, and thus can be used for an optimization analysis applicable to most actinides (e.g. $^{235}$U, $^{238}$U, and most plutonium isotopes).  Also, both tungsten and silver are easily available in the necessary quantities for performing follow-up experimental studies.  The shielding design included cadmium filters in the beam path, as a way of eliminating the thermal flux exiting the moderator and impinging the detector,  via the $^{113}$Cd(n,$\gamma)^{114}$Cd reaction "threshold" at 0.5 eV. See Table \ref{table:parameters} for the dimensions of the system model.

\begin{figure}[h]
\centering\includegraphics[width=0.8\linewidth]{figures/Experimental-Setup-New.pdf}
\centering\includegraphics[angle=270,width=0.8\linewidth]{figures/DT-rendering.pdf}
\caption{Top: schematic of a compact NRTA platform using a DT source (not to scale).  Bottom: a Geant4 rendering of the simulation.  Due to its insignificant contribution to the process of moderation the DT generator is not included in the actual simulation.}
\label{fig:portgeo}
\end{figure}

\subsection{Moderator Optimization} \label{sec:moderator_optimization}

The process of neutron moderation introduces uncertainty into the TOF reconstruction.  An optimal moderator maximizes the epithermalization efficiency and minimizes the spread in TOF.  The epithermalization efficiency was defined as the fraction of incident 14.1 MeV neutrons that were slowed down into the energy range 1-30 eV and exited the rear surface of the moderator.  Figure \ref{fig:efficiencycurve} shows a plot of overall epithermalization efficiency vs. moderator thickness.  The ideal moderator thickness was found to be approximately 10 cm, with minimal change between 7 and 11 cm.  The initial increase in efficiency with increasing moderator thickness is primarily a result of additional scattering of 14.1 MeV neutrons.  The decrease in efficiency beyond 11 cm occurs as down-scattering from the epithermal to thermal ranges begins to dominate. The maximum epithermalization efficiency was determined to be $(9.3\pm 0.1)\times 10^{-5}$.  Even for the optimized moderator thickness, the thermal flux at the rear surface of the moderator is two orders of magnitude higher than the epithermal flux since increasing cross sections at lower energies lead to a build-up of thermal neutrons in the system.  Most of these thermal neutrons can be removed using the cadmium filters.

\begin{figure}[h]
\centering\includegraphics[angle=270,width=0.8\linewidth]{figures/efficiencycurve.pdf}
\caption{Epithermalization efficiency of various moderator thicknesses.  The maximum efficiency occurs at approximately 9-10 cm with a stable range between 7 and 11 cm.  Competing processes of epithermalization and thermalization lead to the local maximum.}
\label{fig:efficiencycurve}
\end{figure}

The simulation output was also analyzed to understand the effect of moderation on the TOF energy reconstruction.  This analysis showed that changes in moderator thickness had little effect on the accuracy of TOF reconstruction in the epithermal regime in the range of flight distances considered.  Studies of moderation times indicates that, depending on energy, the spread in times are in the range of 0.1-1 \textmu s.  Thus, using appropriate digitizer or time-to-digital coincidence electronics with the time resolution of $\ll$ \textmu s is necessary. CAEN V1720 digitizers, as an example, have sampling at the period of 4 ns.  Such sampling will have an insignificant effect on TOF resolution.

\begin{comment}
As can be seen in \ref{fig:moddistribution}, the variation in moderation time increases as neutron energy decreases, likely due to the increase in potential combinations of scattering events required for further epithermalization.
\begin{figure}[h]
\centering\includegraphics[width=0.8\linewidth]{figures/moderation_time_dist.pdf}
\caption{Distributions of simulated moderation times for neutrons from 14.1 MeV to the stated final neutron energy groups.  In the range 0.1-10 eV, as neutron energy decreases, the variation in moderation time increases roughly proportional to the increase in time of flight. Note that the x-axis is on log scale.}
\label{fig:moddistribution}
\end{figure}

The overall effect of thermalization time on the energy reconstruction accuracy can be further decreased by increasing the distance between the moderator and the detector apparatus.  This is also true for other effects which lead to uncertainty in the time-of-flight of epithermal neutrons and hence uncertainty in their energy.  It can be shown that the uncertainty in TOF relates to uncertainty in energy via $\delta E/ E = 2\delta t_0/(t-t_0)$, where $E$ is energy, $t$ is the TOF, and $t_0$ is the time-zero, and $\delta t_0$ is the uncertainty in time-zero~\cite[Suplemental Note 1]{ref:hecla2018nuclear}.  Increasing the length of the flight path increases the energy resolution of the system by the relation $t-t_0=d\sqrt{2E/m}$, where $d$ is the flight path and $m$ is neutron energy.  This however also decreases the overall counting efficiency due to solid angle losses by $\sim 1/d^2$.  A more detailed discussion of pulse width and flight path distance optimization analysis is discussed in Section~\ref{sec:pulsewidth}.
\end{comment}

\subsection{System Geometry}

The configuration described in Figure~\ref{fig:portgeo} succeeded in both minimizing the impacts of room return and in-scattering on the transmission spectrum.  Designs that incorporated less High Density Polyethylene (HDPE) or eliminated the shell suffered significant leakage from the outer surface of the absorber-moderator sphere, resulting in high intensities of room return and epithermal background.  We thus omit these configurations from this discussion.

\begin{table}[h]
\centering
\begin{tabular}{l l}
\hline
\hline
\textbf{Parameter} & \textbf{Dimension} \\
\hline
Distance from DT generator to moderator surface & 5 cm \\
Moderator thickness & 10 cm \\
Distance from surface of moderator to detector & 1.845 m \\
Collimator diameter & 20 cm \\
Outer radius of borated polyethylene  & 75 cm \\
Thickness of boron carbide plating & 1 cm \\
Thickness of detector shielding & 9 cm \\
Length of detector shielding  & 80 cm \\
\hline
\hline
\end{tabular}
\caption{Key dimensions of portable, epithermal DT source. Neutron flight path can be increased or decreased depending on needs of measurement.}
\label{table:parameters}
\end{table}

The detector collimator was also plated with boron carbide to minimize the probability that a neutron scattered off or through the surface of the collimator would be measured by the detector.  With this shielding, any neutron that does not enter the collimator is absorbed by the borated absorber-shield.  These off-axis neutrons have a negligible effect on the room return background of the system, as can be seen in the discussion below.
The geometric model of the detector collimator/shield played an important role in the accuracy of the TOF reconstructions.  While no epithermal or thermal neutrons scattered from the concrete walls were observed to reach the detector, a simple cylindrical collimator geometry resulted in a significant number of detected epithermal neutrons which underwent in-scatter from the inner walls of the collimator.  For this simple geometry the line-of-sight (LoS) neutrons, i.e. neutrons that did not undergo inscatter and reached the detector directly from the source constituted only 85\% of the total.  By comparison when using a polycone geometry, described schematically in Fig.~\ref{fig:portgeo}, this number was increased to $(95.4 \pm 0.2)$\%, thus significantly improving the accuracy of the TOF technique.  The fraction of LoS neutrons can be further improved via future optimizations of the inner surfaces of the shielding.

\subsection{Target} \label{sec:target}
Simulated targets consisted of metallic cylinders 10 cm in diameter and between 5 and 10 mm in width.  Simulations were performed for five targets consisting of mixes of silver, molybdenum, and tungsten with varied molar ratios.  Mixes of silver and tungsten were simulated at 10\%, 20\%, and 30\% silver molar abundance.  Mixes of molybdenum and tungsten were simulated at 50\% and 75\% molybdenum molar abundance.  For brevity, isotopic compositions will be referenced using their composite elements followed by their atomic percent.  Using this convention, the mixes are referenced as Ag10:W90, Ag20:W80, Ag30:W70, Mo50:W50, and Mo75:W25, respectively. All materials were simulated with their natural isotopic abundances.

\section{Results} \label{sec:results}
The accuracy of the energy reconstruction and the isotopic sensitivity were both characterized using a $\chi^2$ test.  The TOF energy reconstructions were determined for the energy range of 1-15 eV using a histogram bin width of 0.1 eV.
The goals of this study were to characterize the following:
\begin{enumerate}
    \item Accuracy of spectral reconstruction - determined by comparing the spectra of the TOF energy reconstruction and the actual neutron energy in an idealized detector.  Lower $\chi^2$ values indicated better reconstruction with an ideal reconstruction having a $\chi^2$ value of zero.
    \item Isotopic sensitivity - determined by directly comparing the TOF energy reconstructions for targets of differing isotopic composition.  Higher $\chi^2$ values indicated higher sensitivity.
\end{enumerate}

\subsection{TOF Energy Reconstruction}

In all MC simulations, the distance from the surface of the moderator to the detector was set to 1.845~m as a balance between TOF reconstruction precision and geometrical efficiency.  TOF energy calculations used a flight path length of 1.87~m because it resulted in more accurate reconstructions.  Further analysis confirmed that the mean flight path for incident LoS epithermal neutrons was approximately 1.87 meters, since most epithermal neutrons are epithermalized between 2 and 3 cm before the rear surface of the moderator.

Figure \ref{fig:TOFrec} shows a plot of the energy spectrum of the detected neutrons and the TOF reconstruction for a silver-tungsten mix (at a 10\% silver molar abundance).  The absorption lines due to the silver resonance (5.2 eV) and tungsten resonances (4.1, 7.7 eV) are clearly observed in both the actual and TOF-reconstructed energy spectra.  Table \ref{table:tof-rec} shows the $\chi^2$ for comparisons for various targets.  The value of $\chi^2/$DoF does not exceed unity for any of the materials or isotopic compositions indicating a reasonably good, but not ideal, energy reconstruction.

\begin{figure}[h]
\centering
\includegraphics[angle=0,width=1.02\linewidth]{figures/TOF-fig3.pdf}
\includegraphics[width=\linewidth]{figures/spectrum.pdf}
\caption{Top:  TOF distribution.  Bottom: actual energy and TOF-reconstructed energy comparison. The simulation is for a Ag10:W90 target.  Dashed and solid gray lines denote the expected locations for the tungsten and silver resonances, respectively.}
\label{fig:TOFrec}
\end{figure}

\begin{table}[h]
\centering
\begin{tabular}{l l l}
\hline
Target & \makecell[l]{Target thicknesses\\(mm)} & $\chi^2/$DoF \\
\hline
Ag10:W90 & 5 & 73 / 132 \\
Ag25:W75 & 5 & 39 / 129 \\
Ag30:W70 & 5 & 39 / 133 \\
Mo50:W50 & 10 & 52 / 135 \\
Mo75:W25 & 10 & 70 / 139 \\
\hline
\end{tabular}
\caption{$\chi^2$ values for comparison of actual energy to TOF-reconstructed energy for various targets.  The reconstructed energy spectra for all targets gives a $\chi^2/$DoF below 1.0, demonstrating reasonably accurate energy reconstruction.  The inaccuracy of the reconstruction is primarily caused by the smearing of the neutron pulse time distribution due to the process of moderation.  The DoF changes for different spectral comparisons due to the exclusion of bins with zero counts.}
\label{table:tof-rec}
\end{table}

\subsection{Isotopic Sensitivity}
One of the main advantages of NRTA is its ability to identify an object's material composition and to enable isotopic-geometric comparisons between two objects, such as in a warhead verification exercise~\cite{ref:hecla2018nuclear,ref:engel2019}.  The system's performance in the latter scenario is used as the basis for its characterization.  The spectral transmission comparisons are performed between objects of different isotopic compositions, enabling characterization of the technique's material discrimination capability.

The isotopic sensitivity of the system was characterized by comparing the TOF reconstructions of varied materials (Ag10:W90, Ag20:W80, Ag30:W70, Mo50:W50, and Mo75:W25).  Molar abundances and thicknesses were chosen to optimize transmission probabilities and resonance reconstruction, as well as to allow for practical follow-up experiments.  Figure \ref{fig:spectracomp} shows the plots of spectral comparisons for a subset of test scenarios.  To evaluate the system's ability to distinguish between targets of different isotopic composition, the transmitted spectra were compared using a $\chi^2$ test.  Table \ref{table:spectra-sep} shows the $\chi^2/$DoF for the TOF-based spectral comparisons of various isotopic configurations.  A comparison of two simulations for the 10\%/90\% Silver-Tungsten target (with different random seeds) was used as a control.  Materials of various isotopic composition were compared with each other directly.  In the last two columns of the table we also present the results of simulations where the neutron pulse was smeared with a 5 $\mu$s square function, to emulate the typical duty of the most commonly available neutron generators, such as P325~\cite{ref:p325}.  All simulation results were scaled to $(1.9\pm0.3)\times10^{10}$ total 14.1 MeV neutron generation events, enabling direct comparison.  The comparisons show a clear differentiation between objects of different isotopic composition, even when the neutron pulse is smeared by 5 $\mu$s.  For a typical DT source, the above sample population would correspond to a total measurement time of approximately 30 minutes, assuming an isotropic source output of approximately $10^9$~neutrons/s.

\begin{table}[h]
\centering
\begin{tabular}{l l l l l l l}
\hline \hline
Material & Ratio1 &  Ratio2 & $\chi^2/$DoF & P-value & \makecell[l]{$\chi^2/$DoF \\ 5$\mu$s}&\makecell[l]{P-value\\ 5$\mu$s } \\
\hline
Ag-W & 10:90 & 10:90 & 102 / 85 & 0.10 & 82 / 85 & 0.57\\
Ag-W & 10:90 & 20:80 & 189 / 83 & $3\times10^{-10}$  & 152 / 85 & 0.00\\
Ag-W & 10:90 & 30:70 & 324 / 81 & $3\times10^{-30}$ & 301 / 83 & 0.00\\
Ag-W & 20:80 & 30:70 & 128 / 81 & $7\times10^{-4}$  & 113 / 82 & 0.01\\
Mo-W & 50:50 & 75:25 & 216 / 151 & $4\times10^{-4}$ & 225 / 154 & 0.00\\
\hline \hline
\end{tabular}
\caption{$\chi^2$ values for spectral separation analysis.  Column 1 denotes the isotopes present.  Columns 2 and 3 denote the molar ratios being compared.  Columns 4 and 5 denote the $\chi^2$ test statistic and the corresponding p-value for an instantaneous pulse and a detector at a distance of 1.87~m.  Columns 6 and 7 denote the same values for a DT pulse of 5$\mu$s and a detector placed at 5~m.  Separation increases monotonically as the molar abundance of silver in the second mix increases.  The control comparison has a $\chi^2$ statistic of approximately order unity.  The silver-tungsten and molybdenum-tungsten comparisons are restricted to energy ranges of 1.5-10 eV and 1.5-15 eV, respectively.}
\label{table:spectra-sep}
\end{table}

The control comparison (Ag10:W90 vs. Ag10:W90) has a reduced $\chi^2$ statistic of approximately order unity which is consistent with two samples drawn from the same distribution.  As the molar ratio of the silver is increased in the mix, the $\chi^2$ value compared to Ag10:W90 increases monotonically.  Figure \ref{fig:spectracomp} shows the spectra for the first four comparisons between mixed compositions of silver-tungsten.  The differences between the reconstructions are clearly visible in the spectra.  The separation between the 20\% and 30\% samples deteriorates due to limited statistics within the silver and tungsten resonance structures.  This is reflected by the relatively low $\chi^2$ value.  The Mo50:W50 and Mo75:W25 separation is also limited, and this is due primarily to low statistics within the resonances combined with a similar overall cross section for both configurations.  In both cases however, the spectra can be clearly differentiated.  In all of these simulations $2.5\times10^9$ neutrons were generated.  With a source that produces neutrons at the rate of $10^9$s$^{-1}$, individual runs would require 25 seconds.  Thus the compact configuration described in this study can distinguish between different isotopic compositions with total measurement times of approximately one minute. This time could be further reduced as new, more intense sources of neutrons become available.

\begin{figure}[h]
    \centering
    \includegraphics[angle=270, width=\linewidth]{figures/spectra_comp.pdf}
    \caption{Spectral Comparisons of Various Targets.  Comparisons of TOF energy reconstructions for transmitted neutrons through different targets are presented. Visual inspection shows clear differentiation between varied isotopic mixes.  Dashed lines represent tungsten and silver resonances.  Table~\ref{table:spectra-sep} lists the results of the $\chi^2$ tests for these comparisons, showing a clear quantitative difference.}
    \label{fig:spectracomp}
\end{figure}

\subsection{Effects of Neutron Generation Pulse Width} \label{sec:pulsewidth}

\begin{figure}[h]
    \centering
    \includegraphics[width=\linewidth]{figures/AgW19_5us_5m_semilog_4610.pdf}
    \caption{Effect of pulse width on TOF Reconstruction Quality. Comparison of Ag10:W90 transmission spectra for an instantaneous neutron generation pulse (blue) and a 5.0 $\mu$s pulse width (red) at 5.0~m TOF distance. Note the slight smearing of the narrow resonance at 7.7 eV and the difference in spectra at higher energy.  These effects are not seen in the spectra in Figures \ref{fig:TOFrec} and \ref{fig:spectracomp} produced without a finite pulse width.}
    \label{fig:pulsewidthspectrum}
\end{figure}

Up to this point, it was assumed that the neutron source generates a neutron beam instantaneously. While fast, sub-microsecond DT systems exist, the most commonly commercially available sources typically produce a $\sim5$ $\mu$s pulse.  To account for this circumstance the output of the simulation was modified by adding a uniformly random 0-5 $\mu$s offset to the TOF of every neutron event.  This pulse time introduces uncertainty into the $t_0$ of any given neutron, and as a result introduces an energy-dependent error in the TOF-reconstructed neutron energy (see Figure \ref{fig:pulsewidthspectrum}).  As with error contributions from the moderation time, this spectral smearing effect is more prominent at higher energies where the time of flight is shorter. This effect could be mediated by moving the detector further from the source and thus increasing the time of flight.  For this analysis the detector was placed at a TOF distance of 5~meters, with an equivalent neutron flux incident on the detector.  The last two columns of Table \ref{table:spectra-sep} show the reduced $\chi^2$ values for the various targets for a time-smeared pulse.  While a significant and expected deterioration of the statistical separation is observed, all spectra can be distinguished from one another.  This indicates that while the system's analyzing power is reduced, it is still capable of performing material differentiation.  For a realistic system, the optimal TOF distance will depend on the isotopic composition and thickness of the targets.

\section{Conclusions}

The proposed neutron resonance transmission analysis (NRTA) system is capable of highly accurate TOF energy reconstruction and is highly sensitive to changes in isotopic composition of the target.  Errors in TOF reconstruction are small and due primarily to pulse smearing within the moderator and in-scattering from the detector shielding.  The polycone design significantly reduces these effects.  The reconstruction might be further improved by increasing the flight path and further optimizing the polycone design. The sensitivity to isotopic composition makes this methodology a powerful tool for on-site measurements and characterizations of nuclear and non-nuclear materials alike.  With further work the system may be able to characterize and identify unknown samples of material with epithermal resonance structures.

The use of compact neutron sources has the potential to significantly expand the utility and applicability of NRTA in a variety of fields.
While this study has focused on the optimization of DT-based sources alone, DD-based and linac-based photodisintegration sources may also provide viable methods for compact neutron generation.  Such sources may require less shielding and enable more precise TOF determination via narrower pulses.  For example, new generation of high intensity DD sources can achieve ${\lesssim}$2$\mu$s pulses~\cite{ref:starfire}, allowing for high TOF precision and significantly reduced shielding.  Furthermore, a 3 MeV linac can produce fast-pulsed neutrons in the $\sim$200 keV range, requiring minimal neutron shielding and collimation.  Future research should explore these alternatives.

Finally, this study did not take into the account the background contributions of gammas from neutron capture on hydrogen and boron.  Past work by \citeauthor{wang2016improved}\cite{wang2016improved} using $^6$Li glass detectors have indicated that suppression of gamma counts by $\mathcal{O}(10^6)$ is feasible.  Furthermore, the pulsed mode of the source operation can allow for further reduction of gamma backgrounds: as capture gammas come primarily from thermal neutron capture their arrival times are more spread out in time than that of the 1-20 eV neutrons.  Nevertheless, the careful choice of the detector and the data analysis algorithms is important for controlling the background and should be part of future experimental efforts.

\section*{Acknowledgements}

This work is supported in part by Department of Energy Award DE-NA0002534, as part of the NNSA Consortium of Verification Technologies (CVT). The authors would like to thank their colleagues in CVT and MIT's Nuclear Science and Engineering for encouragement and advice.