Modeling the Ongoing Dynamics of Short and Long-Range Temporal Correlations in Broadband EEG During Movement

1. Introduction

Movement is the primary mode of interaction with the environment and hence studying the neuronal processes involved in movement generation is interesting. The temporal and spectral changes occur in the neuronal processes during voluntary movement. Detecting the movement intention by identifying these changes in the neuronal processes observed in electroencephalography (EEG), not only helps to understand motor command generation but also has applications in Brain Computer Interfaces (BCI). Traditionally, spectral power changes in the narrowband sensorimotor oscillations in EEG such as Event-Related (De)Synchronization (ERD/S) (Pfurtscheller and Lopes da Silva, 1999) are used to determine movement. Along with these narrowband spectral processes, changes also occur in the temporal processes in EEG. EEG shows Long-Range Temporal Correlation (LRTC) (Linkenkaer-Hansen et al., 2001; Hardstone et al., 2012) because of the power-law decay of its autocorrelation. The alpha band amplitude envelope shows LRTC, which decreases during movement (Linkenkaer-Hansen et al., 2004). Despite both the temporal and spectral changes in EEG, narrowband spectral features such as ERD are explored more commonly (Yuan and He, 2014; He et al., 2015), especially for movement detection and its applications in BCI. The LRTC during movement is also primarily obtained on the narrowband alpha amplitude envelope of the longer segments of EEG (Linkenkaer-Hansen et al., 2001; Zhigalov et al., 2016). Such LRTCs require a repeated movement and cannot be used for detecting single movement in real-time from the ongoing EEG. There is a paucity of research studying the ongoing temporal changes in the broadband EEG to detect movement on the single trial basis. In our previous study, we found that the autocorrelation of the broadband EEG decayed slower during movement (Wairagkar et al., 2018). We, therefore, investigate the temporal dependencies in the broadband EEG here in detail and also study their relationship with the ERD.

The temporal dynamics in the brain processes can be assessed by studying the temporal dependencies in EEG. These temporal dependencies can be directly observed from the autocorrelation function (ACF) of EEG. If the autocorrelation becomes zero after finite time lags or decays exponentially, then the time series exhibits a Short-Range Dependence (SRD), otherwise, if the autocorrelation decays slower than the exponential, then it has a Long-Range Dependence (Torre et al., 2007; Dette et al., 2017). A specific case of a long-range dependence is an LRTC which is characterized by the power-law decay of the autocorrelation function (Botcharova, 2014). LRTCs are widely observed in neuronal processes recorded at different levels that show power-law scaling such as neuronal firings (Hu et al., 2013), neuronal avalanches (Benayoun et al., 2010; Palva et al., 2013), local field potentials (Benayoun et al., 2010), electrocorticography, and non-invasive EEG and magnetoencephalography (Nikulin and Brismar, 2005; Benayoun et al., 2010). The spontaneous EEG spectrum is of the form 1/f, which shows power-law scaling and LRTC (Nikulin and Brismar, 2005; Berthouze et al., 2010). LRTCs have been ubiquitously observed in both oscillatory and non-oscillatory neuronal processes in EEG such as various narrowband oscillation amplitude envelope fluctuations (Linkenkaer-Hansen et al., 2001; Hardstone et al., 2012), phase of alpha oscillation (Botcharova, 2014), broadband phase synchrony (Kitzbichler et al., 2009), avalanches (Benayoun et al., 2010; Palva et al., 2013), and the energy profile of EEG (Parish et al., 2004; Benayoun et al., 2010).

The LRTC in EEG is typically obtained in the alpha band amplitude envelope fluctuations (Linkenkaer-Hansen et al., 2004; Zhigalov et al., 2016). During movement, the LRTC in the alpha amplitude fluctuations decreases, possibly due to the disruption caused in the long-memory process by movement (Linkenkaer-Hansen et al., 2004; Botcharova et al., 2014, 2015). However, these alpha amplitude LRTCs completely disregard the LRTCs in the broadband EEG. The brain rhythms are non-stationary and not strictly restricted to the selected narrow sinusoidal frequency bands (Cole and Voytek, 2018), and hence the LRTCs computed on such bands can overlook important features present in the entire power spectrum. The arrhythmic broadband processes and oscillatory processes coexist in neuronal activity (He, 2014). The arrhythmic broadband brain activity was previously considered background noise; however, recent studies demonstrate that it is physiologically and functionally relevant (He, 2014). The dynamics of arrhythmic broadband EEG change with task demands and cognitive states, and it has been associated with the excitation/inhibition balance of neuronal populations (Haller et al., 2018). It is interesting to determine the unexplored changes in such temporal dynamics of the arrhythmic broadband EEG during voluntary movement. The traditional narrowband alpha amplitude fluctuation LRTC requires long EEG segments (Linkenkaer-Hansen et al., 2001), and LRTC is considered an invariant property of brain dynamics over several scales. This approach does not facilitate observation of the ongoing instantaneous changes in LRTC. Berthouze and Farmer (2012) characterized changes in LRTC using a Kalman filter, but their timescales were several seconds long. Here we focus on shorter timescales by computing instantaneous changes in the LRTC using 2 s broadband EEG sliding windows to detect movement intention as opposed to a single LRTC measure for the entire duration of EEG. We call these continuous, instantaneous LRTC measures on short sliding EEG windows as ongoing LRTC. This allows us to observe the real-time changes occurring in the broadband LRTC, which is not possible with the traditional approach.

The LRTC can be quantified by the exponent of the power-law decay (α) of the ACF (ρ(t)) following the relation ρ(t) = Ct^−α, where C is a constant. Equivalently, LRTC can also be quantified by the exponent (β) of the power spectrum of EEG (S(f)) obtained by taking the Fourier transform of ACF which follows the relation S(f) = Bf^−β, where B is a constant. This gives an equality β = 1 − α (Botcharova, 2014). More conveniently, the Hurst exponent (H) can be used to quantify LRTC reliably because using the above methods is challenging in practice (Rangarajan and Ding, 2000; Delignieres et al., 2006). The relation between H and β is $H=\frac{1-\beta }{2}$ (Rangarajan and Ding, 2000). Detrended Fluctuation Analysis (DFA) (Peng et al., 1995) is commonly used to estimate the Hurst exponent because it is effective on a non-stationary time series (Kantelhardt et al., 2001; Delignieres et al., 2006; Hardstone et al., 2012).

In this study, we not only identify the LRTCs in the ongoing dynamics of the broadband EEG during movement but also model them. One of the models for modeling the long-range dependency is the Autoregressive Fractionally Integrated Moving Average (ARFIMA) model. ARFIMA(p, d, q) contains three components: the Autoregressive (AR) process of order p, the Moving Average (MA) process of order q, and the long-range dependence parameter d (d = H−0.5) (Wagenmakers et al., 2004). ARFIMA is a generalization of the Autoregressive Integrated Moving Average (ARIMA) model which is used to model a non-stationary time series. An integrated model such as ARIMA(p, d, q) makes the time series stationary by differencing it d times, after which the remaining stationary parameters AR and MA can be estimated (Box and Jenkins, 2015). In the case of time series such as EEG which has long-range dependencies and 1/f power spectrum, the time series must be fractionally differenced to make it stationary. Hence the ARFIMA(p, d, q) model with a fractional differencing parameter d is more suitable for EEG (Wagenmakers et al., 2004). We chose the ARFIMA model because it allows modeling both short- and long-range dependencies simultaneously, which enabled us to investigate the changes in the ongoing dynamics of both the types of dependencies during voluntary movement intention and execution. ARFIMA is useful for modeling broadband EEG because these types of parametric models can describe the second-order statistics of time series completely (Schlögl et al., 2000). To our knowledge, ongoing changes in the dynamics of LRTC and SRD in the broadband EEG on short windows were not investigated during motor command generation and used for detection of movement on a single trial basis.

In the literature, parametric models for time series analysis such as autoregressive, adaptive autoregressive, multivariate adaptive autoregressive (Hettiarachchi et al., 2015) models were used for movement detection from EEG. However, these attempted to model only selected frequency band amplitudes such as alpha and beta. None of the studies modeled the broadband EEG along with its short- and long-range components.

We also investigated the relationship between these temporal dependencies in the broadband EEG with ERD during movement, which remains unexplored in the literature. We hypothesize that the ARFIMA parameters related to broadband LRTC and SRD can provide additional information about the movement which is complementary to the commonly used narrowband ERD and provide deeper insights into the processes involved in the motor command generation. We hypothesize that these new neural correlates introduced in this paper can also be used in motor-based BCIs where it is important to detect movement in real-time and identify changes in EEG associated with movement on short timescales. Hence, we focus on the 2 s sliding window to estimate ARFIMA features on a single trial basis to detect movement intention with high accuracies and explore the possibility of using these features for BCI.

The focus of this paper is to investigate the fundamental changes in the temporal dynamics of EEG during movement intention from all the available information present in broadband EEG without restricting to a particular narrowband. Second, to test whether these new broadband characteristics can be used to detect movement intention on single trials for potential applications in BCI. Thus, in this paper, we aim (1) to provide a complete characterization of temporal dependencies (SRD/LRTC) in broadband EEG during movement; (2) to accomplish this using the ARFIMA model and study the effect of movement on its parameters; (3) to investigate the complementarity of narrowband ERD and the broadband temporal dependencies LRTC and SRD; (4) to use the three independent streams of information provided by LRTC, SRD, and ERD to classify movement intention for applications in BCI.

2. Materials and Methods

2.1. Participants

EEG was recorded from 14 healthy participants (8 female, age 26 ± 4 years, 12 right-handed) with no prior EEG and BCI experience. This study was carried out following the recommendations of the human experimentation guidelines of the University of Reading. The ethical approval for the EEG experiment was obtained from the ethics committee of the School of Systems Engineering, University of Reading, UK. Informed written consent was obtained from all participants before the EEG recording in accordance with the Declaration of Helsinki.

2.2. Experimental Paradigm

We chose a self-paced, asynchronous single index finger tapping task for EEG recording. Each EEG trial started with a fixation cross at the center of the screen for 2 s followed by an instruction for a single right index finger tap, left index finger tap and resting state in random order. Our aim was to investigate changes in broadband EEG related to right- and left-hand movement intention. The participants were asked to perform the task at any random time of their choice within a window of 10 s following the instruction. To avoid cue effects, we instructed participants not to tap immediately after the display of the instruction. There was a break of 1–1.5 s at the end of each trial. The experimental paradigm was developed in MATLAB Simulink R2014a (The MathWorks, Inc., Natick, Massachusetts, United States) using the BioSig toolbox (Vidaurre et al., 2011).

We recorded 40 EEG trials per condition at the sampling rate of 1,024 Hz. EEG was later downsampled to 128 Hz for further processing. The impedances of all the electrodes were kept below 7 kΩ. We recorded EEG using a Deymed TruScan 32 EEG amplifier (Deymed Diagnostic s.r.o., Hronov, Czech Republic) and EASYCAP EEG cap (EASYCAP GmbH, Herrsching, Germany). In this study, we used channels C3, Cz, and C4 over the motor cortex according to the international 10–20 system out of the 19 EEG channels recorded with a referential montage with reference on FCz and ground on AFz to study movement-related changes in EEG.

We recorded the onset of a finger tap with a bespoke microcontroller tapping device developed using an 8-bit Microchip PICDEM2 Plus demo board (Microchip Technology Inc., Arizona, USA) at 1,000 Hz. The participants placed both the index fingers in the corresponding finger caps of the tapping device. The two channels of binary tapping signals, capturing the onset and duration of each finger tap were co-registered with EEG using TOBI SignalServer 2.0 protocol (Breitwieser et al., 2012). EEG data is available from http://dx.doi.org/10.17864/1947.117 (Wairagkar, 2017). The details of the experimental paradigm and artifact removal are given in (Wairagkar et al., 2018).

2.2.1. Pre-processing and Artifacts Removal

We performed all EEG analyses offline in MATLAB. EEG was filtered using a fourth-order zero-phase non-causal Butterworth filter to avoid phase distortions. The power-line noise was removed with a notch filter at 50 Hz. We removed the DC offset and the high-frequency noise from EEG by band-pass filtering between 0.5 and 45 Hz.

We performed artifact removal using independent component analysis (ICA) (Jung et al., 2000) with EEGLAB toolbox for MATLAB (Delorme and Makeig, 2004), which uses an automated version of infomax ICA (Makeig et al., 1997). We manually identified and removed the independent components with artifacts. We visually inspected the reconstructed uncontaminated EEG again and eliminated any undesirable trials containing large residual artifacts. EEG was segmented into time-locked trials of 6 s (−3 s to +3 s from the onset of the finger tap). These trials were divided into 2 s sliding windows from time t − 2 s to t shifted by 100 ms. The instantaneous indices obtained at time t for a single trial were computed independently on a 2 s broadband (0.5–45 Hz) EEG causal window from t−2 s to t. Thus, we computed new features every 100 ms to avoid causing a lag in BCI. We chose window length of 2 s because it provided a sufficient number of timescales to robustly compute the Hurst exponent using DFA analysis as detailed in the next section without compromising its instantaneous characterization which was essential for studying ongoing temporal dynamics.

2.3. Identifying Long-Range Temporal Correlation (LRTC) in Broadband EEG Using Detrended Fluctuation Analysis (DFA)

Whether a time series has a long- or short-range dependence can be identified from the decay of its autocorrelation function (ACF) and log-log plot of its power spectrum. If the ACF decays according to the power-law, and the log-log power spectrum shows a straight line, then the time series has LRTC. If the ACF decays exponentially, and the log-log power spectrum is not linear, then it has SRD (Wagenmakers et al., 2004). To identify the temporal dependencies in the 2 s broadband EEG, we computed the ACF and log-log power spectrum obtained by squaring the Fourier transform. We plotted the grand average ACF and power spectrum of all the windows in all the trials in all the participants. The LRTCs can be quantified by Hurst exponent using DFA (Peng et al., 1995) in a non-stationary time series such as EEG by avoiding the artifactual dependencies caused by non-stationarity and trends (Peng et al., 1995; Kantelhardt et al., 2001). We performed the DFA on each 2 s sliding window of the broadband EEG which gave an estimate of the Hurst exponent H every 100 ms. The Hanning window was applied to each 2 s EEG segment before the DFA to avoid edge effects. The DFA was performed as follows:

1. The 2 s EEG window X of length N (256 samples) was integrated according to Equation (1).

$\begin{array}{ll}{Y}_k={\displaystyle \sum _{i=1}^k}{X}_i-\overline{X}& (1)\end{array}$

where, k = 1, …, N and Y is the integrated time series. $\overline{X}$ is the mean of X.

2. The integrated time series Y was divided into N/n non-overlapping boxes of length n, where n is an individual timescale at which we computed the Root Mean Square (RMS) fluctuations. We chose the timescales of n = [10, N/4] samples (i.e., [78 ms, 0.5 s]). The box sizes of n = [10, N/4] are commonly used to get good estimate of RMS fluctuations at each timescale (Delignieres et al., 2006; Botcharova et al., 2013). We used 25 box sizes between n = [10, N/4] equidistant on the log₂ scale as our number of samples was a power of 2. Delignieres et al. (2006) have shown that we can obtain correct estimates of Hurst exponent using a short time series of length 256 with DFA. Our box sizes were within the range suggested by Li et al. (2018) $[max(k+2,\frac{F_s}{F_{max}}),min(\frac{N}{4},\frac{F_s}{F_{min}})]$ where k = 1 (linear detrending in DFA) for filtered data between F_min (0.5 Hz) and F_max (45 Hz).

3. At each scale n, for every non-overlapping segment of Y of length n, a trend was obtained by the least square linear fit. Y_n is a concatenation of trends at a scale n for all the N/n boxes and the RMS fluctuations were computed according to Equation (2) for each n = [10, N/4].

$\begin{array}{ll}{F}_n=\sqrt{\frac{1}{N}{\displaystyle \sum _{i=1}^N{\left(Y_i-Y_{n,i}\right)}^2}}& (2)\end{array}$

4. N was not fully divisible by n for each box size. Hence, we obtained the final RMS fluctuations by averaging the F_n computed using the steps 1–3 from the forward and backward direction of each EEG window X (Kantelhardt et al., 2001).

5. We computed the Hurst exponent by obtaining the slope of the linear fit to the log-log plot of RMS fluctuations at each timescale n (log₂F(n) vs. log₂n).

After obtaining H for each EEG window, an exponential smoothing filter was used to smooth H in the consecutive sliding windows in the single trials to avoid noisy estimates.

The Hurst exponent estimated using DFA is valid and suggests the presence of the power-law in the fluctuations at different timescales only if the log-log DFA plot is linear. We validated the Hurst exponents H by assessing the linearity of the log-log DFA plot by comparing the fit of the linear, polynomial, logarithmic and exponential models to it using the maximum likelihood-DFA (ML-DFA) method detailed in Botcharova et al. (2013). The LRTC is present in the time series if the Hurst exponent is between 0.5 and 1 (Kantelhardt et al., 2001; Linkenkaer-Hansen et al., 2001).

2.4. Modeling the Broadband EEG Using Autoregressive Fractionally Integrated Moving Average (ARFIMA) Model

AFRIMA allows simultaneous estimation of both LRTC and SRD in a time series (Wagenmakers et al., 2004). ARFIMA(p,d,q) incorporates SRD processes through the AR parameters p and the MA parameters q, and the LRTC process through the fractional differencing parameter d. The ARFIMA(p,d,q) model is given by Equation (3) (Botcharova, 2014).

$\begin{array}{ll}\left(1-{\displaystyle \sum _{i=1}^p}{\varphi }_i{B}^i\right){(1-B)}^d{X}_t=\left(1+{\displaystyle \sum _{i=1}^q}{\theta }_i{B}^i\right){ϵ}_t& (3)\end{array}$

B is the backshift operator, such that BX_t = X_t−1 and $B^n{X}_t=X_{t-n}$, ϕ_p are the AR coefficients of the order p, θ_q are the MA coefficients of the order q and ϵ_t is innovation at time t drawn from a normal distribution. For ARFIMA, d can have a fractional value. We estimated the ARFIMA(p,d,q) for each 2 s sliding window of the single trial broadband EEG by firstly fractionally differencing the series with d and then estimating the parameters of ARMA(p,q) as described in the following sections.

2.4.1. Removing LRTC From EEG With Fractional Differencing

The parameter d accounts for the LRTC in the ARFIMA process. The AR and MA parameters can only be estimated accurately for a stationary SRD process. The first step of fitting ARFIMA was to fractionally difference each 2 s EEG window by its corresponding fractional differencing parameter d to remove LRTC and make it stationary. The parameter d = H − 0.5, where H was estimated by the DFA method described in the previous section 2.3. The fractional differencing can be performed using a binomial series expansion as given in the Equation (4) with a Gamma function (Granger and Joyeux, 1980; Baillie, 1996; Liu et al., 2017). We used the Matlab fast fractional difference algorithm provided in Jensen and Nielsen (2014) for fractionally differencing each EEG window.

$\begin{array}{ll}{(1-B)}^d={\displaystyle \sum _{k=0}^{\infty }}\left(\begin{array}{c}d\\ k\end{array}\right){(-B)}^k& (4)\end{array}$

2.4.2. Identification of the Order of the ARMA(p,q)

The AR and MA parameters p and q of the ARFIMA were estimated by fitting the ARMA(p,q) model to the fractionally differenced EEG window. The stability of the model was assessed by confirming the stationarity of the time series using the augmented Dicky-Fuller test for unit root (Im et al., 2003). We then identified the model order of ARMA for each 2 s EEG window using Akaike Information Criterion (AIC) (Akaike, 1974) by comparing the models with p ranging from 1 to 10 and q ranging from 1 to p − 1. The order of the best fitting ARMA(p,q) model which gave the least AIC for the maximum number of EEG windows in all the participants was chosen as the order of the AR and MA parameters of ARFIMA.

2.4.3. Estimation of ARFIMA(p,d,q)

Having identified the order of the ARMA part of the ARFIMA, we then estimated the AR and MA parameters p and q, respectively, using the Matlab functions arima() and estimate() from Econometrics toolbox (The MathWorks Inc., 2018) for each 2 s sliding broadband EEG window. The residual analysis was performed on the estimated model using the Ljung-Box Q test (Ljung and Box, 1978) to assess whether the residuals have any significant autocorrelation and one-sample Kolmogorov–Smirnov test (Massey, 1951) to evaluate whether the residuals have a normal distribution. All the parameters d, p, and q obtained for each window were then plotted to assess the ongoing change during the trial.

2.5. Event-Related Desynchronization (ERD) on Single Trial EEG

The ERD analysis was also done on individual 2 s EEG sliding windows. We used the band power method for characterizing the ERD (Pfurtscheller and Lopes da Silva, 1999). Each EEG window was band-pass filtered between 8 and 13 Hz (alpha band). The mean of each window was subtracted from itself. The analytic signal obtained from the absolute value of the Hilbert transform of EEG was used to get the amplitude of the alpha band. The band power was computed by taking the mean of the squared alpha band amplitude. ERD is the relative change in the alpha band power from the baseline during movement. We computed the baseline alpha band power R for individual participants by averaging the band powers of all the 2 s windows in all the resting state trials. We then computed the percent ERD at each time t by subtracting the baseline R from the band power of 2 s sliding EEG window A_t using $ER{D}_t=\left(\frac{A_t-R}{R}\right)\times 100$.

We also evaluated ERD as above on the fractionally differenced EEG windows to observe the effect of removal of LRTC from EEG on the ERD in the alpha band.

2.6. Hybrid Classifier for Movement Intention Detection

From the temporal and spectral EEG analysis described in the previous sections, we computed three types of features: broadband LRTC obtained from DFA, broadband SRD obtained from the parameters of ARFIMA and well-known narrowband ERD which could be used for detecting movement from EEG. We performed the classification using various combinations of the LRTC, SRD, and ERD features to identify right tap vs. resting state and left tap vs. resting state independently using binary Linear Discriminant Analysis (LDA). Our goal was to compare the performance of the classifier for movement intention detection using these features independently and with hybrid combinations of these features.

The classification was done on each participant independently. We trained a separate LDA classifier for each sliding window with the feature vectors from corresponding windows in all the movement trials and the same number of features vectors randomly selected from the resting state trials of that participant. Each LDA had 40 data samples (equivalent to the number of trials) with the selected number of features for each class. We used 10 × 10 fold cross-validation to assess the performance of the classifier by obtaining the classification accuracies, sensitivities and specificities at each time point in the trial given by 2 s sliding windows. The 95% confidence level for binary classification (movement vs. rest) was obtained from the binomial distribution with n = number of EEG trials and p = 0.05. We also noted the time at which classifier crossed this threshold as the time of significantly identifying the movement intention.

3. Results

3.1. Autocorrelation Function and Power Spectrum of Broadband EEG During Voluntary Movement

We plotted the ACF and the log-log power spectrum of a 2 s EEG window from −1 s to +1 s of movement onset such that it contained the information regarding movement intention and execution. We also plotted the ACF and the log-log power spectrum of the corresponding resting state EEG 2 s window. Figure 1A shows that the ACF decays slowly (slower than exponential) in all the conditions indicating the presence of long-range dependence. The ACF decays slower for right and left finger tapping, indicating an increase in long-range dependence. Figure 1B shows the log-log power spectrum, which is linear in all the three conditions with a peak at 10 Hz as expected. The linear log-log power spectrum suggests that the dependence is long-ranged. If the slope of the power spectrum is between −0.5 to −1.5, then it indicates LRTC in the time series (Wagenmakers et al., 2004). The slopes that we observe are in the valid range for LRTC with a slightly increased slope for right and left finger tapping. We explore the increase in LRTC during the voluntary movement further in the next sections because it is not practical to identify the exponent for LRTC robustly from either autocorrelation or power spectrum of a non-stationary process such as EEG, especially on single trials. Hence, we use DFA analysis in the next section to compute Hurst exponent robustly on single EEG windows. The 10 Hz peaks for right and left finger tapping have lower power than that of the resting state peak, especially in channels C3 and C4; this represents the ERD during voluntary movement.

FIGURE 1

Figure 1. Grand average autocorrelation function (ACF) and power spectrum of EEG during movement. (A) The grand average ACF of 2 s broadband EEG from −1 s to +1 s for all the participants for right finger tapping (red), left finger tapping (blue), and resting state (black) in channels C3, Cz, and C4. The ACF decays slowly indicating the LRD. (B) The grand average log-log power spectrum in the channels C3, C4, and Cz for all the three conditions. The log-log power spectrum is linear in all the cases with slopes between −0.5 and −1.5 suggesting the presence of LRTC.

3.2. LRTC in the Broadband EEG Using DFA

We obtained valid Hurst exponents for the LRTCs in 2 s broadband EEG windows in the range of 0.5 to 1 using DFA. The ML-DFA validated the DFA scaling exponents by selecting the linear model as the best fitting model to the log-log DFA fluctuation plots confirming the presence of the LRTCs in the broadband EEG. The LRTC was present in EEG in the movement state as well as in the resting state. Figure 2 shows the grand average Hurst exponent throughout EEG trial of all the participants during movement and resting states. We can see a clear increase in the LRTC during movement intention and execution. The Hurst exponent H increased significantly from the resting state between the vertical gray dotted bars (p < 0.05, Mann–Whitney U-test, n = 14; number of participants) approximately between −1 and 3 s of the onset of movement. The H was used to obtain the fractional differencing parameter d for ARFIMA by subtracting 0.5 from it. We did not observe any lateralization in the LRTC during right- and left-finger taps. Channels C3, Cz, and C4 showed maximum changes in the Hurst exponent. It would be interesting to study the localization of LRTC in the future using spatial filters.

FIGURE 2

Figure 2. The grand average time evolution of Hurst exponents of broadband EEG (H quantifying LRTC). The progressions of grand average H in all the participants in channels C3, Cz, and C4 during right finger tap (red), left finger tap (blue), and resting state (black). The LRTC increases during movement intention and execution. The movement onset is at 0 s, marked by a solid vertical line. The H of movement trials is significantly different from the H of resting state trials in the time region between the dotted gray vertical lines (p < 0.05). The shaded areas show the standard deviation.

3.3. Modeling the Broadband EEG Using ARFIMA(p,d,q)

3.3.1. Removal of LRTC From EEG With Fractional Differencing

We examined the ACF and the log-log power spectrum of the fractionally differenced EEG windows with d to identify whether the LRTC has been removed so that the ARMA(p, q) model could be fitted. Figure 3 shows the grand average ACF, Partial Autocorrelation Function (PACF) and log-log power spectrum of the fractionally differenced 2 s broadband EEG windows from −1 s to +1 s of movement onset in all the three conditions. Figure 3A shows that the ACF decays fast (after four lags) after fractional differencing and Figure 3B shows that the PACF cuts off after nine lags suggesting the presence of the SRD in the residual EEG.

FIGURE 3

Figure 3. Grand average autocorrelation function (ACF), partial autocorrelation function (PACF), and power spectrum of the fractionally differenced EEG. (A) The grand average ACF of 2 s broadband fractionally differenced EEG from −1 s to +1 s for all the participants for right finger tapping (red), left finger tapping (blue), and resting state (black) in channels C3, Cz, and C4. The ACF decays fast indicating SRD. (B) The grand average of PACF in all the three channels in movement and resting state conditions, which cuts off after nine lags. (C) The grand average log-log power spectrum in the channels C3, C4, and Cz for all the three conditions. The log-log power spectrum is no more linear with slopes outside the range of LRTC (−0.5 and −1.5).

The log-log power spectrum is also no longer linear. Figure 3C shows that the lower frequencies have been flattened creating a bend in the power spectrum, which again confirms the removal of LRTC. The slopes of the power spectra are also outside the range of −0.5 and −1.5 for the LRTC. The decrease in alpha power in the right and left movement conditions is still visible in the channels C3 and C4 as smaller 10 Hz peaks than that of the resting state, which is consistent with the ERD.

3.3.2. Identification of the Order of the ARMA(p,q) Model

The fractionally differenced EEG was confirmed to be stationary by the augmented Dicky-Fuller test for unit roots (p < 0.05). EEG windows did not have unit roots indicating the stability of the ARMA model to be fitted to the short-range dependent process.

Figure 3B shows that the PACF cuts off after nine lags, which can suggest the initial estimate of the AR order could be 9. According to the AIC, the order of the best fitting model was ARMA(10,0) for 85% of times of all the windows in all the trials in all the channels and all the participants. In some cases where the order of MA (q) was greater than zero, the roots were non-invertible, and hence we set q = 0 since ARMA(p, q) can also be represented by AR(p) with higher order. The distribution of the selected orders of the models by AIC remained the same for all channels, all participants, all conditions before movement, during movement and after the movement. Hence we selected ARMA(10,0) for modeling the SRD in fractionally differenced EEG.

3.3.3. Estimation of ARFIMA(10,d,0)

The AR parameters of the order 10 were estimated, and the residual analysis was performed on the residue of the model ARMA(10,0). The residual analysis using the Ljung-Box Q test showed that for 95% of all EEG windows, the estimated ARMA(10,0) fitted well (p < 0.05) to the SRD process in the fractionally differenced EEG and there was no more information in the residuals left to be modeled. The Kolmogorov–Smirnov test confirmed that the residuals of all EEG windows had a normal distribution (p < 0.05). Thus, we modeled 2 s broadband EEG windows successfully using ARFIMA(10,d,0) and estimated the eleven model parameters (d and ten parameters for AR).

3.4. Changes in the Long-Range and Short-Range Temporal Dependence Identified From the Ongoing ARFIMA(10,d,0) During Movement

ARFIMA incorporated the LRTC through the parameter d and the SRD through ten AR parameters. Figure 4 shows the grand average time progressions of these ARFIMA parameters throughout the movement trial in all the participants. A clear increase in the parameter d was observed during movement intention and execution. The parameter d was significantly different in movement trials and resting state trials (p < 0.05, Mann–Whitney U-test, n = 203688; individual windows on which parameters were estimated in all the participants). The first six of the ten AR parameters showed change during movement in all the three channels. Though there is a change during movement, there is no significant difference in the individual parameter in the movement trials vs. resting state trial for all the participants together. Thus, on the grand average, during movement, there was a change in the SRD in broadband EEG though not significant, whereas, the LRTC in the broadband EEG increased significantly.

FIGURE 4

Figure 4. Time progression of the grand average ARFIMA(10,d,0) estimated parameters for the broadband EEG. The time progression of 10 grand average AR parameters reflecting the dynamics of short-range dependence and fractional differencing parameter d reflecting the dynamics of long-range dependence (LRTC) throughout the trial in C3, Cz, and C4 for right finger tapping, left finger tapping and resting state. The movement onset is at 0 s marked by a solid vertical line. There were significant changes in the AR parameters for SRD during the movement for 10 individual participants. However, the grand average AR parameters did not show significant change due to inter-participant variability. The parameter d increased significantly (p < 0.05) during the movement.

The same results were obtained on individual participants for d on single trials. However, there was larger variability in the AR parameters in the individual participants. Ten out of fourteen participants showed a significant change in the AR parameters between resting state trials and movement trials (p < 0.05, Mann–Whitney U-test, n = 40; number of trials, note that the p-value was obtained at each time point during the trial). Out of the 10 participants that showed a significant change, 80% showed a change in all the three channels and 90% showed a significant change in at least the first six AR parameters. These parameters had a higher magnitude. None of the participants showed a significant change in AR parameters 8, 9, and 10, which had a lower magnitude. The absolute values of the AR parameters decreased gradually in the higher order parameters. The consecutive order AR parameter values alternated between positive and negative, as seen in Figure 4. Even though ten individual participants show significant differences in their SRD dynamics, averaging over population in obscuring the information in individuals due to larger population variability than individual variability.

3.5. Complementarity of Long-Range Temporal Correlation and ERD

We established the complementarity of narrowband ERD and broadband LRTC by examining the effect of LRTC on estimating ERD by removing LRTC and by examining the effect of ERD on estimating LRTC by removing the ERD from EEG.

First, we compared the ERD on raw EEG having LRTC and EEG after removing the LRTC using fractional differencing. Figure 5A shows the ERD in EEG with LRTC (solid line) and ERD in the same EEG after removing LRTC (dashed line) in right finger tap, left finger tap and resting state. Both the ERDs follow similar traces and show a significant decrease in the alpha band power during movement as compared to the resting state (p < 0.5, Mann–Whitney U-test, n = 14; number of participants) as expected. The difference between ERD when LRTC was present and ERD when the LRTC was absent was not significant in all the three conditions. Same results were obtained for the individual participants' ERD on single trials. Thus, the ERD is unaffected by the presence of LRTC in EEG.

FIGURE 5

Figure 5. Effect of LRTC on ERD and effect of ERD on LRTC. (A) Solid lines show the time progression of grand average ERD calculated from EEG containing LRTCs in the right finger tap (red), left finger tap (blue), and resting state (black) trials from all the participants. The dotted lines represent the ERD after removing LRTCs from the broadband EEG by fractionally differencing it by d. The movement onset is at 0 s shown by a solid vertical line. Both the ERDs during movement are significantly different from the resting state (p < 0.05). ERDs with and without LRTCs are not significantly different. (B) Solid lines show the time progression of grand average LRTC calculated from EEG containing ERD in right tap, left tap and resting state. The dotted lines represent the LRTC after removing ERD from the broadband EEG by fixing the 8–13 Hz power. Both the LRTCs during movement are significantly different from the resting state (p < 0.05). LRTCs with and without ERD are not significantly different from each other. LRTC and ERD do not have a significant effect on each other because removing one does not change the values of the other.

We then tested the effect of ERD on LRTC by comparing the LRTC on raw EEG and LRTC on EEG with ERD removed. We suppressed the changes in ERD by fixing the magnitude of alpha power in all the windows of a trial. This was done by replacing the magnitude of 8–13 Hz in the Fourier spectrum of all the 2 s windows in a trial with the magnitude of 8–13 Hz from a randomly selected window from that trial and then reconstructing EEG time series by inverse Fourier transform. This ensured that the alpha power did not change throughout the trial resulting in a constant value of ERD at all the time points in a trial. Figure 5B shows the LRTC in EEG with ERD (solid line) and LRTC in EEG after suppressing the ERD (dashed line). There is no significant change in the two LRTCs. Thus, the LRTC is unaffected by ERD. This indicates that the broadband LRTC and narrowband ERD are distinct and complementary processes. Knowledge of ERD does not inform us about the values of LRTC and vice versa.

The ERD decreased, whereas the LRTC increased during the movement. The scatter plot of the grand average ERD and LRTC in Figure 6A shows a strong inverse correlation between them with the high correlation coefficients for right and left finger tap conditions. There was no correlation between ERD and LRTC during the resting state. During movement, there is a switch in the dynamics of ERD and LRTC, and they become coupled from their uncorrelated state during the resting condition.

FIGURE 6

Figure 6. The correlation between ERD and LRTC. (A) Scatter plot shows a strong inverse correlation between ERD and LRTC during right (red) and left (blue) finger tapping in channels C3, Cz, and C4. There is no correlation between ERD and LRTC during resting state (black). The correlation coefficients are shown beside the fitted straight lines. (B) The distribution of lags with maximum cross-correlation between LRTC and ERD in single trials of all the participants in all the three conditions. The LRTC and ERD have maximum correlation at lag 0.

Figure 6B shows the distribution of lags at which there is a maximum correlation between LRTC and ERD in single trials in all the participants in all the three conditions. The correlation is maximum at lag 0, which indicates that LRTC or ERD processes do not precede one another and occur at the same time.

3.6. Hybrid Classifier for Movement Intention Detection

We compared the classification accuracies for movement intention detection using hybrid features of LRTC, ARFIMA and ERD. For the LRTC features, we used H from all the three channels C3, Cz and C4 as LRTCs were observed in all the three channels with equal strength (Figure 2). For the ERD features, we used channel C3 and C4 as Cz showed relatively less distinction between movement and rest ERD (Figure 5). For the ARFIMA features, we used parameters d and AR1 to AR6 from all the three channels since these AR parameters significantly changed in most participants during movement (p < 0.05) leading to 21 features. The average classification accuracies, sensitivities, and specificities of all the participants using all the classifiers are shown in Table 1. The classification accuracies of individual participants are included in Supplementary Tables 1A,B.

TABLE 1

Table 1. The average of peak LDA classification accuracies for the right and left finger movement vs. resting state of all the participants using LRTC, ARFIMA, ERD, hybrid LRTC+ERD and hybrid ARFIMA+ERD features.

The LRTC features showed higher classification accuracies than the ERD features. There was no significant increase in the classification accuracies of the hybrid LRTC and ERD classifier and was similar to the LRTC classifier. Combining the LRTC and ERD features improved the classification accuracies marginally. The classifier with ARFIMA features (which contained features representing both short- and long-range dependence) gave significantly higher accuracies of 87.30±4.87% than LRTC and ERD and their combination (p < 10⁻⁵, Wilcoxon signed-rank test, n = 28; right and left hand accuracies of 14 participants). The hybrid classifier with ARFIMA and ERD gave the highest mean classification accuracy of 88.71±4.12% which was only marginally higher (by 1.4%) than ARFIMA. The ERD, LRTC and ARFIMA accuracies were also correlated, i.e., the participant with lower classification accuracy with LRTC and ARFIMA also showed lower classification accuracies with ERD.

The LRTC classification accuracies, sensitivities and specificities were all similar, which shows that this classifier is robust and reliable. The hybrid classifier of LRTC and ERD and ARFIMA and ERD also showed the same robustness. Interestingly, the ERD classifier had high sensitivity but low specificity leading to lower accuracy. This indicates that by using ERD, we are likely to identify movement with higher accuracy but also get more false positives during the resting state. Considering these results, including the LRTC and SRD temporal dependency features improves the movement intention detection instead of using only ERD. The LRTC and SRD provide complementary information to ERD about motor command generation.

3.6.1. Timing of Movement Intention Detection

Figure 7 shows the time at which the classification accuracies crossed the significance threshold (chance level) in three types of classifiers. The LRTC, ERD and ARFIMA all were able to detect movement intention before its onset in most of the participants. On average, ERD detected movement −0.25 s of the actual movement onset, LRTC detected it at −0.5 s, and ARFIMA detected movement earliest at −1 s. Supplementary Figure 1 shows timings of all hybrid classifiers. All temporal dependency features detected movement earlier than ERD, which shows their suitability for application in BCI.

FIGURE 7

Figure 7. The timings of movement intention detection. The timings of movement intention detection when the classification accuracy crossed the significance threshold for right finger tap and left finger tap are shown for all the 14 participants. Timings obtained from classifiers using LRTC, ARFIMA, and ERD features are shown.

4. Discussion

We completely characterized the broadband temporal dependencies in EEG during movement by estimating the short and long-range temporal correlations. We have not only identified the ongoing changes in the dynamics of these temporal dependencies in the broadband EEG but also modeled them successfully using ARFIMA(10,d,0). Several other studies model EEG with ARMA models (Pfurtscheller et al., 2006; Hu, 2009; Chae et al., 2012; Resalat and Saba, 2016; Antelis et al., 2017), however, these models are suitable for processes with short-range dependencies. We have shown that the broadband EEG contains LRTCs and hence the ARFIMA model is more suitable, which incorporates both LRTC and SRD simultaneously. Monitoring the ongoing changes in both the LRTC and SRD continuously has resulted in the detection of movement intention with high accuracy, which has not been done before.

The LRTC increased significantly during the movement (p < 0.5) consistently in all the participants indicating that it is a robust correlate of movement. The SRD parameters also changed significantly (p < 0.5) in 71.4% participants, but there was a higher variability among participants. ARFIMA is a detailed mathematical model for time series, but it does not explain the underlying physiological processes and interpretation of EEG (Wagenmakers et al., 2004). The order of AR parameters of ARFIMA that we identified using AIC is in a similar range as the orders used by studies modeling EEG with AR which were between 9 and 13 by Florian and Pfurtscheller (1995), 10 by Fabiani et al. (2004), 8.67 by Tseng et al. (1995) and 16 by McFarland and Wolpaw (2008).

Most of the studies use AR and MA models for merely extracting the features from EEG, mostly from sensorimotor rhythms for BCI classification; however, these studies do not describe the nature of the ongoing changes in the model parameters during motor activity (Pfurtscheller et al., 2006; Hu, 2009; Chae et al., 2012; Resalat and Saba, 2016). The parametric ARMA models in the literature were used to characterize spectral features of specific narrowbands, especially of the sensorimotor rhythms for movement (Pfurtscheller et al., 2000; Fabiani et al., 2004; Chae et al., 2012; McFarland et al., 2015; Antelis et al., 2017). The ongoing changes in the broadband LRTCs were not studied previously and used for identification of movement intention. The broadband LRTC might even have implications in self-organized criticality in the brain (Linkenkaer-Hansen et al., 2001). LRTCs are observed in several brain processes (Botcharova et al., 2014). In EEG, LRTCs are typically observed in the envelope of alpha oscillations (Linkenkaer-Hansen et al., 2001). LRTCs are typically attributed to the theory of self-organized criticality in the brain facilitating the neural networks to quickly reorganized during varying processing demands (Linkenkaer-Hansen et al., 2001). However, further work will be required to study the mechanisms of broadband LRTC, and their relationship with self-organized criticality remains to be seen. These broadband temporal dynamics are more robust than ERD because they do not require selection of participant-specific frequency bands for better performance (Pfurtscheller and Lopes da Silva, 1999; Durka et al., 2001).

We discovered that the arrhythmic broadband LRTC and the rhythmic narrowband ERD are complementary processes and provide different information about movement intention. Removing the broadband LRTC from EEG did not affect the ERD and removing ERD from EEG did not affect LRTC significantly (Figure 5). The LRTC represented by broadband 1/f process and ERD represented by the alpha peak that resides over this 1/f spectrum are complementary movement-related neuronal process. Becker et al. (2018) suggested that the alpha power causes the change in LRTC in lower frequency band by observing a negative lag of maximum correlation between alpha power and LRTC in the spontaneous EEG. However, negative lag is not sufficient to conclude the causal effect between the alpha power and LRTC. Our results show that the broadband LRTC and ERD have inverse correlation during movement and the maximum correlation magnitude was at lag zero (see Figure 6), and thus both occur at the same time. The co-evolution of ERD and broadband LRTC may be because of a common input that drives both the processes during motor command generation.

Another independent neuronal process known to occur during movement intention is Motor-Related Cortical Potential (MRCP) or Bereitschafts potential (Shibasaki and Hallett, 2006; Bai et al., 2011). Burke et al. (2005) used an AR model with Bereitschafts potential as an exogenous input for movement detection and observed increased classification accuracy by combining the features from complementary processes. We have previously shown that the changes in the temporal dependencies are also complementary to the MRCPs (Wairagkar et al., 2015) because MRCP which is commonly observed in low frequencies (< 2 Hz) characterizes slow trend in EEG by eliminating the fluctuations, and on the contrary, the broadband LRTC characterizes the dynamics in the fluctuations by eliminating the slow trends using DFA. Thus, we have independent processes of LRTC, ERD, and MRCP containing complementary information about motor command generation.

We have confirmed that the LRTCs can be detected on short windows in broadband with rigorous analysis and ML-DFA (Botcharova et al., 2013) for validating the Hurst exponent. ARFIMA showed that the broadband EEG contained both short- and long-range temporal correlations. We obtained new model parameters every 100 ms describing the ongoing changes in the LRTC and SRD using an online type of processing pipeline that allowed us to detect the movement intention on single trials. The sliding window approach for model parameter estimation has been commonly used for non-stationary time series (Antelis et al., 2017). However, the model estimation is computationally expensive and will need efficient execution for BCI applications. The LRTC estimation, on the other hand, is achievable within the available time window for online BCI applications.

The hybrid classifier with ARFIMA parameters and ERD has a feature vector with high dimensionality and gives high classification accuracies. It was clear that increasing the feature vector dimensionality improved the classification accuracy significantly (p < 10⁻⁵) than just ERD, LRTC or their combination, which had a maximum of five features. We used 10 × 10 fold cross-validation scheme, and hence, the improved accuracy was not due to overfitting. The ARFIMA parameters reinforce the difference between the resting state and movement as multiple parameters show a clear distinction between the two conditions. Higher dimensionality made the feature space sparse, making it easier for LDA to find an optimum classification boundary. It is common in BCI to have feature vectors of high dimensionality, as generally several features are extracted from several channels and used for classification (Lotte et al., 2007). For example, the study in Hettiarachchi et al. (2015) used 52 features with 70 EEG trials per class. The training sets in BCI are often relatively small because of the time-consuming EEG recording processes, which is also exhausting for the participants (Lotte et al., 2007). Our hybrid classifier thus follows standard practices in BCI, giving high classification accuracies. However, it would be interesting to observe whether the classification accuracies using ARFIMA features remain the same after increasing the number of trials significantly. It would be interesting to investigate in future whether these LRTCs provide more information about the movement such as its kinematics. Based on our previous study (Wairagkar et al., 2018) of the autocorrelation of EEG during movement, Robinson et al. (2017) observed that changes in autocorrelation decay of EEG can also determine the speed of movement.

The ERD has higher sensitivity, and it is likely to give more false positives during the resting state leading to lower overall accuracy (Table 1). The LRTC based classifiers (including LRTC+ERD, ARFIMA and ARFIMA+ERD features) do not have bias toward movement unlike ERD and have similar sensitivity and specificity. Hence, these classifiers are robust and have higher classification accuracy. The movement can be predicted earliest using all the hybrid features containing ARFIMA and ERD up to 1 s before its onset. Thus, using complementary neural correlates of broadband LRTC and narrowband ERD could be beneficial for use in BCI applications.

5. Conclusion

We demonstrated broadband LRTCs as novel neural correlates of movement intention in EEG. We propose that the broadband LRTC is the third fundamental neuronal process related to voluntary movement, complementary to the well-known MRCP and narrowband processes such as ERD, providing complementary information about motor command generation.

We identified that the short- and long-range temporal correlations of the broadband EEG estimated from 2 s windows change significantly (p < 0.05) during movement intention and execution. Consistently with this, we used ARFIMA(10,d,0) to identify ongoing short- and long-range temporal correlations.

Our approach enabled us to predict the movement 1 s before its onset, significantly earlier than the prediction based on ERD. We also obtained significantly higher classification accuracies (p < 0.05) than those obtained from ERD. The best classification accuracy was obtained for hybrid ARFIMA (accounting for both short- and long-range temporal correlations) and ERD features (88.3 ± 4.2%). The resulting method offers a robust movement intention detection that might be useful for online BCI.

Data Availability Statement

The datasets generated for this study can be found in the University of Reading Research Data Archive (http://dx.doi.org/10.17864/1947.117).

Ethics Statement

The studies involving human participants were reviewed and approved by the ethics committee of the School of Systems Engineering, University of Reading, UK. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

MW conceptualized the study, designed and conducted the experiments, developed the methodology, analyzed and interpreted the results, created the artwork, and wrote the original manuscript. YH and SN supervised the study, conceptualized the methodology, validated and interpreted the results, and reviewed the manuscript.

Funding

This research was funded by the University of Reading International Research Studentship awarded to MW (reference number: GS14-46, URL: http://www.reading.ac.uk/). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnsys.2019.00066/full#supplementary-material

References

Akaike, H. (1974). A new look at the statistical model identification. IEEE Trans. Automat. Control 19, 716–723. doi: 10.1109/TAC.1974.1100705

CrossRef Full Text | Google Scholar

Antelis, J. M., Montesano, L., Ramos-Murguialday, A., Birbaumer, N. P., and Minguez, J. (2017). Decoding upper limb movement attempt from EEG measurements of the contralesional motor cortex in chronic stroke patients. IEEE Trans. Biomed. Eng. 64, 99–111. doi: 10.1109/TBME.2016.2541084

PubMed Abstract | CrossRef Full Text | Google Scholar

Bai, O., Rathi, V., Lin, P., Huang, D., Battapady, H., Fei, D. Y., et al. (2011). Prediction of human voluntary movement before it occurs. Clin. Neurophysiol. 122, 364–372. doi: 10.1016/j.clinph.2010.07.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Baillie, R. T. (1996). Long memory processes and fractional integration in econometrics. J. Econometr. 73, 5–59. doi: 10.1016/0304-4076(95)01732-1

CrossRef Full Text | Google Scholar

Becker, R., Van De Ville, D., and Kleinschmidt, A. (2018). Alpha oscillations reduce temporal long-range dependence in spontaneous human brain activity. J. Neurosci. 38, 755–764. doi: 10.1523/JNEUROSCI.0831-17.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Benayoun, M., Kohrman, M., Cowan, J., and van Drongelen, W. (2010). EEG, temporal correlations, and avalanches. J. Clin. Neurophysiol. 27, 458–464. doi: 10.1097/WNP.0b013e3181fdf8e5

PubMed Abstract | CrossRef Full Text | Google Scholar

Berthouze, L., and Farmer, S. F. (2012). Adaptive time-varying detrended fluctuation analysis. J. Neurosci. Methods 209, 178–188. doi: 10.1016/j.jneumeth.2012.05.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Berthouze, L., James, L. M., and Farmer, S. F. (2010). Human EEG shows long-range temporal correlations of oscillation amplitude in Theta, Alpha and Beta bands across a wide age range. Clin. Neurophysiol. 121, 1187–1197. doi: 10.1016/j.clinph.2010.02.163

PubMed Abstract | CrossRef Full Text | Google Scholar

Botcharova, M. (2014). Modelling and analysis of amplitude, phase and synchrony in human brain activity patterns (Ph.D. thesis). University College London, London, England.

Google Scholar

Botcharova, M., Berthouze, L., Brookes, M. J., Barnes, G. R., and Farmer, S. F. (2015). Resting state MEG oscillations show long-range temporal correlations of phase synchrony that break down during finger movement. Front. Physiol. 6:183. doi: 10.3389/fphys.2015.00183

PubMed Abstract | CrossRef Full Text | Google Scholar

Botcharova, M., Farmer, S. F., and Berthouze, L. (2013). A maximum likelihood based technique for validating detrended fluctuation analysis (ML-DFA). arXiv (preprint) arXiv:1306.5075. Available online at: https://arxiv.org/abs/1306.5075

Google Scholar

Botcharova, M., Farmer, S. F., and Berthouze, L. (2014). Markers of criticality in phase synchronisation. Front. Syst. Neurosci. 8:176. doi: 10.3389/fnsys.2014.00176

PubMed Abstract | CrossRef Full Text | Google Scholar

Box, G. E. P., and Jenkins, G. M. (2015). Time Series Analysis: Forecasting and Control, 2nd Edn. Hoboken, NJ: John Wiley & Sons.

Google Scholar

Breitwieser, C., Daly, I., Neuper, C., and Müller-Putz, G. R. (2012). Proposing a standardized protocol for raw biosignal transmission. IEEE Trans. Biomed. Eng. 59, 852–859. doi: 10.1109/TBME.2011.2174637

PubMed Abstract | CrossRef Full Text | Google Scholar

Burke, D. P., Kelly, S. P., De Chazal, P., Reilly, R. B., and Finucane, C. (2005). A parametric feature extraction and classification strategy for brain-computer interfacing. IEEE Trans. Neural Syst. Rehab. Eng. 13, 12–17. doi: 10.1109/TNSRE.2004.841881

PubMed Abstract | CrossRef Full Text | Google Scholar

Chae, Y., Jeong, J., and Jo, S. (2012). Toward brain-actuated humanoid robots: asynchronous direct control using an EEG-Based BCI. IEEE Trans. Robot. 28, 1131–1144. doi: 10.1109/TRO.2012.2201310

CrossRef Full Text | Google Scholar

Cole, S. R., and Voytek, B. (2018). Cycle-by-cycle analysis of neural oscillations. bioRxiv 302000. doi: 10.1101/302000

PubMed Abstract | CrossRef Full Text | Google Scholar

Delignieres, D., Ramdani, S., Lemoine, L., Torre, K., Fortes, M., and Ninot, G. (2006). Fractal analyses for ‘short' time series: a re-assessment of classical methods. J. Math. Psychol. 50, 525–544. doi: 10.1016/j.jmp.2006.07.004

CrossRef Full Text | Google Scholar

Delorme, A., and Makeig, S. (2004). EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis. J. Neurosci. Methods 134, 9–21. doi: 10.1016/j.jneumeth.2003.10.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Dette, H., Preuss, P., and Sen, K. (2017). Detecting long-range dependence in non-stationary time series. Electron. J. Stat. 11, 1600–1659. doi: 10.1214/17-EJS1262

CrossRef Full Text | Google Scholar

Durka, P. J., Ircha, D., Neuper, C., and Pfurtscheller, G. (2001). Time-frequency microstructure of event-related electro-encephalogram desynchronization and synchronisation. Med. Biol. Eng. Comput. 39, 315–321. doi: 10.1007/BF02345286

CrossRef Full Text | Google Scholar

Fabiani, G. E., McFarland, D. J., Wolpaw, J. R., and Pfurtscheller, G. (2004). Conversion of EEG activity into cursor movement by a brain-computer interface (BCI). IEEE Trans. Neural Syst. Rehab. Eng. 12, 331–338. doi: 10.1109/TNSRE.2004.834627

PubMed Abstract | CrossRef Full Text | Google Scholar

Florian, G., and Pfurtscheller, G. (1995). Dynamic spectral analysis of event-related EEG data. Electroencephalogr. Clin. Neurophysiol. 95, 393–396. doi: 10.1016/0013-4694(95)00198-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Granger, C. W., and Joyeux, R. (1980). An introduction to long memory time series models and fractional differencing. J. Time Ser. Anal. 1, 15–29. doi: 10.1111/j.1467-9892.1980.tb00297.x

CrossRef Full Text | Google Scholar

Haller, M., Donoghue, T., Peterson, E., Varma, P., Sebastian, P., Gao, R., et al. (2018). Parameterizing neural power spectra. bioRxiv 299859. doi: 10.1101/299859

CrossRef Full Text | Google Scholar

Hardstone, R., Poil, S.-S. S., Schiavone, G., Jansen, R., Nikulin, V. V., Mansvelder, H. D., et al. (2012). Detrended fluctuation analysis: a scale-free view on neuronal oscillations. Front. Physiol. 3:450. doi: 10.3389/fphys.2012.00450

PubMed Abstract | CrossRef Full Text | Google Scholar

He, B., Baxter, B., Edelman, B. J., Cline, C. C., and Ye, W. W. (2015). Noninvasive brain-computer interfaces based on sensorimotor rhythms. Proc. IEEE 103, 907–925. doi: 10.1109/JPROC.2015.2407272

CrossRef Full Text | Google Scholar

He, B. J. (2014). Scale-free brain activity: past, present, and future. Trends Cogn. Sci. 18, 480–487. doi: 10.1016/j.tics.2014.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Hettiarachchi, I. T., Nguyen, T. T., and Nahavandi, S. (2015). “Multivariate adaptive autoregressive modeling and Kalman filtering for motor imagery BCI,” in Proceedings - 2015 IEEE International Conference on Systems, Man, and Cybernetics, SMC 2015 (Kowloon), 3164–3168. doi: 10.1109/SMC.2015.549

CrossRef Full Text | Google Scholar

Hu, J., Zheng, Y., and Gao, J. (2013). Long-range temporal correlations, multifractality, and the causal relation between neural inputs and movements. Front. Neurol. 4:158. doi: 10.3389/fneur.2013.00158

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, J. F. (2009). “New biometric approach based on motor imagery EEG signals,” in FBIE 2009 - 2009 International Conference on Future BioMedical Information Engineering (Sanya), 94–97. doi: 10.1109/FBIE.2009.5405787

CrossRef Full Text | Google Scholar

Im, K. S., Pesaran, M. H., and Shin, Y. (2003). Testing for unit roots in heterogeneous panels. J. Econometr. 115, 53–74. doi: 10.1016/S0304-4076(03)00092-7

CrossRef Full Text | Google Scholar

Jensen, A. N., and Nielsen, M. Ø. (2014). A fast fractional difference algorithm. J. Time Ser. Anal. 35, 428–436. doi: 10.1111/jtsa.12074

CrossRef Full Text | Google Scholar

Jung, T. P., Makeig, S., Humphries, C., Lee, T. W., McKeown, M. J., Iragui, V., et al. (2000). Removing electroencephalographic artifacts by blind source separation. Psychophysiology 37, 163–178. doi: 10.1111/1469-8986.3720163

PubMed Abstract | CrossRef Full Text | Google Scholar

Kantelhardt, J. W., Koscielny-Bunde, E., Rego, H. H., Havlin, S., and Bunde, A. (2001). Detecting long-range correlations with detrended fluctuation analysis. Phys. A 295, 441–454. doi: 10.1016/S0378-4371(01)00144-3

CrossRef Full Text | Google Scholar

Kitzbichler, M. G., Smith, M. L., Christensen, S. R., and Bullmore, E. (2009). Broadband criticality of human brain network synchronization. PLoS Comput. Biol. 5:e1000314. doi: 10.1371/journal.pcbi.1000314

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, R., Wang, J., and Chen, Y. (2018). Effect of the signal filtering on detrended fluctuation analysis. Phys. A 494, 446–453. doi: 10.1016/j.physa.2017.12.011

CrossRef Full Text | Google Scholar

Linkenkaer-Hansen, K., Nikulin, V. V., Palva, J. M., and Ilmoniemi, R. J. (2001). Long-range temporal correlations and scaling behavior in human brain oscillations. J. Neurosci. 21, 1370–7. doi: 10.1002/anie.201106423

PubMed Abstract | CrossRef Full Text | Google Scholar

Linkenkaer-Hansen, K., Nikulin, V. V., Palva, J. M., Kaila, K., and Ilmoniemi, R. J. (2004). Stimulus-induced change in long-range temporal correlations and scaling behaviour of sensorimotor oscillations. Eur. J. Neurosci. 19, 203–211. doi: 10.1111/j.1460-9568.2004.03116.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, K., Chen, Y., and Zhang, X. (2017). An evaluation of ARFIMA (autoregressive fractional integral moving average) programs. Axioms 6:16. doi: 10.3390/axioms6020016

CrossRef Full Text | Google Scholar

Ljung, G. M., and Box, G. E. P. (1978). On a measure of lack of fit in time series models. Biometrika 65, 297–303. doi: 10.1093/biomet/65.2.297

CrossRef Full Text | Google Scholar

Lotte, F., Congedo, M., Lécuyer, A., Lamarche, F., and Arnaldi, B. (2007). A review of classification algorithms for EEG-based brain-computer interfaces. J. Neural Eng. 4, R1–R13. doi: 10.1088/1741-2560/4/2/R01

PubMed Abstract | CrossRef Full Text | Google Scholar

Makeig, S., Jung, T. P., Bell, A. J., Ghahremani, D., and Sejnowski, T. J. (1997). Blind separation of auditory event-related brain responses into independent components. Proc. Natl. Acad. Sci. U.S.A. 94, 10979–84. doi: 10.1073/pnas.94.20.10979

PubMed Abstract | CrossRef Full Text | Google Scholar

Massey, F. J. (1951). The Kolmogorov-Smirnov test for goodness of fit. J. Am. Stat. Assoc. 46, 68–78. doi: 10.1080/01621459.1951.10500769

CrossRef Full Text | Google Scholar

McFarland, D. J., Sarnacki, W. A., and Wolpaw, J. R. (2015). Effects of training pre-movement sensorimotor rhythms on behavioral performance. J. Neural Eng. 12:066021. doi: 10.1088/1741-2560/12/6/066021

PubMed Abstract | CrossRef Full Text | Google Scholar

McFarland, D. J., and Wolpaw, J. R. (2008). Sensorimotor rhythm-based brain-computer interface (BCI): Model order selection for autoregressive spectral analysis. J. Neural Eng. 5, 155–162. doi: 10.1088/1741-2560/5/2/006

PubMed Abstract | CrossRef Full Text | Google Scholar

Nikulin, V. V., and Brismar, T. (2005). Long-range temporal correlations in electroencephalographic oscillations: relation to topography, frequency band, age and gender. Neuroscience 130, 549–558. doi: 10.1016/j.neuroscience.2004.10.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Palva, J. M., Zhigalov, A., Hirvonen, J., Korhonen, O., Linkenkaer-Hansen, K., and Palva, S. (2013). Neuronal long-range temporal correlations and avalanche dynamics are correlated with behavioral scaling laws. Proc. Natl. Acad. Sci. U.S.A. 110, 3585–3590. doi: 10.1073/pnas.1216855110

PubMed Abstract | CrossRef Full Text | Google Scholar

Parish, L. M., Worrell, G. A., Cranstoun, S. D., Stead, S. M., Pennell, P., and Litt, B. (2004). Long-range temporal correlations in epileptogenic and non-epileptogenic human hippocampus. Neuroscience 125, 1069–1076. doi: 10.1016/j.neuroscience.2004.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Peng, C. K., Havlin, S., Stanley, H. E., and Goldberger, A. L. (1995). Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series. Chaos 5, 82–87. doi: 10.1063/1.166141

PubMed Abstract | CrossRef Full Text | Google Scholar

Pfurtscheller, G., and Lopes da Silva, F. H. (1999). Event-related EEG/MEG synchronization and desynchronization: Basic principles. Clin. Neurophysiol. 110, 1842–1857. doi: 10.1016/S1388-2457(99)00141-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Pfurtscheller, G., Müller-Putz, G. R., Schlögl, A., Graimann, B., Scherer, R., Leeb, R., et al. (2006). 15 Years of BCI research at Graz University of Technology: Current projects. IEEE Trans. Neural Syst. Rehab. Eng. 14, 205–210. doi: 10.1109/TNSRE.2006.875528

PubMed Abstract | CrossRef Full Text | Google Scholar

Pfurtscheller, G., Neuper, C., Guger, C., Harkam, W., Ramoser, H., Schlögl, A., et al. (2000). Current trends in Graz Brain-Computer Interface (BCI) research. IEEE Trans. Rehab. Eng. 8, 216–219. doi: 10.1109/86.847821

PubMed Abstract | CrossRef Full Text | Google Scholar

Rangarajan, G., and Ding, M. (2000). Integrated approach to the assessment of long range correlation in time series data. Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics 61, 4991–5001. doi: 10.1103/PhysRevE.61.4991

PubMed Abstract | CrossRef Full Text | Google Scholar

Resalat, S. N., and Saba, V. (2016). A study of various feature extraction methods on a motor imagery based brain computer interface system. Basic Clin. Neurosci. J. 7, 13–20.

PubMed Abstract | Google Scholar

Robinson, N., and Vinod, A. (2017). “Decoding speed of hand movement execution using temporal features of EEG,” in 8th International IEEE/EMBS Conference on Neural Engineering (NER) (Shanghai), 572–575. doi: 10.1109/NER.2017.8008416

CrossRef Full Text | Google Scholar

Schlögl, A., Universit, T., Wissenschaften, T., and Technik, B. (2000). The electroencephalogram and the adaptive autoregressive model: theory and applications (dissertation) (Graz: Technical University of Graz), 68.

Google Scholar

Shibasaki, H., and Hallett, M. (2006). What is the Bereitschaftspotential? Clin. Neurophysiol. 117, 2341–2356. doi: 10.1016/j.clinph.2006.04.025

PubMed Abstract | CrossRef Full Text | Google Scholar

The MathWorks Inc. (2018). Econometrics Toolbox Documentation - MathWorks United Kingdom (R2018a).

Google Scholar

Torre, K., Delignieres, D., and Lemoine, L. (2007). Detection of long-range dependence and estimation of fractal exponents through ARFIMA modelling. Brit. J. Math. Stat. Psychol. 60, 85–106. doi: 10.1348/000711005X89513

PubMed Abstract | CrossRef Full Text | Google Scholar

Tseng, S.-Y. S.-Y., Chen, R.-C., Chong, F.-C., and Kuo, T.-S. (1995). Evaluation of parametric methods in EEG signal analysis. Med. Eng. Phyxs 17, 71–78. doi: 10.1016/1350-4533(95)90380-T

PubMed Abstract | CrossRef Full Text | Google Scholar

Vidaurre, C., Sander, T. H., and Schlögl, A. (2011). BioSig: the Free and Open Source Software Library for Biomedical Signal Processing. Comput. Intell. Neurosci. 2011, 1–12. doi: 10.1155/2011/935364

PubMed Abstract | CrossRef Full Text | Google Scholar

Wagenmakers, E.-J., Farrell, S., and Ratcliff, R. (2004). Estimation and interpretation of 1/fα noise in human cognition. Psychon. Bull. Rev. 11, 579–615. doi: 10.3758/BF03196615

PubMed Abstract | CrossRef Full Text | Google Scholar

Wairagkar, M. (2017). EEG Data for Voluntary Finger Tapping Movement. University of Reading Dataset. doi: 10.17864/1947.117

CrossRef Full Text | Google Scholar

Wairagkar, M., Hayashi, Y., and Nasuto, S. J. (2015). “Movement intention detection from autocorrelation of EEG for BCI,” in Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), Vol. 9250. eds Y. Guo, K. Friston, F. Aldo, S. Hill, and H. Peng (London: Springer), 212–221. doi: 10.1007/978-3-319-23344-4_21

CrossRef Full Text | Google Scholar

Wairagkar, M., Hayashi, Y., and Nasuto, S. J. (2018). Exploration of neural correlates of movement intention based on characterisation of temporal dependencies in electroencephalography. PLoS ONE 13:e0193722. doi: 10.1371/journal.pone.0193722

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuan, H., and He, B. (2014). Brain-computer interfaces using sensorimotor rhythms: current state and future perspectives. IEEE Trans. Biomed. Eng. 61, 1425–1435. doi: 10.1109/TBME.2014.2312397

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhigalov, A., Kaplan, A. Y., and Palva, J. M. (2016). Modulation of critical brain dynamics using closed-loop neurofeedback stimulation. Clin. Neurophysiol. 127, 2882–2889. doi: 10.1016/j.clinph.2016.04.028

PubMed Abstract | CrossRef Full Text | Google Scholar