Fluctuations of spontaneous EEG topographies predict disease state in relapsing-remitting multiple sclerosis.

Spontaneous fluctuations of neuronal activity in large-scale distributed networks are a hallmark of the resting brain. In relapsing-remitting multiple sclerosis (RRMS) several fMRI studies have suggested altered resting-state connectivity patterns. Topographical EEG analysis reveals much faster temporal fluctuations in the tens of milliseconds time range (termed "microstates"), which showed altered properties in a number of neuropsychiatric conditions. We investigated whether these microstates were altered in patients with RRMS, and if the microstates' temporal properties reflected a link to the patients' clinical features. We acquired 256-channel EEG in 53 patients (mean age 37.6 years, 45 females, mean disease duration 9.99 years, Expanded Disability Status Scale ≤ 4, mean 2.2) and 49 healthy controls (mean age 36.4 years, 33 females). We analyzed segments of a total of 5 min of EEG during resting wakefulness and determined for both groups the four predominant microstates using established clustering methods. We found significant differences in the temporal dynamics of two of the four microstates between healthy controls and patients with RRMS in terms of increased appearance and prolonged duration. Using stepwise multiple linear regression models with 8-fold cross-validation, we found evidence that these electrophysiological measures predicted a patient's total disease duration, annual relapse rate, disability score, as well as depression score, and cognitive fatigue measure. In RRMS patients, microstate analysis captured altered fluctuations of EEG topographies in the sub-second range. This measure of high temporal resolution provided potentially powerful markers of disease activity and neuropsychiatric co-morbidities in RRMS.

Introduction

Multiple sclerosis (MS) is characterized by recurrent inflammatory events and progressive neurodegeneration leading to disseminated demyelination and axonal loss in the central nervous system (Compston and Coles, 2008). The search for additional disease markers is needed (Filippi and Agosta, 2010) because conventional structural magnetic resonance imaging (MRI) of the brain and spinal cord (Rovira et al., 2015) does not capture the “hidden” axonal damage known to occur early in the normal-appearing brain tissue (Fu et al., 1998) and because there ispoor correlation between lesion load and clinical impairment (Barkhof, 2002).

MS has received much less attention from the EEG and MEG field than from the MRI community. Very few studies have looked at electrophysiological markers of the clinical course (Fuhr et al., 2001, Schlaeger et al., 2014), and only recently has interest grown in the functional connectivity and resting-state network analysis by EEG, with the aim to characterize altered neuronal conduction known to occur in MS pathology (Lenne et al., 2013, Leocani et al., 2000, Van Schependom et al., 2014). A few MEG studies investigated functional connectivity at rest in patients with MS and found altered connectivity patterns in specific frequency bands (Cover et al., 2006, Hardmeier et al., 2012, Schoonheim et al., 2013b, Tewarie et al., 2013, Van der Meer et al., 2013).

Alternative to frequency-band specific analyses with quantitative EEG or functional connectivity measures, brain activity can be described through scalp potential fields. During rest the broad-band EEG displays topographical fluctuations, among which meta-stable periods of some tens of milliseconds can be detected when a topographical clustering approach is used. These quasi-stable states have been termed “microstates” (Lehmann et al., 1987), and most studies have consistently obtained four prototypical topographies along which the fluctuations are described (Koenig et al., 2002; for reviews see Khanna et al., 2015, Lehmann et al., 2009). Both the spatial and temporal orders of these topographical fluctuations are not random, but follow a complex temporal structure, including long-range dependencies (Gschwind et al., 2015, Van De Ville et al., 2010). Simultaneous EEG/fMRI recordings in healthy subjects showed that the different microstates could be associated with the blood-oxygen-level dependent (BOLD) pattern of established resting-state networks (Britz et al., 2010, Yuan et al., 2012). This supported the idea that EEG microstates represent elementary short-lasting periods of coordinated synchronized communication within large-scale brain networks. Consequently, it was assumed that the EEG microstates could be the electrophysiological correlates of the periods of stable spatial patterns proposed in the global workspace theory (Baars, 2002, Baars, 2005, Changeux and Michel, 2004, Dehaene and Changeux, 2011).

Several studies over the last 20 years have shown that the temporal dynamics of EEG microstates are influenced by different states of consciousness such as hypnosis (Katayama et al., 2007), meditation (Kopal et al., 2014), and sleep (Brodbeck et al., 2012), and that they are altered in diseases such as schizophrenia (Kindler et al., 2011, Lehmann et al., 2014, Strelets et al., 2003), risk for schizophrenia (Andreou et al., 2014, Tomescu et al., 2014, Tomescu et al., 2015), and dementia (Dierks et al., 1997, Nishida et al., 2013, Strik et al., 1997) in a disease-specific way.

In the present study, we compared high-density EEG topographies (microstates) during rest between patients with RRMS and healthy subjects. Given that RRMS is characterized by relapsing episodes of distributed inflammation, leading to demyelination and axonal loss, we expected that a focal alteration of microstate topographies would not be observed, but instead a temporal alteration of the microstate switch, either as an increased speed due to compensatory or adaptive hyperactivation (decrease of microstate durations) or as decreased speed in the context of the slowdown of the neuronal transmission (increase of microstate durations). We hypothesized that such altered temporal dynamics could be related to a specific dysfunction of clinical or neuropsychological parameters (Schoonheim et al., 2013a, Schoonheim et al., 2015b). We asked whether such a temporal alteration would be found in all microstates classes as a general effect of the disease, or whether it would be specfic to certain microstate classes only, and therefore specific to certain brain networks related to these states. If the latter was the case, EEG microstate changes might serve as a possible surrogate marker of MS pathology.

Methods

Subjects

The study was approved by the local ethics committee, and all participants gave written informed consent for their participation in accordance with the Declaration of Helsinki. From a cohort of 86 patients with MS, in which 256-channel high-density EEG was recorded at the Outpatient Clinic of Neurology, University Hospital Basel, a total of 53 patients with confirmed relapsing-remitting multiple sclerosis (RRMS) according to McDonald's diagnostic criteria (Polman et al., 2011) were selected. The following criteria were used: (1) an Expanded Disability Status Scale (EDSS) ≤ 4 (Kurtzke, 1983); (2) right-handedness; (3) neither clinical relapse nor corticosteroid therapy for at least 6 weeks prior to inclusion in the study; and (4) no other neurological diagnoses or major psychiatric illnesess according to the DSM-IV-TR criteria. Matched to this patient sample, we also selected 49 right-handed healthy control subjects with no history of neurological or psychiatric condition. The demographic details are summarized in Table 1. There were more males in the control group, but the difference was not significant.

Table 1

Demographics (mean and standard deviation).

	Patients		Controls
N = 53		N = 49
Females/males	45/8		33/16
Age [y]	37.69	± 7.10	36.35	± 8.20
Education [y]	14.64	± 2.69	15.27	± 2.22
Disease duration [y]	9.99	± 6.09
Annualized relapse rate (2y-ARR)	0.58	± 0.68
Expanded Disability Status Scale (EDSS)	2.12	± 0.97
Kurtzke Functional System Score – visual (FSS-vis)	0.32	± 0.83
Kurtzke Functional System Score – mental (FSS-ment)	0.91	± 0.84

For detailed patient specification, see Supplementary Table 1.

For detailed patient specification, see Supplementary Table 1.

Neuropsychological assessment

Patients with RRMS and healthy controls who participated in the study underwent neurological examination and neuropsychological evaluation, including several established tests known to best describe the cognitive dysfunction in the context of MS (Papadopoulou et al., 2013): Laterality Index of Handedness (LIH) was assessed by the Edinburgh Handedness Inventory (Oldfield, 1971); right-handedness was defined by a LIH > 30. Visuoperceptual processing and psychomotor speed and working memory were tested with the Symbol Digit Modalities Test (SDMT; Benedict et al., 2008), giving the number of items processed in 90 s. Fatigue was evaluated using the Fatigue Scale for Motor and Cognitive Functions (FSMC), which was composed of both a cognitive sub-score (cog) and a motor sub-score (mot) (Penner et al., 2009); each sub-score summed the patient's rating on a 10-item scale. Depression was scored with the German version (Allgemeine Depressions Skala-L [ADS-L]) of the Center for Epidemiologic Studies Depression Scale (CES-D) (Hautzinger and Bailer, 1991), which resulted in a score from the patient's rating on a 20-item scale.

EEG data acquisition and preprocessing

During data acquisition, the participants sat calmly in a chair with their eyes closed, without falling asleep, and scalp EEG was recorded over 12 min with a high-density, 256-channel EEG system (Netstation 200, using a HydroCel Geodesic Sensor Net, Electrical Geodesics, Inc., Eugene, OR). The electrode net was placed relative to the preauricular points and Fz, Cz, and Oz as landmarks. Electrode impedances were kept below 40 kΩ in order to ensure good quality of data acquisition (Ferree et al., 2001). Recording band-pass was 0.1 to 100 Hz; the sampling frequency 1 kHz; and the vertex was used as the recording reference. A subset of electrodes was monitored online during recording in order to check for vigilance fluctuations and ensure that the participant did not fall asleep. Offline, the EEG was band-pass filtered between 1 and 40 Hz; electrodes on the cheeks and neck were excluded, resulting in 204 electrodes that were kept for further analysis. Artifact-free EEG epochs were then selected from each 12-min recording. Epoch length varied depending on the occurrence of visually identified artifacts, which included eye blinks, muscle activity and DC shifts, or any other intermittent high-amplitude deflections on any electrode. The total of artifact-free segments covered at least 5 min of resting EEG. Independent Component Analysis (ICA, infomax based) was applied to these epochs in order to remove cardiac and oculomotor artifacts based on the waveform, topography, and time course of the ICA component (Jung et al., 2000). The data were then down-sampled to 125 Hz, disturbed electrodes were interpolated using a 3-D spherical spline (Perrin et al., 1989), and the data were recomputed to the common average reference. EEG data processing was realized with the freely available Cartool software package (release 3.51 [2268]) (Brunet et al., 2011) and custom MATLAB scripts (release 2012a, Mathworks Inc., Natick, MA).

Topographical analysis

Microstate analysis consists of determining the topographical template maps with subsequent labelling of the whole EEG recording. In order to determine these most dominant topographical template maps, a standard approach applied in previous studies by different groups was used and consisted of running a spatial k-means clustering procedure on every subject separately (Brunet et al., 2011) using the time frames of peaks of Global Field Power (GFP) (Britz et al., 2010, Brodbeck et al., 2012, Brunet et al., 2011, Khanna et al., 2014, Koenig et al., 2002, Koenig et al., 2005, Michel et al., 2009, Nishida et al., 2013). The GFP is calculated as the standard deviation of all electrodes after subtracting the mean potential across the electrodes (Lehmann and Skrandies, 1980); therefore, it is a reference-free measure of total field strength. GFP maxima are considered the time points of maximal synchronized neuronal activity (Skrandies, 2007), and the maps tend to be topographically stable around the GFP maxima (Michel et al., 2009), the GFP maxima represent therefore the moments of best signal-to-noise ratio. The formula is as follows (n = number of channels of the montage, including the reference; ui = the average-referenced potential of the ith electrode at a given time point t):

For each participant, the spatial k-means clustering procedure identified between 2 and 10 different, most representative topographical maps, which were then entered into a second k-means clustering across all individuals (patients and controls). At both steps, the optimal cluster number was chosen based on the Krzanowski-Lai cluster index (Charrad et al., 2014, Krzanowski and Lai, 1988). The Krzanowski-Lai cluster index is determined by the L-corner of the dispersion curve, which is a quality measure of the clustering. The optimal clustering is set when an additional cluster does not lead to a significant gain of the global quality. The criterion had been proposed for EEG microstate clustering in Murray et al. (2008) and was explained in more detail therein. This procedure resulted in one set of four different “grand” template topographies that, at a minimum number, explained a maximum of the variance of the whole data set across patients and controls, as previously demonstrated (Britz et al., 2010, Khanna et al., 2014, Khanna et al., 2015, Koenig et al., 2002). (Note that only the relative topographical configuration and not its polarity was considered). We then fitted these four topographical template maps back to the entire preprocessed EEG recording of each participant, using the spatial correlation between the instantaneous EEG topography and each of the four grand template maps after normalizing both maps by dividing by the GFP. Every time point in the resulting four time-dependent (absolute-valued) correlation curves was then assigned to the template that showed the highest correlation value (winner-takes-all fashion), resulting in the fine-grained temporal sequence of four microstates labelled according to the convention as classes A, B, C, and D (Van De Ville et al., 2010) (Fig. 1).

Fig. 1

Method of microstate analysis. A. The patient's high-density EEG (256 channels) at rest (eyes closed), after standard preprocessing, is displayed as a time series of global field topographies, showing quasi-stable periods between irregular changes (B.). C. The peaks of global field power (GFP) were determined and their specific topographies were selected and submitted to a k-means clustering procedure, for each individual participant, and in a second step across all participants (D.). E. This resulted in a set of four most representative topographies for all participants in both groups equally, the four microstate classes. Only the topography's relative configuration but not its polarity is considered. These four template topographies are then fitted back to the original EEG recording of each individual participant resulting in a labelling of the whole recording according to predominating microstate class. G. The resulting fine-grained time sequence of the labels is called the microstate sequence and used for statistical analysis.

To determine the degree of topographical fluctuations and their relevance in RRMS, the following conventional parameters were calculated from the microstate sequences (Tomescu et al., 2014): (1) the total time coverage giving the percent of total time of the whole EEG recording, which was covered by each of the four microstates; (2) the global explained variance (GEV) for each microstate, calculated as the sum of the squared spatial correlation between the template and each assigned time point's topography and weighted by the GFP at each time point. The GEV therefore gives a ratio of how well each template described the assigned time points through the whole data (Brunet et al., 2011, Koenig and Giannotti, 2009). The GEV formula is as follows (from Murray et al. (2008); where GFPu(t) is the GFP of the data for microstate class U at time point t labelled; Tt is the template map assigned by the segmentation for microstate class U at time point t labelled; and Cu,Tt is the spatial correlation between data of microstate class U at time point t labelled and the template map of microstate class U):

Finally, we calculated (3) the microstate duration (Koenig et al., 2002), given as median, geometric mean (geomean), and geometric standard deviation (Limpert and Stahel, 2011) of each microstate's duration in milliseconds. As demonstrated elsewhere, temporal organization is a key property of EEG microstates (Gschwind et al., 2015).

Statistical analysis

The four microstate template maps were created for both patients and controls together, and both groups were directly compared to these common template maps. The microstate GEV, total time coverage, and duration measures (median, geomean, and geometric standard deviation) of each of the four microstate classes (A, B, C, and D) were entered into a multivariate analysis of variance (MANOVA) across both groups. In order to test how the patient with RRMS's clinical and neuropsychological characteristics related to the neurophysiologic parameters, stepwise multiple linear regression models were calculated (inclusion/exclusion probability levels for the stepwise procedure at < 0.05/>0.1). A factorial analysis was performed on the different clinical characteristics using standard procedures and varimax rotation. Statistical analyses were run on SPSS (version 20, IBM Corporation, Armonk, NY), and the 8-fold cross-validation (McLachlan et al., 2005) was performed using the package DAAG (Maindonald and Braun, 2010) in R version 3.2.3 (The R Foundation for Statistical Computing, 2015, cran.r-project.org). Significance was set globally to an alpha level at 5%.

Results

The results of the neuropsychological assessment of patients with RRMS and controls are displayed in Table 2. Patients differed significantly from controls (two-sample t-tests) in cognitive and motor fatigue (FSMC, patient's rating on a 10 + 10-item scale), processing speed (SDMT, raw score counting the number of items processed in 90 s), as well as depression scores (CES-D, raw score from patient's rating on a 20-item scale).

Table 2

Results of the neuropsychological assessment of patients with RRMS and healthy controls.

	Patients	Controls	t-Test
	Mean ± SD	Mean ± SD	p-Value
Laterality Index of Handedness	88.44 ± 15.70	88.94 ± 17.37	0.88
SDMT	55.88 ± 10.78	63.60 ± 11.49	0.0009
FSMC-cog	24.62 ± 11.42	15.56 ± 4.45	0.0000
FSMC-mot	25.79 ± 11.46	15.98 ± 4.23	0.0000
CES-D	9.84 ± 8.98	6.21 ± 5.51	0.02

SDMT: Symbol Digit Modalities Test; FSMC-cog: cognitive sub-score of the Fatigue Scale for Motor and Cognitive Functions; FSMC-mot: physical sub-score of the Fatigue Scale for Motor and Cognitive Functions; CES-D: Center for Epidemiologic Studies Depression Scale.

Thus, SDMT, FSMC-cog, FSMC-mot, and CES-D scores were looked at together with the disease-characterizing parameters in order to build stepwise multiple linear regression models between these clinical measures and the EEG microstate parameters within patients (see section below).

An explorative factorial analysis across all clinical parameters in patients showed mutual relationships between them. The first component grouped handicap scores and fatigue-related parameters, whereas time-related variables fell in the second component, which was composed of disease duration, age, and relapse rate. The depression scores, cognitive information processing speed, and sex were grouped together in the third component (Table 3).

Table 3

Factorial analysis of clinical and neuropsychological characteristics of patients with RRMS.

	Component
1	2	3	4
Variance explained (%)	31.955	14.536	10.587	9.366
Eigenvalues	3.835	1.744	1.27	1.124
	r	r	r	r
FSMC-cog	0.857
FSMC-mot	0.851
FSS-ment	0.800
EDSS	0.790
FSS-vis	0.608
Age		0.771
2y-ARR		− 0.766
Disease duration		0.697
SDMT			− 0.721
Sex			0.674
CES-D			0.592
Education years				0.684

N = 53. Factor contribution > 0.5 taken into account. FSMC-cog: cognitive sub-score of the Fatigue Scale for Motor and Cognitive Functions; FSMC-mot: physical sub-score of the Fatigue Scale for Motor and Cognitive Functions; FSS-ment: Kurtzke Functional System Score – mental; EDSS: Expanded Disability Status Scale; FSS-vis: Kurtzke Functional System Score – visual; 2y-ARR: Annualized relapse rate; SDMT: Symbol Digit Modalities Test; CES-D: Center for Epidemiologic Studies Depression Scale.

EEG microstate analysis

The two-step cluster analysis (first performed at an individual level and then across the whole population) resulted in four clusters as the optimal number of maps needed to explain the data according to the Krzanowski-Lai cluster index (k). We also performed a separate cluster anaylsis for each group in order to check whether four maps were the best solution for each group. The similarities of the topographies are striking as well as the similarities with the four canonical maps described in the literature (Koenig et al., 2002). Indeed, the first peak of the k-values was at four cluster maps in both groups and this was the highest value in the control group. In the patient group, a second peak with a slightly higher k-value was seen for seven clusters, but the additional clusters contributed very little to the total explained variance. Fig. 2 shows the four mean cluster topographies for the whole population, as well as separately for patients and controls.

Fig. 2

Grand template topographies as results of the two-step k-means clustering for all (N = 102), and separately for patients with RRMS (N = 53) and healthy controls (N = 49).

Given the high degree of similarity between the four cluster maps between the groups, we fitted the four template maps from the whole-population clustering (see Fig. 1E) to each subject and analyzed the temporal microstate dynamics according to commonly used, previously described parameters. Fig. 3 shows this step of the analysis.

Fig. 3

Upper part: Results of microstate analysis separated for each microstate class A, B, C and D for patients with RRMS (red) and healthy controls (blue). The group mean and S.E.M.is shown. Lower part: The histogram of the durations of all 4 microstates (controls above, patients below) details the infra-second time resolution of the microstate temporal sequence. It shows a long right tail, here arbitrarily cut at 400 ms (note, that some rare durations may reach > 2000 ms). It shows differences between maps and between patients and controls as measured by the MANOVA. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

All investigated parameters were submitted to a MANOVA, which revealed a group difference between patients with RRMS and controls (Pillai's trace = 0.32, F(19,82) = 2.07, p = 0.01, partial η2 = 0.32), a significant main effect between the different microstates (Pillai's trace = 1.00, F(19,82) = 11,612.56, p < 0.00000), as well as a significant interaction between group and microstate (F = 1.84, p = 0.04). Follow-up univariate ANOVAs (corrected for multiple comparisons using false discovery rate (Benjamini and Hochberg, 1995), q-value = 0.017) revealed that patients with RRMS had significantly increased GEV compared to controls for microstate classes A and B (GEV class A: F(1,100) = 6.05, p = 0.02; GEV class B: F(1,100) = 5.91, p = 0.02) and significantly increased total time coverage for class B (F(1,100) = 6.69, p = 0.01), while time coverage for class A did not survive multiple comparison (p = 0.047). GEV and total time coverage were not different between groups for classes C and D (p > 0.1 for each of them). Patients with RRMS showed increased measures of microstate duration compared to controls for classes A and B (median duration class A: F(1,100) = 9.18, p = 0.003; median duration class B: F(1,100) = 9.23, p = 0.003; geomean of class A: F(1,100) = 7.47, p = 0.007; geomean of class B: F(1,100) = 8.04, p = 0.006), but not for classes C and D (p > 0.43 for median and p > 0.45 for geomean). The geometric standard deviation did not show any difference between groups.

We further investigated the specific differences of microstate duration between patients and controls and plotted the probability density of durations, expressed as the subtraction of patients from controls (delta) for each microstate (Fig. 4). Again, differences were only observed for microstate classes A and B: There were fewer short duration states (< 120 ms) and more frequent long durations states (> 120 ms) in the RRMS group, indicating an imbalanced processing of microstate classes A and B. Distributions in classes C and D showed no particular differential pattern.

Fig. 4

The difference of probability of microstate duration between patients with RRMS and controls (delta) is shown for the four microstate classes. The delta is negative for durations < 120 ms for microstates classes A and B, indicating that patients display fewer short duration states, but it is positive for durations between 120 ms and 300 ms, indicating more longer durations in patients. This pattern is consistent with the patient's higher geomean and median in microstate class A and B. For classes C and D no such difference between patients and controls was found.

In order to control for basic frequency differences of the raw EEG recordings between groups, notably a diffuse slowing typically seen in patients with more advanced MS (Beach et al., 2011), and which theoretically could cause a temporal slowing also of the here described microstates fluctuations), we estimated the spectral power across the whole EEG recording in every patient with RRMS and healthy control subject using Welch's method (window length of 2500 time frames, 50% overlap, range 0.5 to 45 Hz in 0.2-Hz steps) averaged across all channels. None of the two-sample t-tests on each of the standard frequency bands were significantly different between the two groups (δ-band (1 to 4 Hz), p = 0.98; θ-band (4 to 8 Hz), p = 0.56; α-band (8 to 12 Hz), p = 0.56; β-band (12 to 20 Hz), p = 0.09; γ-band (20 to 40 Hz), p = 0.39).

Prediction of MS-related clinical measures by EEG microstate parameters

In order to determine the predictive value of the altered dynamics of the topographical field fluctuations in RRMS patients for their clinical and neuropsychological characteristics (Table 1), we fitted stepwise multiple linear regression models. RRMS-related clinical variables (see Section 2.1), such as the disease duration, the EDSS, and the 2-year annual relapse rate (2y-ARR), were entered as dependent variables, as well as those neuropsychological scores that were significantly different between patients with RRMS and controls (see Section 3.1), such as SDMT, FSMC-cog, FSMC-mot, and CES-D. All EEG microstate parameters that had survived the initial feature selection by the MANOVA (see Section 3.2) were entered as potential predictors in each of the models. Possible confounding factors such as age, sex, and education were equally included. The stepwise procedure selected the best predicting variables on a purely data-driven basis.

We then tested for multicollinearity in these models by means of the variance inflation factor (VIF), the most widely used diagnostic for multicollinearity, which estimates linear dependence between predictors and, therefore, how much the variance of a coefficient is “inflated”. None of the five models suffered from multicollinearity (disease duration: VIF 1.001; 2y-ARR: VIF 1.062; DES-D, FSMC-cog, and EDSS: VIF 1.000).

We obtained the following results: (1) The strongest model predicted disease duration (p < 0.0001) with a higher median duration of class A. Prediction was improved with age as a covariable, reaching almost 32% of explained variance. (2) High 2y-ARR was significantly predicted by shorter duration of class B and young age (p = 0.007). (3) High depression scores (CES-D) were predicted by short duration of class A (p = 0.009). (4) Cognitive fatigue (FSMC-cog, p = 0.02) and neurological impairment (EDSS, p = 0.03) were also significantly predicted by short duration of class B and high GEV of class A, respectively. For the SDMT and FSMC-mot models, no predictor variables were identified by the stepwise procedure.

Even when conducting a hierarchical stepwise multiple regression and forcing the covariates age and sex into each model, this result pattern did not change substantially: only for disease duration the forced model with age and sex was retained, and age played a significant role in disease duration and 2y-ARR, exactly as it also resulted from the pure stepwise analysis. In none of the models, sex played a significant role (Supplementary Table 2).

Validation of model consistency by 8-fold cross-validation

We performed a validation of these five predictor models (Table 4) using the consistency measure obtained by an 8-fold cross-validation, a standard procedure used in automatic classification approaches (McLachlan et al., 2005). For this purpose, the patient group (N = 53) was divided into eight groups of approximately the same size. The model was then calculated on the patients of seven groups, while the remaining eighth group was left out and used as a test group. This was done eight times with each of the eight non-overlapping groups once. The test error was then estimated by averaging the eight resulting mean square error estimates (its square root is the standard error of estimates,S.E.est). The results of the prediction for the remaining patient group was then compared to the full model calculated on all patients; the 8-fold cross-validation showed a very high model consistency with values of R2 from 0.891 to 0.980 for every one of the five different predictor models (Fig. 5).

Table 4

Stepwise linear regression models, on patient's clinical characteristics as target variables, reveal significant predictors from EEG microstate parametres.

Target	Model	Predictor	p-Value*)	R2	Adj. R2	β	sr	8-fold cross-validation
S.E.est	R2	Adj. R2
Disease duration	1	Median duration Class A	0.000	0.21	0.20	0.46⁎⁎⁎	0.46
	2	Median duration Class A	0.000	0.32	0.29	0.38⁎⁎	0.41	5.29	0.97	0.97
		Age				0.33⁎⁎	0.37
2y-ARR	1	Age	0.02	0.10	0.08	− 0.31⁎	− 0.31
	2	Age	0.007	0.18	0.15	− 0.32⁎	− 0.33	0.67	0.98	0.98
		Geomean duration Class B				− 0.29⁎	− 0.30
CES-D	1	Geomean duration Class A	0.009	0.13	0.11	− 0.36⁎	− 0.36	8.70	0.89	0.89
FSMC-cog	1	Median duration Class B	0.02	0.10	0.08	− 31⁎	− 0.31	11.09	0.95	0.95
EDSS	1	GEV Class A	0.03	0.09	0.07	0.30⁎	0.30	0.94	0.97	0.97

*): Corrected for multiple comparisons using false discovery rate (Benjamini and Hochberg, 1995), q-value = 0.03. Only data of patients were used N = 53. For the models SDMT and FSMC-mot, no predictor variables were selected by the stepwise procedure. 2y-ARR: 2-year annual relapse rate; CES-D: Center for Epidemiologic Studies Depression Scale FSMC-cog: cognitive subscore of the Fatigue Scale for Motor and Cognitive Functions; EDSS: Expanded Disability Status Scale; GEV: global explained variance; sr: partial correlation; S.E.est.: Standard error of estimate.

p < 0.05.

p < 0.01.

p < 0.001.

Fig. 5

Consistency measures for the established predictive models usign 8-fold cross-validation. Each figure compares the predicted model for all N = 53 patients to a model using the same predictor variables but calculated on 7 out of 8 folds only, each time leaving the remaining fold as unknown test set. This procedure is repeated 8 times on the 8 non overlapping folds, so that every patient was once a test case. All models (A.–E.) show a remarkably high coincidence with the predicting model for every of the 53 patients (R2 between 0.89 and 0.98). 2y-ARR: two-years annual relapse rate; EDSS: Expanded Disability Status Scale; CES-D: Center for Epidemiologic Studies Depression Scale; FSMC-cog: cognitive subscore of the Fatigue Scale for Motor and Cognitive Functions.

Discussion

In this study, we investigated the temporal fluctuations of EEG topographies measured with high-density EEG in 53 right-handed patients with RRMS and compared them to 49 matched healthy controls. To the best of our knowledge, this is the first study examining topographical fluctuations in patients with RRMS using the microstate analysis approach, which consists of segmenting the flow of field potential topographies into subsecond-periods of meta-stable canonical maps and determining changes of presence and duration of these segments for each of the maps. The key finding of this study is that in patients suffering from relapsing-remitting multiple sclerosis the temporal characteristics of two of the four canonical microstates were altered. In particular, patients showed prolonged durations of microstate classes A and B compared to age-matched controls, as well as longer total time coverage and higher GEV. These alterations were strongly predictive of patients' clinical parameters such as disease duration, annual relapse rate, and EDSS, as well as the neuropsychological scores of self-rated depression (CES-D) and cognitive fatigue (FSMC-cog).

Microstates as a temporal signature of electrophysiological processes

Our patients were thoroughly stratified into a homogenous patient group using a multitude of established clinical and neuropsychological measurements, which was therefore comparable to other studies with similar patient cohorts (Hasan et al., 2012, Schoonheim et al., 2015a). The EEG microstate analysis, in contrast to functional connectivity studies using fMRI resting state analysis, focuses on the temporal dynamics of large-scale brain networks. During a microstate period with stable topography, a set of sources in the brain operates synchronously (“The observation of microstates thus suggests the existence of subsecond global functional networks in the brain, established through a synchronized firing of the neural elements constituting the network” (Koenig et al., 2005)). Consequently, alterations in the temporal parameters of the EEG microstates due to disease indicates changes in the temporal dynamics of resting state networks in contrast to changes in connectivity between networks or network nodes as studied with fMRI (e.g., Pinter et al., 2016).

The question of the anatomical distribution of the networks underlying the EEG microstate has been the subject of several combined EEG/fMRI studies, which looked at BOLD fluctuations that go along with the (temporally smoothed) microstate fluctuations (Baker et al., 2014, Britz et al., 2010, Yuan et al., 2012). The different temporal dimensions of both methods are not contradictory and several models of temporal integration have been put forward (Baker et al., 2014, Gschwind et al., 2015, Van De Ville et al., 2010). It was suggested that microstate class A is closely related to BOLD activation changes of the bilateral superior temporal and parietal cortex, interpreted as sensorimotor (Yuan et al., 2012) or auditory network (Britz et al., 2010, Yuan et al., 2012). Microstate class B reflected the visual system via BOLD changes in the striate and extrastriate cortex (Britz et al., 2010, Yuan et al., 2012). Microstate class C recruited the bilateral insular and anterior cingulate “saliency network” (Britz et al., 2010) and microstate class D reflected the fronto-parietal attentional network of the nondominant hemisphere (Britz et al., 2010). The question of the relationship between the default-mode network (DMN) and the EEG microstates is still open. None of the combined EEG/fMRI studies were able to clearly identify the DMN as one of the EEG microstate networks, even if some nodes of the DMN contributed to microstate class C as in Britz et al. (2010). The only existing study on direct source imaging of the EEG microstates (Pascual-Marqui et al., 2014) suggested that different areas belonging to the DMN were contributing to each microstate, which indicated that the DMN can be split apart at high-time resolution. Further studies are needed to establish the relationship between EEG microstates and the DMN.

With respect to the putative relationship of the microstates to certain fMRI resting state networks, our modification of the temporal dynamics of microstate classes A and B could be interpreted as aberrant temporal functioning of sensory-related networks. This would be in agreement with structural studies that indicated the brain regions most affected by MS lesions to be in the corpus callosum, the periventricular white matter (especially around the posterior horns), the optic radiation, as well as the fiber tracts in the parietal and temporal lobes (Compston and Coles, 2008, Filippi et al., 2012, Giorgio et al., 2013, Hasan et al., 2012).

Impaired dynamics of microstate fluctuations predict disease state of MS patients with RRMS

The finding of less frequent short duration segments (< 120 ms) and more frequent long duration segments (> 120 ms) of microstate classes A and B in the patient group means an impaired dynamic of the switch from classes A and B to other classes. This observation of slowed topographical fluctuations points to a loss of flexibility of neuronal dynamics in the MS-affected brain, adding complementary information to the current knowledge of MS-related changes. Many studies have reported increased BOLD activity in task-related regions with additionally recruited adjacent regions (Pantano, 2002, Rocca et al., 2005). This was interpreted as a compensatory mechanism around existing structural brain lesions (Bonavita et al., 2011, Roosendaal et al., 2010). An alternative view, however, suggested that this apparent increased (static) connectivity might actually be explained by reduced neuronal diversity due to early axonal loss, rather than compensatory excitation (Hawellek et al., 2011). Our observation of loss of rapid changes in functional interactions would be in line with this interpretation, as well as with the recent finding of an increased low-frequency power in certain dynamic connectivity patterns in patients with RRMS (Leonardi et al., 2013).

The use of graph analytical methods (Bullmore and Sporns, 2009) has also lead to interesting knowledge about the modifications in MS-affected brains (Filippi et al., 2013). In particular, in EEG, MEG and fMRI, these connectivity measures showed a strong decrease in network efficiency, by segregation into increased network modularity (Gamboa et al., 2014) and decreased centrality (Hardmeier et al., 2012, Schoonheim et al., 2013a), with changes in clustering coefficient and path length (Schoonheim et al., 2013b, Schoonheim et al., 2012, Van Schependom et al., 2014), leading to impaired integration of network information (Rocca et al., 2014) and clustering (Helekar et al., 2010). This approach was able to discriminate patients from controls (Richiardi et al., 2012); however, it did not allow to uncover connections with clinical parameters.

The prediction of disease duration by prolonged microstate class A in the present study was a very stable finding, with an 8-fold cross-validation consistency measure of R2 of 0.970, meaning that the prediction of disease duration is highly valid for any new unknown patient. Age as a covariable improved the explained variance from 21.3% to 31.9%, meaning not only did (a) age add supplementary information to the model, which was not composed of disease duration alone (although the factorial analysis grouped age and disease duration together in the same component), but (b) the prediction of disease duration was improved by only 10% when including age, whereas two thirds were already explained by the duration of microstate class A. High 2y-ARR was predicted by shorter durations of microstate class B together with young age and, therefore, holds the potential to be an interesting surrogate marker for disease activity. Also here was age part of the model, and our results are consistent with a recent study that convincingly showed that age was a more important determinant of decline in relapse rate than disease duration and that decline in relapse activity over time was more closely related to patient age than was disease duration (Kalincik et al., 2013).

Interestingly, we found a dissociation between cognitive fatigue (FSMC-cog), which was predicted by altered EEG fluctuations, and motor fatigue (FSMC-mot), for which we did not find any predictors. A high FSMC-cog score was predicted by shorter microstate durations, similar to the depression score (CES-D). This was only true for patients with RRMS, whereas there was no correlation in healthy controls. Note that formally diagnosed depression was an exclusion criterion in our study; therefore,the CES-D score only served as a screening tool for self-reported mood disturbances. However, both CES-D and FSMC were independent of disease duration or relapse rate, and FSMC-cog was associated with clinical disability measures (EDSS), whereas CES-D was correlated to information speed processing (SDMT) and sex. Also, we did not find any microstate parameters predicting information processing speed (SDMT).

The results of the other very few EEG and MEG resting state studies in patients with MS converge with this most frequently observed MS-related brain pathology of axonal loss and atrophy in the corpus callosum, with their finding of described decreased interhemispheric coherence (Cover et al., 2006, Hardmeier et al., 2012, Lenne et al., 2013, Leocani et al., 2000, Schoonheim et al., 2013b, Van Schependom et al., 2014). These studies showed focally higher synchronization likelihood in temporo-occipital and parieto-occipital connections and lower inter-temporal connectivity (Schoonheim et al., 2013b) in patients with RRMS, with a shift from temporal to less efficient parietal processing, correlating with a decrease in cognitive performance (Hardmeier et al., 2012). Abnormalities in oscillatory brain dynamics in the alpha band, especially in the occipital, temporal, and parietal regions, were demonstrated (Van der Meer et al., 2013), and functional connectivity in the alpha band was lower in the DMN and the visual network (correlating with thalamic atrophy), whereas functional connectivity in the beta band was higher in the DMN and in the temporo-parietal network, both correlating with clinical disability (Tewarie et al., 2013). Decreased interhemispheric connectivity is one of the most robust findings among fMRI studies on resting state in patients with MS (Bonavita et al., 2011, Hawellek et al., 2011, Leonardi et al., 2013, Richiardi et al., 2012). Whereas the above-mentioned coherence studies (Cover et al., 2006, Hardmeier et al., 2012, Lenne et al., 2013, Leocani et al., 2000, Schoonheim et al., 2013b, Van Schependom et al., 2014) concentrated on spatial variations of average correlations between regions of interest, they investigate average values across a time window and are primarily blind to the temporal variability of the dynamic connectivity (Allen et al., 2014, Leonardi et al., 2013), which was detected with the microstate approach described here. Note that the basic spectral analysis of the EEG signal in our patients did not reveal any difference between patients with RRMS and controls; in particular, there was no diffuse slowing in the signal as it had been described in patients with advanced MS (Beach et al., 2011, Striano et al., 2003). It is, therefore, unlikely that the found slowing in the topographical fluctuations originated from a basic difference of EEG frequency power. Independence of microstates from frequency bands and absolute spectral power has repeatedly been shown (Britz et al., 2010, Khanna et al., 2014). The EEG microstate analysis is thus a complementary approach to the classical EEG frequency analysis. It specifically looks at spatial patterns of the scalp electric field instead of variations of oscillations in specific frequency bands at certain electrode locations. Since changes in scalp electric field spatial patterns directly indicate changes in the configuration of neuronal generators, the method looks at the switching behavior of different neuronal networks rather than at changes in the strength of oscillatory activities in these networks. Since microstate analysis is usually performed on broad-band activity, it does not give information about the specific frequencies underlying the aberrant temporal dynamics in microstate fluctuations. Further studies using frequency-resolved spatial analysis as proposed by Koenig et al. (2005) are needed to link microstate fluctuations to changes in synchronization of neuronal populations and changes in functional connectivity between the nodes of the large-scale networks thereof.

In patients with MS, thalamic atrophy occurs early and has frequently been described to lead to thalamic dysfunction and thalamo-cortical dysconnections (Cifelli et al., 2002, Houtchens et al., 2007, Schoonheim et al., 2015a, Wylezinska et al., 2003). Therefore, it is possible that the altered microstate parameters we found in our study were reflecting alterations in the thalamo-cortical loop. Because scalp EEG does not record thalamic activity, such alterations cannot be confirmed by EEG recordings alone. In a recent combined EEG/fMRI study, Schwab et al. (2015) demonstrated selective BOLD activity in specific parts of the thalamus with respect to different microstate classes in specific frequency bands, meaning that thalamic activity indeed participates in the microstate fluctuations. However, in this study, thalamic contributions to the microstate fluctuations were found for all microstate classes except class A. Therefore, this study does not directly support the hypothesis that our effects on microstate classes A and B are related to thalamic dysfunctions.

Conclusion

In the present study, we used high-density EEG to map the topographical fluctuations of the scalp potential field in 53 patients with RRMS in order to compare them with 49 healthy controls. We found that distinct microstates displayed altered durations in the patient group, which powerfully predicted the disease duration, annual relapse rate, as well as depression score, and cognitive fatigue, which is typically experienced by patients with MS. Topographical fluctuations (microstates) analysis might therefore have the potential to serve as a tool to characterize evolution in RRMS (i.e., as a prognostic marker of disease progression). Our patients were, on average, minimally disabled, with a mean disease duration of 9 years.

Given the promising results regarding the prediction of the clinical evolution, it would be interesting to evaluate the approach presented here at the time of initial diagnosis in order to determine if microstate classes A and B already display the alteration observed here at the very onset of the disease. Longitudinal prospective designs are mandatory in patients with clinically isolated syndrome in order to characterize possible indicators for the conversion into RRMS and to better determine the value of EEG resting state topographical fluctuations as surrogate markers for disease evolution.

The following are the supplementary data related to this article.

Supplementary Table 1

Detailed Patient specifications.

Supplementary Table 2

Hierarchical multiple regression controlling for age and sex in each model (corresponding to Table 4).

Author contributions

MG designed the study, analyzed the data, and wrote the paper. MH designed the study, recruited the patients, recorded the data, and revised the paper. DVDV contributed analysis tools and revised the paper. MIT preprocessed the data. IKP analyzed the data and revised the paper; YN classified and recruted the patients; and PF, CM, and MS designed the study, obtained the funding, and revised the paper.

Conflict of interest

None of the authors report a conflict of interest in relation to this work.