Diurnal variations of resting-state fMRI data: A graph-based analysis

Circadian rhythms (lasting approximately 24 h) control and entrain various physiological processes, ranging from neural activity and hormone secretion to sleep cycles and eating habits. Several studies have shown that time of day (TOD) is associated with human cognition and brain functions. In this study, utilizing a chronotype-based paradigm, we applied a graph theory approach on resting-state functional MRI (rs-fMRI) data to compare whole-brain functional network topology between morning and evening sessions and between morning-type (MT) and evening-type (ET) participants. Sixty-two individuals (31 MT and 31 ET) underwent two fMRI sessions, approximately 1 hour (morning) and 10 h (evening) after their wake-up time, according to their declared habitual sleep-wake pattern on a regular working day. In the global analysis, the findings revealed the effect of TOD on functional connectivity (FC) patterns, including increased small-worldness, assortativity, and synchronization across the day. However, we identified no significant differences based on chronotype categories. The study of the modular structure of the brain at mesoscale showed that functional networks tended to be more integrated with one another in the evening session than in the morning session. Local/regional changes were affected by both factors (i.e., TOD and chronotype), mostly in areas associated with somatomotor, attention, frontoparietal, and default networks. Furthermore, connectivity and hub analyses revealed that the somatomotor, ventral attention, and visual networks covered the most highly connected areas in the morning and evening sessions: the latter two were more active in the morning sessions, and the first was identified as being more active in the evening. Finally, we performed a correlation analysis to determine whether global and nodal measures were associated with subjective assessments across participants. Collectively, these findings contribute to an increased understanding of diurnal fluctuations in resting brain activity and highlight the role of TOD in future studies on brain function and the design of fMRI experiments.


Introduction
Circadian rhythms are endogenous oscillations with a periodicity of approximately 24 h in most living organisms. They have an impact on a variety of physiological phenomena, including the sleep-wake cycle (Borbély, 1982;Dijk and Lockley, 2002;Schmidt et al., 2012), body temperature (Refinetti and Menaker, 1992), endocrine and metabolic rhythms (Hastings et al., 2007), gene expression (Storch et al., 2002), and musculoskeletal activity (Aoyama and Shibata, 2017), as well as a wide range of brain functions (Dibner et al., 2010;Schmidt et al., 2007). Studies of brain function in humans have shown that circadian variations also have an impact on a wide variety of abilities, such as attention (Valdez et al., 2005), working memory (Ramírez et al., 2006), motor (Edwards et al., 2007), and visual detection (Tassi et al., 2000). These studies investigated function using multiple scales of brain organization, from the level of individual cells and synapses (Gilestro et al., 2009;Kuhn et al., 2016;Vyazovskiy et al., 2008) to brain regions and large-scale functional connectivity (Blautzik et al., 2013;Hodkinson et al., 2014;Orban et al., 2020;Shannon et al., 2013;Steel et al., 2019).
The chronotype-based paradigm is susceptible to circadian and homeostatic rhythms, which provides a suitable way to measure the effects of sleep-wake regulation on cerebral mechanisms (Schmidt et al., 2007). A chronotype-a biologically driven circadian typology -refers to individual differences in sleep-wake cycles, diurnal preferences, and alertness throughout the day (Roenneberg et al., 2003;Susman et al., 2007). Questionnaires have reliably confirmed the differences in chronotypes (Adan and Almirall, 1991;Horne and Östberg, 1976), which have been shown to be strongly correlated with the physiological properties of circadian rhythms, including melatonin levels, core body temperature, rest/ activity cycles, midsleep point, heart rate, blood pressure, and physical activity (Adan et al., 2012;Anderson et al., 2017;Facer-Childs and Brandstaetter, 2015;Roeser et al., 2012). Genetic factors have been shown to associate with diurnal preference and homeostatic regulation of sleep. One of the most extensively studied clock genes is the primate-specific gene PERIOD3 (PER3). In humans, the variable number tandem repeat (VNTR) in the PER3 gene, consisting of either 4 or 5 repeated 54 base pair sequences encoding 18 amino acids, affects circadian typology and sleep homeostatic drive (for a review, see Archer et al., 2018).
Traditionally, individuals fall into the morning ("early larks") or evening ("night owls") chronotypes. Morning chronotypes typically have physiological and mental peaks that are shifted toward the earlier hours of the day, while evening chronotypes are more active later in the day (Bailey and Heitkemper, 2001;Kerkhof and Dongen, 1996). There is some evidence that polymorphism of the PER3 gene is associated with individual differences in circadian and sleep phenotypes. The PER3 4-repeat allele has been associated with "eveningness," whereas PER3 5/5 is linked with "morningness" and greater homeostatic sleep pressure (e.g., Archer et al., 2018;Liberman et al., 2017;Viola et al., 2007). Various studies have shown that people with different chronotypes have significantly different diurnal profiles of cognition and behavior (Horne et al., 1980;Norbury, 2020;Schmidt et al., 2007;Valdez et al., 2012). Circadian variations in performance-related neural activity have been reported in studies utilizing chronotype-based paradigms (e.g., Facer-Childs et al., 2019b;Fafrowicz et al., 2009;Gorfine et al., 2007;Peres et al., 2011;Schmidt et al., 2009Schmidt et al., , 2012Schmidt et al., , 2015Vandewalle et al., 2009Vandewalle et al., , 2011. Although extensive research has been carried out on the effects of circadian rhythms on behavior, few studies have investigated the impact of time of day (TOD) and chronotype on functional magnetic resonance imaging (fMRI) activity (Blautzik et al., 2013;Cordani et al., 2018;Fafrowicz et al., 2019;Gorfine and Zisapel, 2009;Hodkinson et al., 2014;Jiang et al., 2016;Marek et al., 2010;Peres et al., 2011;Schmidt et al., 2012;Shannon et al., 2013;Steel et al., 2019). The studies that have been performed often yielded contradictory or even ambiguous findings. Also, most fMRI studies assume that diurnal fluctuations of brain connectivity patterns and human chronotypes are relatively in-significant and unlikely to lead to a substantial systematic bias when performing group analyses (Orban et al., 2020).
Here, we tried to bridge the gap using a graph-based approach. In recent years, the application of graph theory in neuroimaging studies to analyze the human brain connectome has received much attention (Bassett and Bullmore, 2009;Sporns, 2012, 2009;Farahani et al., 2019b;Rubinov and Sporns, 2010). In terms of circadian rhythms, the few studies that have been conducted using graph theory (Anderson et al., 2017;Farahani et al., 2019a) mainly focused on task-based data and did not consider chronotypes. However, many other cognitive and behavioral studies have performed analyses of network properties that feature graph theory, including those focused on human intelligence (Hilger et al., 2017), lifetime trajectory (Finotelli et al., 2018;Gozdas et al., 2018), working memory performance (Markett et al., 2018), perception (Sadaghiani et al., 2015), and fatigue (Petruo et al., 2018). Graph theory has also been applied to the study of a wide range of neurological and psychiatric disorders, including epilepsy (Evangselisti et al., 2018), Alzheimer's disease (Hojjati et al., 2017), multiple sclerosis (Eijlers et al., 2017), autism (Sadeghi et al., 2017), and attention-deficit/hyperactivity disorder (dos Santos Siqueira et al., 2014). This paper is an extension of the work published by Farahani et al. (2021), and it utilizes the same experimental dataset, though it introduces significant methodological improvements. Here, using concepts from network neuroscience, we examined the effect of TOD on resting-state fMRI (rs-fMRI) functional connectivity while taking into account subject chronotypes. We compared the brain network properties-at different topological scales including local, meso and global-between the morning and evening sessions, as well as "early larks" and "night owls." We found greater small-worldness, assortativity, and synchronization as the waking time increased, although there was no chronotype effect. In the mesoscale analysis, we found that systems were more inclined to integrate in the evening than in the morning. The local/regional analysis revealed significant changes in both factors under study-TOD and chronotype-which were primarily associated with the somatomotor, attention, control, and default-mode networks. Furthermore, due to the homogeneity inefficiency in the parcellation we used previously (i.e., the automated anatomical labeling atlas; Tzourio-Mazoyer et al., 2002)-which may not represent the structure of resting-state functional connectivity (FC) well (Craddock et al., 2012;Gordon et al., 2016;Shen et al., 2013)-in this study, we identified areas of interest using cortical Schaefer/Yeo parcellation (Schaefer et al., 2018), in which each node is preassigned to a functional network. As a result, we noticed significant changes, particularly in the local and mesoscale analyses, that were not detected in the previous study. Within the context of our findings, it is clear that TOD may influence connectivity patterns in resting-state fMRI data, making this variable a potentially vital factor to consider in future rs-fMRI experiments.
wake-up and bedtimes). Individuals reporting excessive daytime sleepiness were excluded from the study, as determined by a cutoff ESS score of 10 or fewer points. Four hundred fifty-one participants were divided into morning or evening chronotypes. All participants underwent genotyping to identify PER3 VNTR polymorphisms in DNA isolated from buccal swabs using a DNA Gen-eMATRIX Swab-Extract DNA Purification Kit (EURx, Gdańsk, Poland) according to the manufacturer's protocol. Only individuals who were homozygous for the PER3 4 (evening-type [ET] circadian typology) and PER3 5 (morningtype [MT] circadian typology) alleles were included in the study. Other selection criteria included an age between 20 and 35 years, right-handedness-as indicated by the Edinburgh Handedness Inventory (Oldfield, 1971)-a regular TOD schedule with no sleep debt, no neurological or psychiatric disorders, no addiction, normal or corrected-to-normal vision, and no contraindications to magnetic resonance imaging (MRI). We identified 73 young, healthy participants (39 women; age: 23.97 ± 3.26 years) who met these criteria and selected them for the study. Demographic information and the results of the questionnaires are provided in Table 1.
Resting-state fMRI was performed twice, in a morning and evening session, about 1 and 10 h, respectively, after awakening from nighttime sleep. The order of the sessions was counterbalanced across the study sample. The participants were asked to maintain a regular sleep-wake schedule for one week before the study, which was monitored using Motion Watch 8 actigraphs. These actigraphs were also worn to supervise participants' sleep length and quality during the study days. Actigraphy results are provided in Table 2. The night before the morning sessions, the participants slept in the same building as the MRI scanner. All participants abstained from alcohol (48 h) and caffeine (24 h) before the MRI scanning sessions and were only allowed to engage in non-strenuous activities during study days. The study was approved by the Institute of Applied Psychology Ethics Committee of the Jagiellonian University (Krakow, Poland). Written informed consent was provided by all participants in accordance with the Declaration of Helsinki.

Data acquisition
MRI studies were conducted using a 3T Siemens Skyra MR system equipped with a 64-channel head coil. Anatomical data were acquired through a sagittal 3-dimensional T1-weighted MPRAGE sequence. Ten minutes of resting-state blood oxygenation leveldependent (BOLD) images were scanned using a gradient-echo single-shot echo-planar imaging sequence with the following parameters: repetition time (TR) = 1800 ms; echo time (TE) = 27 ms; field of view (FOV) = 256 × 256 mm 2 ; slice thickness = 4 mm; and voxel size = 4 × 4 × 4 mm 3 with no gap. A total of 34 interleaved transverse slices and 335 vol were collected for each participant. The subjects were instructed to remain awake with their eyes open and to avoid thinking of anything deliberately throughout the scanning session. The participants were monitored using an eye-tracking system to ensure that they remained awake throughout the scan (Eyelink 1000, SR Research, Mississauga, ON, Canada).
MATLAB platform (MathWorks, Inc., Natick, MA). Due to signal equilibration, the first ten volumes were discarded. This was followed by slice timing and realignment with an appraisal of voxel-specific head motion. The head motion parameters were determined for each participant; those with movements above 3 mm translation and 3° rotation were excluded from further examination. A total of four participants were excluded due to excessive head movements. Next, functional scans were registered using T1 images and normalized to the Montreal Neurological Institute (MNI) template using DARTEL (Ashburner, 2007) at a resolution of 3 × 3 × 3 mm 3 . In total, seven participants were excluded due to low-quality image registration. Functional data were spatially smoothed using a 4-mm fullwidth half maximum (FWHM) Gaussian kernel to increase the signal-tonoise ratio. The signal was band-pass filtered (0.01-0.1 Hz). Finally, the nuisance signals (24 motion parameters, cerebrospinal fluid, and white matter signals) were removed from the time course of each voxel (Behzadi et al., 2007). We did not regress out the global signal to keep any additional information (Liu et al., 2017).

Brain network construction
Our analysis pipeline is shown in Fig. 1. A large-scale brain network consists of a finite set of nodes (e.g., single neurons or anatomical brain regions) that are connected by edges (e.g., structural or functional connections between nodes). We specified the nodes by parceling the brain into seven systems/networks consisting of 200 cortical regions of interest (ROIs) from the Schaefer-Yeo atlas (Schaefer et al., 2018;Fig. 1C). In this atlas, each node is preassigned to one of the following systems/networks: visual, somatomotor, dorsal attention, salience/ventral attention, limbic, frontoparietal, or default mode. Then, the average BOLD time series across all voxels within each ROI were extracted separately (Fig. 1D). The connectivity between each pair of ROIs was then computed using Pearson's correlation coefficient. The correlation values were converted into z-values using Fisher's r-to-z transformation to improve the normality. At this stage, a symmetrical weighted connectivity matrix (adjacency matrix) with a size of 200 × 200 was constructed for each participant (Fig. 1E).
To reduce the number of spurious connections in the fully weighted matrices (Power et al., 2011), we adopted a thresholding procedure based on network density to preserve the ratios of the strongest connections and remove the weaker links (van den Heuvel et al., 2017). This procedure leads to equal network density across all participants (i.e., an equal number of edges), crucial for comparing network topology within or between participants (Gamboa et al., 2014). The sparsity threshold used in this study ranged from 0.05 to 0.5 with an interval of 0.05 to prevent the creation of either disconnected or densely connected networks (Wang et al., 2020). This step was followed by binarizing the thresholded matrices to render the computational complexity more tractable (Fig. 1F). We used the absolute value of all correlations in weighted matrices for binarization. In this study, we used both the weighted and binary matrices depending on the type of analysis (i.e., global, mesoscale, local, or hub analysis) described in the following sections. Fig. 2 displays the average weighted and binary matrices across participants in the morning and evening sessions. This figure indicates the allocation of nodes to each of the Schaefer-Yeo systems by colored rectangle patches as tick labels along axes.

Computation of graph measures
2.5.1. Global and local metrics-Using binary undirected matrices, we examined the topological features of functional brain networks for each subject across a range of cost thresholds at the global and local levels with the brain connectivity toolbox (BCT; Rubinov and Sporns, 2010) and the GRETNA toolkit (Wang et al., 2015). Table 3 provides mathematical definitions and descriptive explanations of each network statistic. Global metrics principally measure the functional segregation and integration of brain networks. Thus, we calculated global efficiency, mean clustering coefficient, characteristic path length, small-worldness, assortativity, and synchronization. Local network measures were calculated separately for each node (region) by examining the nodal centrality and density of network hubs (i.e., nodes with more than the average number of links). Hubs can be classified as either provincial or connector, which contain mostly local connections within a module or both local and long-range links that connect nodes in different modules, respectively. We calculated the most common local properties, including degree centrality, betweenness centrality, nodal clustering coefficient, nodal efficiency, and participant coefficient (Rubinov and Sporns, 2010).

Mesoscale metrics-
We utilized a multi-layer (or multi-slice) community detection algorithm (Mucha et al., 2010) to explore the modular structure of the resting-state brain networks across individuals. Each layer corresponds to an individual's functional connectivity matrix (weighted matrix). This algorithm ensures that the community assignments (labels) are preserved across layers, thus making them comparable to each other. The multi-layer modularity function Q is generally initialized with two crucial parameters, structural resolution γ and interlayer coupling ω (see Table 3 for mathematical definitions). They tune the size of modules within each layer and the number of modules across layers, respectively. Although finding the optimal ω and γ is not straightforward, strategies for achieving reasonable solutions have been proposed (Tardiff et al., 2021). In this study, we first formed a 2D discrete parameter space inspired by previous studies (γ ∈ [0.5, 1.5] with a step size of 0.05; ω ∈ [0, 1] with a step size of 0.05). Then, we performed the modularity maximization procedure for all (γ, ω) combinations in the space and selected the corresponding parameter values with the highest Q, which resulted in γ = 1.2 and ω = 0.1.
In each run, the multi-layer modularity function outputs the community labels in addition to Q, which are used to create the module allegiance matrix.
The module allegiance matrix (Figs. 4A and 4B) represents the fraction of layers in which two nodes belong to the same community . To construct an allegiance matrix, we created a co-occurrence matrix (200×200) for each layer, wherein the ij th element is equal to 1 if the nodes i and j have a shared community label, and 0 otherwise. The average of all co-occurrence matrices across layers (62 layers per session) forms the allegiance matrix; thus, its elements range from 0 to 1. Based on the module allegiance matrix, we computed two network coefficients at the mesoscale called recruitment and integration  to compare community structure between the target populations (see Table 3 for mathematical definitions).

Statistical tests
For global and local analyses, we applied a non-parametric permutation test (p-values were estimated from 30,000 permutations of group labels) in a mixed-design ANOVA to compare the mean differences between groups (considering the interaction effect) in which the within-subjects factor represented TOD, and the between-subjects factor represented chronotype (Anderson, 2001). For mesoscale and correlation analyses, we used one-way permutation tests to examine the difference between sessions (i.e., morning and evening) and detect correlations' significance, respectively. This approach does not require distributional assumptions (Nichols and Holmes, 2002). A false discovery rate (FDR) correction was applied to all statistical tests (Benjamini and Hochberg, 1995).

Global analysis
We compared FC between morning and evening sessions and found significant differences in small-worldness, network synchronization, and assortativity ( Fig. 3). No statistical evidence of change was found in any other global measures. Small-worldness (Fig. 3A) decreased with higher network sparsity in both sessions. Compared with the morning session, results from the evening session revealed higher small-worldness at a sparsity of 0.05 and 0.1 (p < 0.05, FDR corrected), which did not differ between chronotypes. Assortativity (Fig. 3B ) and network synchronization (Fig. 3C) increased with higher network sparsity in both sessions. Our analysis revealed that assortativity and synchronization were significantly higher during the evening session than the morning session at a sparsity of 0.3 to 0.5 (p < 0.05, FDR corrected), with no differences between chronotypes.

Mesoscale analysis
We used a multi-layer (multi-subject) modularity framework (Mucha et al., 2010) to compare the community structure of the functional connectivity between sessions. Although most multi-layer modularity studies are about uncovering time-varying patterns, other applications such as the study of communities across subjects or task states have been performed (Betzel et al., 2019;Zamani Esfahlani et al., 2021). So, we assumed each layer as the weighted connectivity matrix of each individual in this study. By operating the multi-layer modularity function and taking the community labels, we created a module allegiance matrix  per session (Figs. 4A and 4B). The module allegiance matrix displays how 200 brain regions and 7 Schaefer-Yeo networks/systems are cohesively engaged across individuals (Mattar et al., 2015). Then we extracted the recruitment and integration coefficients from allegiance matrices  to compare the modular structure between the morning and evening groups. See Methods on how to build these metrics.
These coefficients allow the functional interplay among brain regions and predefined/static functional systems to be quantified. Recruitment measures how a region is recruited to its own system across individuals, and integration measures the extent to which a region is integrated with other systems across individuals. Each row/column of the allegiance matrix corresponds to a brain region whose average values inside and outside of its static system yield the recruitment and integration coefficients of that region, respectively. According to Figs. 4A and 4B, the warm block-like patterns along the diagonal of each quadrant in allegiance matrices confirm that predefined systems generally tend to be recruited than integrated with other systems across individuals (Cole et al., 2014). Fig. 4C compares the recruitment coefficients of the 200 brain regions between the morning and evening sessions as a scatterplot with the linear regression fit (red line). We applied a permutation test by shuffling the group labels to see if the regression line differed significantly from the identity line. As can be seen, the recruitment results were consistent with the null (p > 0.05, FDR corrected). The recruitment coefficients are also plotted on top of the brain glass schematics in Fig. 4D for both sessions. This plot shows how brain regions differed in recruitment from their predefined/static system across participants. A similar scatterplot and brain glass for the integration coefficient are also shown in Figs. 4E and 4F, respectively. Contrary to recruitment results, these plots confirm that the difference in integration coefficients between the morning and evening sessions was inconsistent with the null (p < 0.05, FDR corrected). Notably, we found that the brain areas in the evening session were more integrated with the regions of other systems than in the morning.
In addition to regionally studying the brain as shown in Fig. 4 (i.e., thoroughly investigating the role of each region, both within its network and into other networks), a coarser granularity of brain interactions at the systemic scale across individuals could be explored by smoothing out regional information (Mattar et al., 2015). Therefore, we defined the systemlevel module allegiance matrix so that its kl th element was computed as the average of the values of all pairs of regions between systems k and l (including k = l) from the regionallevel module allegiance matrix (i.e., the mean of squares gridded with white lines). We then merged the information from both hemispheres, resulting in a 7 × 7 allegiance matrix that exhibits how regions from large-scale systems are engaged in the functional brain network (Fig. 5A). Next, using these allegiance matrices derived from both morning and evening sessions, we calculated and compared their recruitment and integration coefficients for systems instead of regions. For a given system, the recruitment coefficient is the probability that any region of that system has the same community label as the other regions within that system. Simply put, the diagonal elements of the allegiance matrix correspond to the recruitment coefficients of the large-scale systems. Also, the integration coefficient between each pair of systems is the average probability that regions in one system share the same community label as regions in another system (off-diagonal elements of the allegiance matrix). See Table 3 for detailed mathematical definitions. Fig. 5A shows that the systems differed in their strength of network recruitment and integration in both morning and evening sessions. Some systems were more consistently recruited across participants, such as the visual and somatomotor networks. At the same time, some were less recruited, such as the dorsal attention, ventral attention, and frontoparietal networks (see diagonal elements). For a more transparent illustration of the integration measure, we also created a chord diagram per session (Fig. 5B) using off-diagonal elements of the allegiance matrices, in which the edges represent the network integration between the brain systems. Among all systems in both sessions, the visual network was less integrated with other networks (also less in the morning than in the afternoon). We also found that the limbic, frontoparietal, and default networks tended to be well integrated among themselves. Finally, we observed that the somatomotor and ventral attention systems were well integrated. This graphical representation helps to understand better the complex patterns of integrations in a heterogeneous set of large-scale systems.
We also compared each system's recruitment and average integration (to all other systems) between the morning and evening sessions in Fig. 5C. Thus, for each coefficient and in each system, we performed a permutation test in which the group labels were shuffled repeatedly to check whether the difference in sessions was in line with the null or not. Significant differences are marked with an asterisk above each pair of bars (p < 0.05, FDR corrected). We found that all systems were more integrated with most other systems in the evening than in the morning. While there were significant differences for the recruitment coefficient only in the dorsal attention and frontoparietal systems, the former decreased and the latter increased during the day. Despite all these variations, the systems generally tended to be more recruited than integrated with other systems across individuals, reflecting the cohesive nature of large-scale systems in the brain (Mattar et al., 2015).

Local analysis
Table 4 summarizes brain regions that exhibited significant differences between the morning and evening sessions (first factor) and between the "lark" and "owl" participants (second factor), based on their nodal/local properties in more than half of the network sparsity. As shown in Table 4, the measures of degree, betweenness centrality, clustering coefficient, and nodal efficiency were computed for group comparisons. Most of these differences involved regions and their homotopic partners in the opposite hemisphere. No significant differences in participation coefficient and nodal shortest path were found for either factor (p > 0.05, FDR corrected). The results of the area under the curve (AUC) analysis for the degree, betweenness centrality, clustering coefficient, and nodal efficiency for all 200 brain regions are presented in Fig. 6. The AUC was calculated for each metric to provide a scalar not contingent on a particular threshold value (Wang et al., 2009;Zhang et al., 2011).
According to Table 4 and Fig. 6, compared with the morning session, the evening session showed a significantly higher nodal degree in the somatomotor network and in areas such as the bilateral superior temporal gyrus, and postcentral gyrus as well as a decreased degree of centrality in the left ventral attention network and regions such as the supramarginal and middle frontal (p < 0.05, FDR corrected). Similar results were obtained for nodal efficiency, together with a significant reduction in FC throughout the day in the left angular gyrus (p < 0.05, FDR corrected). Betweenness centrality analysis also showed a substantial increase in the evening session compared with the morning session in areas such as the bilateral precuneus as well as the right angular and right supramarginal gyri (p < 0.05, FDR corrected). Finally, the nodal clustering coefficient was higher in the evening session than in the morning session in the right superior temporal gyrus, while this value was lower in the left superior frontal gyrus and right angular gyrus later in the day (p < 0.05, FDR corrected).
In chronotype analysis (Fig. 7), degree centrality and nodal efficiency underwent a significant decrease in the bilateral dorsal anterior cingulate cortex and left insular cortex in the ET participants compared with the MT (p < 0.05, FDR corrected). Also, a comparison of the clustering coefficient and nodal efficiency characteristic of the left superior frontal gyrus showed significantly lower values in the ET than in the MT (p < 0.05, FDR corrected). Finally, the betweenness centrality in areas that included the bilateral dorsal anterior cingulate gyrus and right precentral gyrus showed significantly lower values in FC among the ET subjects compared to the results from the MT group (p < 0.05, FDR corrected).

Hub analysis
In this subsection, using a pre-determined modular classification that includes the visual, somatomotor, dorsal attention, ventral attention, limbic, frontoparietal, and default mode network (Yeo et al., 2011), we identified network hubs (connector or provincial) for the morning and evening sessions, as well as for the MT and ET groups. In this regard, we found changes in hub organization across the TOD (Fig. 8). No significant differences were detected with respect to chronotype. The results presented in Fig. 8 are based on the average connectivity matrix (across all individuals for each scanning session). For ease of visualization, we selected a network density of 0.05. As can be seen, the identified hubs nearly overlap with one another in the morning and evening sessions, except for changes in the left frontal operculum insula (LH_SalVentAttn_FrOperIns_3, 4), right superior parietal gyrus (RH_DorsAttn_Post_5), right precentral ventral gyrus (RH_DorsAttn_PrCv_1), left precuneus (LH_Cont_pCun_1), and bilateral posterior cingulate cortex/precuneus (Default_pCunPCC_1, 2). We found that the somatomotor network contains more hubs than any of the other networks in both morning and evening sessions. Notably, its hubs are both provincial (i.e., within modular connections) and connector (i.e., between modular connections), at nearly the same ratio. In contrast, hubs in the ventral attention network are connector, while hubs identified in the visual network were mostly provincial. To examine connections among all regions and to identify hub types (i.e., connector or provincial), connectograms of both sessions were created using Circos software (Krzywinski et al., 2009). The results are illustrated in Fig. 9.

Correlation analysis
A correlation analysis was performed to determine whether global and nodal measures throughout the day were significantly associated with variables of interest (e.g., ME scale, AM scale, and ESS) across participants while controlling for the differences among the covariates of no interest (e.g., age, sex, and clinical variables). The ME, AM, and ESS scores represent an individual's chronotype preference, the strength of this preference, and the degree of sleepiness during the day, respectively. Overall, the number of significant associations was greater in the morning session than in the evening session. From a global perspective (Table 5 and Fig. 10), correlation analysis revealed significant negative associations between AM scores and both small-worldness and modularity in the morning session (p < 0.05, FDR corrected). We also found significant positive correlations between ESS and average path length and assortativity in the morning session and positive correlations between AM score and path length and assortativity in the evening session.
From the nodal perspective, we found significant correlations between the degree centrality of various brain regions in both hemispheres and these subjective indicators (ME scale, AM scale, and ESS), mostly across the morning scanning session (Table 6 and Fig.  11). In the morning session, significant negative correlations were found between AM scores and areas within the default network, including the left rostral anterior cingulate gyrus (Default_PFC_6), left precuneus (Default_pCunPCC_2, 4), right medial prefrontal cortex (Default_PFCm_4), and right posterior cingulate cortex (Default_pCunPCC_2); between ESS and the left pole of the superior temporal gyrus (Limbic_TempPole_3), the right lateral fronto-orbital gyrus (Cont_PFCl_1), and the right pole of the middle temporal gyrus (Default_Temp_1); and between ME scores and the left postcentral gyrus (SomMot_4) and the left pole of the superior temporal gyrus (Limbic_TempPole_4). Significant positive associations were found between AM scores and the bilateral precentral gyrus (left DorsAttn_FEF_1 and right SomMot_11) and between ME scores and the bilateral insula (SalVentAttn_FrOper_2 and Cont_PFCv_1), the left anterior cingulate gyrus (Default_PFC_8), and the right precentral gyrus (DorsAttn_FEF_1).
During the evening session, we identified significant positive associations between AM scores and the left middle occipital gyrus (Vis_11), right fusiform gyrus (Vis_2), and right superior occipital gyrus (Vis_14), as well as negative correlations between AM scores and the left precentral gyrus (DorsAttn_PrCv_1), bilateral lateral fronto-orbital gyrus (Limbic_OFC_1, 2), and right middle temporal gyrus (Default_Temp_4). No significant correlations were found between ESS and ME scores and the degree centrality of these brain regions.

Discussion
The current study, used rs-fMRI, a chronotype-based paradigm, and graph theory to examine the diurnal fluctuations of whole-brain connectivity architecture in 62 young, healthy participants. The study revealed meaningful information regarding the topological variations of the brain network during the day and organizational differences between the "lark" and "owl" groups, as well as associations of graph theory metrics with selected variables of interest (i.e., ME, AM, and ESS scores).
The main results can be summarized as follows: (1) Among the global measures, there was a significant increase in small-worldness, assortativity, and network synchronization in the evening session over the morning session (p < 0.05, FDR corrected). However, there was no compelling evidence of changes in any of the global metrics in chronotype (i.e., between MT and ET participants). (2) Mesoscale analysis showed higher brain systems/regions integration among themselves in the evening session. (3) Local graph measures varied during the day and between the two chronotypes, predominantly across the somatomotor, attention, and default-mode networks. (4) Analysis of the hubs revealed that the somatomotor network was the densest area of the brain in both sessions, but more so during the evening session, including both provincial and connector types, whereas hubs in the ventral attention network and visual network were primarily connector and provincial, respectively. (5) Correlation analysis revealed significant associations between the variables derived from the questionnaires (ME, AM, and ESS) and the nodal characteristics of several brain regions in both scanning sessions, most of which were associated with the morning session.

Diurnal variations in the brain network as a whole (global properties)
A small-world network is an intermediary between a random and a regular grid that contains many short-range connections alongside a few long-range shortcuts (Watts and Strogatz, 1998). There are other ways to define a small-world network (e.g., by considering the physical length of connections), all of which suggest mathematically that small-world networks share relatively high transitivity and small mean geodesic distance (i.e., shortest path) between nodes. This property strikes an optimal balance between network integration and segregation (Bassett and Bullmore, 2006;Rubinov and Sporns, 2010). In our analysis of rs-fMRI data, we found high values of small-worldness for both scanning sessions (small-world networks generally have an σ of at least 1), albeit with significant superiority at highly sparse networks in the evening compared with the morning session. In another study, Anderson et al. (2017) explored how TOD affects functional brain networks in older healthy adults, and they found no topological changes in small-worldness during the resting state. It is noteworthy that the significant changes that we found occurred only at lower densities, and-as in Anderson et al. (2017), who considered a sparsity range of 0.1 to 0.7-we found no significant difference at rest during the day at densities above 0.1.
However, the small-worldness provides little information on the actual network organization, and comparing its values is not straightforward because this measure does not follow a linear relationship. A review study (Mark D. Humphries and Gurney, 2008) shows that when the wiring structure of a regular network is gradually moving towards an utterly random architecture, the small-worldness, σ, first increases to reach a peak, and then decreases. The nonlinear nature of the small-worldness complicates the assessment of two given networks as they may have similar values falling on either side of the peak, while their topological structure is different (i.e., one side could be more like a regular network while the other side could be more like a random network). To further elucidate the small-worldness, we performed a system-level analysis to see how communities/modules were organized in the network for both sessions (see mesoscale analysis in Results). Therefore, we found that the daily increase in small-worldness, primarily at lower densities (global analysis), was due to the higher tendency of large-scale networks to be integrated with one another in the evening session (mesoscale analysis). Collectively, these findings reflect a more efficient topology of the information flow-due to the slight addition of randomness-in the evening session than in the morning session. Furthermore, our findings revealed an increase in network assortativity over the course of the day, which is a finding that overlaps with our previous results (Farahani et al., 2019a). Increased assortativity is related to a higher propensity for a node to connect to other nodes with the same or a similar degree (Newman, 2003;Foster et al., 2010), thereby increasing the likelihood that a nearby hub will be capable of supporting a faulty node. Finally, the results of network synchronization, which is a measurement used to assess how well all nodes oscillate in the same wave pattern, were not consistent with the findings reported by Barahona and Pecora (2002). They showed that in networks of low redundancy, small-worldness results in higher synchrony than what is found in standard deterministic graphs, random graphs, and ideal constructive schemes. We discovered that small-worldness decreased and synchronization increased in morning and evening sessions as we moved from low to high sparsity. Notably, network synchronization was significantly higher during the evening scanning session than in the morning session. Overall, the findings of all global measures indicate that brain network organization varied throughout the day, a determination that might be associated with increased brain function and interaction from morning to evening.

Diurnal variations in the community structures (mesoscale properties)
Higher integration of DMN, limbic system, and frontoparietal regions in the evening session, revealed in the current study, has been previously proven to be associated with deliberate (not spontaneous) mind-wandering (Golchert et al., 2017). The above state of mind is thought to reflect cognitive control. Other studies indicate that individuals with higher control are better at modulating mind-wandering under task demands, presumably making it more efficient (Smallwood and Andrews-Hanna, 2013;Smallwood and Schooler, 2015). Interestingly, according to our different results, the stronger the chronotype, the more integrated neural networks are in the evening. The above results suggest that those evening hours benefit people with extreme chronotypes. The high predominance of spontaneous over deliberate mind-wandering is usually related to attention deficits (Seli et al., 2015b) and higher reactivity to inner experiences (Seli et al., 2015a). Moreover, the frontoparietal control network has been proven to include two separate subsystems where FPNa is interconnected with DMN and FPNb with DAN. As a result, the former part is highly associated with introspective processes, while the latter is related to the overall processing of the external stimuli (Dixon et al., 2018). Furthermore, somatomotor and ventral attention networks, which also showed higher integration in our study, are believed to be associated with sustained attention (Mitko et al., 2019). Overall, our mesoscale findings on network integration are congruent with other results in the study, pointing out a more prominent role of the frontoparietal regions in executive control during the evening hours (Dixon et al., 2018).
We also observed significant variation in the strength of network recruitment associated with the time of the day. The most striking example was the increased recruitment of frontoparietal systems during the course of the day. The above structure is believed to be associated with cognitive control and subsequent working-memory storage (Dormal et al., 2012); therefore, results from the current study may be related to the accumulation of information during the day. Notably, FPNb, which is interconnected with DAN, is associated with processing abstraction, monitoring, and manipulating sensorimotor contingencies (Dixon et al., 2018). Current results revealing increased recruitment in the FPN and decreased recruitment in the DAN are fascinating and deserve further investigation in the future.

Nodal/local changes affected by time of day
In this study, the somatomotor, attention, and default mode networks experienced the largest quantity of topological variations among brain networks during the day. Moreover, the study of hubs in the exact same networks was proved to host several densely connected nodes, which affect brain functional integration and segregation.
In the somatomotor network, resting-state findings from both morning and evening sessions revealed higher within-network connectivity than other brain networks and higher betweennetwork connectivity, particularly with the ventral attention network. Moreover, both connectivity and hub analyses indicated that the somatomotor network contained the most highly connected areas in the brain, mainly during the evening session. We also found significantly higher FC in the bilateral superior temporal and left postcentral gyri as waking time increased. Consistent with our results, similar changes in the neural response across TOD were previously reported in other rs-fMRI studies (Anderson et al., 2014;Fafrowicz et al., 2019;Jiang et al., 2016;Maire et al., 2018), as well as in magnetoencephalography (MEG) studies evaluating oscillatory activity at rest and during a finger-tapping task (Wilson et al., 2014). On the contrary, conflicting results were presented in a morphometric study by Trefler et al. (2016), who discovered a significant decrease in cortical thickness as a function of TOD across the lateral surfaces of the left frontal, temporal, and parietal lobes. Nevertheless, findings from the current study provide insight into how changing functional activity in sensorimotor networks is associated with the course of the day. The increase of FC in somatomotor regions indicates that neural synchronization is enhanced in these areas.
Ventral and dorsal attention networks are believed to be involved in stimulus-driven and goal-directed attention, respectively (Vossel et al., 2013). Both connectivity and hub analyses revealed the ventral attention system to be the second most densely connected network, after the somatomotor areas. The above results were primarily observed during the morning session, focusing on the frontal operculum insula and dorsal anterior cingulate gyrus. Moreover, the study confirmed decreasing FC throughout the day within ventral areas such as the left supramarginal gyrus and left middle frontal gyrus (dorsal prefrontal cortex) -a finding that is in line with previous studies (Anderson et al., 2017;Jiang et al., 2016;Vandewalle et al., 2009). Other changes were also observed in the right dorsal attention areas, such as the angular and supramarginal gyri.
In addition, the current study revealed the default mode network (DMN) to be involved in relatively high neural activity at the resting state during both sessions, particularly in the posterior cingulate cortex, where it was identified as a hub. However, this activity was higher in the evening than in the morning session. The DMN consists of functionally-connected and specialized neural units, contributing to many cognitive functions (Andrews-Hanna et al., 2007;Mayer et al., 2010;Raichle et al., 2001). Moreover, the posterior cingulate cortex plays a crucial role in meditating intrinsic activity through the DMN (Fransson and Marrelec, 2008). Notably, some areas within the DMN were found to be prone to variations in their connectivity profiles across the day, such as the bilateral posterior cingulate cortex, precuneus, angular gyrus, superior temporal gyrus, and left superior frontal gyrus, a finding that is consistent with the results from previous rs-fMRI studies (Facer-Childs et al., 2019a;Fransson, 2005;Jiang et al., 2016;Ku et al., 2018;Lunsford-Avery et al., 2020;Orban et al., 2020;Raichle et al., 2001;Shannon et al., 2013).
Further studies reported a rhythmic FC pattern of the DMN during the day that peaked in the morning and declined during the afternoon (Blautzik et al., 2013;Hodkinson et al., 2014). It is believed that diurnal changes in these areas indicate the functional coordination of spatially disparate gyri (Jiang et al., 2016). Considering the findings in this study, total DMN neural activity decreased during the day; however, the FC increased in hub regions such as the posterior cingulate cortex and precuneus. Both the posterior cingulate cortex and precuneus are believed to be highly affected by decreased consciousness, which could explain the lower FC in the morning hours (Luppi et al., 2019). However, the results from the whole network are consistent with reports on sleep inertia, which is associated with increased DMN functional connectivity shortly after awakening (Vallat et al., 2019). This finding may suggest compensatory mechanisms of the mentioned areas to balance neuronal interactions, such as coupling or decoupling within DMN subregions. It could perhaps redouble efforts to increase the adaptability of the network under the continued wakeful condition throughout a day.

Nodal/local changes affected by chronotype differences
The current findings prove that the frontoparietal and dorsal attention networks underwent the most topological changes between earlier and later chronotypes. However, hub analysis did not show a statistically significant difference between the two groups. In the present study, participants with ET chronotypes revealed less neural activity and less network integration of the frontoparietal network (or central executive network) than the MT group, particularly in the dACC parcel. Consistent with our results, Horne and Norbury (2018) reported a significant reduction in FC between the dACC and amygdala in the latter chronotype, which they believed led to impaired emotion regulation. Overall, these findings could explain the role of chronotype in the interaction between the alerting functions and executive control networks (Martínez-Pérez et al., 2020).

Correlation between the network properties and subjective variables
The strength of the chronotype preferences was negatively correlated with small-worldness and modularity in the morning session. The graph measures presented are believed to be responsible for the integration of large-scale brain activities (Chavez et al., 2010). Furthermore, small-worldness has been reported to be characteristic of a healthy brainbecause lower levels signify abnormal brain functioning (Brier et al., 2014;Liu et al., 2008). Moreover, strong chronotype preference was positively correlated with path length and assortativity in the evening session. Higher assortativity is another sign of a well-functioning network, confirming that individuals with a stronger chronotype preference present higher integration of the neuronal networks in the evenings compared to the mornings. As a result, the more extreme the chronotype, the more integrated the neuronal network in the evening, which was a finding that was not dependent on the chronotype itself. Sleepiness during the day was positively correlated with the average path length and assortativity during the morning hours.
Local graph measures revealed significant differences in resting-state activity associated with the preferred chronotype. People with the ET chronotype revealed less centrality in the somatomotor and limbic networks during the morning hours, which is believed to be associated with less efficiency. These results also confirm previous findings concerning the effect of sleep inertia on later chronotypes (Ritchie et al., 2017).
Furthermore, declared sleepiness during the day negatively correlated with the centrality of the DMN structure during the morning session. Interestingly, according to Tian et al. (2020), the DMN mediates the association between chronotype and sleep quality. The same study displayed higher precuneus and medial prefrontal cortex connectivity in late chronotypes. Moreover, during the evening session, the strength of the chronotype preferences positively correlated with the degree of centrality in several regions within the visual network, such as the middle and superior occipital gyrus and fusiform gyrus. These regions are believed to be associated with visual attention.

Limitations and future directions
Certain limitations associated with this study should be considered in future research. Firstly, the relatively small number of subjects in this study might constrain the translational value of our results. Future studies with larger samples are needed to confirm our findings and increase the reproducibility of the research.
Secondly, brain nodes were derived from the cortical Schaefer-Yeo atlas (200-parcel/ 7-network parcellation; Schaefer et al., 2018). Schaefer's parcellations are available at multiple resolutions (100 to 1000 parcels). Further studies on appraising network topology using finer parcellation schemes are warranted. Also, because the Schaefer atlas considers only cortical areas, whole-brain studies could be carried out by adding subcortical regions with other atlases or segmentation algorithms. Although each node in this atlas is preassigned to a system/network (which is an advantage over many atlases), the Schaefer system labels do not allow for individual variation in the topography of brain systems; that is, they force everyone to have the same systems, which may not be appropriate and it is an issue that requires more attention (Gordon et al., 2017;Kong et al., 2019).
A third limitation concerns the measures used to compare the group-representative functional brain networks with one another. These measures tend to be correlated with one another; for example, a brain network with high efficiency must necessarily have a shorter path length (Betzel et al., 2018). Thus, if we find significant differences in one measure, they will probably be found in others. Our analyses could therefore be extended in future work to determine which of these measures might be driving the others and how an exhaustive set of metrics might be designed to fit the study specification from a neurobiological perspective. Another issue about network measures is that global statistics are often non-specific (i.e., they are not entirely informative and revealing). For example, the meaning of the phrase "the patient group has a lower efficiency than the control" may not be apparent to a neurosurgeon. Future work should be directed to better interpretation of such metrics.
Another limitation concerns the applicability of the small-world property. In real systems, the early definitions of small-worldness initiated by Watts and Strogatz (1998) are ineffective because they confuse regular networks with small-world structures and neglect the weight and physical length of connections and the network density. Most of the definitions present the network on the border of a circle; however, real systems are not embedded in this way. There are many ways for a network to be small-world other than starting from a regular grid and adding random links to reduce path length. Numerous researchers have addressed these constraints by introducing several practical metrics (Bolaños et al., 2013;Muldoon et al., 2016;Rubinov and Sporns, 2010;Telesford et al., 2011). Applying these modified metrics in future work would bring the study of the smallworld brain closer to reality.
Yet another possible extension of this work involves studying the recent theory-driven techniques that emphasize the importance of machine learning, algorithmic optimization, and parallel computing in functional neuroimaging (Cohen et al., 2017;Douglas et al., 2013). For example, various algorithms known as graph neural networks, including graph convolutional networks (GCNs), have been proposed to show how graph theory can be used to train deep learning models (Kipf and Welling, 2016;Wu et al., 2020) and to discover neurological biomarkers using fMRI data . As another example, a growing trend has developed in a family of algorithms known as hyperalignment (or functional alignment) that permit a projection of individuals' data into a shared space across participants based on how voxels respond to stimuli (Guntupalli et al., 2016;Haxby et al., 2011) or how they are connected to other voxels (Guntupalli et al., 2018;Haxby et al., 2020). Combining these techniques with network neuroscience will open a new generation of studies to transform our knowledge of neural representations in complex brain networks.

Conclusion
This study, presents evidence for topological changes in functional brain networks throughout the day (TOD effect; morning and evening sessions) using rs-fMRI data and graph theory analysis. We also consider inter-individual differences in diurnal preferences (chronotype effect; "lark" or "owl" types) in addition to the impact of TOD. In summary, the results from the global examination showed more efficient functional topology in the evening session, regardless of the chronotype. Moreover, the mesoscale results represented how different systems/regions interacted with one another at both sessions, providing an intuitive assessment of modular organization using measures of recruitment and integration. To counterbalance, the local analysis revealed chronotype-specific modulation of diurnal fluctuation prominently across the somatomotor, ventral attention, and default-mode networks. These findings provide insight into diurnal variations in resting-brain networks, reflecting the universal effect of TOD on neural functional architecture when designing experiments. The findings also indicate the need to control for circadian typology, which could influence experimental results in neuroimaging studies.

Funding statement
This study was funded in part by the Polish National Science centre through a Harmonia grant (no. 2013/08/M/HS6/00,042) and supported by the Foundation for Polish Science project, "Bio-inspired Artificial Neural Networks" (POIR.04.04.00-00-14DE/18-00). Schematic representation of our graph-based analysis. After preprocessing (B) the raw rs-fMRI data (A) and parcellating the brain into 200 regions of interest using Schaefer-Yeo atlas (C), corresponding time courses were extracted from each region (D) to compute the weighted correlation matrix (E). To reduce the complexity, a binary correlation matrix (F) and the corresponding functional brain network (G) were constructed. A set of global and local graph theory measures were then derived from these connectivity matrices (H). Finally, non-parametric statistics were applied to identify significant group means and correlations (I).  Recruitment and integration coefficients across brain regions. Panels (A) and (B) show the module allegiance matrices at the regional level for the morning and evening sessions, respectively, each representing the probability that two brain regions are assigned to the same community across individuals. The regions are arranged according to their belonging to 7 Schaefer-Yeo networks and for each hemisphere separately. In panels (C) and (D), we compared the recruitment coefficient between the sessions across brain regions, which are displayed by a scatterplot with the linear fit (p > 0.05, FDR corrected) and brain glasses, respectively. Similar plots, as for the recruitment, are shown for the integration coefficient in panels (E) and (F), where there is a significant increase in the evening compared to the morning (p < 0.05, FDR corrected). In general, as shown in (D) and (F), neighboring areas tend to have similar recruitment and integration coefficients, confirming the presence of cohesive large-scale structures in the brain. Area under the curve for the "lark" (morning-type) and "owl" (evening-type) participants for degree centrality (A), betweenness centrality (B), clustering coefficient (C), and nodal efficiency (D) in all 200 brain regions of interest. Each node in either the left or right hemisphere is labeled with a color that is matched to Schaefer-Yeo 7 network parcellation. Significant diurnal changes are represented by dashed red lines. Hubs (highly connected regions) in the morning (A) and evening (B) sessions were determined at a sparsity of 0.05. Regional information is provided in Appendix  The mean connectogram across all participants in the morning (A) and evening (B) sessions at a thresholding value of 0.05. Parcellated elements within the outermost circle include the 200 Schaefer/Yeo brain areas marked with a unique RGB code that has been associated with one of the predefined modules in each hemisphere. The outer circle circumscribes a set of five inner circular heatmaps that were created to represent the values of five different centrality measures. The range of colors for each metric represents the minimum to the maximum values. Toward the center, these measures are degree centrality, participation coefficient, K-coreness centrality, eigenvector centrality, and PageRank. The values of all measurements and the functional connections in each connectogram are derived from the mean of all participants in the corresponding session. The red and black curves indicate the functional connections between and within modules, respectively. An unambiguous abbreviation scheme was created to label each parcellation, as summarized in Appendix   Correlation analysis between the nodal centrality of brain regions and the questionnaire variables: AM (A), ESS (B), and ME (C). Nodes are colored according to the magnitude of the correlation.   Clustering coefficient: the extent to which the neighbors of a given node are interconnected (i.e., the fraction of triangles around a node)

Appendix
Clustering coefficient of the network (Watts and Strogatz, 1998): Efficiency: the efficiency of parallel information transfer of a given node determined as the average of the reciprocal shortest path length from a node to all other nodes Efficiency for a given node i (Latora and Marchiori, 2001): Betweenness centrality: the ratio of all shortest paths in the graph that contain a given node Betweenness centrality of node i (Freeman, 1978): Assortativity: the extent to which a network can resist failures in its main components. If r ≥ 0, the nodes with a high degree are more likely to connect to others that are similar in degree (an assortative network), while r < 0 reflects a tendency for high-degree nodes to attach to nodes with a low degree (a disassortative network).
Assortativity coefficient of the network (Newman, 2002 Table 4 List of brain regions of interest (ROIs) that differed significantly between morning and evening sessions (t: time) and between the "lark" and "owl" participants (c: chronotype). P-values were computed using 30,000 permutations followed by Benjamini-Hochberg correction in a two-way ANOVA; FDR was set to 0.05 and the asterisks indicate statistically significant corrected p-values).   Table 5 Global analysis: significant correlations between global metrics and AM, ESS, and ME scores (n = 62) for the morning and evening sessions. The significance level was set at p < 0.05 with non-parametric permutations.  Table 6 Nodal analysis: significant correlations between degree centrality and AM, ESS, and ME scores (n = 62) for the morning and evening sessions. The significance level was set at p < 0.05 with non-parametric permutations.