Functional brain networks reveal the existence of cognitive reserve and the interplay between network topology and dynamics

We investigated how the organization of functional brain networks was related to cognitive reserve (CR) during a memory task in healthy aging. We obtained the magnetoencephalographic functional networks of 20 elders with a high or low CR level to analyse the differences at network features. We reported a negative correlation between synchronization of the whole network and CR, and observed differences both at the node and at the network level in: the average shortest path and the network outreach. Individuals with high CR required functional networks with lower links to successfully carry out the memory task. These results may indicate that those individuals with low CR level exhibited a dual pattern of compensation and network impairment, since their functioning was more energetically costly to perform the task as the high CR group. Additionally, we evaluated how the dynamical properties of the different brain regions were correlated to the network parameters obtaining that entropy was positively correlated with the strength and clustering coefficient, while complexity behaved conversely. Consequently, highly connected nodes of the functional networks showed a more stochastic and less complex signal. We consider that network approach may be a relevant tool to better understand brain functioning in aging.

"protective" factor by reducing the incidence of dementia 5 , and its pathological consequences in AD patients 7,8 . However, recent studies did not find this effect 9,10 , suggesting that CR actually may have a "masking" effect.
One of the new approaches to explore how brain functioning may be modulated by CR is functional network organization (based on functional connectivity; FC), although there is still scarcity of studies focused on this field. Solé-Padullés et al. 11 described that healthy elders with higher CR showed reduced functional magnetic resonance imaging (fMRI) activation during cognitive processing. In the same line, a previous study from our group 12 observed that healthy elders with a high CR score in comparison with those with low CR score, exhibited less magnetoencephalographic (MEG) functional connectivity during the performance of a memory task. These results may suggest that higher CR is related to more effective use of cerebral resources although this hypothesis has not been tested yet.
In addition, FC may be used to explore brain's organization in terms of structural and functional networks, by the application of different methodologies coming from Network Science 13 . This kind of analysis can introduce a different perspective to better understand how CR modules functional resilience in healthy aging. Recently, Yoo et al. 14 and Marques et al. 15 have proposed the analysis of the topology of the associated functional networks as an alternative way of quantifying CR. Interestingly, functional brain networks obtained during resting state show a positive correlation between CR and network efficiency 16,17 . Network efficiency, first introduced by Latora et al. 18 is a metric evaluating the inverse of the topological distance between nodes of a network. In the context of functional brain networks, the topological distance between two brain regions is normally defined as the inverse of their level of coordination. In this way, nodes that are not functionally connected (i.e., with zero coordination) have infinite distance between them. However, information between not connected regions (at an infinite topological distance) can travel indirectly through intermediate nodes. Network efficiency measures how close or distant the nodes of a network are in terms of their topological distance using the functional connections through intermediate nodes 18 . Under this framework, the positive correlation between CR and network efficiency during resting state relies on fMRI recordings, which report higher levels of synchronization between brain regions in individuals with higher CR 15,16 . However, during a cognitive process, the correlation between CR and synchronization turns to be negative, as shown in 11 . With the aim of clarifying whether CR is related to the topology of the functional networks during a memory task, we are concerned about the biomagnetic network profiles of two groups of healthy elders with different CR levels. Specifically, we obtained the MEG recordings of both groups during the performance of a modified Stenberg's test 17 and we explored the topological properties of the corresponding functional networks. Among the diversity of network metrics that one may extract from the functional networks we focused on the most extended ones: the network strength, which is a measure of the average level of synchronization along the functional network 19 ; the weighted clustering coefficient, measuring the existence of triangles in the network and related to the local resilience 20 ; the eigenvector centrality a measure of node importance obtained from the eigenvector associated to the largest eigenvalue of the adjacency matrix 21 ; the average shortest path d(i), measuring the average number of steps to travel from one node to another through the shortest topological distance 22 ; the strength of the nearest neighbors snn(i), indicating what is the average strength (i.e., sum of the weights of a node) of the neighbors (i.e., nodes directly connected) of a given node 19 ; and the outreach o(i), a network metric obtained by multiplying the weight of the connections by their Euclidean distance 23 (see Supp. Info. for a detailed definition of all network metrics).
We also quantified the dynamical properties of the brain regions in terms of entropy and complexity to compare both high and low CR groups. To the best of our knowledge, the present MEG network's approach, where the study of network metrics is combined with the dynamical properties of the nodes, has never been used within this field of research. Combining both approaches we have been able to detect differences between low and high CR individuals at the network level, showing that the network outreach and, consequently, the energetic expenses, is higher in individuals with low CR. At the same time, we show how the topological properties of the nodes (i.e., brain regions) are related to the complexity and entropy of their corresponding time series.

Materials and Methods
Participants. The sample consisted of 20 healthy elders recruited from the "Geriatric Unit of the University Hospital San Carlos" (Madrid, Spain). All of them were right-handed 24 , native Spanish speakers, and between 65 and 85 years old. Regarding the neuropsychological assessment, all of them met the following criteria: (1) a score above 28 in MMSE 25 ; (2) no memory impairment as evidenced by delayed recall from Logical Memory II subtest of the Wechsler Memory Scale Revised (WMS-III-R 26 ) and (3) normal daily living activities measured by the Spanish version of the Functional Assessment Scale 27 . Exclusion criteria included: (1) a history of psychiatric or neurological disease; (2) psychoactive drugs consumption; and (3) severe sensory or comprehension deficits.
The 20 participants did not differ in physical, cognitive and social activities during late-life, and were further subdivided into two subgroups according to their CR level. The CR index (CRI) was estimated for each individual by adding the following two categories: (1) formal educational level (from 1 -illiterate/functional illiterate-to 5superior studies-) and (2) occupational attainment (from 1 -housewife-to 5 -business management/research-). According to their CRI score, 12 participants were classified into the low CR group if their CRI score was between 1 and 5; or into the high CR group (8 subjects) if their CRI score was between 6 and 10. In addition, both groups did not differ in age and in their performance during the cognitive task (see Table 1).
Ethics statement: Methods were carried out in "accordance" with the approved guidelines. The investigation was approved by the local Ethics Committee (Madrid) and all participants signed a written informed consent before the MEG recordings.

MEG task.
A modified version of the Sternberg's letter probe task was employed 17 (see Fig. 1). Each subject was asked to remember a single set of 5 letters (i.e. 'SMAQE'). After that, there were presented a series of single letters, one at a time (500 ms in duration with a random ISI between 2 and 3 seconds). They were asked to make a match/non-match decision by pressing a button with their right hand if the presented letter was one of the five initially memorized. The experiment consisted of 250 letters, in which half were initially presented and the rest were distracters, meaning they were not in the initial 5 letter set. All of the participants completed a training session until demonstrating they had correctly memorized the initial letter set. The task was projected through a LCD video projector (SONY VPLX600E) located outside the shielded MEG room through a series of mirrors. The screen was suspended approximately 1 meter above the participant's face. The letters subtended 1, 8 and 3 degrees of horizontal and vertical visual angle respectively. MEG recordings. MEG data was recorded during the execution of the modified Sternberg's task with a 254 Hz sampling rate and a band pass of 0.5 to 50 Hz, using a 148 channel whole-head magnetometer (MAGNES 2500 WH, 4D Neuroimaging) placed in a magnetically shielded room. In order to reduce external noise we employed an algorithm using reference channels at a distance from the MEG sensors. After trial segmentation, artifacts were visually inspected by an expert and epochs containing muscular artifacts or eye movements were discarded for further analysis. Only hits were considered, removing false alarms, correct rejections and omissions. We randomly selected a set of 35 trials for each individual, since it was the lowest number of successful epochs for all subjects in the study. Each trial consisted of each 230 time points for each of the 148 channels. Trials are not  Table 1. Mean values ± standard deviation for high and low CR groups. Significant p values are p < 0.05 and marked with an asterisk (*). Both groups differ in their educational level, occupation attainment and CRI scores, but not in the other variables. Acc. Smaqe, accuracy in the modified Stenberg's memory task.

Figure 1.
Representation of the modified version of the Sternberg's letter probe task. In the encoding phase, participants are instructed to memorize 5 letters (i.e.: "SMAQE"). In the recognition phase, participants are instructed to make a match/non-match button-press to indicate that the presented letter matched any of the encoded ones. necessary consecutive, however, this issue is not crucial for the analysis of the coordination between sensors since each trial is analyzed individually and averaged later.

Functional connectivity analysis.
To calculate the connectivity networks for each subject we used Synchronization Likelihood algorithm (SL) 28 which is a nonlinear measure of the synchronized activity that has been proven to be a suitable quantifier for datasets obtained from magnetoencephalographic recordings 29 . Specifically, the SL algorithm detects windows of repeated patterns within the time series of a channel A and, next, checks whether a channel B also shows a repeated pattern at the same time windows, no matter if it is the same or different to that observed in channel A. The range of values of SL is 0 ≤ SL ≤ 1, being 0 when the time series of both channels are uncorrelated, and 1 for maximal synchronization 28 . SL was calculated for each pair of sensors since, in our study, sensors were the nodes of the functional networks. All pairs of SL are then included into a correlation matrix W{w ij } where w ij have values comprised between ~0.05 and ~0.5. We apply a linear normalization that leads to a probability matrix P{p ij }, where the values of p ij are obtained as In this way, the probability matrix P{p ij }, whose values are within the interval [0, 1], reflects the probability of the presence of a link between nodes (brain regions) i and j.
Graph metrics. We calculated a series of metrics related with the role of the nodes within the network: the strength s(i), the weighted clustering cw(i), the eigenvector centrality ec(i), the average shortest path d(i), the strength of the nearest neighbours snn(i) and the outreach o(i). We used the inverse of the values of the probability matrix P{p ij } for obtaining the "topological distances" between nodes to quantify the shortest paths. Another parameters such as the within-module degree z(i) (also kwon as the z-score) and the participation coefficient p(i) 23,30 were also computed using the community affiliation vector Ccom(i) extracted from the classical partition of the brain into six lobes: left-frontal, right-frontal, central, left-temporal, right-temporal and occipital. A series of global network features such as the global efficiency Eg and the average shortest path (d) were also calculated. The node average (average along the same node for all subjects within the same group) of several metrics was computed to obtain the following mean values: ( ) Network averages of the preceding features were also appraised in order to compare both groups. Finally, we obtained the global efficiency E g of the networks from the harmonic mean of the inverse of the shortest paths.
We constructed a set of randomized versions of the former functional networks, in order to evaluate to what extent the deviation of the network parameters between groups was a consequence of a topological reorganization or, on the contrary, just a matter of the number of links and their weights. With this purpose, we generated a group of 100 networks for each of the functional network. The randomization maintains the value of the links' weights by reshuffling the components of the weighted probability matrix P{p ij } (see 23 for a similar normalization). In this way, we guaranteed that the average strength of the network was maintained. Next, we calculated the network parameters for each of the randomized versions and obtained an average value of each metric. Finally, we normalized all network parameters with the average of the set of surrogate matrices, i.e., for a parameter X its normalized value would be: Note that this normalization allows focusing on changes related to the structure of the functional networks, since they exclude variations of the average weights of the networks (and their influence on the network metrics).
Evaluation of the dynamical properties of the nodes. In order to calculate the entropy and complexity brain signal, we used for each of the nodes in the network their corresponding time series comprising the 35 trials of the whole experiment. Since each trial has 230 time steps, we obtained times series of M = 8050 points for each node of each functional network.
Next, we applied the methodology introduced by 31 to quantify the dynamical properties of nonlinear dynamical systems. We calculated the normalized permutation entropy of the signal (S) as described in 32 to capture the level of uncertainty of a signal.
We calculated the complexity (H) of the signal for each subject following the procedure described in 33 . This measure is effective to quantify the complexity of different dynamical systems 34 .
Statistical analysis. We performed permutation tests as well as non-parametric ones for all statistical hypotheses of this work. The permutation test is based on a t-statistical distance between two samples. For robustness, we constructed a distribution of 5000 different p-values associated to the number of permutations in each hypothesis. In each distribution, we found the percentile associated to the original p-value, thus dividing this percentile by the number of permutations as a correction of the resampling. We considered the rejection of the null hypothesis when the aforementioned ratio is bellow or equal to α = 0.05 level of significance. Results of permutation tests are presented throughout the paper. Complimentary, we performed a nonparametric Mann-Whitney U-test (rank-sum test) to show how the results do not change. The latter results are shown in the Supplementary Information. . Figure 2 shows that both groups exhibited the most influencing nodes placed at the same cortical regions (i.e., at the occipital lobe). Nevertheless, we observed signs of a displacement of the node centrality toward the central and left-temporal lobes.
In order to quantify the differences between groups, we obtained the average difference of the node centrality ∆ev i ( ( )), within-module degree ∆z i ( ) and participation coefficient ∆p i ( ) of each node for both groups. Figure 3 shows the significant differences obtained in node centrality (A), within-module degree (B) and participation coefficient (C) between groups ≤ . p ( 005) val . The high CR group showed higher values of the eigenvector centrality in nodes over the central lobe, near the parieto-occipital sulcus, as well as some tiny regions in the frontal lobe (i.e., > ev i e v i ( ) ( ) high low and subsequently, ∆ > ev i ( ) 0). Conversely, cortical tissue placed at the lefttemporal and occipital lobes have higher centrality in the low CR population (i.e., ∆ev i ( ) < 0). High CR group presented higher within-module degree over both temporal brain areas and in the occipital lobe, while low CR group showed higher within-module degree in some regions located at the central lobe. Finally, there were few statistical differences between groups in the participation coefficient, showing the high CR group higher values than the low CR one, especially at regions placed at the frontal lobe. Statistical analysis based on a Man-Whitney U test shows similar results (see Supp. Info. for details).
It is worth mentioning that both the eigenvector centrality and the within-module degree are quantifying the percentage of importance a node has with regard to the other nodes of the network. Therefore, they should be interpreted as a way of evaluating the increase/decrease of importance with respect to all nodes, in the case of  eigenvector centrality, or the nodes belonging to the same lobe, in the case of the within-module degree. Note that, if nodes belonging to the temporal and occipital lobes have higher centrality, it should be at expenses of other nodes of the functional network (in this case, those located at the central lobe), since the sum of the centralities of all nodes is constant (due to the fact that the modulus of the eigenvector is one). In other words, when comparing both groups, it would not be possible to observe higher centrality for all nodes of a given group (High CR or Low CR).

Macro-scale:
Analysing the topology of the whole network. We analysed how the average network parameters were modified according to the level of CR of the participants. First, we calculated the mean and standard deviation of six network parameters for each subject (network strength S, outreach O, weighted clustering coefficient C w , average neighbour strength S nn , global efficiency E g and average shortest path d ). Second, we carried out a permutation test and a nonparametric Mann-Whitney U-test (see Section Statistical analysis for details) for each network metric. Finally, those metrics with p-values < 0.05 were accepted to have statistically significant differences.
In our results, the low CR group exhibited higher averages values than the high CR group in all parameters calculated, except for d (see Fig. 4 for details). However, only the network outreach O showed statistical significant differences both in the permutation (p = 0.0232) and the rank-sum (p = 0.039) tests. Interestingly, the low CR group showed a higher O (O low = 13.57 and O high = 12.29). It is also remarkable that the average shortest path length d exhibit significant differences (p = 0.0491) at the rank sum test, despite not passing the permutation test. In this case the average values of shortest path in low CR are lower than in the high CR group (d low = 12.67 and d high = 14.13) (see also

Dynamical Analysis of MEG Time Series.
We found a group of nodes showing statistical differences between the high and low CR groups, both for H(i) and C(i) (p val < 0.05). In accordance with previous sections, the occipital lobe is the brain region that concentrates the highest amount of nodes with statistical differences (see Figs 5 and 6). We obtained a group of 23 nodes whose entropy was higher in the low CR group (except for one node) (see Fig. 5). On the contrary, we obtained lower complexity of the signals in the low CR group (except for two nodes), but in a very similar brain localization that in the case of entropy (only two nodes did not overlap) (see Fig. 6). These results indicate that a higher level of CR could be associated to a brain dynamics with lower entropy and higher complexity.
Correlation between Entropy and Complexity. Then, we explored the relationship between these two dynamical properties. Figure 7 shows a complexity-entropy scatter plot for all nodes, with the inset plotting only those nodes with statistical differences. In the figure, each point corresponds to the average of each node over the subjects of a given group. We can observe a significant negative correlation between H(i) and C(i) for low (r = 0.9874; p val = 0.0010) and high CR groups (r = 0.8288; p val r = 0.8288; p val = 0.0010), revealing that those nodes with less entropy are, at the same time, those nodes with higher complexity. Correlations between Topology and Dynamics. We explored the relationship between the dynamic properties of the nodes (i.e., H i ( ), C i ( )) and two of their topological parameters, namely the node strength S i ( ) and the weighted clustering coefficient C i ( ) w for each CR group. Figure 8 shows all possible combinations. For the low CR group, we found a positive correlation between H i ( ) and S i ( ) (r 2 = 0.5833; p val < 0.001), which indicates that nodes with higher strength, i.e. the hubs of the network, are in turn those nodes with higher entropy at their dynamics (see Fig. 8A). On the contrary, the correlation between the node complexity and its strength was negative (r 2 = 0.7124; p val < 0.001) (see Fig. 8B). The same pattern was found between H i ( ) and C i ( ) w (positive correlation) (r 2 = 0.5713; p val < 0.001), and between C i ( ) and C i ( ) w (negative correlation) (r 2 = 0.7055; p val < 0.001) (see Fig. 8C,D).
Correlations found for the high CR group are qualitatively similar to the low CR but the values of the correlation parameter are systematically lower (See Table S2 of Supp. Info. for details).

Discussion
In the present study, we explored if CR might play a role in both the topology and the dynamics of the functional brain networks in healthy elders while performing a memory task. To this end, we classified the participants according to a CR index that combined educational level and occupational attainment into two groups: high or low CR level.
We firstly focused on the topological characteristics of the two groups, finding differences at microscale (i.e. nodes) and macroscale (i.e. average over all nodes) levels. At the microscale level (i.e., when nodes are analyzed), we observed that the most influential nodes for both groups were mainly placed over the occipital region, although there was a slight displacement towards the central lobe in the group with high CR. Specifically, we found that the low CR group exhibited higher eigenvector centrality over both temporal and occipital areas, while the high CR group presented higher values in central areas. These observations indicate that the importance of regions placed at the occipital and temporal regions is higher at the low CR group. At the same time, since the sum of the eigenvector centrality of all nodes of the network is constant, other nodes unavoidably have lower  centrality values to compensate the higher values reported at the occipital and temporal regions. These nodes are placed at the central lobe, where the high CR group has higher centrality when compared with the low CR one. Interestingly, the within-module degree behaves in the opposite way. Comparing Fig. 3A,B, we can observe that, in general, regions that have higher eigenvector centrality in the high CR group reduce their within-module degree. The explanation of this, in principle, counter-intuitive observation is that nodes have higher eigenvector centrality as a consequence of having higher weights at their long-range connections (i.e., in weight of the functional connections with nodes laying outside of their module).
Besides, high CR group showed higher participation coefficient values than the low CR group, especially in regions located at the frontal lobe. These results pointed out that the importance of a node and its connections with their neighbours inside and outside their community (lobe, in our case) may be differently located during the execution of a cognitive task, according to the level of CR. Red squares correspond to low CR group and blue circles to the high CR group. The inset shows only nodes with statistical significant differences using the permutation test. Specifically, these nodes are the channels that have statistical differences when comparing both the mean entropies and complexities of high vs. low CR groups at node level. See Supp. Info. for similar results obtained with the ranksum test. At macroscale level, the low CR group exhibited higher values in different networks parameters compared to the high CR group (network strength, weighted clustering coefficient, average neighbour strength, global efficiency), being this increment statistically significant only in the outreach parameter both for the permutation and the rank-sum tests. Although strength results did not reach statistical significance, the low CR group presented a higher average in this network parameter suggesting that these participants may require higher brain synchronization to perform the memory task as well as the participants with high CR. As a result of this rise of functional connections, the average shortest path showed lower values, however, only the rank-sum test leaded to statistically significant differences between low and high CR groups. Along with the average shortest path, the outreach was the parameter more altered in the low CR group (in this case passing the two statistical tests). The higher outreach observed in the low CR group indicated the existence of a higher amount of long-range connections, implying a greater energetic cost to accomplish the task. Interestingly, this pattern was also described in a study carried out by our group 23 , in which MCI and healthy subjects were compared performing the same Sternberg's letter-probe task. Network's results showed that MCI subjects exhibited an increment of strength and outreach values, and a decrement in the average shortest path, as formerly mentioned. Taken together all these findings, and especially the fact that the network outreach (and as a consequence, the energy consumption) is higher in the low CR group we may consider that brain functioning and network organization of those subjects with low CR have some similarities with MCI subjects during the execution of a memory task, and therefore we suggest that CR should be investigated as a possible factor with a remarkable impact on the cognitive status within the healthy aging. However, it is worth mentioning that, in our case, none of the normalized parameters showed significant statistical differences between low and high CR groups, while in MCI subjects showed differences in the normalized parameters indicating that their topology was more random than networks of the control group 23 . In our case, since the normalization is made by equivalent random networks, the lack of statistical differences between groups indicates that we can not claim that functional networks of the low CR group have a higher/lower random organization.
It is worth mentioning how the synchronized activity along the whole functional networks was negatively correlated with CR during the cognitive task 11 . In this sense, individuals with higher CR required less synchronization between brain regions to successfully carry out the memory test. On the contrary, fMRI recordings obtained during resting state reveal negative correlations, i.e., individuals with higher CR maintain higher synchronization in resting state 15,16 . This discrepancy suggests different mechanisms behind the creation and organization of the functional networks during resting state and cognitive tasks. Future studies should investigate this phenomenon with resting and task data with MEG.
The interpretation of the brain network profile described in the low CR group may be explained as the result of a neural compensatory effect 4 . In this case, as both groups did not differ in their memory task accuracy, we may consider it as a successful neural compensation since the low CR group had to engage additional resources to maintain or improve their performance 12 . This effect may be related to the Compensation-Related Utilization of Neural Circuits Hypothesis (CRUNCH) 35 , which proposes that there is a loss of efficiency in brain networks of aging, when compared with young adults, since additional brain areas are needed to perform cognitive tasks. In the present study, we used a low demanding memory task, but if we would had modified its difficulty, we would had probably found an unsuccessful compensation attempt in the low CR group, as previously described in MCI subjects 36 . The cause of this compensatory effect may be the beginning of a network malfunctioning since individuals with low CR exhibited a lower energetic efficiency to perform the memory task, as previously described in MCI subjects 23 . In this line, these individuals may exhibit clinical symptoms of pathological aging (i.e. MCI) earlier or more severe than those subjects with a higher CR.
Another noticeable novelty explored in the present work was the dynamical properties of the cortical brain regions while performing the memory task. In line with the microscale results, the main differences between both groups were located over the parieto-occipital lobe. In fact, a lower CR level was related with lower complexity and higher entropy. This means that the brain signals of the low CR group presented a pattern more disorganized than the high CR group within the occipital lobe. In addition, the low CR group exhibited the opposite pattern in the right temporal lobe (i.e. less entropy and more complexity). We then may hypothesize that the recruitment of posterior brain areas was not providing an efficient contribution to adequately perform the task, and therefore the individuals with low CR needed to involve right temporal regions to achieved it. We then verified that these two dynamical properties were inversely correlated, independently of the CR level. Finally, one of the most original approaches performed in this study was the exploration of the relationship between the topology (i.e. node strength and weighted clustering coefficient) and the dynamic characteristics of the nodes (i.e. entropy and complexity). We found that in both CR groups, greater values in both topological measures were related with a higher entropy and a lower complexity. These results demonstrated that the most important nodes or hubs within the functional brain networks were those that exhibited more random dynamics, no matter what the level of CR of the individuals.
The results obtained in the current study are limited due to the lack of follow-up of the participants and the small sample size, especially in the case of the high CR group. Notwithstanding, it still demonstrated that CR plays a relevant role in the topological and dynamical network's properties in healthy aging. It should be pointed that one of the advantages of conducting the analysis in the sensor space was that MEG signal was not manipulated, affecting neither the complexity nor the entropy. In addition, these results corroborated that network analysis may be a suitable approach for studying CR, and that MEG may be an appropriate tool to detect functional brain changes in healthy aging. Thus, the profiles described in this work for the low CR group indicate a dual pattern of compensation and network impairment. As a consequence, if one of these subjects suffers brain damage or starts a neurodegenerative disease process, their closeness to random network organization would probably be a risk factor for developing a more severe cognitive impairment, although some compensation mechanisms are still present. Thus, this study highlights the importance of preventive interventions, including cognitive training, which could induce a more protective brain network organization in the case of brain damage in elderly subjects.