**2.1 Participants**

A total of 492 volunteers (307 females, aged 19 to 80 years) were recruited from Southwest University (SWU, China) 17. All participants reported no psychiatric disorder, substance abuse, or MRI contraindications. The project was approved by the Research Ethics Committee of the Brain Imaging Center of Southwest University, following the Declaration of Helsinki. Written informed consent was obtained from each participant.

**2.2 Imaging acquisition and preprocessing**

All rs-fMRI data were obtained from a 3T Siemens Trio MRI scanner (Siemens Medical, Erlangen, Germany) at the Brain Imaging Center of SWU. Subjects were asked to close their eyes, rest without thinking about any in particular, but refrain from falling asleep. Two hundred and forty-two volumes were acquired for each subject using the T2-weighted gradient echo planar imaging (EPI) sequence: 32 slices of 3 mm, slice gap = 1 mm, TR/TE = 2000/30 ms, flip angle = 90°, field of view = 220 mm × 220 mm, resulting in a voxel with 3.4 × 3.4 × 4 mm3. The MRI data used in this study are available to the public from the International Data-sharing Initiative (http://fcon_1000.projects.nitrc.org/indi/retro/sald.html).

Image preprocessing was conducted using the Data Processing Assistant for Resting-State fMRI package (DPARSF, http://www.restfmri.net) 18 according to steps in previous studies 5,19: removing the first 12 volumes, slice timing and realignment. Subjects whose translational and rotational displacement exceeded 2.0 mm or 2.0° or mean frame-wise displacement (FD) exceeded 0.2 were excluded. The remaining sample included 322 subjects (194 females; mean age = 41.48, SD = 17.36). Images were then normalized to the standard EPI template, resampled to a 3 × 3 × 3 mm3 cube, and spatially smoothed (6-mm FWHM Gaussian kernel). Linear detrend, white matter, cerebrospinal fluid signals, and Friston 24 motion parameters were used as regressors to reduce head movement and non-neuronal information 20.

**2.3 Power spectrum analysis of GS fluctuations**

The GS was obtained by averaging signals over all gray matter voxels constrained by the binary automated anatomical labeling (AAL) 90 mask 21,22. The Welch method with hamming window (window width 0.031 Hz, overlap rate 50 %) was applied to transform time series into frequency domain 23. Data were cutoff within 0.007 ~ 0.25 Hz for de-noising 24. The power-law function *y* = *a* × *x**b* was applied to separate the fractal trend from oscillations because the original BOLD signal consisted of a scale-free trend and multiple oscillations 25. Frequency boundaries of oscillations were determined by the local minima on the mean power density curve of all subjects 16. The spectral centroid (SC) of each oscillation was calculated with Eq. (1), representing the center of gravity of the power spectrum within the given range of oscillation 26.

Where *f* = 0.25/256 Hz, representing the width between two successive frequency points, *P*(*i*) indicates the power at the *i*th frequency point within *i*1 ~ *i*2 Hz.

**2.4 Hemodynamic response function (HRF) de-convolution**

The basic hypothesis underlying the BOLD signal is the convolution of neural events and neurovascular coupling 27. In order to determine whether the relationship between the power of GS and age is resulted from neural activity, the blind hemodynamic response function (HRF) de-convolution approach was performed. According to our previous studies 22,28, the following steps were conducted. After noise regression, the point process analysis was adopted to detect spontaneous neural events 29. BOLD signals larger than mean plus one SD were detected and the onsets of neural events were extracted for HRF reconstruction 30. The HRF in each voxel was evaluated by matching BOLD signal with the canonical HRF and its time derivative. After that, neural level signals were recovered by Wiener de-convolution (http://users.ugent.be/,dmarinaz/HRF_deconvolution.html) 31.

**2.5 Contributor detection for the relationship between GS fluctuations and age**

The same analysis as the original data was performed for de-convolved data. The only difference was that the linear function *y* = *ax* + *b* is used to separate the trend from oscillations because the power-law trend disappeared after HRF de-convolution (see Fig. 1).

Using Pearson’s correlation, we evaluated the relationship between age and relative indices, including the mean and SD of GS, GS power at each frequency point, SCs of two oscillations, coefficients (*a*, *b*) of power-law and linear functions. Paired-samples t tests on SCs were performed to examine the effect of neurovascular coupling on the frequencies of two oscillations. Lastly, the correlation between the FD and age was calculated to evaluate the contribution of head motion to our results. Multiple comparisons were corrected with the false discovery rate (FDR) method (q < 0.05).