Participants

Fourteen healthy male volunteers were included in this study (mean age = 23.8, range = 20–27). Four additional participants were recruited for this study, but had to be excluded from analysis because of non-completion of the entire study protocol or non-correctable imaging artifacts. To enable a quantification of potential habituation effects, participants came in on two days, and were assessed multiple times throughout the day. Because of this resource- and cost-intense study design, sample size was limited. No participant had a history of neuropsychiatric or vascular disease, took medication or psychoactive substances, had recent or actual sleep disturbances, worked shifts, or consumed caffeine excessively (>3 cups of coffee or other caffeinated beverages per day). We restricted this study to male participants because of reported gender variations in resting-state dynamics36 and possible interactions between the menstrual cycle and circadian rhythms37. All participants gave their informed consent and received payment for their participation. This study has been approved by the ethics commission of the Medical Faculty of the Goethe University Frankfurt (GZ 244/09) in accordance with the Declaration of Helsinki.

We instructed all participants to follow regular bedtimes and wake-up-times (between 23:00–00:00 h and 7:00–8:00 h, respectively) for five days before scanning until completion of the experimental phase. Compliance was monitored with a diary of activities and by continuous wrist-actigraphy recordings38. The actigraph (Actiwatch Mini; CamNtech Ltd., Papworth Everard, UK) was worn at all times on the non-dominant wrist except during scanning sessions. Individual differences in the phase of entrainment – i.e., the chronotype39 – and increasing homeostatic sleep pressure during time spent awake, but also increasing sleep debt as a marker of sleep deprivation, are known to modulate the ToD-dependent profile of neurobehavioral variables[7].Therefore, the participants’ individual chronotype was assessed with the Munich Chronotype Questionnaire and characterized with a single value, the midpoint of sleep on free days, corrected for sleep debt accumulated during work days (MSFsc)39,40,41. We a priori excluded extreme chronotypes. The mean ± SD MSFsc of the participants was 4.87 ± 0.92 (range = 3.89–6.71). The MSFsc could not be assessed in two participants due to non-completion of the questionnaire.

Individual wake-up time on scan days was extracted from actigraphy and used as an indicator for sleep pressure in each participant at each ToD38,42. Since participants were always scanned at the same times, an earlier wake-up time indicates higher sleep pressure at each ToD compared to later wake-up times. Mean ± SD wake-up time of the participants was 06:10 h ± 00:39 h (range = 05:00–07:10 h). Wake-up time could not be calculated in one participant, because actigraphy could not be obtained on scan days. We estimated the participants’ sleep debt, indexed by the difference between mean weekly sleep duration from the Munich Chronotype Questionnaire and sleep duration during scanning days extracted from actigraphy40,41,42. Because of non-completion of the questionnaire and missing actigraphy data (see above), sleep debt could not be assessed in three participants. Mean ± SD sleep debt of the participants was 1.81 ± 1.00 (range = 0.17–3.09). Higher values indicate higher sleep debt. Because actigraphy-based sleep estimates often under-estimate the duration of sleep, as compared to self-reports43, this quantification possibly overestimates the true extent of sleep pressure in the population.

Before each fMRI session, body temperature was measured with an ear thermometer and the participants' sleep propensity was assessed with a self-report questionnaire based on a modified version of the Epworth Sleepiness Scale (ESS) questionnaire[44], in which participants were asked to rate their sleep propensity by projecting their actual vigilance state in the situations provided by the questionnaire to quantify diurnal sleepiness. Participants did not nap or consume stimulants during scanning days. Time in this paper is expressed in Central European Time (CET), the German local time. Times of twilight were obtained from the German Meteorological Office.

Resting-state fMRI

The participants underwent six fMRI scanning sessions at fixed times on two subsequent days (8:00, 11:00, 14:00, 17:00, 20:00, and 23:00 h on each day) to account for habituation effects to the experimental setting. To minimize interferences with the participants’ sleep, no measurements were scheduled between 24:00 and 07:00 h. Each fMRI session included a 7 min resting-state measurement in dim light. The participants were instructed to lie still and comfortably in a supine position, to keep their eyes open and fixate a white crosshair on a black background (< 1 lux, maximum illuminance measured at corneal level using a VOLTCRAFT LX-1108 lux meter; Voltcraft, Hirschau, Germany), not to fall asleep, and not to think of anything in particular and let their mind wander freely. A minimum of five minutes elapsed between switching off the main lights in the scanning room and the beginning of the functional measurement, allowing participants to adapt to the dim light conditions. Participants wore headphones and earplugs to protect against the scanner noise, and their heads were immobilized with foam cushions to minimize movement artifacts. If visual correction was necessary, participants were provided with appropriate MRI compatible corrective lenses. A coil-mounted mirror allowed viewing the projection screen and the white crosshair on a black background was displayed using Presentation software (Neurobehavioral Systems Inc., Berkely, CA, USA; http://www.neurobs.com).

During fMRI scans, we monitored the participants visually and recorded their cardiorespiratory and movement parameters. Prior to the experimental phase, participants underwent an ‘acclimatization session’ in the scanner to acclimatize to the unfamiliar environment and experimental procedures.

Visual detection task fMRI

Nine of the original participants additionally performed a visual detection task in dim light while still being scanned for additional 4 min 52 s at each ToD. They were instructed to press a button with their right index finger as soon and as quickly as possible when they saw a low-contrast orange crosshair (< 1 lux) flashing shortly (500 ms) on the center of a black screen inside the scanner (background light in the dimly lit scanner room was <0.1 lux). Three one minute task blocks with 11 crosshair presentations each (randomized interstimulus interval from 2–10 s) alternated with four 25 s fixation blocks (white crosshair, <1 lux, as during resting-state measurements).

Data acquisition and preprocessing of resting-state fMRI data

Resting-state data were collected using a 3 Tesla MRI-scanner (Siemens MAGNETOM Allegra, syngo MR 2004A, Erlangen, Germany). A total of 12 resting-state sessions of a gradient-echo T2*-weighted transverse echo-planar imaging sequence (EPI) were acquired for each participant. A resting-state session comprised 210 EPI volumes (6 sessions on two subsequent days = 12 × 210 volumes for each participant). Each volume contained 30 axial slices acquired in a sequential manner, covering the whole brain (TR/TE/flip angle = 2000ms/30 ms/90°, FOV = 192 mm, matrix size (resolution) = 64 × 64, voxel size of 3 × 3 × 3 mm³, distance factor 25%). The first 4 volumes were discarded to avoid magnetic saturation effects. Spatial preprocessing was performed using standard algorithms implemented in SPM 12 (Statistical Parametrical Mapping 12, Version 6225; Wellcome Trust Center for Neuroimaging, London, UK; http://www.fil.ion.ucl.ac.uk/spm/): head-motion correction via intra-subject spatial realignment (least squares approach and 6 parameter rigid body spatial transformation), spatial normalization to the standard EPI template of the Montreal Neurological Institute (MNI) and resampling to a voxel size of 2 × 2 × 2 mm³ (using a 12 parameter affine transformation), and spatial smoothing with an isotropic 8 mm full-width at half maximum (FWHM) Gaussian kernel45.

Resting-state BOLD data are contaminated by non-neural signals originating from residual motion artifacts, respiration, cardiac beats, and scanner drifts, which makes additional preprocessing steps necessary46. Using in-house scripts in MATLAB 7.7 (The MathWorks, Inc., Natick, MA, USA; https://www.mathworks.com/), data were first high pass filtered with a 0.01 Hz cutoff using a 6th order Butterworth filter. Then, in order to remove the global effect of the various noise sources, the six motion parameters computed in the realignment step and fMRI time series from cerebrospinal fluid (CSF) and white matter (WM), as well as the first derivatives of each, were detrended and regressed out of the resting-state data using least squares linear regression. The CSF and WM time series (CSF: averaged within a 2 mm sphere at MNI coordinates −8, −8, 21, WM: averaged within a 5 mm sphere at MNI coordinates 30, −20, 30) were extracted using the MarsBaR toolbox (Version 0.44; http://marsbar.sourceforge.net). Finally, a mask based on the Automated Anatomical Labeling (AAL) atlas47 was applied to all volumes to restrict further analysis to gray matter-voxels only.

A high-resolution T1-weighted 3D magnetization-prepared rapid acquisition gradient echo (3D MP-RAGE) sequence was also acquired for each participant (144 sagittal slices, 1 slab, TR/TE/TI/flip angle = 2250 ms/2.6 ms/900 ms/9°, FOV = 256 mm, matrix size (resolution) = 224 × 256, voxel size 1 × 1 × 1.10 mm³, distance factor 50%). The high-resolution anatomical T1 sequence was used to exclude any structural anomalies and to map the functional results onto an averaged brain with higher anatomical resolution (see Supplementary Fig. 1).

Single subject analyses of resting-state fMRI data

Three variables of interest were calculated from the resting-state BOLD time series for each session (ToD) separately, using in-house scripts in MATLAB 7.7. Whole brain, voxel-wise standard deviation (SD) of the BOLD signal: We calculated the voxel-wise SD of the BOLD time series for each ToD, resulting in an SD-map for each ToD. We obtained voxel-wise SD with:

$$SD = \sqrt {\frac{{\mathop {\sum }

olimits_{i = 1}^n \left( {x_i - \bar x} \right)^2}}{{n - 1}}}$$

{x 1 , x 2 , …, x n } represents the single values of a time series of length n and \(\bar x\) is the mean value of the whole time series.

Whole brain, voxel-wise amplitude of low-frequency fluctuation (ALFF) of the BOLD signal: The ALFF of the BOLD signal represents the power spectral density within the frequency band of 0.01–0.1 Hz48,49. Since the power spectrum of a time series describes how the variance of the data is distributed over the frequency domain, BOLD ALFF is analogous to BOLD SD in the given frequency range. We first transformed the BOLD time series into the frequency domain with the fast Fourier transform algorithm and then obtained BOLD ALFF by averaging the resulting power spectrum across the frequency band of interest50. Functional connectivity of the thalami: BOLD time series from the left and right thalamus, as defined by the AAL atlas, were extracted and averaged within each cluster for each ToD. The time series was correlated (Pearson’s product-moment correlation coefficient) with each voxel separately, resulting in a functional connectivity map of correlation coefficients for each ToD18.

Group analysis of resting-state fMRI data

Group statistical analyses were performed using SPM 12. First, for each participant, the BOLD SD, ALFF, and thalamic connectivity values acquired on the two consecutive days were averaged voxel-wise at each ToD separately, using the ImCalc tool in SPM, yielding BOLD SD, ALFF and thalamic connectivity whole brain maps for each of the six time points of interest. Then we used whole-brain, random-effects, one-way repeated measures ANOVAs (ToD as independent variable with 6 levels, homoscedasticity and sphericity were assumed) for group inference and tested for a main effect of ToD with an F-test, for BOLD SD, ALFF and thalamic connectivity values, respectively. All statistical brain maps were corrected for multiple comparisons by using cluster-extent based thresholding at P < 0.05, FWE. The cluster-defining threshold was set at P < 0.001 (uncorrected). We also report the MNI coordinates of local maxima in these clusters. The denomination of local maxima was adopted from the SPM Anatomy toolbox51. To identify the ToD-differences driving the main effect of ToD in BOLD SD, we performed additional post-hoc t-tests. First, BOLD SD values from significant clusters were averaged within each cluster and extracted for each ToD and participant separately using MarsBaR. Then, using IBM SPSS Statistics (Version 22; IBM, Armonk, NY, USA; https://www.ibm.com/analytics/us/en/technology/spss/), we performed a repeated measures ANOVA with subsequent pairwise comparisons between all ToD (dependent, paired samples, two-tailed t-tests with Bonferroni correction, significance threshold at P < 0.05) for each cluster. Clusters in right, left and mesial somatosensory regions were considered as one. BOLD ALFF values from significant clusters were extracted in the same way as for BOLD SD and used for planned pairwise comparisons based on the results of the analysis of BOLD SD (difference between 08:00 and 20:00 h and, midday measurements 11:00, 14:00, and 17:00 h). Again, clusters in right, left and mesial somatosensory cortex were considered as one. Clusters in left and right visual cortex were also considered as one. The three clusters in the left temporal cortex were also considered as one. Thalamic connectivity values were not investigated further, because we found no significant main effect of ToD in the whole-brain, repeated measures ANOVA. Cohen’s d statistic for repeated measures was calculated after Dunlap et al.52. Functional results of the BOLD SD analyses were rendered on the ch2better.nii template using MRIcron (Version 4 August 2014; http://www.mccauslandcenter.sc.edu/mricro/mricron/). For visualization purposes, the high-resolution T1 images (only 13 participants, one participant failed to appear for the recording of the MP-RAGE sequence) were segmented and normalized in SPM 12, using the standard MNI templates. Then, using ImCalc, we averaged the T1-images and skull stripped the averaged image by masking with the gray-matter and white-matter images from the segmentation step. Functional results for Supplementary Fig. 1 were rendered on the averaged T1-image using MRIcroGL (Version 1 Jan 2015; http://www.mccauslandcenter.sc.edu/mricrogl/).

Multi-variable adjustment of linear mixed models

Resting-state activity can be affected by a range of potential confounders. To test whether the observed ToD effects in BOLD SD were driven by such effects, we calculated a linear mixed model with resting-state BOLD SD in the visual cortex as dependent variable. This model included ToD and scanning day as repeated measures fixed effects (autoregressive covariance structure AR(1)), as well as heart rate, breathing rate, body temperature, subjective sleepiness, amplitude and variance of head motion parameters, and the number of outliers in head motion parameters (defined as values above the mean ± 2,5 standard deviations), chronotype, sleep pressure, and sleep debt as fixed effects. Subjects were included as a random effects variable. We used the restricted maximum likelihood estimation method and tested whether the main effect of ToD in BOLD SD in the visual cortex remained significant (P < 0.05), when including all these potential confounding factors. To study whether chronotype, sleep pressure, and sleep debt, as well as scanning day explained part of the ToD modulation in BOLD variance, we additionally included multiplicative interaction terms in the multi-variable adjusted, linear mixed model. Interactions were considered significant at P < 0.05. A significant interaction between ToD and scanning day could potentially reveal habituation effects due to repetitive measurements on two subsequent scanning days. Yet, the interaction was not significant (F(5, 41.829) = 0.871, P = 0.509, type III F-test, n = 14).

Region of interest analysis in the suprachiasmatic region

We additionally investigated ToD effects on resting-state BOLD SD in a region of interest in the suprachiasmatic region, because this area contains the circadian pacemaker, the suprachiasmatic nucleus. The region of interest was defined using a published sphere20 and served for a small volume correction in the BOLD SD voxel-wise group analysis. We tested for a main effect of ToD using an F-contrast with a statistical threshold of P < 0.05, small volume corrected.

Analysis of visual detection task fMRI data

The same data acquisition, preprocessing and single-subject statistical modeling parameters as for the resting-state BOLD SD analyses were applied to obtain voxel-wise BOLD SD for the whole brain, in each session (ToD) of the visual detection task blocks (see above). The first 5 volumes at the beginning of each block were removed to avoid carry-over effects from the previous block. Then, the three task blocks were concatenated and BOLD SD calculated separately for each ToD which resulted in a voxel-wise BOLD SD map for each ToD. Fixation blocks were not further analyzed, because the experimental conditions were similar to resting-state. Using SPM 12 for group statistical analyses, task-related BOLD SD in the visual cortex was probed for the ToD-dependent effects observed during resting-state. As with the resting-state data, task data on the two consecutive days were averaged voxel-wise at each ToD for each participant. We tested for BOLD SD reductions at 08:00 and 20:00 h compared to midday measurements (11:00, 14:00, and 17:00 h) using a one-sample t-test. The analysis was restricted to the visual cortex cluster that showed a main effect of ToD during resting-state and results were thus corrected within this mask at P < 0.05, FWE cluster-corrected (cluster-forming threshold P < 0.001, uncorrected). Values from significant clusters were averaged within each cluster and extracted for each ToD and participant separately using MarsBaR and used for subsequent analyses (see below). We ascertained the consistency of the effect by using the global conjunction hypothesis53 of the contrasts (08:00 h separately against 11:00, 14:00, and 17:00 h) and (20:00 h separately against 11:00, 14:00, and 17:00 h).

Comparing resting-state and task-related BOLD SD

We tested for significant overall decreases in BOLD variance during visual detection compared to rest using dependent, two-tailed, paired samples t-tests in SPSS (significance threshold at P < 0.05, Bonferroni corrected). The peak-voxel BOLD SD values of significant clusters from the visual detection task were extracted for each ToD and participant using MarsBar, then averaged over all ToD and compared with the ToD-averaged resting-state BOLD SD values at the same coordinates.

Visual detection task performance

Omission errors (i.e., lapses) and reaction times for each participant and ToD were analyzed using SPSS Statistics. To test for ToD-dependent differences in visual detection mirroring the ToD-dependent changes of BOLD SD in sensory cortices, we performed a repeated measures ANOVA (ToD as independent variable with 6 levels, Mauchly’s test of sphericity was P ≥ 0.05, unless otherwise stated) with subsequent planned pairwise comparisons using dependent, two-tailed t-tests, testing the difference between 08:00 and 20:00 h and, midday measurements (11:00, 14:00, and 17:00 h), respectively (significance threshold set at P < 0.05).

Correlation analysis of omission errors and resting-state BOLD SD

A Pearson product-moment correlation analysis was performed to test for a linear relationship between resting-state BOLD SD in the visual cortex cluster and omission errors in visual detection at the same times. Since an ordinary correlation is not appropriate for repeated measures54, a Pearson correlation coefficient accounting for the dependence among the repeated measures was calculated, which removes the variance between subjects55. Therefore, the resulting correlation coefficient represents the association between BOLD SD and omission errors over all ToD for all participants, corrected for interindividual differences. The significance threshold was set at P < 0.05.

Correlation analysis of omission errors and task-related BOLD SD

Using the same procedure as above, we also tested for a linear relationship between task-related BOLD SD in visual cortex and omission errors in visual detection at the same times of day. Since the analysis of task-related BOLD SD resulted in two significant clusters in the visual cortex, correlations were calculated for BOLD SD in the two clusters separately, significance threshold P < 0.05.

Power analysis for thalamic resting-state functional connectivity

To ascertain that the negative finding in ToD-dependent thalamocortical connectivity as an objective marker of vigilance did not result from insufficient statistical power, we performed a power analysis using a bootstrapping method on a published dataset in which vigilance changes have been investigated18. For this re-analysis we used exactly the same parameters that we applied in our study here and investigated whether vigilance changes (N1 sleep stage vs. wakefulness) could be detected in 14 subjects in the previous dataset. In detail, a total of 93 continuous (≈3.5 min) epochs of wakefulness and 54 of N1 sleep were selected from Tagliazucchi and Laufs18. Two sets of 14 epochs were randomly selected (with replacement) from the wakefulness and N1 sleep sets and the voxel-wise thalamic functional connectivity (as determined by the AAL atlas47) was computed. We tested for significant group differences using a mass-univariate two-sample and two-tailed Student’s t-test, thresholded at P < 0.001 cluster-forming threshold and P < 0.05, FWE corrected on the cluster level. For each iteration of the randomly selected 14 wakefulness/N1 sleep epochs, voxels were flagged as significant if their associated P-value fulfilled these criteria. This was repeated for 1000 iterations, and a voxel-wise map, showing the proportion of times a voxel was deemed significant, was obtained. Since significant clusters were observed, the group size of 14 participants as used in our study can be judged sufficiently large to be sensitive to detect vigilance changes.

Data availability

The imaging data that support the findings of this study are available at G-NODE with identifier https://doi.org/10.12751/g-node.16591b. The code to calculate voxel-wise variables of interest in fMRI data is available upon demand.