Stochastic non-linear oscillator models of EEG: the Alzheimer's disease case

1. Introduction

Quantitative analysis of human brain electroencephalography (EEG) recordings aimed at enhancing our understanding of brain injuries and disorders is currently an important research area. In addition to being useful in diagnosis, such analysis can provide insights into the underlying neurophysiology of the injury or disorder, thereby leading to better treatment and preventive strategies. Alzheimer's disease (AD) is the most common form of dementia and is the subject of intense research. While no known cure exists, certain medications have shown promise in delaying the symptoms (Dauwels et al., 2010) prompting researchers to seek early diagnosis and intervention strategies. In this context, analysis of the EEG is a potential non-invasive tool that may aid early diagnosis of AD. However, the use of EEG signal analysis in order to improve the diagnosis of AD is a complex problem where, despite significant advances, a number of fundamental questions remain open (Elgendi et al., 2011).

Considering now the characteristics of the EEG, since the non-stationary nature of the signal is generally well-recognized (see, for instance Akin, 2002), decomposition using a fast Fourier transform (FFT) with sliding windows and the wavelet transforms have been the most popular techniques employed to capture the spectral properties of EEG (Darvishi and Al-Ani, 2007; Dauwels et al., 2010). However, linear transformation methods fail to address the non-linear characteristics of the EEG signal (Stam, 2005). Therefore, non-linear dynamic approaches have been attempted as well, mostly involving computationally complex time series analysis (Jeong, 2004). Several other aspects of non-linear modeling and analysis in this context have also been studied in the literature (see, for instance, Stam, 2005 for a review). These include frameworks based on a neural mass model (Valdes et al., 1999; Huang et al., 2011), coupled oscillators (Baier et al., 2005; Leistritz et al., 2007), continuum models (Kim et al., 2007), non-linear non-stationary models (Celka and Colditz, 2002; Rankine et al., 2007), random neural networks (Acedo and Morano, 2013), and chaotic phenomena and stability aspects (Rodrigues et al., 2007; Dafilis et al., 2009). Stochastic approaches based on Markov chain Monte Carlo methods (Hettiarachchi et al., 2012) and Markov process amplitude (Wang et al., 2011) that take into account the inherent randomness of the EEG signal have also been reported. In the same vein, limit cycle oscillators (Hernandez et al., 1996; Burke and Paor, 2004) as well as stochastic synchronization (Bressloff and Lai, 2011) and stochastic approximation (Fell et al., 2000; Sun et al., 2008) methods have been considered in EEG modeling. Notably, limit cycle behavior at each of the brain frequency bands appears to provide a more accurate representation of the EEG signal than one based on chaotic phenomena.

Some of the most important features in non-linear dynamic and stochastic approaches are signal information content and complexity as measured using various forms of information entropy. Measures such as Shannon entropy (Shannon, 1948) characterize the information content in a signal and higher entropy corresponds to increased randomness and chaotic behavior (Abasolo et al., 2006). Importantly, one observes that, with respect to the EEG signal, higher information content correlates with better brain function (Shin et al., 2006). Furthermore, it has been reported that variations in information entropic measures may be used to detect functional abnormalities in the brain caused by disorders or injuries (Slobounov et al., 2009). Hence, information content of the EEG signal, characterized by information-entropic measures, may be expected to be important in identifying distinct states of the brain. This is further reinforced by the recent results of McBride and colleagues on the role of information entropic and spectral analysis in the study of the early stages of Alzheimer's disease and mild Traumatic Brain Injury McBride et al. (2013a,b, 2014).

Entropy may also be utilized to measure signal complexity. For instance, embedding entropy provides information about how the EEG signal fluctuates in time by comparing the time series with a delayed version of itself (Abasolo et al., 2006). Moreover, the concept of approximate entropy was introduced as a measure of system complexity (Pincus, 1991) and has been applied to brain wave signals (Quiroga et al., 2001). However, the approximate entropy measure suffers from drawbacks such as bias and inconsistency (Xu et al., 2010). Hence, the notion of sample entropy was introduced (Richman and Moorman, 2000) as an improvement over approximate entropy.

In recent work, the authors proposed a phenomenological model of the EEG signal based on the dynamics of a stochastic, coupled, Duffing- van der Pol oscillator network (Ghorbanian et al., 2015). An optimization scheme was adopted to match model output with actual EEG data obtained from healthy subjects in the two distinct resting eyes-open (EO) and eyes-closed (EC) conditions and it was shown that the actual EEG signals in both cases were distinct realizations of the model with qualitatively different non-linear dynamic characteristics. Moreover, the model output and the actual EEG data were shown to be in good agreement in terms of both the power spectra (frequency content) and Shannon entropy (information content).

In the present effort, we improve the model introduced in Ghorbanian et al. (2015) by matching the sample entropy of the model output and EEG signal to capture its complexity. A global optimization routine is employed in order to match the output of with EEG recordings in terms of power spectrum, Shannon entropy, and sample entropy. The EEG signals were recorded under resting EC and EO conditions in an earlier pilot study of Alzheimer's disease (AD) patients vs. age-matched healthy control (CTL) subjects (Ghorbanian et al., 2013). The model parameters obtained for the oscillators representing EC and EO EEG signals for CTL and AD patients are compared in order to establish statistically significant, distinct models for AD and CTL subjects under each condition. In addition, we present new results from a phase space reconstruction analysis of the model output to match the actual EEG signal. The results indicate that the analytical model effectively captures the frequency spectrum and non-linear characteristics of the EEG signal in terms of complexity and information content. Furthermore, it is shown that the addition of sample entropy significantly enhances the model performance in terms of complexity and non-linear dynamic characteristics, as demonstrated by phase space reconstruction. The results suggest exciting new pathways to develop better tools for distinguishing pathological and normal brain states in AD and perhaps other neurological diseases and disorders.

The rest of the article is set as follows. Details of the EEG recordings, the analytical model, the optimization scheme and the phase space reconstruction technique are provided in Section 2. The results are presented in Section 3 and discussed in Section 4. The articles concludes with comments on further research in Section 5.

2. Materials and Methods

2.1. EEG Recording Blocks

Twenty six AD patients and healthy age-matched CTL subjects were selected for this study (“A Brain-Computer Interface for Diagnosing Brain Function,” Aspire IRB, Human Subject Protocol Number PDMC-001, approved on October 7, 2010). Of the 26 subjects selected, one withdrew and one did not qualify as AD or CTL. Subjects were asked to relax and wear an EEG recording headset during alternating blocks of EC and EO followed by a variety of cognitive and auditory tasks and a final EC-EO resting period.

The EEG signals were recorded through a single-dry electrode device at position Fp1 (based on a 10–20 electrode placement system) with a Bluetooth enabled telemetric headset. The headset's effective sample rate is 125 Hz. Frequencies below 1 Hz and above 60 Hz (near Nyquist frequency) were filtered out by the device hardware. On comparison of the EEG recordings by the device with those from other widely accepted devices, frequencies within 2–30 Hz were deemed to be very accurate.

The recording device eliminated frequently observed artifacts including line noise. Other artifacts were mainly due to eye- and muscle-movements, which are common at Fp1 location and can be clearly identified by their high amplitudes compared to true EEG signal recordings during resting states. These artifacts were removed using a simple artifact detection that eliminated any part of the signal greater than 4.5σ (standard deviation). The algorithm also reconstructed the nulled samples using FFT interpolation of the trailing and subsequent recorded data (Ghorbanian et al., 2013).

The EEG recordings in this study were obtained from subjects in an AD pilot study with 14 control (CTL) subjects and 10 Alzheimer's Disease (AD) patients presented in our earlier work (Ghorbanian et al., 2013). Recording blocks of 40-s duration (approximately 5000 sample signals) from resting eyes-closed (EC) and eyes-open (EO) conditions were selected. In all, 60 random blocks were selected from the pilot study: 40 blocks from control CTL subjects (20 EC and 20 EO) and 20 blocks from AD subjects (10 EC and 10 EO). Note that, the smaller number of AD patients along with smaller number of AD patient recording sessions that were were not dominated by artifacts resulted in the selection of smaller AD sample size.

2.2. EEG Features

The time-varying power spectrum in each of the major brain EEG frequency bands was calculated using short time fast Fourier transform (FFT) with sliding window, since a good model must produce signals that can match EEG's frequency content. Specifically, the power spectrum was computed in seven ranges: lower δ (1–2 Hz), upper δ (2–4 Hz), θ (4–8 Hz), α (8–13 Hz), lower β (13–20 Hz), upper β (20–30 Hz), and γ (30–60 Hz). However, lower δ and γ band powers, which happen to have little power, were ignored due to unreliability of the device in those frequency ranges.

Shannon entropy was measured based on a sliding temporal window technique. A temporal window was defined to slide along the signal time representation with a sliding step (interval or bin) to sample a part of the signal. A discrete entropy estimator was applied, in which 10 uniform intervals equally divided the range of the normalized observed signal. Then the probability that the sampled signal belongs to the interval is the ratio between the number of the samples found within each interval and the total number of samples of the signal. The Shannon entropy is then calculated based on these probabilities (Shin et al., 2006), separately for each 40-s EEG recording block (5000 samples).

Sample entropy (SE) is the negative natural logarithm of the conditional probability that two sequences of a time series, similar for m points, remain similar at the next point. For given N data points from a time series, [x(1), x(2), …, x(N)], we calculated SE of each 40-s EEG recording block (5000 samples) by the statistic (Abasolo et al., 2006):

$\begin{array}{cc}\text{SE}(m,r,N)=\left\{-\ln \left[\frac{U^{m+1}(r)}{U^m(r)}\right]\right\},& \left(1\right)\end{array}$

where m is the run length, r is the tolerance window size, and

$\begin{array}{cc}{U}^m(r)=\frac{1}{(N-m)(N-m-1)}{\displaystyle \sum _{i=1}^{N-m}{U}_i}.& \left(2\right)\end{array}$

In the above equation, U_i indicates the number of k's (1 ≤ k ≤ N − m) such that the Euclidean distance between X_m(i) and X_m(k), k ≠ i, is less than or equal r and X_m(i) = [x(i), x(i + 1), …, x(i + m − 1)].

Generally, large m or small r values result in number of matches being too small for confident estimation of the conditional probability and vice versa (Lake and Moorman, 2011). In this study, we used m = 2 and r = 0.25σ based on the consistency of the results and recommended ranges in previous studies (Richman and Moorman, 2000; Xu et al., 2010).

2.3. Stochastic Coupled Non-linear Oscillators

We recall that the EEG has been modeled in the literature taking into account characteristics including non-linearity (both chaotic and non-chaotic), non-stationarity, and randomness of the signal (Fell et al., 2000; Rankine et al., 2007; Sun et al., 2008). The EEG has also been studied as the manifestation of underlying limit cycle oscillations at a given frequency and other such periodic solutions (Hernandez et al., 1996; Burke and Paor, 2004). While inspired by the above, the authors were fundamentally motivated to develop models that can better reproduce the significant linear and non-linear characteristics of actual EEG signals. Hence, we proposed a phenomenological model of the EEG based on a coupled system of Duffing—van der Pol oscillators subjected to white noise excitation (Ghorbanian et al., 2015). This particular oscillator was selected because the Duffing non-linearity allows a system with only two oscillators capture the major brain frequency spectra and van der Pol non-linearity provides self-excited limit cycle behavior which have been previously reported for each major brain frequency bands (Burke and Paor, 2004).

We consider a phenomenological model of the EEG based on a coupled system of Duffing—van der Pol oscillators subject to white noise excitation, as shown in Figure 1. The equations for the model may be written as:

$\begin{array}{cc}\{\begin{array}{l}{\ddot{x}}_1+(k_1+k_2)x_1-k_2{x}_2=-b_1{x}_1^3-b_2{(x_1-x_2)}^3\hfill \\ +{ϵ}_1{\dot{x}}_1(1-x_1^2),\hfill \\ {\ddot{x}}_2-k_2{x}_1+k_2{x}_2=b_2{(x_1-x_2)}^3\hfill \\ +{ϵ}_2{\dot{x}}_2(1-x_2^2)+\mu dW,\hfill \end{array}& \left(3\right)\end{array}$

where x_i, $\dot{x}$_i, $\ddot{x}$_i, i = 1, 2 are positions, velocities, and accelerations of the two oscillators, respectively. Parameters k_i, b_i, ϵ_i, i = 1, 2 are, respectively, linear stiffness, cubic stiffness, and van der Pol damping coefficient of the two oscillators. Parameters b_is indicate the strength of the Duffing non-linearity resulting in multiple resonant frequencies while ϵ_is indicate the strength of van der Pol non-linearity and determine the extent of self-excitation and the shape of the resulting limit cycle. Parameter μ represents the intensity of white noise and dW is a Wiener process (Gardiner, 1985; Higham, 2001) representing the additive noise in the stochastic differential system. The input excitation to the system is provided through μdW. The output may be selected as any combination of the positions and velocities to mimic an EEG signal. Note that, the Euler-Maruyama method (Higham, 2001) was selected to integrate the stochastic differential equations in Equation (3) since standard numerical integration methods are not applicable.

FIGURE 1

Figure 1. Schematic of the stochastic coupled Duffing—van der Pol oscillators.

2.4. Optimization Formulation

We have selected the velocity of the second oscillator as the system output approximating the EEG signal since it is directly impacted by the noise. A global optimization search method based on a multi-start algorithm (Ugray et al., 2007) was adopted to determine the oscillator model parameters that can produce the output matching various EEG signals. The optimization objective function was selected as the root mean squared of the errors in power spectrum of each selected brain frequency bands plus weighted values of the errors in absolute Shannon and sample entropies. Hence, the optimization goal is error minimization:

$\begin{array}{cc}\underset{p}{\min }J=\sqrt{{\displaystyle \sum _{j=1}^m{(P_E{}_j-P_O{}_j)}^2}}+w_1|S_E-S_O|+w_2|S{P}_E-S{P}_O|,& \left(4\right)\end{array}$

where J is the objective function, p = [k₁, k₂, b₁, b₂, ϵ₁, ϵ₂, μ] the decision variables, P_Ej and P_Oj the powers in the major brain frequency bands for the normalized EEG signal and the model output, respectively, m is number of frequency bands (m = 7), S_E and S_O the Shannon entropies of the EEG signal and the model output, respectively, SP_E and SP_O the sample entropies of the EEG signal and the model output, respectively, w₁ and w₂ are weighting factor for absolute Shannon and sample entropies, respectively, and | | represents absolute value. The weighting factors w₁ and w₂ are required to give equal importance to power spectrum and entropy characteristics of the signal. Note that the magnitude of the output signals are matched through normalization of both the model output and the EEG signal with respect to their standard deviations.

The objective function minimization is subject to equality constraints represented by the state (Equation 3) and inequality constraints represented by the decision variable lower and upper bounds:

$\begin{array}{cc}\begin{array}{l}0<k_i\le 1e4,0<b_i\le \frac{1}{2}{k}_i,0<{ϵ}_i\le \frac{1}{3}{k}_i,\\ i=1,2,0\le \mu \le 2.\end{array}& \left(5\right)\end{array}$

The constraints for b_i's and ϵ_i's were imposed to avoid the chaotic regime (Li et al., 2006) and provide a periodic stochastic response. Noise intensity is also constrained to avoid a response dominated by random noise. The initial guesses for the global optimization search are randomly generated within the bounds defined in Equation (5).

The stochastic component was introduced as white noise, which was generated through a normally distributed random variable and applied to the model via Wiener process. A new random process was generated and applied to the model during integration of the equations, at each iteration of the optimization algorithm.

2.5. Statistical Analysis

A key objective of the phenomenological modeling in this work is the ability to establish a correspondence between variations in model parameters and the variations in the data obtained from different physiological conditions. Hence, the parametric unpaired t-test and non-parametric Wilcoxon rank sum statistical testing methods were employed to determine the relative significance of the model parameters. Furthermore, Bonferroni correction was applied due to multiple comparisons problem and adequacy of sample sizes for statistical tests were established using power analysis.

2.6. Phase Space Reconstruction

In addition to matching Shannon and sample entropies of the model output and EEG signal through the optimization process, it is of interest to investigate matching other features such as the phase plot which plays a significant role in non-linear time series analysis (Kantz and Schreiber, 2004). It is known that any dynamic system can be completely recovered in the phase space, which maybe reconstructed from the measured time domain response of the system (Nie et al., 2013). While phase space consists of velocity and position variables for a mechanical system, in the case where just the time representation of a signal is available, a phase space reconstruction technique based on the method of delays is used (Kantz and Schreiber, 2004).

The main idea is that one does not need the derivatives to form a coordinate system in which to capture the structure of phase space, but instead one could directly use the lagged variables:

$\begin{array}{cc}x(n+T)=x(t_0+(n+T)\Delta {\tau }_s),& \left(6\right)\end{array}$

where x(n) is the nth sample of the time series, Δτ_s the time step, and T the delay integer to be determined. Then, a vector of (embedding) dimension d may be constructed using the time lags as:

$\begin{array}{cc}[x(n),x(n+T),x(n+2T),\cdots ,x(n+(d-1)T)].& \left(7\right)\end{array}$

Time-delay embedding is probably one of the best systematic methods for converting scalar data to multidimensional phase space (Abarbanel et al., 1993; Burke and Paor, 2004; Nie et al., 2013). An appropriate and successful reconstruction depends on the choice of both time delay T and the embedding dimension d (Nie et al., 2013).

In this study, the appropriate value of the time lag was determined using the average mutual information method applied to each EEG recording block. The idea behind mutual information is to identify the amount of information that can be learned about a measurement at one time from a measurement taken at another time. Consider the time series nth sample x(n) and its value after time delay T with the associated probability distributions of P(x(n)) and P(x(n + T)), respectively. The average information which can be obtained about x(n + T) from x(n) is given by the mutual information of the two measurements (Abarbanel et al., 1993; Mizrach, 1996):

$\begin{array}{cc}I(x(n),x(n+T))={\log }_2\left[\frac{P(x(n),x(n+T))}{P(x(n))P(x(n+T))}\right],& \left(8\right)\end{array}$

where P(x(n), x(n + T)) is the joint probability of the measurements x(n) and x(n + T) calculated using a binning-based method, in which 20 uniform intervals divided the range of the measurements equally. The average mutual information between measurements of any value x(n) and x(n + T) is the average over all possible measurements of I(x(n), x(n + T)) (Abarbanel et al., 1993):

$\begin{array}{cc}I(T)={\displaystyle \sum _{x(n),x(n+T)}P}(x(n),x(n+T))I(x(n),x(n+T)).& \left(9\right)\end{array}$

If T is too small, the measurements x(n) and x(n + T) will have too much overlap. However, if T is too large, then I(T) will approach zero and nothing relates x(n) to x(n + T). It is suggested that the proper T can be chosen as the first minimum of I(T) which is not necessarily optimal but has been shown to work well (Abarbanel et al., 1993; Nie et al., 2013). If in a case, no minima exists for I(T), the choice of T = 1 or 2 has been suggested (Abarbanel et al., 1993).

After specifying the correct time delay T, an appropriate embedding dimension, d, should also be found for the phase space reconstruction. If d is too small, the trajectories will not be unique. On the other hand, too large a d will result in additional computational cost by requiring extra dimensions (Nie et al., 2013).

3. Results

The optimization algorithm was separately applied to determine the model parameters (decision variables) for each of the 60 selected EEG signals using the weighting factors w₁ = w₂ = 0.35. These weighting factors give equal importance to the entropy measures and power spectrum. We then categorized the resulting 60 set of model parameters into four groups based recording conditions and subject diagnosis: EC-CTL, EO-CTL, EC-AD, and EO-AD.

3.1. Healthy Eyes-Closed and Eyes-Open Results

Initially, we studied the models derived for the EC and EO EEG signals of CTL subjects for validation purposes. The means and standard deviations of the optimal values of the model parameters for EC-CTL and EO-CTL EEG signals are listed in Table 1. The p-values from the two statistical tests and non-parametric method after Bonferroni corrections indicate that the differences between of all parameters of the two models are strongly statistically significant with the exception of noise intensity. Note that, μ is also found statistically significant using t-test but is slightly off when the non-parametric method is used.

TABLE 1

Table 1. Optimal parameters of the Duffing—van der Pol oscillator for EC and EO of CTL subjects (N = 40) and the p-values from unpaired t-test, Wilcoxon rank sum test, and Bonferonni correction.

In order to ensure that adequate sample sizes are used, the minimum required difference between means of two groups of data for each parameter are computed. As expected due to very small p-values, the sample size for statistical testing is found to be sufficient with more than 99.9% power for all parameters except noise intensity μ, which was not found to be statistically significant using the non-parametric method.

Power spectrums of the optimal stochastic oscillator model output and EEG signals for the EC and EO cases of CTL subjects are presented in Figure 2 where θ, α, and β band powers show excellent agreement. The comparison revealed that, as expected, the optimal model is closely following the α-band dominance in the EC cases. While, in the EO cases, the optimal model follows a more flat frequency distribution from upper δ to lower β frequency bands. Furthermore, Shannon and SE values of the EEG signals and the model outputs for the EC and EO cases show close agreement. Shannon entropy values were 1.80 ± 0.08 and 1.92 ± 0.08 for EC EEG and model output, respectively, and 1.71 ± 0.11 and 1.57 ± 0.15 for EO EEG and model output, respectively. While, SE values were 1.04 ± 0.20 and 1.17 ± 0.22 for EC EEG and model output, respectively, and 0.97 ± 0.20 and 1.20 ± 0.18 for EO EEG and model output, respectively. These results show a significant improvement over our previous model where only Shannon entropy was used (Ghorbanian et al., 2015). The improvement is clearly observed in the the power spectra of sample EC and EO EEG signals and their corresponding optimal model outputs, respectively shown in Figures 3, 4. Both figures demonstrate more distributed spectra of the model outputs with similar noise complexities to the actual EEG signals when SE is added to the objective function; i.e., power spectra of the signals without matching of SE have very discrete peaks unlike the EEG.

FIGURE 2

Figure 2. Comparison of major brain frequency band mean powers of CTL EEG signals and optimal oscillator model output; EC (top), EO (bottom).

FIGURE 3

Figure 3. Power spectrum of a sample CTL EC (top) EEG signal; (middle) output of stochastic oscillator model using Shannon entropy; (bottom) output of stochastic oscillator model using Shannon and sample entropies.

FIGURE 4

Figure 4. Power spectrum of a sample CTL EO (top) EEG signal; (middle) output of stochastic oscillator model using Shannon entropy; (bottom) output of stochastic oscillator model using Shannon and sample entropies.

The impact of SE to match signal complexity is further demonstrated through phase plot reconstruction of the time series. Average mutual information for a sample EC EEG signal and outputs of the optimal stochastic oscillator models are shown in Figure 5 as a function of lag time. The first minimum occurs at T = 5 lag samples for both the EEG signal and the optimal model derived with both Shannon and sample entropies while T = 3 for the output of the model derived solely based on Shannon entropy. The resulting reconstructed phase plots of the EC EEG signal and the outputs of the two optimal models are presented in Figure 6. Clearly, the reconstructed phase plots of the EEG and the output of the model derived using both Shannon and sample entropies, display similar behavior. While the output of the model derived using only Shannon entropy is qualitatively different form the EEG signal in terms of complexity and noise. Indeed this result provides further affirmation that the stochastic Duffing—van der Pol model yields an output that matches the actual EEG data in terms of non-linear characteristics observed in the phase space.

FIGURE 5

Figure 5. Average mutual information for a sample CTL EC (top) EEG signal; (middle) output of stochastic oscillator model using Shannon entropy; (bottom) output of stochastic oscillator model using Shannon and sample entropies.

FIGURE 6

Figure 6. Reconstructed phase plot of a sample CTL EC (top) EEG signal; (middle) output of stochastic oscillator model using Shannon entropy; (bottom) output of stochastic oscillator model using Shannon and sample entropies.

3.2. Alzheimer's Disease vs. Control Results

Next, we studied the models derived for the EC and EO EEG signals of AD subjects. The mean and standard deviation of the optimal values of the model parameters for EC-AD and EO-AD EEG signals are listed in Table 2 along with the p-values from the two statistical tests and the non-parametric test after Bonferroni corrections indicating that the differences between only the first four parameters of the two models are statistically significant. Next, we separately compared the model parameters of EC and EO EEG signals of CTL subjects with those AD patients.

TABLE 2

Table 2. Optimal parameters of the Duffing—van der Pol oscillator model for EC and EO of AD subjects (N = 20) and the p-values from unpaired t-test, Wilcoxon rank sum test, and Bonferonni correction.

Table 3 lists the p-values from the two statistical testing methods and the non-parametric method after Bonferroni corrections comparing CTL vs. AD subjects under separate EC and EO conditions. The results indicate that the difference between all model parameters of CTL and AD subjects under EC condition are strongly statistically significant except for noise intensity. Again, μ is also found statistically significant using t-test but is slightly off when non-parametric method is used. The difference between the model parameters of CTL and AD subjects under EO condition are not, however, as strong, though they are still mostly statistically significant. In the EO case, parameter μ is not statistically significant using either method and t-test does not find b₁ to be statistically significant either.

TABLE 3

Table 3. The p-values from unpaired t-test, Wilcoxon rank sum test, and Bonferonni correction for comparison of model parameters between AD (N = 20) and CTL (N = 40) subjects.

The power analysis results for 90%, 95%, 99%, and 99.9% for two statistical are listed in Tables 4, 5 for EC and EO cases, respectively. The actual difference between means are given within parentheses following each parameter. The results indicated that our sample size for statistical testing in EC case between AD and CTL subjects was sufficient for all parameters except μ with more than 99.9% power. However, in the EO case, only sample size for parameters b₂ and ϵ₂ has more than 99.9% confidence and k₂ shows a 90% confidence. The sample size for the remaining parameters did not provide sufficient confidence.

TABLE 4

Table 4. Minimum required difference between model parameter mean values of EC AD vs. EC CTL for various desired powers of statistical tests.

TABLE 5

Table 5. Minimum required difference between model parameter mean values of EO AD vs. EO CTL for various desired powers of statistical tests.

Power spectrums of the optimal stochastic oscillator model output and EEG signals for the EC and EO cases of AD subjects are presented in Figure 7 where again θ, α, and β band powers show excellent agreement. The comparison revealed that the optimal model was closely and correctly slightly θ-band dominated in the EC cases for AD subjects (Ghorbanian et al., 2013). While, in the EO cases, the optimal model followed the more flat frequency distribution. Again, it should be noted that the higher error rates are related to those frequency bands with lower powers. Furthermore, Shannon and SE values of the EEG signals and the model outputs for the EC and EO cases show close agreement. Shannon entropy values were 1.78 ± 0.04 and 1.70 ± 0.10 for EC EEG and model output, respectively, and 1.63 ± 0.32 and 1.62 ± 0.27 for EO EEG and model output, respectively. While, SE values were 1.06 ± 0.19 and 1.17 ± 0.21 for EC EEG and model output, respectively, and 1.02 ± 0.39 and 1.29 ± 0.24 for EO EEG and model output, respectively.

FIGURE 7

Figure 7. Comparison of major brain frequency band mean powers of AD EEG signals and optimal oscillator model output; EC (top), EO (bottom).

Power spectra of outputs of the optimal stochastic oscillator models and EEG signals for sample EC and EO cases of AD subjects are presented in Figures 8, 9. Again, it is clear that the addition of SE to the objective function results in output signals with power spectra patterns which are much more similar to the EEG signal in terms of distribution and noise complexity. As expected, the power spectrum plots demonstrated that the EC EEG signals from AD subjects were slightly θ band dominated unlike α band dominance of EC EEG recordings from CTL subjects.

FIGURE 8

Figure 8. Power spectrum of a sample AD EC (top) EEG signal; (middle) output of stochastic oscillator model using Shannon entropy; (bottom) output of stochastic oscillator model using Shannon and sample entropies.

FIGURE 9

Figure 9. Power spectrum of a sample AD EO (top) EEG signal; (middle) output of stochastic oscillator model using Shannon entropy; (bottom) output of stochastic oscillator model using Shannon and sample entropies.

4. Discussion

Power spectra of the optimal stochastic oscillator model output and EEG signals show excellent agreement in the brain's major frequency bands. The comparison revealed that the optimal model is closely following the α-band dominance in EC recordings for the control subjects. Furthermore, the model for EC recordings of AD patients closely followed θ-band power dominance indicating the slowing of the EEG signal for these patients. In the EO cases, the optimal model, as expected, followed a more flat frequency distribution from upper δ to lower β frequency bands for both AD and CTL subjects. Further evidence of robustness of the the models derived in this study is that the models derived for healthy subject EC and EO EEG signals in our earlier study (Ghorbanian et al., 2015) fall within the same distributions obtained for the CTL subjects in the clinical study.

Moreover, Shannon and SE values of the EEG signals and the model outputs for the EC and EO cases show close agreement for both CTL and AD subjects. However, the difference between the entropy values of the CTL subjects and AD patients were not statistically significant for neither the EEG signal nor the model output. This aspect needs to be further studied since EEG signals from AD patients may be expected to have lower complexity and thus lower entropy values.

The contributions of the article are as follows. Firstly, the objective function of the optimization scheme that yields model parameters based on comparison with actual EEG data in our previous work was extended to include both Shannon and sample entropies, with the latter being a measure of signal complexity. The procedure yielded model outputs that were in agreement with the actual EEG signals with respect to the frequency content (power spectra), information content (Shannon entropy) and complexity (sample entropy). It was shown that the addition of SE significantly enhances the performance of the optimal model in terms of both power spectrum and non-linear characteristics displayed through phase space reconstruction. The results demonstrate the feasibility of stochastic non-linear oscillator models which can be further studied for greater insight into EEG signal dynamic characteristics.

Secondly, the model parameter differences for EC and EO EEG recordings were statistically significant leading to qualitatively and quantitatively distinct realizations of the underlying models for the cases considered. This is a key result of the work since it verifies that distinct models represent the EEG signals recorded under different brain states. Potentially, this could lead to unique models for different brain disorders and injuries.

Thirdly, the study provided unique models for EC and EO EEG recordings from AD patients. The results showed that almost all of the model parameters were statistically significant for the EC and EO cases when comparing the AD and CTL subjects. Moreover, the power spectrum plots showed a good match between the generated signal from the stochastic model and the actual EEG signal from AD patients. However, the results for the EC case of AD were more accurate and reasonable than the results of EO cases mainly due to the ability of the optimization scheme to provide a better match in EC cases. The important conclusion here is that unique stochastic non-linear oscillator models can be developed to represent EEG signals from patients with a brain disorder.

Of particular interest is the potential connection between our model and the neural mass models studied in the literature. For instance, characterization of functional connectivity between remote cortical areas has been studied using neural mass models (David and Friston, 2003; David et al., 2004). These and other efforts (Sotero et al., 2007) represent intriguing attempts to capture actual neural dynamics using coupled oscillator models and suggest that, after all, models such as the one discussed in this article may be of broader scope than being purely phenomenological. Extrapolating further, it would then be of immense interest to understand the manifestation of phenomena such as synchronization (Mirollo and Strogatz, 1990) within the framework of our model and the implications for EEG characterization.

5. Conclusions

In this article, we presented results that further develop our recent work on modeling the EEG signal as the response of a stochastic, coupled Duffing—van der Pol system of two oscillators. The results presented verify that unique and statistically significant stochastic Duffing—van der Pol oscillator models represent EEG recorded from AD patients vs. health controls. Overall, the results presented in this article further affirm the efficacy of a stochastic Duffing—van der Pol oscillator network model in capturing the key characteristics of actual EEG data under different brain states as well as brain conditions in terms of healthy controls vs. patients with a brain disorder. The validation provided by the results certainly motivates further research toward improving the analytical model and testing it against larger data sets. Furthermore, the results suggest that the modeling approach could potentially help develop novel diagnostic and interventional tools for neurological diseases and disorders.

Conflict of Interest Statement

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.

References

Abarbanel, H., Brown, R., Sidorowich, J., and Tsimring, L. S. (1993). Analysis of observed chaotic data. Rev. Mod. Phys. 65, 1331–1391.

Google Scholar

Abasolo, D., Hornero, R., Espino, P., Alvarez, D., and Poza, J. (2006). Entropy analysis of the EEG background activity in Alzheimer's disease patients. Physiol. Meas. 27, 241–253. doi: 10.1088/0967-3334/27/3/003

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Acedo, L., and Morano, J. (2013). Brain oscillations in a random neural network. Math. Comput. Model. 57, 1768–1772. doi: 10.1016/j.mcm.2011.11.028

CrossRef Full Text | Google Scholar

Akin, M. (2002). Comparison of wavelet transform and FFT methods in the analysis of EEG signals. J. Med. Syst. 26, 241–247. doi: 10.1023/A:1015075101937

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Baier, G., Hermann, T., and Muller, M. (2005). “Polyrhythmic organization of coupled nonlinear oscillators,” in International Conference on Information Visualization (Washington, DC), 5–10.

Bressloff, P., and Lai, Y. (2011). Stochastic synchronization of neuronal populations with intrinsic and extrinsic noise. J. Math. Neurosci. 1, 1–28. doi: 10.1186/2190-8567-1-2

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Burke, D., and Paor, A. D. (2004). A stochastic limit cycle oscillator model of the EEG. Biol. Cyber. 91, 221–230. doi: 10.1007/s00422-004-0509-z

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Celka, P., and Colditz, P. (2002). Nonlinear nonstationary wiener model of infant EEG seizures. IEEE Trans. Biomed. Eng. 49, 556–564. doi: 10.1109/TBME.2002.1001970

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Dafilis, M., Frascoli, F., Cadusch, P., and Liley, D. (2009). Chaos and generalised multistability in a mesoscopic model of the electroencephalogram. Physica D 238, 1056–1060. doi: 10.1016/j.physd.2009.03.003

CrossRef Full Text | Google Scholar

Darvishi, S., and Al-Ani, A. (2007). “Brain-computer interface analysis using continuous wavelet transform and adaptive neuro-fuzzy classifier,” in International Conference of the IEEE Engineering in Medicine and Biology Society (Lyon), 3220–3223.

Google Scholar

Dauwels, J., Vialatte, F., and Cichocki, A. (2010). Diagnosis of Alzheimer's disease from EEG signals: where are we standing. Curr. Alzheimer Res. 7, 487–505. doi: 10.2174/156720510792231720

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

David, O., Cosmelli, D., and Friston, K. J. (2004). Evaluation of different measures of functional connectivity using a neural mass model. Neuroimage 21, 659–673. doi: 10.1016/j.neuroimage.2003.10.006

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

David, O., and Friston, K. J. (2003). A neural mass model for meg/eeg: coupling and neuronal dynamics. Neuroimage 20, 1743–1755. doi: 10.1016/j.neuroimage.2003.07.015

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Elgendi, M., Vialatte, F., Cichocki, A., Latchoumane, C., Jeong, J., and Dauwels, J. (2011). “Optimization of EEG frequency bands for improved diagnosis of Alzheimer disease,” in International Conference of the IEEE Engineering in Medicine and Biology Society (Boston: MA), 6087–6091.

Google Scholar

Fell, J., Kaplan, A., Darkhovsky, B., and Roschke, J. (2000). EEG analysis with nonlinear deterministic and stochastic methods: a combined strategy. Acta Neurobiol. Exp. 60, 87–108.

PubMed Abstract | Full Text | Google Scholar

Gardiner, C. (1985). Handbook of stochastic methods for Physics, Chemistry and the Natural Sciences. New York, NY: Springer.

Ghorbanian, P., Devilbiss, D., Verma, A., Bernstein, A., Hess, T., Simon, A., et al. (2013). Identification of resting and active state EEG features of Alzheimer's disease using discrete wavelet transform. Ann. Biomed. Eng. 41, 1243–1257. doi: 10.1007/s10439-013-0795-5

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Ghorbanian, P., Ramakrishnan, S., Whitman, A., and Ashrafiuon, H. (2015). A phenomenological model of EEG based on the dynamics of a stochastic Duffing–van der Pol oscillator network. Biomed. Signal Process. Control 15, 1–10. doi: 10.1016/j.bspc.2014.08.013

CrossRef Full Text | Google Scholar

Hernandez, J., Valdes, P., and Vila, P. (1996). EEG spike and wave modeled by a stochastic limit cycle. Neuroreport 7, 2246–2250.

PubMed Abstract | Full Text | Google Scholar

Hettiarachchi, I., Mohamed, S., and Nahavandi, S. (2012). “A marginalised Markov chain Monte Carlo approach for model based analysis of EEG data,” in IEEE International Symposium on Biomedical Imaging (Barcelona), 1539–1542.

Google Scholar

Higham, D. J. (2001). An algorithmic introduction to numerical simulation of stochastic differential equations. SIAM Rev. 43, 525–546. doi: 10.1137/S0036144500378302

CrossRef Full Text | Google Scholar

Huang, G., Zhang, D., Meng, J., and Zhu, X. (2011). Interactions between two neural populations: a mechanism of chaos and oscillation in neural mass model. Neurocomputing 74, 1026–1034. doi: 10.1016/j.neucom.2010.11.019

CrossRef Full Text | Google Scholar

Jeong, J. (2004). EEG dynamics in patients with Alzheimer's disease. Clin. Neurophysiol. 15, 1490–1505. doi: 10.1016/j.clinph.2004.01.001

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Kantz, H., and Schreiber, T. (2004). Nonlinear Time Series Analysis. New York, NY: Cambridge University Press.

Google Scholar

Kim, J., Shin, H., and Robinson, P. (2007). Compact continuum brain model for human Electroencephalogram. Complex Syst. II 6802, T1–T8. doi: 10.1117/12.759005

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Lake, D., and Moorman, J. (2011). Accurate estimation of entropy in very short physiological time series: the problem of atrial fibrillation detection in implanted ventricular devices. Am. J. Physiol. Heart Circ. Physiol. 300, H319–H325. doi: 10.1152/ajpheart.00561.2010

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Leistritz, L., Putsche, P., Schwab, K., Hesse, W., Susse, T., Haueisen, J., et al. (2007). Coupled oscillators for modeling and analysis of EEG/MEG oscillation. Biomed. Tech. 52, 83–89. doi: 10.1515/BMT.2007.016

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Li, Z., Xu, W., and Zhang, X. (2006). Analysis of chaotic behavior in the extended duffing–van der Pol system subject to additive non-symmetry biharmonical excitation. Appl. Math. Comput. 183, 858–871. doi: 10.1016/j.amc.2006.06.033

CrossRef Full Text | Google Scholar

McBride, J., Zhao, X., Munro, N., Smith, C., Jicha, G., Hively, L., et al. (2014). Spectral and complexity analysis of scalp EEG characteristics for mild cognitive impairment and early Alzheimer's disease. Comput. Methods Programs Biomed. 114, 153–163. doi: 10.1016/j.cmpb.2014.01.019

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

McBride, J., Zhao, X., Munro, N., Smith, C., Jicha, G., and Jiang, Y. (2013a). Resting EEG discrimination of early stage Alzheimer's disease from normal aging using inter-channel coherence network graphs. Ann. Biomed. Eng. 41, 1233–1242. doi: 10.1007/s10439-013-0788-4

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

McBride, J., Zhao, X., Nichols, T., Vagnini, V., Munro, N., Berry, D., et al. (2013b). Scalp eeg-based discrimination of cognitive deficits after traumatic brain injury using event-related tsallis entropy analysis. IEEE Trans. Biomed. Eng. 60, 90–96. doi: 10.1109/TBME.2012.2223698

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Mirollo, R. E., and Strogatz, S. H. (1990). Synchronization of pulse-coupled biological oscillators. SIAM J. Appl. Math. 50, 1645–1662.

Google Scholar

Mizrach, B. (1996). Determining delay times for phase space reconstruction with application to the FF/DM exchange rate. J. Econ. Behav. Organ. 30, 369–381.

Google Scholar

Nie, Z., Hao, H., and Ma, H. (2013). Structural damage detection based on the reconstructed phase space for reinforced concrete slab: experimental study. J. Sound Vibrat. 332, 1061–1078. doi: 10.1016/j.jsv.2012.08.024

CrossRef Full Text | Google Scholar

Pincus, S. (1991). Approximate entropy as a measure of system complexity. Proc. Natl. Acad. Sci. U.S.A. 88, 2297–2301.

PubMed Abstract | Full Text | Google Scholar

Quiroga, R., Rosso, O., Basar, E., and Schurmann, M. (2001). Wavelet entropy in event-related potentials: a new method shows ordering of EEG oscillations. Biol. Cyber. 84, 291–299. doi: 10.1007/s004220000212

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Rankine, L., Stevenson, N., Mesbah, M., and Boashash, B. (2007). A nonstationary model of newborn EEG. IEEE Trans. Biomed. Eng. 54, 19–28. doi: 10.1109/TBME.2006.886667

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Richman, J., and Moorman, J. (2000). Physiological time-series analysis using approximate entropy and sample entropy. Am. J. Physiol. Heart Circ. Physiol. 278, H2039–H2049.

PubMed Abstract | Full Text | Google Scholar

Rodrigues, S., Goncalves, J., and Terry, J. (2007). Existence and stability of limit cycles in a macroscopic neuronal population model. Physica D 233, 39–65. doi: 10.1016/j.physd.2007.06.010

CrossRef Full Text | Google Scholar

Shannon, C. (1948). A mathematical theory of communication. Bell Syst. Tech. J. 27, 623–656.

Google Scholar

Shin, H., Tong, S., Yamashita, S., Jia, X., Geocadin, R., and Thakor, N. (2006). Quantitative EEG and effect of hypothermia on brain recovery after cardiac arrest. IEEE Trans. Biomed. Eng. 53, 1016–1023. doi: 10.1109/TBME.2006.873394

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Slobounov, S., Cao, C., and Sebastianelli, W. (2009). Differential effect of first versus second concussive episodes on wavelet information quality of EEG. Clin. Neurophysiol. 120, 862–867. doi: 10.1016/j.clinph.2009.03.009

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Sotero, R. C., Trugillo-Barreto, N. J., Itturia-Medina, Y., Carbonell, F., and Jimenez, J. C. (2007). Realistically coupled neural mass models can generate eeg rythms. Neural Comput. 19, 478–512. doi: 10.1162/neco.2007.19.2.478

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Stam, C. J. (2005). Nonlinear dynamical analysis of EEG and MEG: review of an emerging field. Clin. Neurophysiol. 116, 2266–2301. doi: 10.1016/j.clinph.2005.06.011

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Sun, S., Lan, M., and Lu, Y. (2008). “Adaptive EEG signal classification using stochastic approximation methods,” in International Conference Acoustics, Speech and Signal Process (Las Vegas, NV), 413–416.

Ugray, Z., Lasdon, L., Plummer, J., Glover, F., Kelly, J., and Marti, R. (2007). Scatter search and local nlp solvers: a multistart framework for global optimization. Inf. J. Comput. 19, 328–340. doi: 10.1287/ijoc.1060.0175

CrossRef Full Text | Google Scholar

Valdes, P., Jimenez, J., Riera, J., Biscay, R., and Ozaki, T. (1999). Nonlinear EEG analysis based on a neural mass model. Biol. Cyber. 81, 415–424.

PubMed Abstract | Full Text | Google Scholar

Wang, J., Wang, B., Wang, X., and Nakamura, M. (2011). “MPA EEG model–based vigilance level estimation by artificial neural network,” in International Conference on Biomedical Engineering and Informatics (Shanghai), 795–799.

Google Scholar

Xu, G., Zhang, X., Yu, H., Ho, S., Yang, Q., Fu, W., et al. (2010). Complexity analysis of EEG under magnetic stimulation at acupoints. IEEE Trans. Appl. Superconductivity 20, 1029–1032. doi: 10.1109/TASC.2010.2040726

CrossRef Full Text | Google Scholar