QuPWM: Feature Extraction Method for MEG Epileptic Spike Detection

TL;DR

This paper introduces QuPWM, a feature extraction method combined with SVM classification for automatic epileptic spike detection in MEG signals, achieving high accuracy and efficiency over a balanced dataset.

Contribution

The paper presents a novel feature extraction approach using Position Weight Matrix with quantization for MEG spike detection, improving accuracy and reducing feature size.

Findings

Achieved up to 98% accuracy in spike detection.

Effective feature extraction with PWM improves detection performance.

Method reduces feature vector size for efficient processing.

Abstract

Epilepsy is a neurological disorder classified as the second most serious neurological disease known to humanity, after stroke. Localization of the epileptogenic zone is an important step for epileptic patient treatment, which starts with epileptic spike detection. The common practice for spike detection of brain signals is via visual scanning of the recordings, which is a subjective and a very time-consuming task. Motivated by that, this paper focuses on using machine learning for automatic detection of epileptic spikes in magnetoencephalography (MEG) signals. First, we used the Position Weight Matrix (PWM) method combined with a uniform quantizer to generate useful features. Second, the extracted features are classified using a Support Vector Machine (SVM) for the purpose of epileptic spikes detection. The proposed technique shows great potential in improving the spike detection accuracy and reducing the feature vector size. Specifically, the proposed technique achieved average accuracy up to 98\% in using 5-folds cross-validation applied to a balanced dataset of 3104 samples. These samples are extracted from 16 subjects where eight are healthy and eight are epileptic subjects using a sliding frame of size of 100 samples-points with a step-size of 2 sample-points

Tables

Table 1. Table 1: Example of a position-wight matrix P ​ W ​ M + 𝑃 𝑊 superscript 𝑀 PWM^{+} containing the frequencies of five levels [ q 1 , q 2 , … , q 5 ] subscript 𝑞 1 subscript 𝑞 2 … subscript 𝑞 5 [q_{1},q_{2},\dots,q_{5}] .

$n^{th}$ position	Quantization levels
	$q_1$	$q_2$	$q_3$	$q_4$	$q_5$
1	1104	257	331	384	372
2	1209	213	271	291	324
3	1217	245	298	336	336
4	987	263	309	409	474
5	941	253	383	442	532
6	1273	243	308	317	367
7	1494	208	243	261	245
8	1301	264	299	364	388
9	61	101	306	829	1331
10	144	88	266	599	1133
..	..	..	..	..	..

Table 2. Table 2: Example of position-wight matrix P ​ W ​ M j + 𝑃 𝑊 superscript subscript 𝑀 j PWM_{\textbf{j}}^{+} containing the frequencies of di-mer motif.

$n^{th}$ position	The di-mer motifs (j)
	$q_1{q}_1$	$q_1{q}_2$	$q_1{q}_3$	$q_1{q}_4$	$q_2{q}_1$	$q_2{q}_2$	$q_2{q}_3$	$q_2{q}_4$	..
1	15	13	11	14	17	13	11	8	..
2	10	10	12	5	9	11	16	6	..
3	0	0	0	0	0	0	0	0	..
4	0	0	0	0	0	0	0	0	..
5	8	10	8	12	7	14	3	12	..
6	298	290	310	294	303	293	322	299	..
7	93	93	101	112	94	109	89	116	..
8	0	0	0	0	0	0	0	0	..
..	..	..	..	..	..	..	..	..	..

Table 3. Table 3: Comparison of the feature-size of proposed feature generation methods

Method	Parameters	Classifier	Feature size	Accuracy	Sensitivity	Specificity
Raw data	–	–	2400	60.06	78.27	41.84
PWM	M=8, $𝐫$=0.18	SVM	2	92.10	91.10	93.10
mPWM	M=10, $𝐫$=0.15 k-mers: $𝐣=\{1,2\}$	SVM	220	98.23	98.07	98.39

Table 4. Table 4: Comparison of the proposed methods with existing methods reported recently in [ 9 ]

Method	# Subject	Feature	Accuracy	Sensitivity	Specificity
ICA [7]	4	Independent Components Analysis	-	86.91	81.19
CSP-LDA [8]	20	Common Spatial Patterns and Linear Discriminant Analysis	-	86.14	90.38
AT-DTW [9]	30	Amplitude Thresholding and Dynamic Time Warping	-	92.45	95.81
SCSA-RF [10]	16	Semi-Classical Signal Analysis Random Forest classifer	94.33	93.68	95.08
PWM-SVM	16	Position Weight Matrix SVM classifer	92.10	91.10	93.10
mPWM-SVM		motif-based Position Weight Matrix SVM classifer	98.22	98.06	98.38

Table 5. Table 5: The 5-Folds cross-validation performance with L = 70 𝐿 70 L=70 for different subjects combination.

# Subjects	# Electrodes	size	L	step	Method	Parameters	Classifier	Feature size	Accuracy	Sensitivity	Specificity
16	24	4862	70	2	Row_Data	–	SVM	1680	55.44	59.97	48.92
					PWM	M=8, $𝐫$=0.18		2	54.34	85.12	31.97
					mPWM	M=12, $𝐫$=0.12 , $𝐣=\{1,2\}$		312	66.98	81.52	52.43
14	24	4730	70	2	Row_Data	–	SVM	1680	55.48	63.13	47.82
					PWM	M=12, $𝐫$=0.12		2	58.54	85.12	31.97
					mPWM	M=12, $𝐫$=0.12, $𝐣=\{1,2\}$		312	66.98	81.52	52.43
12	24	4048	70	2	Row_Data	–	SVM	1680	56.62	64.28	48.96
					PWM	M=8, $𝐫$=0.18		2	65.02	74.80	55.23
					mPWM	M=12, $𝐫$=0.12, $𝐣=\{1,2\}$		312	68.36	77.62	59.09
10	24	3464	70	2	Row_Data	–	SVM	1680	59.12	70.95	47.28
					PWM	M=8, $𝐫$=0.18		2	64.72	72.40	57.05
					mPWM	M=12, $𝐫$=0.12, $𝐣=\{1,2\}$		312	70.61	80.89	60.34
8	24	3260	70	2	Row_Data	–	SVM	1680	58.31	72.02	44.60
					PWM	M=10, $𝐫$=0.15		2	63.87	84.29	43.44
					mPWM	M=10, $𝐫$=0.15, $𝐣=\{1,2\}$		220	72.15	84.05	60.25
6	24	2728	70	2	Row_Data	–	SVM	1680	60.67	74.56	46.78
					PWM	M=8, $𝐫$=0.18		2	64.66	81.00	48.30
					mPWM	M=12, $𝐫$=0.12, $𝐣=\{1,2\}$		312	73.46	87.90	59.02
4	24	2644	70	2	Row_Data	–	SVM	1680	60.10	76.63	43.57
					PWM	M=12, $𝐫$=0.12		2	94.82	94.70	94.93
					mPWM	M=8, $𝐫$=0.18, $𝐣=\{1,2\}$		144	99.09	98.26	99.92

Table 6. Table 6: The 5-Folds cross-validation performance with L = 100 𝐿 100 L=100 for different subjects combination

# Subjects	# Electrodes	size	L	step	Method	Parameters	Classifier	Feature size	Accuracy	Sensitivity	Specificity
16	24	3102	100	2	Row_Data	–	SVM	2400	60.06	78.27	41.84
					PWM	M=8, $𝐫$=0.18		2	92.10	91.10	93.10
					mPWM	M=10, $𝐫$=0.15, $𝐣=\{1,2\}$		220	98.23	98.07	98.39
14	24	3079	100	2	Row_Data	–	SVM	2400	59.21	76.56	41.85
					PWM	M=8, $𝐫$=0.18		2	92.76	92.40	93.11
					mPWM	M=10, $𝐫$=0.15, $𝐣=\{1,2\}$		220	98.38	98.31	98.44
12	24	2819	100	2	Row_Data	–	SVM	2400	61.44	78.87	44.00
					PWM	M=8, $𝐫$=0.18		2	91.98	90.92	93.04
					mPWM	M=10, $𝐫$=0.15, $𝐣=\{1,2\}$		220	98.44	98.30	98.58
10	24	2503	100	2	Row_Data	–	SVM	2400	62.44	84.03	40.85
					PWM	M=8, $𝐫$=0.18		2	94.17	94.49	93.84
					mPWM	M=10, $𝐫$=0.15, $𝐣=\{1,2\}$		220	98.52	98.88	98.16
8	24	2443	100	2	Row_Data	–	SVM	2400	61.69	82.98	40.37
					PWM	M=8, $𝐫$=0.18		2	95.13	94.85	95.41
					mPWM	M=10, $𝐫$=0.15, $𝐣=\{1,2\}$		220	99.55	99.10	100.00
6	24	2263	100	2	Row_Data	–	SVM	2400	62.40	81.18	43.59
					PWM	M=8, $𝐫$=0.18		2	95.27	95.14	95.40
					mPWM	M=10, $𝐫$=0.15, $𝐣=\{1,2\}$		220	99.51	99.03	100.00
4	24	2246	100	2	Row_Data	–	SVM	2400	63.00	82.19	43.81
					PWM	M=8, $𝐫$=0.18		2	94.61	94.12	95.10
					mPWM	M=10, $𝐫$=0.15, $𝐣=\{1,2\}$		220	98.04	96.97	99.11

Equations

$$Q(n)=q_i\text{mbox}wherei=1,2,\dots ,M$$

$$|X(n)-\mu |\le 3\sigma$$

$$\mu =\frac{1}{N_s}\sum _{n=1}^{N_s}{\mu }_n\text{mbox}and\sigma =\frac{1}{N_s}\sum _{n=1}^{N_s}{\sigma }_n,$$

$$PW{M}^{+}[i,q_j]=\sum _{s=1}^{N_s^{+}}\xi (Q_s^{+}(i),q_j)$$

$$PW{M}^{-}[i,q_j]=\sum _{s=1}^{N_s^{-}}\xi (Q_s^{-}(i),q_j)$$

$$\xi(a,b)=\begin{cases}1&~{}\mbox{if the $a=b$ }\\ 0&~{}elsewhere\end{cases}$$

$$Scor{e}^{+}(\text{Q})=\sum _{n=1}^N\frac{PW{M}^{+}\left[n,Q(n)\right]}{{\sum }_{j=1}^{M}PW{M}^{+}[n,q_j]}$$

$$Scor{e}^{-}(\text{Q})=\sum _{n=1}^N\frac{PW{M}^{-}\left[n,Q(n)\right]}{{\sum }_{j=1}^{M}PW{M}^{-}[n,q_j]}$$

$$m_{q_{1}}(n)=\begin{cases}1&~{}\mbox{if $Q(n)=q_{1}$ }\\ 0&~{}elsewhere\end{cases}$$

$$m_{q_{1}q_{3}}(n)=\begin{cases}1&~{}\mbox{if $Q(n)=q_{1}$ and $Q(n+1)=q_{3}$ }\\ 0&~{}elsewhere\end{cases}$$

$$PW{M}_{\text{j}}^{+}(n)=\sum _{s=1}^{N_s^{+}}\xi (m_{j,s}^{+}(n),1)$$

$$PW{M}_{\text{j}}^{-}(n)=\sum _{s=1}^{N_s^{-}}\xi (m_{j,s}^{-}(n),0)$$

$$\{\begin{array}{l}{\Omega }_1=\{q_1,q_2,q_3,\dots \}\text{mbox}{}^{\prime }mono-mer{s}^{\prime }\\ {\Omega }_2=\{q_1{q}_1,q_1{q}_2,q_1{q}_3,\dots \}\text{mbox}{}^{\prime }di-mer{s}^{\prime }\\ {\Omega }_3=\{q_1{q}_1{q}_1,q_1{q}_2{q}_3,q_3{q}_3{q}_3,\dots \}\text{mbox}{}^{\prime }tri-mer{s}^{\prime }\\ \dots \end{array}$$

$$Scor{e}^{+}(m_j)=\sum _{n=1}^N\frac{m_j(n)PW{M}_{\text{j}}^{+}(n)}{{\sum }_{n=1}^{N}PW{M}_{\text{j}}^{+}(n)}$$

$$Scor{e}^{-}(m_j)=\sum _{n=1}^N\frac{(1-m_j(n))PW{M}_{\text{j}}^{-}(n)}{{\sum }_{n=1}^{N}PW{M}_{\text{j}}^{-}(n)}$$

Full text

K****eywords magnetoencephalography (MEG)  $\cdot$ Position Weight Matrix (PWM)  $\cdot$ Epileptic spike detection  $\cdot$ machine learning

1 Introduction

Epilepsy is not a singular disease entity but a variety of dysfunctions reflecting brain disorders or abnormal electrical activities that may strike patients of all ages [1]. An epileptic seizure occurs when a burst of electrical impulses in the brain exceeds its normal limits. These pulses spread to neighboring areas in the brain, which may create an uncontrolled storm of electrical activity sent to body organs. The electrical impulses could be transmitted to the muscles, causing twitches or convulsions. In particular, epileptic patients may stare blankly for a few seconds during a seizure, while others have uncontrollable jerking movements of the arms and legs. Among the diagnosis tools for epilepsy is the magnetoencephalography (MEG) [2]. The MEG is a recent functional neuro-imaging technology that measures the magnetic activity of the brain. This new technology uses an array of highly sensitive sensors or magnetometers called superconducting quantum interference devices (SQUIDs). However, as the brain magnetic field is very weak compared to surrounding magnetic sources, the MEG measurement needs a shielded room to isolate the patient from the external magnetic fields such as the magnetic field of the earth or the electronic devices. For Epileptic patients, the MEG is used during two main phases of the treatment [3, 4]: first, localizing the region of the brain which produces the abnormal electrical activities that cause the neurological disorder. If the localized region is not responsible for a vital function in the brain such as speech, this region is usually extracted surgically if the medications were not effective for the treatment. Second, assessing surgeries outcomes.

Neurologists often spend hours to manually read the MEG recordings. For this reason, automated interictal epileptic spike detection in MEG signals has attracted research interest over the last decade. Four methods, to the best of our knowledge, have been proposed for spike detection of the MEG signals in literature [5, 6]. The first method uses independent component analysis (ICA).The ICA method is a multi-channel MEG spikes localization method that decomposes spike-like and background components into separate spatial topographies and associated time series [7]. The detection is performed using a thresholding technique. The second method is based on the common spatial patterns (CSP) method and linear discriminant analysis (LDA) (CSP-LDA) [8]. The CSP-LDA performs eigenvalue decomposition of the covariance matrices of the input data to find the most discriminative features. LDA classifier is then employed to perform the spikes detection. The third method is the Amplitude Thresholding and Dynamic Time Warping (AT-DTW) proposed in [9]. This method first uses amplitude thresholding to determine the most likely spiky segments. For the spike detection, the method employs dynamic time warping. The maximum specificity and sensitivity reported in the aforementioned detection methods are 95.8%,and 92.4%, respectively. A more recent detection method has been proposed in [10] where the decomposition of the signal into squared eigenfunctions of the Schrodinger operator was used for detection[11]. This method uses the largest negative eigenfunction in absolute value of the discrete spectrum of the Schrödinger Operator as a feature with the Random Forest (RF) classifier. With this low dimension feature vector this method could achieve a sensitivity and specificity of 93.68% and 95.08%, respectively.

In this work, we propose a novel feature generation method for multi-channel MEG signal which improves patient-independent interictal spike detection. This method is based on the combination of the Position Weight Matrix (PWM) method and uniform quantization scheme. We name this approach QuPWM. This method takes advantage of the efficiency of the PWM method, which is usually used for DNA sequences classification. QuPWM shows great potential in improving the accuracy of the interictal spike detection models.

The paper is organized as follows. Section II includes a description of the MEG dataset and the proposed features extraction process methods based on the PWM method and the used quantization schemes in Section II. Section III presents the obtained results using eight healthy and eight epileptic patients which are presented in Section III. Section IV presents a discussion of the findings, where a conducted sensitivity analysis on the proposed feature extraction method is analyzed with respect to the frame length and number of subjects. Finally, Sections V summarizes our concluding remarks.

1.1 Classification problem definition

An Epileptic seizure happens when a burst of electrical impulses in the brain exceed their normal limits. They spread to neighboring areas in the brain and might create an uncontrolled storm of electrical activity sent to body organs. The electrical impulses can be transmitted to the muscles, causing twitches or convulsions.

For epileptic patient treatment, it is important to localize the region in the brain that has an abnormal activity to extract this part with surgical intervention. After the surgery, the patient needs to take MEG test to ensure that the abnormal brain activities are absent. Up to now, the visual assessment is the "only" way to assess the MEG signals. Indeed, there is no efficient computer based aid to analyse MEG signals. Therefore, there is a need for developing machine learning based algorithms to assist clinicians in quickly and efficiently detecting and predicting abnormalities in MEG signals.

2 Materials and Methods

The proposed methodology consists of three main steps as shown in Figure 1. First, the MEG recordings are split into segments using framing technique. For instance, the first frames of the 24 channels are concatenated to form the first sample, and so on. Then the samples are converted into sequences using a uniform quantization scheme. Second, the quantified samples are mapped to extract two types of PWM-based features. Finally, the classification model SVM is used for spike detection. These three steps are explained in the following sections 2.2 and 2.3.

2.1 Experimental data acquisition and analysis

MEG data were recorded in a shielded room at National Neural Institute- King Fahad Medical City (NNI-KFMC), Riyadh, Saudi Arabia with an Elekta Neuromag system. As the MEG signals are much weaker than normal environmental magnetic noise, the shielded room blocks the majority of environmental magnetic fields so that the magnetic fields generated by the brain can be accurately detected. Elekta Neuromag head system (helmet) contains 102 magnetometer and 204 gradiometer sensors. These sensors are further categorized according to the different brain regions. Clinically, the brain is divided into eight regions; left temporal (LT), right temporal (RT), left frontal (LF), right frontal (RF), left parietal (LP), right parietal (RP), left occipital (LO), and right occipital (RO). Each element of the Elekta Neuromag system is comprised of three sensors, one magnetometer, and two gradiometers. Magnetic brain activity was recorded at a sampling frequency of 1 kHz. MEG data was filtered by tSSS (Spatiotemporal signal space separation) method [14]. The data were then off-line band-pass filtered 1–50 Hz for visual inspection. A total of 18 MEG data segments, each of 15 minutes duration and 26 channels, were taken from 8 epileptic patients and eight healthy patients. These segments are analyzed by specialized neurologists from NNI, KFMC, Riyadh. The neurologists marked the MEG spikes locations, in different brain regions, by visual inspection. The total number of spikes in these recordings is 166. As mentioned earlier, there are 306 sensors to cover the whole head. These sensors are further marked according to the brain regions.

Written informed consent was signed by each participant or responsible adult before they participated in the study. The study was conducted in accordance with the approval of the Institutional Review Board at KFMC (IRB log number: 15-086, 2015).

2.2 Multi-channels MEG signal pre-processing

Multi-channels MEG pre-processing aims essentially to prepare classification samples. The pre-processing consists of two steps:

• •

signal Framing: build equisized samples that combine sliding widows or frames from all the 24 MEG channels. These combined frames will represent the raw data samples used as inputs to the classification which will be quantified in the next step.

• •

signal Quantization: convert the real-values samples to sequences of letters. This step is needed to make the samples suitable for the used feature generation method explained in Sections 2.3.1 and 2.3.2.

2.2.1 MEG signal Framing

The multi-channel MEG signals are segmented into frames using sliding frames using a sliding window of size $100$ sample-points with a step of $2$ sample-points as shown in Figures 1 and 3. For instance, the first 100 sample-points of all channels are concatenated together to form the first sample of size 2400 sample-points of the input raw-data. A binary label or class is assigned to each sample as follows:

• •

Positive class: represent samples that have epileptic spike event occurring in MEG records of the epileptic patients.

• •

Negative class: represent samples from uniformly distributed time-locations in MEG records of the healthy patients.

After signal framing, the positive and negative samples will be quantified using a uniform quantizer as explained in the next section.

2.2.2 MEG signal quantization

In this work we employed the PWM method for feature extraction as explained in Section 2.3.1. However, as the PWM method deals only with sequences, the input MEG samples should be converted into sequences, as shown in Figure 4. For this reason, we used a uniform quantization scheme. The quantization is utilized to convert the real-valued signal $X$ to a sequence $Q$ of different levels $q_{1},q_{2},...,q_{M}$ defined as follows [15, 16]:

[TABLE]

One of the purposes of using quantization in this classification is to illuminate noise effect by mapping a range of value of the input signal, as presented in Figure 4, to a single level $q_{i}$ based on the probability density of this range of values in the dataset. The quantization scheme depends on two major parameters which are the desired number of levels $M$ and the resolution $\mathbf{r}$. For an appropriate choice of these parameters, the probability distribution of the signal is analyzed for both classes.

Indeed, it is known in probability theory and mathematical statistics that the value of a random variable is considered impossible, in other words abnormal, if it lies $3\sigma$ far from its mean $\mu$. This observation is called the 3-sigma rule [17] [18]. Therefore, the most significant variability and randomness that characterizes a Gaussian random variable $X$, where the other value can be considered as outlier values which can be neglected. This fact is observed if :

[TABLE]

where $\mu$ and $\sigma$ are defined to be the mean and standard deviation of the Gaussian random variable $X$.

The standard rule is used to define quantization scheme as shown in Figure 5. It is important to note that the probability distribution of $X$ varies from subject to subject which affects the quantization scheme. Therefore and in order to unify the quantization scheme, the mean $\mu$ and standard deviation $\sigma$ are defined as the average mean and standard deviation of all subjects as follows:

[TABLE]

where $N_{s}$ is the number of subjects. $\mu_{n}$ and $\sigma_{n}$ are the mean and the standard deviation of the data values of the $n^{th}$ subject.

Since the 3-sigma rule indicates that most of the information of a Gaussian random variable are located within $3\sigma+\mu$, the quantization resolution is set to be equal to the standard deviation $\sigma$ as shown in Figure 5 and explained in detail in Algorithm 1.

2.3 QuPWM-based features generation

The proposed Quantization-based position weight matrix (QuPWM) feature generation method is based on combining the position weight matrix (PWM) method with Quantization. This method uses two approaches to generate the PWM features as explained in the next section.

2.3.1 Standard PWM-based features

A position weight matrix (PWM) , also known as a position-specific weight matrix (PSWM) or position-specific scoring matrix (PSSM), has been presented for the first time by Gary in [19]. This method is widely used technique for motifs characterization and discovery in biological sequences such as DNA/mRNA [20, 21, 22, 23, 24]. It showed high potential in sequences characterization and motifs extracting with remarkable binary classification accuracy. For MEG signals, we adapted the PWM-based method to extract relevant features after signals quantization methods. The PWM is based on building two matrices, $PWM^{+}$ derived from the positive training set and the second $PWM^{-}$ represents the negative set. The PWM matrices indicate the significance of each position along the input sequences for each class. For binary classification, the two matrices $PWM^{+}$ and $PWM^{-}$ are defined as follows:

[TABLE]

and

[TABLE]

where $N^{+}_{s}$ and $N^{-}_{s}$ are the total number of sequences in the positive and negative classes. The $\xi(a,b)$ function is defined as follows:

[TABLE]

Table 1 shows an example of a positive matrix $PWM^{+}$. Each column of the matrix represents the frequency of a specific letter $q_{i}$ in a specific position $n$ among all positive sequences.

Then, these two PWMs are used to generate two scores representing the extracted features. The two scores $Score^{+}(\textbf{Q})$ and $Score^{-}(\textbf{Q})$ represent the two probabilities of a given sequence to be in the positive or negative class. In other words, for a given sequence $Q$ to be positive, the score induced by $PWM^{+}$ should be greater than the score induced by $PWM^{-}$. The two scores $Score^{+}(\textbf{Q})$ and $Score^{-}(\textbf{Q})$ of a sequence $Q$ are defined as as follows:

[TABLE]

and

[TABLE]

where $N$ is the number of sample-points in the sample $Q$.

It is very important to mention that the PWMs should be reconstructed only from the training dataset in order not to violate the classification rules.

2.3.2 Motif-based PWM features

In order to extract more advanced patterns, we adopted a new approach which deals with binary sequences, extracted from the original sequence, that represent the presence of a specific motif of one letter or more in the original sequence, as shown in Figure 6. We called this approach the motif-based Position Weight Matrix (mPWM). The main idea is to decompose every sequence into binary sequences reflecting a specific pattern of levels in this sequence such as :mono-mers ($q_{1}$, $q_{2}$, $q_{3}$, etc), di-mers ($q_{1}q_{1}$ , $q_{1}q_{2}$ , $q_{1}q_{3}$, etc), tri-mers ($q_{1}q_{1}q_{1}$, $q_{1}q_{2}q_{3}$, $q_{3}q_{3}q_{3}$, etc). For instance, a sequence of $3$ levels $q_{1}$, $q_{2}$, $q_{3}$ will give: $3$ possible mono-mers, $3x3=9$ possible di-mers combinations, and $3x3x3=27$ possible tri-mers combinations. Figure 6 shows the different k-mer encodings. For instance, the binary sequence of the mono-mers $m_{q_{1}}$ is defined as follows:

[TABLE]

Similarly, the binary sequence of the di-mers $m_{q_{1}q_{3}}$ is defined as follows:

[TABLE]

The binary mapping of k-mers motif extraction is summarized in Algorithm 2.

Similarly to the standard PWM, the binary sequences $m_{j}$, representing the k-mer motifs $\textbf{j}\in\Omega_{k}$, which is defined in Eq. 13, is used to reconstruct multiple pair of PWM matrices $PWM_{j}^{+}$ and $PWM_{j}^{-}$ defined as follows:

[TABLE]

and

[TABLE]

where $\textbf{j}\in\Omega_{k}$, such that $\Omega_{k}$ represents the set of the possible k-mers combination defined as:

[TABLE]

$m_{j,s}$ denotes the motif $m_{j}$ extracted from the $s^{t}h$ sequence $Q_{s}$ . $N^{+}_{s}$ and $N^{-}_{s}$ are the total number of sequences in the positive and negative classes.

Similarly, the two scores $Score^{+}(m_{j})$ and $Score^{-}(m_{j})$ of a k-mer $m_{j}$, extracted from a sequence $Q$, are defined as as follows:

[TABLE]

and

[TABLE]

where $\textbf{j}\in\Omega_{k}$ is defined in Eq. 13.

2.4 Classification models

The input MEG raw data after the pre-procesing step has 1734 spiky samples and 1734 healthy samples from different MEG test sessions of eight healthy and eight epileptic patients.The quantified samples are used to generate QuPWM features and then fed to different classification models. The Support Vector Machine (SVM) model outperforms the other models in 5-folds cross-validation (CV) process.The performance is measured using the average accuracy, sensitivity, and specificity.

3 Results

3.1 Comparison of the PWM-based features

Each of the utilized PWM approaches captures a different property of the spikes characteristic. In order to study how these features perform and compare their performances, the different PWM approaches have been tested on the same dataset, eight healthy and eight epileptic subjects, using the same classifier which is the SVM. The obtained results are shown in Table 5. The results show that the performance of the PWM method can be improved significantly using the di-mer feature. However, the mPWM method needs larger memory and longer time. This limitation can be addressed using parallelization programming and feature selection to select the most important motif to extract.

3.2 Comparison with existing methods

Several detection methods have been proposed for spikes detection using MEG signals. The proposed features generation methods are compared to the reported performance of some recent works [9]. Table 4 shows the proposed feature improves the detection sensitivity and specificity especially using the mPWM features. Table 3 shows that the proposed methods have great potential in improving the detection accuracy and reducing the feature vector size.

4 Discussion

4.1 Choice of the frame length $L$

As explained in Section 2.1, the training and testing data are extracted from the MEG records using a sliding window or frame of size $L$. The choice of the frame-size is motivated by the statistical properties of the spikes time-duration for the different eight subjects. Most of the spikes last almost $95$ sample-points$~{}=95ms$ as shown in Figure 7. Therefore, the frame should be long enough to capture the common characteristics of the different spikes. In order to study how sensitive are the proposed feature generation methods to the input datasets, we performed a sensitivity analysis to study the effect of the frame-length and step-size on the obtained average performance in 5-folds cross-validation using SVM classifier with two frame-lengths $L=\{100,70\}$. Tables 5 and 6 demonstrate that a frame-length $100$ sample-points$=100ms$ gives better detection accuracy. Therefore, the optimal frame-length should be around the mean spike duration for all subjects.

4.2 Subject-independent classification

For epileptic spikes detection, it is important to build subject-independent models which can provide higher specificity and sensitivity for any number of subjects. In order to study the effect of the number of subjects on the detection performance, we studied the combination of randomly combined subjects using SVM classifier and two frame-lengths $L=\{100,70\}$. The QuPWM-based features show a constant performance for most of the combinations using the optimal frame-size 100 sample-points, as shown in Table 6 . This is because the QuPWM-based features extract the common patterns within the dataset regardless of its size.

5 Conclusions

We developed a feature extraction method, called QuPWM, for epileptic spikes detection in MEG signals. This method is based on combining the position weight matrix (PWM) method with digital quantization. The method shows great potential in improving the spike detection accuracy. Moreover, we achieved up to 98.06% in sensitivity and 98.38% in specificity for a dataset consisting of eight healthy and eight epileptic patients. A cluster of the input sequences can be adopted in the future, to build PWM-specific feature for each cluster of similar subjects, to build a hybrid detection model which might improve the accuracy, especially for outliers samples.

Acknowledgment

Research reported in this publication was supported by King Abdullah University of Science and Technology (KAUST) in collaboration with King Abdulaziz City for Science and Technology (KACST) and King Saud University (KSU).

Funding

This research project has been funded by King Abdullah University of Science and Technology (KAUST) Base Research Fund (BAS/1/1627-01-01), in collaboration with King Abdulaziz City for Science and Technology (KACST) and King Saud University (KSU).