Mild cognitive impairment (MCI) is an intermediate state between normal aging and dementia. Patients suffering from MCI have an increased risk of developing dementia, in particular Alzheimer’s disease (AD; Shah et al., 2000; Farias et al., 2005). Those MCI subjects with positive biomarkers for AD (i.e., amyloid deposition and neural injury markers such as accumulations of intracellular tau or medial temporal lobe atrophy [MTA]) are regarded to be at the symptomatic pre-dementia phase of AD, and are often referred to as “MCI due to AD” (Albert et al., 2011). The pathological processes that these biomarkers indicate are well described in AD (Braak and Braak, 1991) and are known to produce synaptic disruptions. AD has been described as a “disconnection syndrome,” not only at the cellular level, but also at the macroscale, since the connections in the brain networks also seem to be disrupted (Blennow et al., 1996; Selkoe, 2002; Delbeuck et al., 2003; Arendt, 2009; Takahashi et al., 2010). In fact, structural and functional changes have been described in MCI subjects, suggesting that this disconnection of brain networks already begins during the MCI-stage of AD (Pijnenburg et al., 2004; Koenig et al., 2005; Buldú et al., 2011; Wang et al., 2013).

One of the main concepts used to understand how the different brain areas interact is functional connectivity (Friston, 1994), which reflects the statistical interdependencies between two-time series of physiological activity. Several resting-state electroencephalography (EEG) and magnetoencephalography (MEG) studies have found decreases in functional connectivity, especially in higher frequency bands (i.e., alpha and beta bands), in MCI patients compared to healthy controls (Moretti et al., 2008; Gómez et al., 2009; López et al., 2014b; Cuesta et al., 2015). This change in synchronization pattern is quite similar to that described for AD patients (Stam and van Dijk, 2002; Jeong, 2004; Stam et al., 2006), although increased connectivity has also been described for AD involving posterior brain regions (Stam et al., 2006; Alonso et al., 2011). Reductions in functional connectivity in MCI have also been observed between regions of the default mode network (DMN), with a parallel disruption of the anatomical connections (Garcés et al., 2014; Pineda-Pardo et al., 2014). However, some resting-state studies comparing different MCI groups have also detected a specific hypersynchronization pattern in high frequency bands (alpha and beta) in those MCI subjects who finally developed AD (López et al., 2014a), or those who presented abnormal concentration of phospho-tau (p-tau) protein in the cerebrospinal fluid (CSF; Canuet et al., 2015). In a recent multicenter study, this profile of hypersynchronization was used to obtain a high percentage of correct classification of MCI and healthy controls (Maestú et al., 2015).

Based on the estimated functional connectivity between time series, a weighted network can be reconstructed using graph theory (Bullmore and Sporns, 2009). For this purpose, brain systems are described as sets of nodes (i.e., brain regions or sensors) and links (i.e., functional connections between nodes). The topology of these networks can then be characterized, for example, providing information about the local integration of the network (e.g., the “clustering coefficient”) as well as the global integration (e.g., the “path length”; see Rubinov and Sporns, 2010 for a review). A small-world network has a high local connectedness (quantified by a large clustering coefficient) and a high global integration (quantified by a short path length) and has been regarded as a network with an optimal topology for the transfer of information. AD patients exhibit brain networks that appear to have a sub-optimal topology in which the networks have shifted toward a more random configuration. This was mainly characterized by a loss of small-worldness (see Tijms et al., 2013 for a review), supporting the hypothesis of the disconnection syndrome. However, Tijms et al. (2013) also show that the results differ drastically between studies. Only few studies, using different approaches and modalities, have explored the network topology in MCI, reporting a disturbed balance between local and global integration [functional Magnetic Resonance imaging (fMRI); Wang et al., 2013; MEG; Buldú et al., 2011; Pineda-Pardo et al., 2014]. Methodological difficulties make the comparison between networks of different sizes and different edge densities challenging, if not impossible (Lee et al., 2006; van Wijk et al., 2010; Stam, 2014; Tewarie et al., 2015). This might lead to contradicting results that could be due to differences in modalities (Tijms et al., 2013), but also due to methodological biases (van Wijk et al., 2010). One solution is to reconstruct the minimum spanning tree (MST; Stam, 2014). The MST is a sub-graph of the complete network, which forms a backbone of the original graph. It is uniquely defined whilst avoiding the arbitrary choices of traditional approaches, therefore solving the limitations of previous graph studies (van Wijk et al., 2010; Tewarie et al., 2015). Despite its advantages and application in other neurological disorders (Ortega et al., 2008; Fraschini et al., 2014; Olde Dubbelink et al., 2014; Otte et al., 2015; Tewarie et al., 2015), there is only one fMRI study and one EEG study that has used the MST in comparing healthy elders and AD patients (Çiftçi, 2011; Engels et al., 2015), and none that studied the MST in patients with MCI.

For this reason, we performed an MEG study with MCIs and healthy controls with the aim to characterize how the functional network organization in the MCI stage differs from that of controls. To this end, we estimated MEG resting-state functional connectivity between cortical regions and characterized the topology of the reconstructed MST. We expected to find subtle differences in functional connectivity between MCI patients and controls while the graph theoretical measures will show a clear disrupted topological pattern in MCI. We expected that the novel MST measures give more insight in the network changes of MCI than the traditional network measures (normalized Cw and normalized Lw).

Magnetoencephalography recordings were obtained from 58 subjects (29 MCI patients and 29 healthy elderly subjects). The MCI group was recruited from the Hospital Universitario San Carlos (Madrid), and the control group from the Seniors Center of the District of Chamartin (Madrid). All subjects were right handed (Oldfield, 1971) and native Spanish speakers.

All participants were screened by means of standardized diagnostic instruments and also received an exhaustive neuropsychological assessment. To evaluate their global and cognitive functional status there were used: the Spanish version of the Mini-Mental State Examination (MMSE; Lobo et al., 1979), the Global Deterioration Scale (GDS; Reisberg et al., 1982), the Functional Assessment Questionnaire (FAQ; Pfeffer et al., 1982), the Geriatric Depression Scale-Short Form (GDS-SF; Yesavage et al., 1982), the Hachinski Ischemic Score (Rosen et al., 1980), the questionnaire for Instrumental Activities of Daily Living (Lawton and Brody, 1969), and the Functional Assessment Staging (FAST; Auer and Reisberg, 1997).

Additionally, all subjects underwent an extensive neuropsychological assessment to explore their cognitive functioning by using the following tests: direct and inverse digit span test (DDS and IDS, Wechsler Memory Scale III, WMS-III; Wechsler, 1997), immediate and delayed recall (IR and DR, WMS-III; Wechsler, 1997), phonemic and semantic fluency (PhF and SF, controlled oral word association test; Benton and Hamsher, 1989), ideomotor praxis of Barcelona test (IP; Peña-Casanova, 1990), Visual Object and Space Perception Test (VOSP; Warrington and James, 1991), Boston Naming Test (BNT; Kaplan et al., 1983), and Trail-Making Test (TMT), parts A and B (TMT-A and TMT-B; Reitan, 1958).

The MCI diagnosis was established according to the National Institute on Aging- Alzheimer Association (NIA-AA) criteria (Albert et al., 2011), which includes: (i) self- or informant-reported cognitive complaints; (ii) objective evidence of impairment in one or more cognitive domains; (iii) preserved independence in functional abilities; and (iv) not demented (McKhann et al., 2011). Besides meeting the clinical criteria, MCI participants had signs of neuronal injury (hippocampal volume measured by magnetic resonance imaging (MRI). So, they might be considered as “MCI due to AD” with an intermediate likelihood (Albert et al., 2011).

None of the participants had a history of psychiatric or neurological disorders (other than MCI). General inclusion criteria were: an age between 65 and 80, a modified Hachinski score ≤ 4, a short-form Geriatric Depression Scale score ≤ 5, and T1/T2-weighted MRI within 12 months and 2 weeks before MEG screening without indication of infection, infarction, or focal lesions (rated by two independent experienced radiologists; Bai et al., 2012). Patients were off those medications that could affect MEG activity, such as cholinesterase inhibitors, 48 h before recordings.

Methods were carried out in accordance with the approved guidelines. The study was approved by the Hospital Universitario San Carlos Ethics Committee (Madrid), and all participants signed a written informed consent prior to participation.

Magnetoencephalography signals were measured by a 306 channel Vectorview system (Elekta Neuromag Oy) at the Center for Biomedical Technology (Madrid, Spain), inside a magnetically shielded room (VacuumSchmelze GmbH, Hanau, Germany) within several days after the neuropsychological assessment. Brain magnetic fields were recorded in the morning during a task-free 3 min eyes-closed resting-state condition, while subjects sat comfortably. Participants were instructed to move as little as possible, and were monitored during the recording to ensure that they did not fall asleep.

Sampling frequency was 1000 Hz with an online filter with bandwidth 0.1–300 Hz. The position of the head inside the sensor array was determined using a head-position indicator (HPI) with four coils attached to the scalp (two on the mastoids and two on the forehead). These four coils along with the head shape (∼500 points) of each subject (referenced to three anatomical fiducials: nasion, and left-right preauricular points) were acquired using a three-dimensional Fastrak Polhemus system (manufacturer: Polhemus, Inc., USA). Vertical ocular movements were measured by two bipolar electrodes attached above and below the left eye, and a third one to the earlobe, for electrical grounding.

Maxfilter software (version 2.2, Elekta Neuromag Oy) was used to remove noise from the MEG data using the temporal extension of signal space separation (tSSS) with movement compensation (Taulu and Simola, 2006). Flat channels, or those that contained excessive artifacts, were manually discarded after visual inspection of the data by one of the authors (M. E. López) before estimation of the SSS coefficients. The tSSS filter was subsequently used to remove noise signals that SSS failed to discard, typically from noise sources near the head, using a subspace correlation limit of 0.9 (Medvedovsky et al., 2009) and a sliding window of 10 s.

3D T1 weighted anatomical brain MRI scans were collected with a General Electric 1.5T MRI scanner, using a high resolution antenna and a homogenization PURE filter [Fast Spoiled Gradient Echo (FSPGR) sequence with parameters: TR/TE/TI = 11.2/4.2/450 ms; flip angle 12°; 1 mm slice thickness, a 256 × 256 matrix and FOV 25 cm].

FreeSurfer software (version 5.1.0; Fischl et al., 2002) was used to obtain the hippocampal volumes, which were normalized with the overall intracranial volume (ICV) of each subject.

The outline of the scalp, as obtained from the subject’s structural MRI, was used for co-registration with the MEG data using the VUmc Amsterdam co-registration surface matching software, resulting in an estimated co-registration accuracy of approximately 4 mm (Whalen et al., 2008). A single sphere fitted to the scalp surface was used as a volume conductor model for the beamformer analysis. An atlas-based beamformer was used to project the MEG sensor signals to 78 cortical regions-of-interest (ROIs) from the Automatic Anatomical Labeling (AAL) atlas (see Supplementary Material) (Tzourio-Mazoyer et al., 2002; Gong et al., 2009). Based on the broad-band (0.5–48 Hz) beamformer weights, time series of neuronal activity were reconstructed for the voxel with the maximum power within a ROI for each frequency band separately, i.e., a virtual electrode that was representative for that specific ROI was reconstructed. A detailed description of this procedure is given in Hillebrand et al. (2012).

Per subject, five artifact free trials of approximately 16.384 s (four times 4096 samples) were selected after careful visual inspection, giving a total of 20 epochs of 4096 samples for further analysis. Time-series of neuronal activation were computed for the six frequency bands: delta (0.5–4 Hz), theta (4–8 Hz), lower alpha (8–10 Hz), upper alpha (10–13 Hz), beta (13–30 Hz), and gamma (30–48 Hz). Selected epochs were converted to ASCII-files and imported into an in-house developed software package BrainWave version 0.9.125, developed by one of the authors (C. J. Stam) and available at: http://home.kpn.nl/stam7883/brainwave.html.

Functional connectivity was assessed with the PLI, which quantifies the consistency of a phase relationship between two signals while zero-lag (mod π) phase differences are ignored (Stam et al., 2007b). Therefore, the PLI is insensitive to spurious interactions caused by the effects of volume conduction and/or field spread (Stam et al., 2007b; Hillebrand et al., 2012). The PLI ranges between 0 and 1 in which 0 represents no consistent coupling or coupling with zero-lag and one represents consistent phase-lagged coupling. First, the instantaneous phase for each time series is computed by taking the argument of the analytic signal (Stam et al., 2007b) as computed using the Hilbert transform. Second, we calculate the asymmetry of the distribution of instantaneous phase differences between two time series:

where the phase difference ΔΦ_t is defined in the interval [-π, π], <> denotes the mean value, sign stands for signum function, || indicates the absolute value, and t corresponds to time samples 1, …, Ns, where Ns is the number of samples. By calculating the PLI values between all pairs of ROIs, we obtained a 78 × 78 adjacency matrix, which we used for the network analyses (see below).

A low characteristic path length (L) and a high clustering coefficient (C) characterize a small-world network (Watts and Strogatz, 1998). In an unweighted network, the C represents the probability that two nodes are connected when they share a neighboring node and the L represents the average of the shortest distance between pairs of nodes, with distance defined by the number of links between nodes. From the weighted graph, the weighted clustering coefficient (Cw) and weighted characteristic path length (Lw) were calculated as described in Stam et al. (2007a). Fifty random control networks were created by randomly shuffling the PLI values in each adjacency matrix while keeping the matrix symmetry intact. For each ensemble of 50 random networks, the average Cw (random) and Lw (random) were computed. The observed network values were divided by the average values obtained for the random networks in order to create normalized values. The resulting normalized clustering coefficient (normalized Cw) and normalized path length (normalized Lw) were used for further analyses. These measures were computed for each epoch, and then averaged over the epochs for each subject.

We constructed the MST from the weighted adjacency matrix containing the PLI values. The MST is a unique subgraph that connects all nodes in the network by the strongest connections (defined as the network links with the highest PLI values) without forming cycles (Stam et al., 2014), and was reconstructed using Kruskal’s algorithm (Kruskal, 1956). By using 1/PLI as input to the algorithm, strongest connections are likely to be included in the MST, as long as no cycles are formed. The MST was characterized by the following measures: degree, betweenness centrality (BC), eccentricity, degree distribution (κ), the number of leafs, degree correlation (R), tree hierarchy and diameter. MST measures were computed for each epoch, and then averaged over the epochs for each subject. The degree describes how many links each node has. BC is a measure of the importance of a node within the network. The BC of node i is defined as the number of shortest paths in the network that run through a specific node, divided by the total number of shortest paths from any node to all other nodes in the MST. In our calculations, we used the maximum BC across all nodes as well as per node for further testing. The eccentricity of a node is defined as the longest distance between that node and any other node in the network. The degree distribution is formed by the likelihood (P) that a randomly chosen node of the network will have degree κ (the number of connections of a specific node); it is a plot of P(κ) as a function of κ (Stam and Reijneveld, 2007). The degree correlation is an index of how much the degree of a node is correlated to the degree of nodes it is connected to. The leaf number is the number of nodes that have a degree of 1 representing the “leafs” or “extremities” of the network. Tree hierarchy quantifies the trade-off between large scale integration in the MST and the overload of central nodes. The diameter represents the longest path in the MST. For more information about these measures, we refer to Boersma et al. (2013), Stam et al. (2014), Tewarie et al. (2015).

Subject’s characteristics were tested using independent samples t-tests or chi-square tests where appropriate using SPSS (20.0 for windows). For the PLI, permutation testing, based upon t-statistics, was used for each pair of regions [among the 78 studied, (78 × 77)/2 in total], and in each frequency band, with the aim of comparing both groups (Maris and Oostenveld, 2007). To this end, participants were randomly divided into two sets with the same size as the original groups (29 vs. 29 subjects). This procedure was repeated 2000 times (2000 permutations). A new t-test between each pair of regions [(78 × 77)/2 couples] in each frequency band, was then carried out using these two newly created groups, getting a t-value for each pair of regions and each frequency band. After sorting these 2001 t-test results [2000 corresponding to “randomly divided” groups and another one for the original (MCI vs. Control) subject’s distribution], only p-values within the 5% of lower values were considered statistically significant [note that this is for each pair of regions and each frequency band studied: hence it was repeated ((78 × 77)/2) × (number of frequency bands) times]. These analyses were uncorrected for the number of permutations performed and therefore serve as exploratory analyses (see Figures 1, 2 and Tables 2, 3). Afterward, we performed an FDR correction on the data with the goal to examine the statistical significant results that survive a multiple comparison correction. For the network analysis, we focused on those frequency bands in which the connectivity analyses showed (uncorrected) significant differences between the groups. Again, corrected and uncorrected (exploratory) permutation testing was used to compare the groups. All statistical analyses were performed using MATLAB (R2015b, Mathworks).

Subject characteristics are shown in Table 1. Controls and MCI patients did not differ in age, gender, or educational level. As expected, the scores of MMSE and two measures of episodic memory (immediate and delayed recall) were both lower in MCI patients compared to controls. Additionally, hippocampal volumes were both lower in the MCI group.

Significant differences between MCI patients and controls in PLI using the uncorrected permutation tests were obtained in the lower (8–10 Hz) and upper (10–13 Hz) alpha band (Figure 1). In the lower alpha band, compared to the control group, PLI in the MCI group was lower between left superior frontal gyrus -orbital part- and left gyrus rectus, and between left superior temporal gyrus and left insula (p < 0.01). Additionally, in the same frequency band, MCI subjects showed higher PLI values than controls between four regions, namely between right superior frontal gyrus (medial orbital) and right precentral gyrus, right superior frontal gyrus (dorsolateral) and left superior frontal gyrus (dorsolateral), left anterior cingulate and paracingulate gyri and left postcentral gyrus, and between right cuneus and right superior occipital gyrus.

Furthermore, in the upper alpha band, MCI subjects presented lower PLI values than controls between four regions, namely between right gyrus rectus and right olfactory cortex, left insula and left middle temporal gyrus, right parahippocampal gyrus and right fusiform gyrus, and right inferior temporal gyrus and left precuneus.

After FDR-correction for multiple comparisons, none of these significant differences survived.

We focused the network analysis on those frequency bands in which there were differences in connectivity, namely in lower and upper alpha bands.

Using the uncorrected permutation tests, no differences between controls and MCI subjects were found for normalized Cw and normalized Lw.

We found group differences in BC and degree in lower and upper alpha bands using the uncorrected permutation tests (see below).

Betweenness centrality results are shown in Figure 2 and Table 2. In the lower alpha band, there were no differences between controls and MCI patients in the maximum value of BC globally. However, compared to the healthy controls, the MCI group showed higher BC values in eight brain areas: the left pre- and post-central gyri, the left and right superior parietal gyri, the right middle frontal gyrus, the right superior and middle temporal gyrus (temporal pole) and the right insula; and lower values in one brain area: the left anterior cingulated/paracingulate gyri (see Supplementary Material).

In the upper alpha band, we did not find global differences between the groups, but there were differences in BC for specific brain areas. Compared to controls, MCIs exhibited higher values in three brain areas: the left superior frontal gyrus dorsolateral (lSFGdor), the right supramarginal gyrus (rSMG), and the right anterior cingulate/paracingulate gyri (rACG); and lower BC values in three brain areas: the left middle temporal gyrus (lMTG), the right fusiform gyrus (rFFG) and the right parahippocampal gyrus (rPHG).

Finally, the MCI group showed higher MST degree values than controls in the lower alpha band in five brain areas: the left superior frontal gyrus, dorsolateral, the left postcentral gyrus, the left Heschl gyrus, the right middle frontal gyrus and the right calcarine fissure and surrounding cortex; while they exhibited a lower degree value in one brain area: the right olfactory cortex. In the upper alpha band, MCI subjects exhibited higher degree values in four brain areas: the left superior frontal gyrus, dorsolateral, the left superior parietal gyrus, the rSMG and the right anterior cingulate and paracingulate gyri; and lower degree values in three brain areas: the lMTG, the right inferior frontal gyrus, orbital part and the rPHG. MST degree results are shown in Figure 2 and Table 3.

After FDR-correction for multiple comparisons, none of these significant differences survived. Also, we did not find any differences between the controls and the MCI group for any of the other MST measures in these two frequency bands.