Frontoparietal connectivity correlates with working memory performance in multiple sclerosis

Working Memory (WM) impairment is the most common cognitive deficit of patients with Multiple Sclerosis (MS). However, evidence of its neurobiological mechanisms is scarce. Here we recorded electroencephalographic activity of twenty patients with relapsing-remitting MS and minimal cognitive deficit, and 20 healthy control (HC) subjects while they solved a WM task. In spite of similar performance, the HC group demonstrated both a correlation between temporoparietal theta activity and memory load, and a correlation between medial frontal theta activity and successful memory performances. MS patients did not show theses correlations leading significant differences between groups. Moreover, cortical connectivity analyses using granger causality and phase-amplitude coupling between theta and gamma revealed that HC group, but not MS group, presented a load-modulated progression of the frontal-to-parietal connectivity. This connectivity correlated with working memory capacity in MS groups. This early alterations in the oscillatory dynamics underlaying working memory could be useful for plan therapeutic interventions.

time-frequency eeG analysis. We studied the effect of both memory load and successful memory performance on the power of brain oscillatory activity. Regarding the effect of memory load, in the control group we found two effects as expected, a positive modulation in theta activity and a negative modulation in alpha/beta activity 13,23 . In contrast, we did not find modulations of theta band (6)(7)(8)(9)(10)(11) in patients with MS, leading to a significant difference (p = 7 e-7; Cluster Based Permutation (CBP) test) between groups in the initial stage of encoding in this frequency band (between 0.3 and 1 s, Fig. 3D-F). This difference in theta had a topographic distribution located in electrodes of the left hemisphere (Fig. 3I). When analyzing group differences during the maintenance stage, we observed differences in theta activity (5)(6)(7)(8)(9) in the period of 1.3 to 3.3 seconds (p = 0.008), which corresponded to the final part of the encoding stage and the entire maintenance phase. The topographic distribution of this modulation was placed over frontoparietal regions with left predominance (Fig. 3M). Next, we analyzed the oscillatory activity related to successful memory performance (SMP), that is, the specific activity in the trials in which subjects correctly respond to the target stimuli in relation to those trials where subjects make mistakes (Fig. 4). We found that during the maintenance period (2.3 to 2.7 s), patients with MS showed a negative modulation in theta activity (5)(6)(7)(8)(9)(10)(11), while healthy subjects presented a positive modulation (Fig. 4A). This led to a significant difference between groups (p = 0.01) showing a medial frontal topographic distribution (Fig. 4C). Regarding the analysis of alpha activity (10-15 Hz), we observed that healthy subjects presented a decrease in oscillatory activity in this band in the final period of encoding and beginning of maintenance with an occipital distribution, which was not observed in patients with MS ( Fig. 4 right). This led to differences between both groups (p = 0.008) in parieto-occipital electrodes (Fig. 4C right).
Source reconstruction. We carried out source reconstructions of the significant differences between groups in theta activity. For encoding we found that theta modulation was placed in left parietal and temporal cortex (Fig. 5, False Discovery Rate, FDR, q < 0.05 vertex corrected, and p < 0.05 cluster-corrected). Two additional clusters were found in right temporal and medial frontal cortex (p < 0.05 cluster-corrected). For maintenance, we found that the main modulation was placed again in left temporal and parietal lobes (Fig. 5, FDR q < 0.05 vertex corrected, and p < 0.05 cluster-corrected). Additionally, we found a cluster in right orbitofrontal cortex and right parietal cortex (Fig. 5, FDR q < 0.05 vertex corrected, and p < 0.05 cluster-corrected). For the modulation related to successful memory performance, the source of theta was placed over the medial prefrontal cortex (Fig. 5, FDR q < 0.05 vertex corrected, and p < 0.05 cluster corrected).
connectivity. Considering the results of the oscillatory activity, we carried out a connectivity analysis. We selected a frontal electrode (Fz) and a left parietal electrode (Cp3), since modulation in theta for both memory load and successful memory maintenance were found in these electrodes. The following source reconstruction enabled us to infer that they represent frontal and temporo-parietal activity respectively. We first used GC that measures statistic dependency between signals (e.i., if one signal is useful in forescasting another signal). During the maintenance stage, we found that healthy subjects presented an increase in the parietal-to-frontal connectivity in the time domain, as an indicator of successful memory performance (−0.72, p = 0.01). This led to a significant difference in patients with MS, who did not present this modulation (diff, p = 0.009, Fig. 6 and Table 2). Additionally, this modulation changed in the interaction between memory load and successful memory performance, reversing the direction on high memory loads from frontal-to-parietal (p = 0.0089, see Table 2 and Fig. 6). Patients with MS did not show this pattern, demonstrating a significant difference when compared to healthy subjects (p = 0.04, see Table 2). No significant modulations were found for the encoding period in both groups ( Table 2).
We then used a complementary approach to infer the flow f the information. Brain connectivity between cortical regions is regulated principally by the projections between supragranular and infragranular layers 24,25 . Connections originating in superficial and deep layers have distinct spectral fingerprints 26,27 . Thus, the interactions between frequency bands (analyzed by mean of cross-frequency coupling or PAC) can reveal the information flow 28,29 . Thus, it is possible to infer that if the phase of a low frequency of one cortical area correlates with the amplitude of a higher frequency of another cortical area, the information flow could have this direction. Using a similar single trial modeling as in GC analyses, we found similar results. In healthy subjects, we found a modulation between the delta/theta phase of frontal area (Fz) and the gamma amplitude of parietal areas (CP3). This modulation correlated negatively with successful memory performance but positively with the interaction between successful memory performance and memory load (Fig. 6C). In contrast, we found a second modulation between the theta/alpha phase of parietal areas (CP3) and the beta amplitude of frontal areas (Fz), with the opposite patterns. This second modulation correlated positively with successful memory performance and memory load, but negatively with the interaction between successful memory performance and memory load (Fig. 6C). None of these modulations were observed in patients with MS, leading to significant differences between groups. clinical correlations. Finally, we explore whether the frontal to parietal connectivity correlated with clinical assessment. We used the regressor of the interaction between memory load and successful memory performance in Granger Causality analysis and the Neuropsychological test. Interestingly, in the patient group, the frontal-to-parietal connectivity was correlated with PASAT that measure WM capacity (Spearman correlation, rho = 0. 58  www.nature.com/scientificreports www.nature.com/scientificreports/

Discussion
In this study we assessed the cortical circuits that support verbal working memory in patients with relapsing-remitting MS with minimal or no burden of neurological disability (EDSS < 3 and PASAT > −1.5 SD). In the early stage of the disease, patients commonly do not present objective cognitive alterations in the neuropsychological evaluations, however, they manifest a subjective sensation of difficulties in their cognitive performance in daily activities (e.g., occupational or academic tasks) 30 . These patients report a common clinical pattern of difficulties in the performance of daily cognitive tasks involving working memory, but without clear evidence in the clinical tests applied in their evaluation routines [30][31][32] . Accordingly, we observed no significant differences in WM performance, but clear neurophysiological differences between patient and control groups. In the early stage of the disease, the clinical evaluations applied to MS patients without impairment often do not detect a deterioration of this cognitive function [33][34][35] . The RT difference that we found could reflect an increase in the cognitive effort necessary to correctly solve the task and could be an early behavioral marker of WM impairment. In fact, reaction time analyses are a good marker of attention and cognitive control dynamics [36][37][38][39][40] .
In spite of no clear evidence for behavioral impairments in the accuracy in the WM task, patients with MS demonstrated a distributed oscillatory activity reorganization. The ability to maintain a sustained activity in the fronto-parietal network in the absence of sensory stimuli depends on the synchronous structured activity of different frequency ranges 41 . Our working memory task elicited a recognized oscillatory activity in a frontal and parietal network 11,[42][43][44] . The left lateralization of the power of the theta oscillatory activity could be related to the areas involved in language, i.e., the phonological loop, required for the specific information of our stimuli (consonants) 13,45 . Low frequency synchronization between the medial frontal region and temporal cortex has been demonstrated in non-verbal WM tasks, probably reflecting cognitive control 15,23 . Indeed, our connectivity analyses showed specific dynamics, in which frontal to parietal/temporal influences increased in function of memory load when a subsequent successful memory performance occurs. These findings are in line with several recent reports indicating that frontal interactions with other cortical areas are key aspects of successful memory performance 14,46,47 . In contrast, patients with MS demonstrated a loss of the WM oscillatory dynamics. These patterns could represent an inefficient cognitive effort to keep the stimuli information active in WM, which would be in accordance with the longer reaction times for erroneous responses in high memory load. In addition to this, the patient group did not present an increase in frontal-to-parietal connectivity 48,49 . WM requires the synchronization of neural network connections distributed in the prefrontal and parietal regions in order to integrate complex information for generating appropriate responses. These distributed brain networks are especially sensitive to the diffuse damage of white and gray substance found in multiple sclerosis 17,18,50 . www.nature.com/scientificreports www.nature.com/scientificreports/ Thus, our results are in accordance with recent findings in patients with mild or minimal cognitive deficit that show reorganization in electrophysiological activity 49,51 . In spite of the fact that patients could have fewer resources to maintain the stimuli in WM they can achieve a similar performance to healthy subjects, as long as the tasks are not too, by means of redistribution of the remaining and available resources. It has been proposed that the central executive may be the main component of WM that is disrupted in cognitively impaired multiple www.nature.com/scientificreports www.nature.com/scientificreports/ sclerosis patients 52 . This proposition is supported indirectly by evidence that suggests that the dorsolateral prefrontal cortex, which is thought to underlie executive control, is commonly recruited when there is heavy demand placed on WM in individuals with brain injury 53 . Our findings give further evidence to support this interpretation. Indeed, the frontal low frequency influences the gamma power in parietal and temporal areas, and this influence increases in relation to memory load increase. Patients fail to display this influence, and this failure could reflect the loss of the mechanisms by which the control process carried out by frontal areas produces a successful memory performance.
Studies using fMRI have found several changes in brain activity during WM tasks including decrease and increase in both connectivity and activity in patients with MS 54,55 . Indeed, follow-up studies show an inverted-U form in the evolution of the disease in resting-state functional connectivity 56 . Early changes in the connectivity/ functional patterns have been interpreted as compensatory changes in prefrontal cortical regions underlying modulations of executive aspects of WM [57][58][59] . For instance, studies in very early states of MS have found medial PFC activity increases, although later meta-analyses studies have revealed a decrease in medial prefrontal activity. In this context, our study gives several insights into the participation of medial prefrontal activity in working memory and its relationship with the cognitive dysfunction in patients with MS. Medial prefrontal theta activity is related specifically to subsequent successful memory performances. As this oscillatory activity has been generally related to cognitive control in several tasks 60 , it is possible to interpret this activity index as the cognitive effort related to successful performance. Interestingly, patients present more activity for errors and also present an increase of reaction time for incorrect responses. Thus, the oscillatory pattern demonstrated by patients with MS may indicate a decrease in the detection of the need for additional cognitive control 61,62 . In line with this, patients do not present an increase of frontal to parietal influences in relation to successful memory. Indeed, patients with MS and cognitive impairments demonstrate low frequency amplitude decreases during cognitive tasks 63 . Along with this, it has been observed that the nested frequencies between theta and gamma can represent the way in which multiple maintenance items are organized into WM 41 . The nested frequencies or PAC refers to a specific type of oscillatory activity, which would reflect a general mechanism by which the cortical areas organize and structure information [64][65][66][67] . Accordingly, we observed a modulation between delta/theta and gamma frequencies in healthy subjects. This communication channel seems specific to the frontal to parietal communication and shows an increase in relation to memory load. While, for low memory load, in information flow seem to be predominantly from parietal to frontal, using a different nested frequency. Patients did not demonstrate these dynamics, reflecting the poor or inefficient frontal parietal communication. Interestingly, this loss of interplay between the frontal and parietal oscillatory dynamics could be an early marker of working memory deficits in patients with MS. Thus, these patients' performance suggests an impairment in the establishment and maintenance of a fluid oscillatory dialogue between the areas involved in the successful development of a WM task. Indeed, connectivity pattern correlated with clinical test of WM.
In summary, in spite of a high incidence of cognitive deficits multiple sclerosis and a large body of literature investigating cognitive dysfunction, the specific oscillatory features that drive cognitive deterioration in multiple sclerosis remain unclear. Hence, our results demonstrate a specific electrophysiological mechanism underlying the WM deficit in patients with MS. Several investigations have revealed that cerebral oscillatory activity supports different cognitive processes and have indexed their alteration in clinical populations [68][69][70] . Indeed, recent evidence indicates that increasing theta activity by means of transcranial Alternating Current Stimulation can improve WM performance in healthy and aging subjects 14,47 . It has been seen that cognitive training can increase memory capacity and that this increase correlates with changes in connectivity of distant brain areas, specifically between frontal and parietal regions. Therefore, the specific oscillatory features related to WM deficits identified in our study could serve to implement non-pharmacological treatments using non-invasive brain stimulation and cognitive training, in order to contribute to the improvement of the quality of life of these patients.

Methods
Design and participants. Our study is a case-control design that include 40 participants. A sample consisted of 20 patients with relapsing-remitting multiple sclerosis in an early stage, with minimal to no clinical evidence of cognitive alterations (Table 1). According to the 2010 McDonald Criteria, the medical diagnosis was made by a Neurologist 71 . Stable patients without episodes of relapses in the last month, with scores of three or more on the Expanded Disability Status Scale (EDSS), with less than −1.5 z-score of Paced Auditory Serial www.nature.com/scientificreports www.nature.com/scientificreports/ Addition Test (PASAT), with non-correctable visual alterations, with a history of traumatic brain injury, neurological and/or psychiatric pathologies, and abuse or regular consumption of drugs or alcohol were excluded. The patients participated in the study while on their usual disease modifying therapies (i.e., immunomodulation therapy only, without other treatment, such as antidepressants or Benzodiazepines). The control group was composed of 20 healthy volunteers, comparable in age, sex, manual preference, and educational level ( Table 1). As well as in the patient group, healthy subjects with non-correctable visual alterations, with a history of neurological and/ or psychiatric pathologies, traumatic brain injury, and abuse or regular consumption of drugs or alcohol were excluded. All participants were Spanish native speakers and provided signed informed consent prior to participation in the study. Patients underwent neuropsychological assessment during the month previous to the EEG session. This was their first neuropsychological evaluation. These assessments included PASAT that measures cognitive processing speed and working memory, Symbol Digit Modalities Test (SDMT) that measures cognitive processing speed, the Brief Visuospatial Memory Test-Revised (BVMT-R) that measures visuospatial memory, and the World Health Organization-University of California-Los Angeles Auditory Verbal Learning Test (WHO UCLA AVLT) that measures verbal memory. Details are summarized in Table 1.
The experimental protocol and all methods were performance in accordance to institutional guidelines and were approved by the Ethical Committee of the Pontificia Universidad Catoĺica de Chile.
Sample size. For the estimation of the minimum required sample size the following parameters were considered: (a) Effect size for the mixed ANOVA statistical test (2 ×2, with interaction effects), b) Statistical power (1β)=0.95 and c) Significance level =0.05. Considering an effect size 2 = 0.09 (effect size F = 0.3 13,72 ), the sample size amounts to a total of 40 participants (n1 = 20; n2 = 20).
Experimental task. In this study we implement a modified version of Sternberg´s Memory Scanning task [Jensen, 2002]. This task consisted of a list of consonants simultaneously presented and displayed in a circular arrangement with a fixation cross in the center of a computer monitor located 57 cm from the subject. The letter "Y" was not included in the memory set to avoid the generation of words that could be used as clues by the subjects. Each memory set arrangement consisted of groups of two, four, or six consonants generating three levels of WM load. The latter refers to the progressive number of stimuli to be stored and manipulated in WM. All the stimuli were placed foveally, minimizing the effect generated by saccadic movements (see the experimental task outlined in Fig. 1). Each memory set was presented for 1800 ms (encoding period), followed by a black screen with the fixation cross (maintenance period) of 2000 ms and then a recovery period in which the fixation point was replaced by a target stimulus for 1000 ms. The subjects were instructed to memorize the memory set and then report whether the target stimulus was present or absent in the memory set, using the right or left hand alternately. Subjects had 2200 ms to answer. Each subject had to respond 270 trials (90 trials for each memory load set). The trials were presented in two main blocks divided by a pause regulated by the subject. In addition, each main block was constituted by 15 sub-blocks, formed by 9 trials of the same memory load each. The order of the memory load was randomized. Both the presentation of the stimuli and the recording of the test responses were done with the Software Presentation ® (Version 13.0, www.neurobs.com).
Electrophysiological recordings. Continuous EEG recordings were obtained with a 40-electrode EEG System (NuAmps, Neuroscan). All impedances were kept under 5kΩ. Electrode impedance was retested during pauses to ensure stable values throughout the experiment. All electrodes were referenced to averaged mastoids during acquisition and the signal was digitized at 1 kHz. Electro-oculogram was obtained with four electrodes. All recordings were acquired using Scan 4.3 and stored for off-line analysis. At the end of each session, electrode position and head points were digitalized using a 3D tracking system (Polhemus Isotrak).
Electrophysiological data analysis. EEG signals were preprocessed using a 0.1-100 Hz band-pass filter. Eye blinks were identified by a threshold criterion of ±100 μV, and their contribution was removed from each dataset using Independent Component Analysis (ICA). Other remaining artifacts (e.g., muscular artifacts) were detected by visual inspection of both the raw signal and the spectrogram. We thus obtained 243 ± 28 artifact-free trials per subject. All artifact-free trials were transformed into current source density (CSD) that was estimated using the spherical spline surface Laplacian algorithm suggested by Perrin et al. 73 , and implemented by Kayser and Tenke 74,75 . Induced power distribution was computed using Wavelet transform, with a 5-cycle Morlet wavelet, in a − 0.5 to 3.8 s window around the onset of the memory set stimuli. This time-window includes 0.5 seconds of inter stimulus interval, 1.8 seconds of the stimulus of the memory set and 2 seconds of maintenance period. For all analyses, we used the dB of power related to the baseline (15 seconds acquired in the beginning of each block).
Source reconstruction. The neural current density time series at source levels were calculated by applying a weighted minimum norm to estimate inverse solution 76 with unconstrained dipole orientations in single trials as in prior work 68,77,78 . We used a default anatomy of the Montreal Neurological Institute (MNI/Colin27) wrapped to the individual head shape (using ~300 head points per subject). We defined 3 ×4000 sources constrained to the segmented gray cortical volume (3 orthogonal sources at each spatial location) in order to compute a three-layer (scalp, inner skull, outer skull) boundary element conductivity model and the physical forward model 79 . Since the inverse solution is a linear transformation, it does not modify the spectral content of the underlying sources. Therefore, it is possible to undertake time-frequency analyses directly in the source space. Finally, we reduced the number of sources by keeping a single source at each spatial location that pointed into the direction of maximal variance. For this, we applied a principal component analysis to the covariance matrix obtained from the 3 orthogonal time series estimated at each source location. Since we used a small number of electrodes (40) and no individual anatomy for head model calculation, the spatial precision of the source estimations is limited. In order Scientific RepoRtS | (2020) 10:9310 | https://doi.org/10.1038/s41598-020-66279-0 www.nature.com/scientificreports www.nature.com/scientificreports/ to minimize the possibility of erroneous results, we only present source estimations if there are both statistically significant differences at the electrode level and the differences at the source levels survive a multiple comparison correction (cluster-based permutation test and vertex correction using false discovery rate, q = 0.05).

Statistical analysis.
We used the Kolmogorov-Smirnoff to test for normality. When the data did not meet the normal assumption, we used non-parametric tests. We evaluated pair comparisons using Wilcoxon test and Bonferroni correction. For the EEG statistical analysis, we first fitted a General Linear Model (GLM) of the power of the oscillatory activity per trial in each subject (first level analysis, see 13,70,80 ), where b 1 is the intercept, and b 2 is the slope or coefficient for the variable Memory Load (ML, ordinal variable that takes the value 2, 4, or 6 depending on the memory load condition) and b 3 is the slope for the variable Successful Memory Performance (SMP, dummy variable that takes the value 0 if the subject makes a mistake in this trail or 1 if the target stimulus is correctly identified). We thus obtained a 3D matrix of t-values (sensor, time, frequency) for each regressor and subject. We then explored for differences between groups and conditions using the Wilcoxon test (second level analysis). To correct for multiple comparisons in time-frequency charts, we used the Cluster-based Permutation (CBP) test 81 using 1000 permutations. For more detail see prior work 70,82 . causal interactions. To evaluate the influence of frontal regions over temporo-parietal regions, we estimated Granger Causality (GC) 83 between selected electrodes (Fz and CP3). The causality was calculated over time series per trials. See detail in 82,83 . We obtained a GC term per trial that was then used in the modeling analyses (see below).

phase-amplitude coupling (pAc).
To evaluate cross-frequency modulation, we carried out a phase-amplitude coupling analysis (PAC) as described in 64 . Briefly, for a given frequency pair, the raw signal was filtered separately in both frequencies (zero phase shift non-causal finite impulse filter with 0.5 Hz roll-off). The lower frequency ranged from 1 to 35 Hz (0.4 Hz increments, 0.8 Hz filtered bandwidth) and the higher ranged from 5 to 120 Hz (1 Hz increments, 5 Hz filtered bandwidth). The phase of the lower frequency range and the amplitude of the higher frequency range were computed using Hilbert transformation. For each epoch of maintenance per trial (1.8 to 3.8 s), we computed the circular-linear correlation between the phase of lower frequencies and the amplitude of higher frequencies. Thus, we obtained a circular correlation coefficient (CCC) per trial and per combination of pairs between the selected electrodes (Fz-Fz, Fz-Cp3, Cp3-Fz, Cp3-Cp3). Thus, we obtained a CCC per trial and pair of electrodes that were used in the modeling analyses (see below).
Gc and pAc modeling. For both GC and PAC analyses, we used the following modeling. At subject level analysis (first level analysis), we compared if the circular correlation coefficient (CCC) or granger coefficient (GC) variation correlated with the task parameters using multiple linear regression as the following equations depict: At group level analysis per frequency pair, we compared whether t values of each regressor were statistically different from zero using the Wilcoxon signed-rank test. For PAC analyses, we corrected for multiple comparisons using a CBP test. The initial threshold for cluster detection was p < 0.05, and the final threshold for significant cluster was p < 0.01. For between-group comparison we used Wilcoxon rank-sum test and Bonferroni correction.
Software. All behavioral statistical analyses were performed in R. The EEG signal processing was implemented in MATLAB using CSD toolbox, in-house scripts (available online as LANtoolbox http://neuroCICS.udd.cl/ LANtoolbox.html), BrainStorm 84 and open MEEG toolboxes 85 .

Data availability
All Data and script used for data analyses are freely available on http://neuroCICS.udd.cl.