Introduction

Rhythms of the brain are generated by neural oscillations occurring across multiple frequencies (Buzsaki 2006). The natural logarithm linear law (N3L) offers a theoretical framework for parcellating these brain oscillations into multiple frequency intervals linking to distinct physiological roles (Penttonen and Buzsáki 2003). Remarkably, when graphed on the natural logarithm scale, the centers of each frequency interval fall on adjacent integer points. Thus, distances between adjacent center points are isometric on the natural logarithm scale, resulting in a full parcellation of the whole frequency domain where each parcel of the frequencies is fixed in theory, namely frequency intervals. These frequency intervals have been repeatedly observed experimentally (Buzsaki and Draguhn 2004). This characteristic suggests that distinct physiological mechanisms may contribute to distinct intervals. These brain oscillations can be measured by different technologies such as EEG and MEG. Much like EEG/MEG recording signals, functional magnetic resonance imaging (fMRI) represents a non-invasive and safe technique with an acceptable trade-off between spatial and temporal resolution by recording the blood oxygen level dependent (BOLD) as the fMRI time series. FMRI has the potential to contribute to the study of certain neural oscillations in the human brain in vivo.

In early fMRI studies of the human brain, researchers tended to treat oscillations across different frequencies without differentiation. Low-frequency oscillations measured by resting-state fMRI (rfMRI) have been assessed primarily in the frequency range of 0.01 to 0.1 Hz, a range in which spontaneous brain activity has high signal amplitude (Biswal et al. 1995; Lowe et al. 1998). While such efforts have been somewhat informative, treating this broad frequency range in a unitary manner may conceal information carried by different frequency intervals. To address this issue, an early study decomposed the rfMRI signals into multiple frequency intervals using the N3L theory (Slow-5: 0.01 - 0.027 Hz, Slow-4: 0.027 - 0.073 Hz, Slow-3: 0.073 - 0.198 Hz, Slow-2: 0.198 - 0.25 Hz)(Zuo et al. 2010). This demonstrated the feasibility of mapping distributional characteristics of oscillations’ amplitude in both space and time across multiple frequency intervals in the brain.

Since then, an increasing number of rfMRI studies have employed such methods by directly applying these frequency intervals, and have detected frequency-dependent differences in brain oscillations in patients. Specifically, these differences were mostly evident between Slow-4 and Slow-5 amplitudes (Han et al. 2011; Jing et al. 2012; Zhao et al. 2015; Mascali et al. 2015; Li et al. 2017; Ren et al. 2016). Such frequency-dependent phenomena have also been explored using other rfMRI metrics including regional homogeneity detected in the Slow-3 and Slow-5 frequency ranges (Wang et al. 2016). While the lower and upper bounds of the frequency intervals are fixed in theory, their highest and lowest detectable frequencies and frequency resolution are determined by the sampling parameters (e.g., rate and duration) in computational practice. However, the above-mentioned studies applied the frequency intervals from earlier studies (Di Martino et al. 2008; Zuo et al. 2010) rather than to use those matching their actual sampling settings. To address this situation, we developed an easy to use toolbox to decode the frequency intervals by applying the N3L theory. This toolbox, named DREAM, is based on MATLAB with a graphical user interface (GUI). Here, we introduce the N3L algorithm and its DREAM implementation. Neural oscillations reflected by the human brain spontaneous activity measured with resting-state functional MRI and head motion data during mock MRI scans were employed as two worked examples to demonstrate the use of DREAM to perform frequency analyses.

Methods and Algorithms

Neuronal brain signals are temporally continuous but they are almost always measured as discrete data for practical reasons. The characteristics of the sampled data should meet the criterion of the sampling theorem proposed by American electrical engineers Harry Nyquist and Claude Shannon. The core algorithm to determine the frequency boundaries of measured neuronal signals in DREAM is based on the Nyquist-Shannon sampling theorem. Specifically, per the theorem, sampling frequency and sampling time determine the highest and lowest frequencies that can be detected and reconstructed. Sampling data retains most of the information contained in the original signals if the sampling frequency is at least twice the maximum frequency of the continuous signals. As for neuronal signals, the highest frequency that could be detected and reconstructed is determined by the sampling frequency, or by the sampling interval which is equal to the reciprocal of the sampling frequency:

$$ f_{max}=\frac{1}{2T_{R}} $$
(1)

where fmax represents the highest frequency that could be detected in the neuronal signal and TR represents the sampling interval.

The lowest frequency in neuronal signals that could be detected depends on the sampling time. As shown in formula (2), in order to distinguish the lowest frequency in neuronal signals, the sampling time should be equal to or larger than the reciprocal of two times the lowest frequency:

$$ T\geq\frac{1}{2f_{min}} $$
(2)

where T represents the sampling time, and fmin represents the lowest frequency in neuronal signals that could be distinguished.

Since the sampling time is equal to the number of samples multiplied by the sampling interval, the lowest frequency can be calculated by:

$$ f_{min}=\frac{1}{2NT_{R}} $$
(3)

where N represents the number of samples.

According to the N3L theory, neural oscillations in mammalian brain formed a linear hierarchical organization of multiple frequency bands when represented on a natural logarithmic scale. The center of each band would fall on each integer of the natural logarithmic scale (Fig. 1-1). Thus, adjacent bands have constant intervals that equals to one, which correspond to the approximately constant ratios of adjacent bands on the linear scale (Fig. 1-2). With the highest and lowest frequencies reconstructed, N3L can derive the number of decoded frequencies and the boundaries of each frequency interval (Fig. 1-3). Accordingly, when graphed on the natural log scale, the center of each decoded frequency is an integer. Thus, adjacent center points on the natural log scale are equidistant, which corresponds to the same proportion of adjacent center points’ values on the linear scale. Based upon this theorem, after performing a linear regression analysis for the highest and lowest frequencies acquired previously, we can determine the central frequencies, as well as the number frequency intervals that can be decoded.

Fig. 1
figure 1

The flowchart on the DREAM algorithm. (1) N3L theory defines an oscillator with a length-one frequency band centered at n, i.e., OSC(n), in the natural log space. (2) In original frequency space, it expands the frequency band (en− 0.5,en+ 0.5) Hz. (3) This frequency band can be discretized with a sampling procedure with N points and TR rate in terms of the classical signal theory. (4) This computational frequency band is for a band-pass filtering process to extract the OSC(n) from the raw time series

Full size image

Finally, the decoding process integrated in DREAM performs band-pass filtering with the frequency intervals provided by DREAM in the previous steps (Fig. 1-4). This is implemented by the MATLAB built-in function fft and ifft to perform direct and inverse time-frequency transformation on the signals for individual decoded frequency intervals, respectively. All the above steps are illustrated as the flowchart in Fig. 1, and relevant algorithms are presented as Pseudo-codes in Supplementary ??.

Computational Performance

The ability to handle big data is an important aspect of software performance, which has been optimised by DREAM. The frequency decoding procedure described above (Fig. 1-1 to 1-3) needs two input parameters, namely sampling frequency and the number of samples. The operation speed is thus very fast and is almost not affected by the data size. The operation speed of the second procedure, band-pass filtering, is dependent on the data size and operational environment. In order to achieve the optimal speed of its computation, we embedded three algorithms into DREAM: 1) the input data was treated as a whole matrix, 2) the input data was divided into ten chunks with approximately the same number of time series and then processed sequentially, and 3) the input data was divided into ten chunks and then processed in parallel. The first algorithm is suitable for most cases, especially for regular-size data, such as most 1.5T and 3T fMRI data. The second algorithm performs better with large-scale data using a computer with limited memory. It is a strategy that balances time (speed) and space (memory). When the computer has enough memory, the third algorithm works best for processing large-scale data. The toolbox will automatically select the best algorithm based on the data size (e.g., number of time series, number of sampling points) and the user’s hardware operating environment (e.g., memory size, number of CPU cores). We tested the performance of these three algorithms in the MATLAB version of DREAM respectively, using the same data from an individual brain of the https://github.com/zuoxinian/CCNP/tree/master/3R-BRAIN database from the https://github.com/zuoxinian/CCNP (Liu et al. 2020) (image size: 104 × 90 × 60, number of samples: 840, repetition time: 0.72s) and the same hardware environment (memory: 16 GB, CPU: four cores four threads). The speed of the first, second and third algorithms are 706.897 seconds, 744.475 seconds and 641.419 seconds respectively. The third algorithm on a high-performance computer can achieve the maximum speed in dealing with big data (e.g., 7T fMRI data).

Interface and Usage

DREAM has been shared and released with the Connectome Computation System (CCS)(Xu et al. 2015). After downloading the package at https://github.com/zuoxinian/CCS/tree/master/H3/DREAM, users need to add the directory where the package is stored into the MATLAB path. The package can be launched by entering ”DREAM” in the MATLAB command line. DREAM integrates its GUI (two buttons) into its flash screen (Fig. 2).

Fig. 2
figure 2

DREAM’s flash screen

Full size image

Program Interface

DREAM supports CCS data structures by default. Users should enter or organize their data into the predefined directory structure (Fig. 3) before start processing the data. The work directory is where the subject directories are stored (full path). Individual data should be stored in each subject directory or a sub-folder inside (data directory). DREAM has a main interface (Fig. 4) for setting up the structure (the left side) and previewing the plots of time series from the data selected (the right side). DREAM is also compatible with the Brain Imaging Data Structure (https://bids.neuroimaging.io) (Gorgolewski et al. 2016).

Fig. 3
figure 3

DREAM’s directory structure

Full size image
Fig. 4
figure 4

The main interface of DREAM

Full size image

GUI Usage

We introduce how to use the graphical interface step by step in below. The circled numbers in Fig. 4 correspond to the analyzing steps in this section. The following steps are applicable to users using the CCS data format. For BIDS users, just specify the data directory, and the toolbox will automatically identify the directory structure and read in data.

After all the above parameters are set up, data meeting the requirement will appear in the list-box (Figs. 4, 5, 6), from where the user can remove unwanted data by selecting the file name and clicking the Remove button. Finally, by clicking the Divide button, a user can start the decoding program. The outputs contain a set of decoded files and a csv file that records the boundary frequencies of each decoded band. The outcomes can be directly used for subsequent analyses.

Fig. 5
figure 5

A preview of the original FD time series from a participant

Full size image
Fig. 6
figure 6

DREAM decodes FD time series into the five bands

Full size image

DREAM-1: Frequency-dependent Oscillations of In-scanner Head Motion in 3-16 Years-old Children

Kids usually move more than adults, representing a natural characteristic improved across development. Frequency insights of this behavioral trait would benefit the understanding of its neurodevelopmental underpinnings. We thus performed a purely behavioral research of the head motion to unveil the frequency characteristics of head motion as well as their relationships with age and sex. Specifically, we took head motion oscillations as the behavioral trait of interest to study their physiological characteristics. We believe that characterizing head motion properties in multiple frequencies can provide new insights of head motion during fMRI scanning. In-scanner head motion has long been treated as a confounding factor in most fMRI studies, especially in studies of children and patients. Many studies have shown the effects of motion on fMRI results such as increases of short-distance correlations and decreases of long-distance correlations in rfMRI-derived connectivity metrics(Power et al. 2012; Power et al. 2015; Yan et al. 2013). Researchers also proposed various methods to correct the motion effects in fMRI studies. In contrast, studying head motion as a neurobehavioral trait has been overlooked (see an exception in (Zeng et al. 2014)), especially in children. Here, we employed DREAM to quantify head motion data acquired from preschool and school children in a mock scanner using a novel multi-frequency perspective. We hypothesized that: 1) head motion is a behavioral trait associated with age; 2) there are sex differences in head motion in children; and 3) the head motion effects are frequency-dependent.

Participants and Data Acquisition

We recruited 94 participants (47 females) between 3 to 16 years of age as part of the https://github.com/zuoxinian/CCNP(Yang et al. 2017; Zuo et al. 2017; Liu et al. 2020), a long-term (2013-2022) large-scale effort on normative research for lifespan development of mind and brain(Dong et al. 2020). All participants were from groups visiting during the Public Science Open Day of the Chinese Academy of Sciences, with the approval of at least one legal guardian. The experiment was performed in a mock MRI scanner at the site of the MRI Research Center of the Institute of Psychology, Chinese Academy of Sciences. The mock scanner was built by PST (Psychology Software Tools, Inc.) using a 1:1 model of the GE MR750 3T MRI scanner in use at the institute. It is used for training young children to lie still in a scanner before participating the actual MRI scanning session. It is decorated with cartoon stickers to provide a children-friendly atmosphere. Head motion data were acquired with the MoTrack Head Motion Tracking System (PST-100722). The system consists of three components: a MoTrack console, a transmitter and a sensor. The sensor is worn on the participant’s head and provides the position of the head relative to the transmitter. For each participant, head motion is displayed on the computer screen in real-time. The original sampling rate of the system is 103 Hz. The averaging buffering size is 11 samples, which results in a recording sampling rate of 9.285 Hz. The participants were instructed to rest quietly on the bed of the mock scanner for around three and half minutes without moving their heads or bodies. They were watching a cartoon film inside the scanner during the ”scanning” to simulate movie-watching scanning. The data acquisition period was designed to resemble the real MRI scanning environment, with a recording of scanning noises of the real MRI machine played as the background noise. This design can setup much faster sampling rates than that of the slow fMRI sampling rate and thus lead to more accurate and higher frequency bands. Such a behavioral analysis of the head motion would benefit the traditional fMRI analysis if the fast-sampling head motion data is collected and analyzed in tandem with the fMRI data by the brain-motion association studies.

Data Analysis

Head motion data are recorded in text files consisting of six parameters for each time point, three translation (millimeters) and three rotation (degrees) measures. The first three parameters are displacements in the superior, left and posterior directions, respectively. The last three parameters are rotation degrees in the three cardinal rotational directions. We converted the original data into frame-wise displacement (FD), a single parameter scalar quantity representing head motion proposed by Power and colleagues (Power et al. 2012). To correct for spikes caused by sudden movements, which may bias mean FD values, we applied the https://afni.nimh.nih.gov 3dDespike command (version 17.3.06) to the FD time series. Of note, these time seires are not fMRI data and stored as text files. Data without this preprocessing was also analyzed and supported reproducible patterns. Then time-windows were determined and applied before feeding the data into DREAM. We retained 1672 sampling points from the zeroed time point (time point when the original six parameters were set to zero), which equaled a duration of three minutes. After preprocessing, we used DREAM to decode the data. Of note, the original FD values were all positive. After decoding, the time series of decoded bands were demeaned, which means the average values of all decoded time series were very near to zero. Thus, we took the absolute value of decoded frequency intervals to calculate mean FD values, which were used in subsequent statistical analyses. Inspired by many human growth curves modeled by exponential function and the scatter plots on the head motion data, we first converted the head motion data using the natural logarithm transformation and then assess the relationship between FD and age by using linear regression models to fit the FD data in each frequency interval with age. We conducted this regression for boys and girls, respectively, and tested whether the slopes and intercepts are significantly different between boys and girls. Of note, this method is equivalent to an Analysis of Covariance (ANCOVA)(Rosner 2015). These analyses were also applied to the standard deviation of FD time series to test the stability of head motion.

Results

Six participants were excluded from further data analysis due to sampling periods less than three minutes. Another four participants were excluded because their mean FD values were three standard deviations higher than the mean value of the whole group (i.e., outliers). Total 42 boys (age: 3 - 14 years, 8.7 ± 3.0) and 42 girls (age: 4 - 16 years, 8.4 ± 3.1) were included in our final analyses. No significant differences in age were found between males and females. All the findings derived with the head motion data without despike preprocessing are highly similar to those of using despike, which are reported as following. Meanwhile, all the results derived from the linear regression models are replicated by the ANCOVA model.

Frequency Decomposition

Since all the head motion data have the same sampling frequency and sampling period, DREAM decoded all the FD time series into the same six frequency intervals named according to Buzsaki and Draguhn (2004) (Slow-4: 0.033 to 0.083 Hz, Slow-3: 0.083 to 0.22 Hz, Slow-2: 0.222 to 0.605 Hz, Slow-1: 0.605 to 1.650 Hz, Delta: 1.650 to 4.482 Hz, Theta: 4.482 to 4.643 Hz). This theta band is too narrow comparing with its full range (up to 10 Hz) to be reliable for the analyses, and thus not included in our analyses. The full band and the five frequency bands from an individual child are depicted in Figs. 5 and 6.

Age-related Head Motion Changes Across Frequencies

Results from the linear regression analysis yielded significant negative correlations between age and mean FD values across all the five bands for both boys and girls (df = 40, FDR corrected p < 0.05):

  • Slow-4: boys, p = 0.018,R2 = 0.218; girls, p = 0.034,R2 = 0.195

  • Slow-3: boys, p = 0.008,R2 = 0.249; girls, p = 0.027,R2 = 0.203

  • Slow-22: boys, p = 0.001,R2 = 0.314; girls, p = 0.017,R2 = 0.221

  • Slow-1: boys, p < 0.001,R2 = 0.358; girls, p = 0.013,R2 = 0.230

  • Delta: boys, p < 0.001,R2 = 0.380; girls, p = 0.008,R2 = 0.250

The relationship between age and mean FD values are plotted in Fig. 7, indicating that younger children tend to move more than older ones, and this trait correlation held in both boys and girls. We also performed a similar linear regression analysis between the standard deviations of decoded FD values and age, and observed similar outcomes that the standard deviations were significantly negatively correlated with age across frequency bands and sexes. This showed older children are more stable with their head motion than younger children.

Fig. 7
figure 7

Nonlinear age-motion relationship across the five frequency bands. The plots are based upon the log transformed motion data, indicating the exponential growth model \(y_{\text {motion}}=\mathbf {e}^{(ax_{\text {age}}+b)}\). The upper-left panel shows the Mock scanning facility in the Magnetic Resonance Imaging Research Center at the Institute of Psychology, Chinese Academy of Sciences

Full size image

We further tested if the two lines are different between boys and girls. Statistical results revealed no such sex-related effect (df = 80, FDR corrected p > 0.05):

  • Slow-4: slope, p > 0.5,F = 0.383; intercept, p = 0.494,F = 3.979

  • Slow-3: slope, p > 0.5,F = 0.531; intercept, p = 0.385,F = 4.428

  • Slow-2: slope, p > 0.5,F = 1.177; intercept, p = 0.486,F = 4.010

  • Slow-1: slope, p > 0.5,F = 1.326; intercept, p = 0.849,F = 3.042

  • Delta: slope, p > 0.5,F = 1.222; intercept, p = 0.968,F = 2.822

Inspired by the trend that sex-related differences in mean FD are smaller in higher frequency bands, especially evident for early stages, we thus divided all the participants into three age groups (3 to 6 years: 14 boys, 18 girls; 7 to 9 years: 14 boys, 15 girls; 10 to 16 years: 14 boys, 9 girls) and compared mean FD values between males and females in each age group using two-way (sex and frequency band) ANOVA with repeated measures. Figure 8 summarized the results of an increasing pattern of head motion from slow to fast bands for all the age groups (3-6yrs: F(4) = 10.90,p = 1.65 × 10− 7; 7-9yrs: F(4) = 20.62,p = 1.20 × 10− 12; 10-16yrs: F(4) = 23.95,p = 3.06 × 10− 13). Meanwhile, we observed a significant interaction between sex and frequency band in 7 to 9 years old children (F(4,1) = 3.22,p = 0.0154) but not for the other groups (3-6yrs: F(4,1) = 0.195,p = 0.940; 10-16yrs: F(4,1) = 1.065,p = 0.380).

Fig. 8
figure 8

Sex-frequency interactions on head motion across ages. All the participants (3 to 16 years old) are divided into three age groups: 3 to 6 years, 7 to 9 years, 10 to 16 years). A two-way (sex and frequency band) ANOVA with repeated measures compares mean FD values between males and females in each age group

Full size image

DREAM-2: Frequency-Dependent Spatial Ranking and Reliability of Low-Frequency Oscillations

The amplitude of low frequency fluctuation (ALFF) is a common metric used in fMRI studies that reflects regional amplitude of the signal intensity’s fluctuations in a frequency range (Zang et al. 2007). Previous studies revealed variations of ALFF in both spatial and frequency domains in the resting-state brain. From the perspective of spatial distribution, in the typical resting-state frequency range (e.g., 0.01-0.1 Hz), the neural oscillations showed higher ALFF in grey matter than white matter (Biswal et al. 1995; Turner et al. 1993). It is noted that ALFF reaches its peaks in visual areas (Kiviniemi et al. 2003), posterior structures along brain midline (Biswal et al. 1995; Zou et al. 2009) and in cingulate and medial prefrontal cortices (Ghosh et al. 2008). In frequency domain, BOLD oscillations distributed to grey matter were mainly in Slow-4 and Slow-5, while its white matter oscillations were dominated by Slow-3 and Slow-2 (Zuo et al. 2010). Specifically, higher ALFF in Slow-4 was detected in the bilateral thalamus and basal ganglia whereas the slow-5 oscillators exhibited higher ALFF in the ventromedial prefrontal cortex, precuneus and cuneus (replicated in Xue et al. 2014). These findings revealed the frequency-specific characteristics of resting-state ALFF. The previous studies are limited by their sampling precision (TR ≤ 2000ms), and studies on the ALFF distribution across more accurate bands and their reliabilities are still lacking. For examples, the Slow-2 frequency band derived in Zuo et al. (2010) has quite small overlap with its theoretical range and thus may limit both reliability and validity of its findings. Here, we use DREAM to decompose the fast (TR = 720ms) rfMRI data from the Human Connectome Project (HCP) (Van Essen et al. 2013) test-retest dataset, to 1) map the ranks of ALFF values through Slow-1, Slow-2, Slow-3, Slow-4, Slow-5 and Slow-6 bands and 2) evaluate the test-retest reliability of the ALFF metrics in these different frequency bands.

Participants and Data Acquisition

The test-retest dataset from HCP consisting of 45 subjects were used for this analysis. All subjects were scanned with an HCP-customized Siemens 3T scanner at Washington University, using a standard 32-channel receive head coil. Three participants were excluded from the substantial analyses because their resting-state scan durations were shorter than others. Forty-two subjects (aged 30.3 ± 3.4 years, 29 males) were included in the present study. In the dataset, each subject paid two visits. The average interval between the two visits is 4.7 months. In each visit, each subject was scanned two times in two consecutive days and each scan contained structural images (T1w and T2w), two rfMRI, seven runs of task fMRI and high angular resolution diffusion imaging (see details of http://protocols.humanconnectome.org/HCP/3T/imaging-protocols.html). The two visits constitute a long-term test versus retest contrast, while the two scans within a visit constitute a short-term test versus retest contrast. Since one subject lacked a scan session in the first visit, we only included 41 subjects in the short-term analyses for visit 1 and the long-term analyses. In the present work, we only used the rfMRI data, which consisted of 1200 volumes (TR = 720 ms; TE = 33.1 ms; flip angle = 52∘, 72 slices, matrix = 104 × 90; FOV = 208 × 180 mm; acquisition voxel size = 2 × 2 × 2 mm). The data were preprocessed according to the HCP MR preprocessing pipeline (Glasser et al. 2013), resulting in the preprocessed surface time series data fed to the following DREAM analysis. The preprocessing pipeline includes: 1) Gradient distortion correction; 2) Motion correction; 3) EPI image distortion correction; 4) Registration to T1w image; 5) One step spline resampling; 6) Intensity normalization and brain masking; 7) Transfer volume-based timeseries into surface-based timeseries.

Amplitude Analysis

For each rfMRI scan, we first extracted the representative time series for each of the 400 parcels (Schaefer et al. 2018) by averaging all the preprocessed time series within a single parcel. DREAM decomposed the time series into its components across the potential frequency bands. We performed ALFF analysis for all the bands of each run and each subject according to Zuo et al. (2010) implemented by CCS (Xu et al. 2015). Subject-level parcel-wise ALFF maps for each frequency band were standardized into subject-level Z-score maps (i.e., by subtracting the mean parcel-wise ALFF of the entire cortical surface, and dividing by the standard deviation). In a single visit, the two standardized ALFF maps in the same session were then averaged, resulting in two (short-term test versus retest) standardized ALFF maps per frequency band for each subject. For the long-term test-retest contrast, the two standardized ALFF maps in the same visit were averaged, and resulting in two long-term ALFF maps per frequency band per subject. To investigate both short-term and long-term test-retest reliability of ALFF across the five frequency bands, we calculated the parcel-wise intraclass correlation (ICC) based upon short-term and long-term ALFF maps respectively (Zuo et al. 2013; Xing and Zuo 2018). We averaged the two standardized ALFF maps of all the subjects in each visit separately to obtain the group-level standardized ALFF maps per visit. And we also averaged the group-level maps of the two visits to obtain the group-level standardized ALFF maps for the whole dataset. In order to evaluate the spatial distribution of the ALFF values for each parcel, we assigned its rank of ALFF values to the parcel (from 1 to 400). All the above analyses were done for each of the five frequency bands, leading to three ALFF ranking maps for each frequency band (visit 1, visit 2 and whole dataset).

Results

We did the same amplitude analysis and ICC calculation for the two short-term contrasts (visit 1 and visit 2) and the long-term contrast. Results from the two visits are similar. Therefore, in the following sections, we will only display results for one short-term (visit 2), along with the long-term results. DREAM decomposed the rfMRI timeseries into six frequency bands (Slow-6: 0.007 - 0.012 Hz; Slow-5: 0.012-0.030 Hz; Slow-4: 0.030-0.082 Hz; Slow-3: 0.082-0.223 Hz; Slow-2: 0.223-0.607 Hz; Slow-1: 0.607-0.694 Hz). Spatial rankings on ALFF for one visit and the whole dataset are mapped in Figs. 9 and 10 respectively. The ranking trend is similar in the two maps. It is noticed that ALFF spatially ranked from high in ventral-temporal areas to low in ventral-occipital areas when the frequency band increased from low to high, while those in part of parietal and ventral frontal regions were reversed.

Fig. 9
figure 9

Spatially ranking ALFF across six frequency bands for data in one visit. LH: left hemisphere; RH: right hemisphere; Vis: visual network; SomMot: somatomotor network; DorsAttn: dorsal attention network; SalVentAttn: salience ventral attention network; Cont: frontal parietal control network; Default: default network; Limic: limbic network; see details of https://github.com/ThomasYeoLab/CBIG/tree/master/stable_projects/brain_parcellation/Schaefer2018_LocalGlobal

Full size image
Fig. 10
figure 10

Spatially ranking ALFF across six frequency bands for whole dataset

Full size image

Test-retest reliability maps of ALFF are also generated (Figs. 11 and 12) by mapping ICC using the linear mixed models. It is clear that, for both short-term and long-term, the higher frequency bands, the more reliable ALFF measurements. However, the long-term ICC values are lower than short-term across frequencies, which suggests that intra-subject variation grows with the scan interval. In the short-term maps, the slow-2 (0.223-0.607 Hz) demonstrated the highest test-retest reliability of ALFF across the six frequency bands.

Fig. 11
figure 11

Short-term test-retest reliability of ALFF across six frequency bands

Full size image
Fig. 12
figure 12

Long-term test-retest reliability of ALFF across six frequency bands

Full size image

Discussion

DREAM is a free and publicly available software that can decode oscillation data into multiple frequency bands. The simple interface was designed to allow all users to easily perform multi-band frequency analyses. The computational methods employed in DREAM to calculate the numbers and ranges of decoded frequency bands apply the Nyquist-Shannon sampling theorem and the brain oscillation theory (Buzsaki and Draguhn 2004). Such a theory has been proven of great potentials to understand the brain dynamics as well as their behavioral correspondences. From a theoretical perspective, the oscillation theory can be independent of any modalities (e.g., EEG, MEG, ECoG, TMS, fMRI, fNIRS, eye tracking, etc.) for measuring these oscillations as windows into brain waves (Balduzzi et al. 2008). DREAM is thus applicable for multiple forms of discrete sampling data, as long as the data are entered in the supported format. Currently, DREAM can process both NIFTI formatted neuroimaging data and text file formatted behavioral data while more other formats will be supported in its forthcoming releases.

As a demonstration of its utility, the results derived with DREAM for pure behavioral recordings suggest that head motion may be a behavioral feature reflecting both state and trait of individuals. We showed that head movements in the high frequency bands are more evident than those in the low frequency bands. This could be a behavioral reflection of the hierarchical organization of brain oscillations for their synchronization at multiple scales in space. Neural oscillations of the higher frequency-bands are related to more local information processing while the lower frequency-bands are for more distant communications in the brain. Our findings are consistent with the previous observation that the head motion had more impacts on the short-distance brain connectivity. While the head motion during fMRI scanning has been treated an important confounding factor in the neural signal (Power et al. 2015), some recent work also argued its neurobiological components related to individual traits of the motor behaviors (e.g., Zeng et al. 2014, Zhou et al. 2016). The current researh offers data for an alternative explanation on such neurobehavioral trait likely driven by brain systems operating within a multi-band frequency landscape. In the context of development, as we expected, younger children moved more than older children across all the slow frequency bands. The stability of head motion during the experiment also varied with age, with head motion becoming less variable or more stable in older children. This is more evident in higher frequency bands, an implication that more sudden and sharp movements in younger children. Moreover, in a specific age range (7 - 9 years), boys moved more than girls across Slow-6 to Slow-1 bands but such differences vanished in the delta frequencies. This age range is a critical period for developing the ability to apply effective cognitive control (i.e., cognitive flexibility during executive function) (Anderson 2002), and our findings might reflect the sex differences in the cognitive development. In summary, our results demonstrate the necessity to study the frequency-specific characteristics of head motion, especially a perspective on understanding the neurobiological mechanism behind these behavior-related oscillations. This is of great potential to enrich our knowledge on the lifespan development such as children, the elderly and patients with neurologic or psychiatric conditions where both distance-related brain and the head motion measurements have been observed to correlate with each other (Andrewshanna et al. 2007; Satterthwaite et al. 2012; Fair et al. 2007, 2013).

Differences in head motion across ages or between cohorts may reflect differences of certain traits, which may co-vary with detected brain signals and behavioral outcomes. The different properties of head motion in different frequency bands show that there may be different mechanisms associated with different frequencies. Head motion at higher frequencies varies more with age, and this may reflect that cognitive control associated with higher frequencies develops better with age. Of note, interpolation analyses indicated that this observation is not related to an issue of better signal-to-noise ratio at higher frequencies because there are more events per unit time. Within the narrow age range of 7 to 9 years old, boys moved more than girls in most frequency bands, although sex differences were larger at lower frequencies. This may indicate that the development of controlling system associated with lower frequencies may have larger sex-related differences for this age range. The above results lead us to speculate that there may be two control systems that are associated with different frequency bands of head motion which develop differently with age and between boys and girls. More detailed experimental studies deserved to test this postulation in future. The strategies of dealing with head motion issues in human brain mapping may also need updates regarding its measurement reliability and validity in terms of the possible neurobiological correlates (Xing and Zuo 2018; Zuo et al. 2019a, 2019b). One promising direction is to separate various sources of the head movements by using additional recordings or developing novel motion metrics (e.g., the recent progress in Power et al. 2018, 2019, 2020). These efforts identified seven kinds of in-scanner motion in resting-state fMRI scans, and five of them related to respiration. Some pseudomotion occurred only in the phase encode direction and was a function of soft tissue mass, not lung volume. Using the Mock scanning experimental design as in the present work, together with the aforementioned approaches, could be of high value in further understanding neurobiobehavioral underpins of the human head movements.

Using fast fMRI from HCP, at the first time, we revealed the spatially configuring pattern of ALFF ranking gradually from low to high frequency bands. This indicates a trend along the two orthogonal axes. Along the dorsal-ventral axis, higher ALFF ranks were moving from the ventral occipital and the ventral temporal lobe up to regions in the parietal lobe as the frequency increasing. Along the anterior-posterior axis, from lower to higher bands, higher ALFF ranks were walking from the posterior to the anterior regions in the ventral part. This frequency-dependent ALFF pattern is similar to the findings of previous studies on the association between brain structure and gene expression, which also reported orthogonal gradations of brain organization and the associated genetic gradients (Chen et al. 2013; Kremen et al. 2013). The underlying physiological mechanism and functional significance of the frequency-dependent ALFF patterns deserve further investigations. It is interesting that the frequency-dependent pattern of ICC is quite uniform across the brain and as the frequency increased, its reliability increased alongside. This observation illustrated that compared with the low frequency bands, higher frequency bands might be more suitable for detecting individual differences in ALFF. Most of the previous studies have adopted ALFF of the lower frequency bands (i.e., Slow-5 and Slow-4 or around 0.01 to 0.1 Hz) where their ICCs rarely met the reliability requirement (ICC ≥ 0.8) of clinical applications. In contrast, our findings suggest that both Slow-2 and Slow-1 ALFF could be the usable and reliable marker of the brain oscillations for these applications. It is noticed that the reliability of Slow-1 ALFF is slightly lower than those of Slow-2 ALFF in the short-term results, and this may be an indication on the limited Slow-1 band here compared to its theoretical range (around 0.6065 − 1.6487 Hz). Our results further show that long-term test-retest reliability is worse than short-term across all frequency bands. This indicates that the intra-subject variation increase with scan duration. In order to maintain the stability of data, it is necessary to select an appropriate scan interval. While studies of the very fast sampled fMRI signals such as HCP are sparse, it is quite promising for future studies with multiple neuroimaging modalities (e.g., Balduzzi et al. 2008, He et al. 2008) to DREAM as an integrative tool across frequencies. An open toolbox such as DREAM is essential for large-scale projects inspired by the increasing practice of open sciences coming with more and more fMRI and EEG datasets openly shared as well as their applications (e.g., Zuo2020).

Limitations and Future Works

Despite the advantages of DREAM presented in the paper, some limitations should be noted. First, DREAM is based on MATLAB, which is a commercial computing software. In order to generalize the application of DREAM, a standalone version will be developed. For now, a light online version of DREAM based on Python has been opened to public. The bandpass filter is highly amenable to optimization: the individual columns of the matrix can be parallelized for linear speedup, and the fft itself can be offloaded onto a GPU if the data will fit into the GPU’s memory. It may also be possible to reduce the complexity of the bandpass to O(N × T) using the approach of (Pankovski 2016) in the future. Second, in the current version, limited data formats are supported. More data formats will be added in subsequent software updates.

Information Sharing Statement

The DREAM toolbox is fully open to the public by sharing both https://github.com/zuoxinian/CCS/tree/master/H3/DREAMand http://ibraindata.com/tools/dream using Python. To ensure the reproducibility of our findings, all the codes and head motion data for generating the figures and other results in the present work are also shared via DREAM and CCS website. Please credit both DREAM and CCS work with their citations if you use our DREAM.