Sugihara causality analysis of scalp EEG for detection of early Alzheimer's disease.

Recently, Sugihara proposed an innovative causality concept, which, in contrast to statistical predictability in Granger sense, characterizes underlying deterministic causation of the system. This work exploits Sugihara causality analysis to develop novel EEG biomarkers for discriminating normal aging from mild cognitive impairment (MCI) and early Alzheimer's disease (AD). The hypothesis of this work is that scalp EEG based causality measurements have different distributions for different cognitive groups and hence the causality measurements can be used to distinguish between NC, MCI, and AD participants. The current results are based on 30-channel resting EEG records from 48 age-matched participants (mean age 75.7 years) - 15 normal controls (NCs), 16 MCI, and 17 early-stage AD. First, a reconstruction model is developed for each EEG channel, which predicts the signal in the current channel using data of the other 29 channels. The reconstruction model of the target channel is trained using NC, MCI, or AD records to generate an NC-, MCI-, or AD-specific model, respectively. To avoid over fitting, the training is based on the leave-one-out principle. Sugihara causality between the channels is described by a quality score based on comparison between the reconstructed signal and the original signal. The quality scores are studied for their potential as biomarkers to distinguish between the different cognitive groups. First, the dimension of the quality scores is reduced to two principal components. Then, a three-way classification based on the principal components is conducted. Accuracies of 95.8%, 95.8%, and 97.9% are achieved for resting eyes open, counting eyes closed, and resting eyes closed protocols, respectively. This work presents a novel application of Sugihara causality analysis to capture characteristic changes in EEG activity due to cognitive deficits. The developed method has excellent potential as individualized biomarkers in the detection of pathophysiological changes in early-stage AD.

Introduction

Mild cognitive impairment (MCI) is the stage between cognitively normal and dementia. Patients with MCI are at high risk to develop Alzheimer's disease (AD) or other neurological conditions (Petersen, 2003; Petersen et al., 2006). Previous research has also shown that MCI patients progress to AD at a rate of approximately 10–15% of patients per year (Petersen et al., 2001; Wee et al., 2011). Currently there is no effective treatment for older patients with Alzheimer's disease (AD) to restore lost memory ability. Therefore, early detection and interventions to slow the disease progression are vitally important.

Diagnosis of MCI and AD of this cohort was based on clinical assessments, neuropsychological tests and patient history evaluations (Schmitt et al., 2012; Albert et al., 2011; McKhaon et al., 2011). Once abnormal cognitive decline is recognized, physicians employ more quantitative diagnostic tools including cerebrospinal fluid protein analysis, positron emission tomography (PET) scans, and magnetic resonance imaging (MRI) (Anoop et al., 2010; Bartzokis et al., 2004; Ikonomovic et al., 2008; Jack et al., 2008; Schmitt et al., 2012). Unfortunately, the expense and lack of access to these diagnostic methods in a primary care setting often deter physicians from ordering them routinely. More recent research effort has explored the development of more convenient and noninvasive means for screening for MCI/AD, including the use of scalp electroencephalography (EEG) (Dauwels et al., 2010; McBride et al., 2013; McBride et al., 2014), to predict cognitive impairment in older adults. Recent evidence has shown that cognitive event-related potentials of averaged EEG signals can predict AD pathology years before of clinical diagnosis (Olichney et al., 2013).

The effects of cortical neuronal deaths, axonal pathology, cholinergic deficits, and other neural network disconnections as concomitants of MCI/AD are manifested by multiple changes in EEG characteristics (Brenner et al., 1986; Jeong, 2004; Signorino et al., 1995; Soininen et al., 1989). Scalp EEG signals are synchronized oscillations of brain electric potentials from post-synaptic neural activity. Consistent with loss of synaptic contacts in both the neocortex and hippocampus in early AD neuropathology, the UK-ADC longitudinal cohort reported a loss of afferents from the medial temporal cortex as underlying synapse loss seen in MCI and early AD (Scheff et al., 2006). The medial temporal cortex has connections to the inferior temporal cortex where object memory is processed, and the prefrontal cortex that is associated with many cognitive functions. EEG data from MCI and AD patients have been shown to have lower mean levels of channel-to-channel synchronization than those of healthy controls (Koenig et al., 2005; Stam et al., 2005). Decreased synchronization has also been observed in graphical network models of EEG activity with AD (Stam et al., 2003; Stam et al., 2005; Stam et al., 2007; McBride et al., 2013). Analysis of inter-channel coherence in MCI/AD has revealed greater uniformity in alpha (7.5–12.5 Hz) and beta (12.5–25 Hz) band activity in MCI and greater uniformity in alpha and gamma (25–40 Hz) band activity in AD during resting states (McBride et al., 2014).

Granger causality, based on statistical hypothesis tests, has been used by various authors to assess causal connectivity in scalp EEG. For example, Barrett et al. (2012) conducted Granger causality analysis of steady-state EEG during propofol-induced anesthesia. Their work indicates significant increases in bidirectional Granger causality during loss-of-consciousness, especially in the beta and gamma frequency ranges. Bressler and Seth (2011) presented a comprehensive introduction to Granger causality for its merits and limitations for applications in neural time series data.

Recently, Sugihara and collaborators proposed a new causality concept, which, in contrast to Granger causality, characterizes underlying deterministic causation of the system (Deyle and Sugihara, 2011; Sugihara et al., 2012). Sugihara causality is based on the theory of nonlinear state space reconstruction and indicates that nearby points on the manifold of one variable can be used to estimate nearby points on the manifold of another variable. Nonlinear state space reconstruction techniques have been exploited by various authors for signal reconstruction (McBride et al., 2011), noise reduction (Richter et al., 1998), and stability prediction (Petrie and Zhao, 2012). In this article, we investigate the hypothesis that causality measurements of scalp EEG in Sugihara sense have different distributions for NC, MCI, and AD subjects and explore the potential of inter-channel causation as biomarkers for MCI and early AD.

Methods

Participants

The EEG data used in this study were collected at the Aging Brain and Cognition laboratory in the Behavioral Science Department and Sanders-Brown Center on Aging at the University of Kentucky (UK) College of Medicine. Scalp EEG was recorded within a week to 6 months of diagnosis. Participants between the ages of 60 and 90 years were recruited from a study cohort of cognitively normal older adults identified by the Alzheimer's Disease Center (ADC) of the UK College of Medicine. The normal older participants are screened annually and were followed into likely MCI or AD stages until autopsy. When screenings indicated possible cognitive decline they were referred to the ADC's Research Memory Disorders Clinic. MCI and AD participants were diagnosed and recruited by cognitive neurologists Drs. C. Smith and G. Jicha at the UK ADC Clinical Core and from its Research Memory Disorders Clinic (Schmitt et al., 2012; McBride et al., 2013). All MCI participants belonged to the amnestic MCI subtype. In addition, a few of the MCI participants also presented with executive dysfunction. Differences between single- and multiple-domain MCI subtypes are not the focus of the current study. MCI and early AD participants' EEG data were recorded as soon as possible after diagnoses were made. IRB approval and consent were obtained with appropriate measures in place to ensure protection of those with impaired cognition. Table 1 shows the demographics of the participants.

Table 1

Demographics of the participants.

	NC	MCI	AD
Number	15	16	17
Age (std)	75.7 (5.5)	74.6 (8.8)	76.7 (5.2)
Gender	Male = 6, Female = 9	Male = 12, Female = 4	Male = 7, Female = 10
MMSEMedian (range)	n = 1430 (28, 30)	n = 1627.8 (23, 30)	n = 15,24.5 (19, 28)
Logic memory IMedian (range)	n = 1413.5 (9, 18)	n = 169.5 (2, 17.5)	n = 146(2, 14)

Note: std = standard deviation. MMSE = mini-mental state examination. Logic Memory I is from Wechsler Memory Scale.

Participants were screened to exclude active or unstable medical conditions, depression, and other psychiatric disorders, or history of neurological or neurosurgical conditions. No participants were to have any psychoactive medication other than antidepressants. Systolic hypertension, donepezil, and anti-depressants are known to modulate EEG markers including event-related potentials (ERPs) associated with cognitive functions. These conditions are very common in the population of interest, complicating outright exclusion; however, uncontrolled hypertension or use of sedatives, such as the benzodiazepines, was basis for exclusion from the study.

Participants were well matched with regard to age, with normal controls (NC), MCI, and AD participants having mean ages of 75.7 years (SD 5.5 years), 74.6 years (SD 9.0 years), and 76.7 years (SD 5.2 years), respectively. NC and AD participants were also well matched in regard to gender, with NC and AD participants being comprised of 60% females. MCI participants were only 25% female. Difference in MCI gender was likely due to recruiting and does not reflect population trends.

Data

Participants were connected to 64- or 32-channel EEG caps using a Neuroscan II system (10–20 montage). In either case, only the 32 common channels were included in the analysis. EEG data were recorded under three different protocols: (1) resting with eyes open for 5 min, (2) resting with eyes closed while counting backwards by ones for 10 min while tapping a finger, and (3) resting with eyes closed without tapping for 10 min. EEG recordings were performed without interruption at the same appointment for each participant. EEG data were acquired at 500 Hz. The 32 EEG channels included 2 ocular channels that were used to determine the dominant eye blink frequency. A simple 2nd order Butterworth filter was used to attenuate frequencies greater than 200 Hz. Notch filters were used to remove dominant eye blink frequencies and to remove 60 Hz frequencies, which may have been amplified by background electronic devices. In addition, a localized subtraction of baseline for identified muscular and other low frequency artifacts was performed. Data were then downsampled to 125 Hz. More advanced methods for noise reduction can be found in Hoffmann and Falkenstein (2008); Muthukumaraswamy, (2013); and Zhou and Gotmanb (2009).

Analysis

Data analysis consists of three key steps, as shown in the following pseudo code.

Feature generation using causality estimation

For each record j

Assembling all other records into an NC, MCI, or AD training set according to their cognitive states.

For each channel i

, q, and q, respectively.

End

Rearrange the calculated quality scores for record j in the order: q1,…, q30, q1,…, q30, q1,…, q30.

End

Feature reduction using regional average and principal component analysis

Classification using the principal components as features

In Step 1, features based on causality estimation are calculated using the leave-one-out principle (Seymour, 1993; Kohavi, 1995). All records except one target record are assigned into an NC, MCI, or AD training set according their cognitive states. In the Model Creation block, a reconstruction model is developed for each EEG channel, which predicts the signal in the current channel using data of the other 29 channels. The reconstruction model of the target channel is trained using the three training sets to generate an NC-, MCI-, or AD-specific model, respectively. Fig. 1 shows a schematic creation of NC models. MCI and AD models are created in a similar manner. Note that M(NC,i) represents the model for channel i trained using the NC training set. Likewise, M(MCI,i) and M(AD,i) represent the models for channel i trained using the MCI, and AD training sets, respectively. Artificial neural nets (Haykin, 1999) are adopted as the reconstruction models. All neural network simulations in this work are carried out using the neural network toolbox in Matlab™ (Mathworks, 2013). Each network model consists of a single hidden layer of five neurons and a single output node. Neurons of the hidden layer use tansig transfer functions and the output node uses a pureline transfer function. For each network, concurrent data from the remaining 29 channels are used to generate a reconstruction for the given channel. Each model is trained using 2-minute data samples.

Fig. 1

A schematic of the creation of NC models. MCI and AD models are created in a similar manner. Note that M(NC,i) represents the model for channel i trained using the NC training set. Likewise, M(MCI,i) and M(AD,i) represent the models for channel i trained using the MCI, and AD training sets, respectively. Since an EEG record has 32 channels (except 2 ocular channels), there are 30 NC models, 30 MCI models, and 30 AD models.

In the Feature Generation block, the reconstruction models are applied to the target channel i to estimate the causality between channel i and other channels. Fig. 2 shows schematic causality estimation for channel i using an NC model. Causality estimation for MCI and AD models is similarly conducted. The predicted channel is compared to the target channel to compute a quality score. Denote the predicted channel from the reconstruction model by and denote the target channel by y. The quality of the reconstruction is determined by Eq. (1):where σ represents the standard deviation of the data in y and MSE represents the mean squared error (i.e., the average of the squared differences between and y). In other words, the quality score, q, is the larger value between the R-squared value and 0. This quality score is a characterization of the causality of the target channel i (to other channels) (Sugihara et al., 2012). Denote the quality scores for channel i using the NC, MCI, and AD models by , , and , respectively. Since there are 30 channels (2 ocular channels are not included in these calculations) and 3 types of reconstruction models (NC, MCI, and AD), it leads to 90 raw quality scores for each target record. The raw quality scores are rearranged in the following order for later convenience: ,…, , ,…, , ,…, , where the superscript indicates the generating model and the subscript indicates the channel index.

Fig. 2

Schematic of causality estimation for channel i using an NC model. Causality estimation using MCI and AD models is similarly conducted. The leave-one-out principle is adopted to prevent information from leaking into the trained models. Causality of channel i is characterized using a quality score q defined in Eq. (1).

In Step 2, we reduce the raw features to a manageable number using regional average followed by principal component analysis. The 30 channels are grouped into 6 scalp regions based on their arrangement and location on the scalp. The regions include: (1) central (C); (2) frontal (F); (3) left temporal (L); (4) occipital (O); (5) parietal (P); and (6) right frontal (R); see Fig. 3 for regional boundaries. The raw scores are reduced by averaging the scores within the same scalp regions. This results in 18 regional average scores: , where the superscript indicates the generating model and the subscript indicates the scalp regions.

Fig. 3

Regional boundaries. C = central, F = frontal, L = left temporal, O = occipital, P = parietal, and R = right temporal.

Then, the regional average scores are further reduced using principal component analysis. For each protocol condition, a 48 × 18 matrix R of regional average quality scores for all participants is arranged as shown in Eq. (2):where qk is a row vector consisting of 18 regional average quality scores for participant k (k = 1,…,48). The mean is subtracted from the matrix R and the resulting matrix θ is then decomposed via singular value decomposition:where U and V are unitary matrices such that UU = I and VV = I and S is a diagonal matrix, the diagonal of which is comprised of singular values s which indicate the relative variance () in the corresponding principal components p, j = 1,…,18. The principal components p of the regional average quality score matrix R are then defined by Eq. (4),where vj is the jth column of V. Only the dominant principal components are adopted as features in the following classification analyses.

In Step 3, a three-way classification based on support vector machine (SVM) is conducted to predict the cognitive group for each participant. SVM models are binary classifiers that find the hyperplane in the feature space maximizing the minimum (signed) distance between the hyperplane and training points (Cortes and Vapnik, 1995). Consider feature vectors from two classes, denoted by −1 and +1, respectively. Using a quadratic kernel function, the distance from the hyperplane for a given set of features x is given by Eq. (5):where αk is the weight, sk is the support vector (selected subset of vectors of training data features), b is the bias, and N is the number of support vectors. If d < 0, then the classifier classifies x as class −1, otherwise as class +1.

For ease of description, we denote the three classes NC, MCI, and AD by classes 1, 2, and 3, respectively. First, we construct three binary SVM classifiers (i.e., MCI vs. NC, AD vs. NC, and MCI vs. AD) using leave-one-out cross validation (Seymour, 1993; Kohavi, 1995). For a given record, if at least two out of the three classifiers predict this record as belonging to class i, then the final decision of the three-way classifier is to classify the record as belonging to class i. Otherwise, the probability that a record belongs to each class, P, i = 1, 2, 3, is estimated using the conditional probabilities of the binary classifiers via pairwise coupling as proposed by Hastie and Tibshirani (Hastie and Tibshirani, 1998; Hastie et al., 2001). We denote by c = Prob(i|i or j) the conditional probability that a given record will be classified as class i by the ij classifier. Then, c can be estimated as follows. Define z = d − b and let m and m be the means of z for all participants in classes i and j, respectively. Furthermore, define θ as the standard deviation of z − (m + m)/2. The conditional probability c for a given record can then be defined by Eq. (6):

where φ(μ,σ) denotes the Gaussian density with mean μ and standard deviation σ. In order to determine probabilities P, we define a model for the true conditional probability as

We then use an iterative procedure to estimate the probabilities p: (1) Start with an initial guess for each class assuming equal probability, with ; (2) compute the corresponding ; (3) update using Eq. (8):where Nij is the number of training examples for binary classifier ij; (4) renormalize such that ; and (5) repeat steps (2)–(4) until convergence is achieved.

The final decision of the three-way classifier is then to choose the class corresponding to the largest probability, argmax(P).

Results

We explore scalp EEG collected under three conditions: resting eyes open, counting eyes closed, and resting eyes closed. Fig. 4 shows the relative variances of the principal components for all three conditions. The first two principal components account for more than 37% of the total variances for all conditions. The variances of the remaining components are significantly lower than those of the first two principal components. Based on such observations, we keep only the first two principal components in the following classification analyses. Fig. 5 shows contributions to the first two principal components from each regional average. Note that the first principal components (for all EEG conditions) contain contributions mainly from the MCI and AD models, each with opposite signs. This pattern suggests that the first principal component reflects the differences between AD-based regional averages and MCI-based regional averages. For the second principal component, the signs for NC models are opposite to those for MCI and AD models, suggesting that the second principal component reflects the differences between NC-based regional averages and MCI- or AD-based regional averages.

Fig. 4

Energies of principal components of regional averages. Note that for all three protocol conditions the first two principal components account for approximately 37% of the variance. Blue circles correspond to NC participants, green diamonds to MCI participants, and red stars to AD participants.

Fig. 5

Contributions to the principal components (PC1: circles and PC2: diamonds) from each regional average. See Fig. 3 for layouts of the regions. The superscript NC indicates causality scores computed using the NC model, the superscript MCI indicates the MCI model, and the superscript AD indicates the AD model. Note that PC1 contains contributions mainly from the MCI and AD models each with opposite signs. In PC2, the signs for AD models are opposite to those for MCI and AD models.

Plots of the first principal component versus the second principal component are shown in Figs. 6–8 for the resting eyes open, counting eyes closed, and resting eyes closed protocol conditions, respectively. These plots clearly show that the NC, MCI, and AD groups are separable using only the first two principal components under all EEG protocol conditions.

Three-way classification results for the resting eyes open condition are presented in Table 2. Of those records classified as NC, 93.3% were NC; of those records classified as MCI, 93.8% were MCI; and of those classified as AD, 100% were AD. Of all the normal participants, 93.3% of them were correctly identified; of all the MCI participants, 93.8% were correctly identified; and of all the AD participants, 100% were correctly identified. Overall, the prediction accuracy was 95.8%.

Table 2

Confusion table of 3-way SVM discrimination for resting eyes open.

		Predicted classes
		NC	MCI	AD
True classes	NC	14	1	0	93.3%
MCI	1	15	0	93.8%
AD	0	0	17	100%
		93.3%	93.8%	100%	Overall Acc.:95.8%

Three-way classification results for the counting eyes closed condition are presented in Table 3. Of those records classified as normal, 93.3% were normal; of those classified as MCI, 100% were MCI; and of those classified as AD, 94.4% were AD. Moreover, 93.3% of normal participants, 93.8% of MCI participants, and 100% of AD participants were correctly identified. Overall, the prediction accuracy was 95.8%.

Table 3

Confusion table of 3-way SVM discrimination for counting eyes closed.

		Predicted classes
		NC	MCI	AD
True classes	NC	14	0	1	93.3%
MCI	1	15	0	93.8%
AD	0	0	17	100%
		93.3%	100%	94.4%	Overall Acc.:95.8%

Three-way classification results for the resting eyes closed condition are presented in Table 4. Of those records classified as normal, 100% were normal; of those classified as MCI, 100% were MCI; and of those classified as AD, 94.4% were AD. Moreover, 93.3% of normal participants were classified as normal, and 100% of MCI and AD participants were classified as MCI or AD correctly. The overall accuracy was 97.9%.

Table 4

Confusion table of 3-way SVM discrimination for resting eyes closed.

		Predicted classes
		NC	MCI	AD
True classes	NC	14	0	1	93.3%
MCI	0	16	0	100%
AD	0	0	17	100%
		100%	100%	94.1%	Overall Acc.:97.9%

The resting with eyes open and eyes closed with counting tasks' three-way classifiers performed equally well in overall accuracy while the resting with eyes closed classifier performed the best. Misclassified NC participants were classified as MCI or AD, misclassified MCI participants were classified as NC, and no AD participants were misclassified.

Conclusions and discussion

We have developed and illustrated a novel method for assessing causal connectivity from scalp EEG in Sugihara sense. In contrast to Granger causality, Sugihara causality is based on the theory of nonlinear state space reconstruction to study the causation that characterizes deterministic dynamics of the system (Sugihara et al., 2012). Numerical results show great potential of Sugihara causality of EEG for capturing alterations in EEG patterns associated with changes in neurological functional organization due to MCI and AD. Our analyses demonstrate excellent accuracies in distinguishing the different cognitive groups using EEG records collected under three different protocol conditions. Results from principal component analysis show that participants from different cognitive groups are clearly clustered into different groups. The method based on Sugihara causality analysis captures characteristic changes in EEG activity due to cognitive deficits, and demonstrate excellent potential as biomarkers in the detection of pathophysiological changes in early-stage AD. We envision that risk scores can be developed using Sugihara causality analysis for cognitive diagnosis at the individual patient level. To the authors' best knowledge, this work is the first application of Sugihara causality to the analysis of scalp EEG.

While information for the EEG comes from synchronous sources, each electrode contains different information and all the electrodes together construct a state space. Causality in this work is defined in Sugihara sense based on nonlinear state space reconstruction; that is, a project of the state space (pseudo-state space) can be constructed using a selected number of the electrodes so that the other electrodes can be estimated from the pseudo-state space. The high discrimination accuracy achieved here is due to the fact that, on average, models performed more poorly in reconstructing EEG records of participants belonging to groups other than the generating model group. For example, the NC-based model produced significantly more accurate reconstructions, on average, for NC participants than for MCI or early AD participants; see Fig. 9. While the reason for the poorer quality of reconstructions for MCI and AD participants is unclear, it may be hypothesized that the causality relationships between EEG channels for normal individuals are less applicable to MCI and AD participants. That is, EEG nonlinear inter-channel reconstruction patterns appear to exhibit greater within-group similarities than between-group similarities. The differences between groups' relationships between EEG channels may reflect differences in the functional organization of information processing in the brain. This work presents for the first time the potential of Sugihara causality analysis as a diagnostic tool for detecting dementia.

Fig. 9

Regional average quality scores. Blue circles correspond to average of NC participants, green diamonds correspond to average of MCI participants, and red stars correspond to average of AD participants. Each row of plots presents the regional average quality scores for one of the three protocol conditions. Each column of plots corresponds to a different group used to generate the reconstruction models. For example, the first column of plots presents regional averages obtained using a reconstruction model trained on data from normal participants. Note that, on average, the qualities of reconstructions for participants other than those belonging to the group used to train the model are significantly lower. Recall that models were trained in a leave-one-out manner.

Signs and symptoms of dementia are caused by the loss of cortical neurons and cholinergic deficits (among other concomitants) which subsequently cause a loss of local and global neuronal connectivity in the brain (Brenner et al., 1986; Jeong, 2004; Signorino, 1995; Soininen et al., 1989). The loss of neuronal connectivity can lead to changes in the functional organization of the brain. One example of this phenomenon is compensatory mechanisms, in which one region of the brain takes over tasks that were previously processed in another region of the brain. Neuroimaging tools such as structural MRI can detect changes in the brain structure critical to cognitive ability long before clinical diagnosis of cognitive decline. Other biomarkers of the AD pathology cascade can signal early progression of the disease before MRI detection (Jack et al., 2013). EEG technology is sensitive to cellular changes in brain functions that are likely earlier biomarkers compared to structural MRI indicators. The results of this work suggest the possible use of Sugihara causality analysis of scalp EEG as a means for objectively discriminating cognitive states among normal, MCI, and AD participants and assessing an individual's risk profile for cognitive impairment. Interestingly, published work (Scheff et al., 2006) in our cohort has demonstrated that the number of synapses showed a significant correlation with the subject's mini-mental state score and with cognitive tests involving delayed recall. However, the synaptic loss showed no relationship to AD pathology indexed by Braak stage and apoE genotype. Thus, the current EEG method might be sensitive to general cognitive decline before measurable biomarker of AD amyloid deposits measured by PET scan and Cerebrospinal fluid in humans. Our results suggest that a simple discrimination model utilizing SVM and the features presented here may be a viable basis for future development of a diagnostic screening tool for MCI and early AD with applicability in the primary care setting. Such a rapid, simple, and cost-effective tool could also prove useful in the drug discovery process.

One limitation of this study is the small number of participants. Future work will increase the sample size to test the robustness and generality of the results here. Future research should also include other nondegenerative disease controls (e.g. vascular dementia or mixed-AD and other dementia) to test the specificities of the results. We shall validate the proposed method with known AD biomarkers to test the potential of the method as a diagnostic tool. An interesting future study is to directly compare the results of Sugihara causality to those of Granger causality. Future research should also consider the effect of gender differences and the presence of antidepressants.