{'fraction_non_alphanumeric': 0.030391292683861572, 'fraction_numerical': 0.019737094239834437, 'mean_word_length': 4.927760109041345, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Steady-State BOLD Response Modulates Low Frequency Neural Oscillations

Neural oscillations are the intrinsic characteristics of brain activities. Traditional electrophysiological techniques (e.g., the steady-state evoked potential, SSEP) have provided important insights into the mechanisms of neural oscillations in the high frequency ranges (>1 Hz). However, the neural oscillations within the low frequency ranges (<1 Hz) and deep brain areas are rarely examined. Based on the advantages of the low frequency blood oxygen level dependent (BOLD) fluctuations, we expected that the steady-state BOLD responses (SSBRs) would be elicited and modulate low frequency neural oscillations. Twenty six participants completed a simple reaction time task with the constant stimuli frequencies of 0.0625 Hz and 0.125 Hz. Power analysis and hemodynamic response function deconvolution method were used to extract SSBRs and recover neural level signals. The SSEP-like waveforms were observed at the whole brain level and at several task-related brain regions. Specifically, the harmonic phenomenon of SSBR was task-related and independent of the neurovascular coupling. These findings suggested that the SSBRs represent non-linear neural oscillations but not brain activations. In comparison with the conventional general linear model, the SSBRs provide us novel insights into the non-linear brain activities, low frequency neural oscillations, and neuroplasticity of brain training and cognitive activities.

N eural oscillations, the intrinsic characteristics of dynamic brain activities, are the patterns of a particular temporal beat (e.g., the occipital alpha wave around 10 Hz 1 ). Using electrophysiological techniques, previous studies have revealed a number of important processes of neural oscillations in the high frequency ranges (.1 Hz), such as the spectral fingerprints 2 , the phase-amplitude coupling 3 , and the multiple generative mechanisms 4 at various brain regions and frequency bands. However, the mechanisms of low frequency (,1 Hz) neural oscillations are rarely examined 5 , although their physiological significance has been emphasized by the local field potentials (LFP), electroencephalogram (EEG), functional magnetic resonance imaging (fMRI) and behavioural studies [6][7][8] .
The lack of sufficient knowledge of low frequency neural oscillations is largely due to the deficiency of conventional approaches or methods. In the domain of electrophysiological techniques, the local field potentials (LFP) are invasive and spatially constrained 8 . The electroencephalography (EEG) and magnetoencephalography (MEG) are limited in recording signal bandwidth (.0.5 Hz) and spatial resolution 9 . In contrast, the fMRI is a promising technique to explore low frequency neural oscillations because of its low frequency blood oxygen level dependent (BOLD) fluctuations (,0.25 Hz for 2 sec repetition time) and high spatial resolution. However, the BOLD responses measure the downstream effects regulated by the neurovascular coupling rather than the neural activity itself 7 . Giving the fact that traditional general linear model (GLM) heavily relies on the neurovascular coupling 1 , the question of whether the low frequency BOLD fluctuations reflect the underlying processes of neural oscillations is still under debate. Therefore, an appropriate tool is necessary to explore the low frequency neural oscillations.
In the high frequency ranges, an effective index of modulation for neural oscillations, the steady-state evoked potential (SSEP), is widely used for investigating cognitive activities in cognitive or clinical neuroscience, and brain-computer interface 10 . The SSEP is an evoked periodic response that peaks at the same frequency and harmonics of a regularly repetitive stimulus after fast Fourier transform 10 . With the advantages of stable spectrum 11 , high signal-to-noise ratio (SNR), and sensitivity to sensory stimuli 10,12,13 , the SSEP is a critical and unique means in examining the mechanisms of basic cognition (e.g., visual attention 14 and working memory 15 ). However, because of the low spatial resolution of EEG and the limited response areas in the brain, using high frequency stimulus presentations would result in an unproductive application of the SSEP when exploring higher-level cognition and deep brain activities.
In the low frequency ranges, taking the advantages of the fMRI technique, we examine (1) whether the steady-state BOLD response (SSBR) could be evoked at the fundamental frequency of stimuli and harmonics, and (2) whether the SSBR could modulate the underlying responses of neural oscillations by two steps: imitating the form of SSEP and using the method of blind hemodynamic response function (HRF) deconvolution to obtain the neural level signals from BOLD signals 16 . For the first step, we hypothesized that SSBR, a SSEP-like phenomenon, could be induced because of the close relationship between BOLD signal and neural activities and between low and high frequency bands. First, the low or infra-slow frequency BOLD fluctuations are closely associated with the signal oscillations of EEG and LFP 6,8 , indicating that BOLD fluctuations measure the substantial neural activities. Second, the distinct frequency ranges of BOLD and EEG signals may also be linked by the nested frequencies 5,9 . This phenomenon characterizes the power of a high frequency component phase-locking to a low frequency component 3 . Furthermore, the dominant frequencies of different brain regions or functional networks are also found in low frequency BOLD fluctuations 17 . For the second step, the blind HRF deconvolution is established on the idea that BOLD spikes are derived from the point events with non-random patterns 18 . By matching BOLD signals with canonical HRF and its time derivative 16 , the blind HRF deconvolution could eliminate the effect of HR in a maximal extent 19 . Therefore, we can compare the SSBR before and after HRF deconvolution to determine whether the SSBR reflects neural oscillations.
In summary, we aimed to devise a reasonable tool to investigate low frequency neural oscillations. First, we evoked the SSBR based on the theory of steady-state neural response in high frequency ranges and the close relationship between high and low frequency neural activities. Second, we adopted the recently developed blind HRF deconvolution method to recover neural level signals form BOLD level signals to inspect the neural basis of SSBR.

Results
SSBRs of the whole brain. A SSEP-like waveform was shown in Figure 1A. The HRF deconvolution changed the energy distribu-tion along frequency ranges, but did not eliminate SSBRs ( Figure 1B). This suggests that the SSBRs are independent of neurovascular coupling. Comparing to the resting state, the LF condition significantly increased the power within 0. These results suggest that SSBRs can be evoked at the fundamental frequency of stimulus and its harmonics. This observation is similar to those of the SSEPs. On the other hand, the HRF did not markedly exert influence on the SSBRs, suggesting that the SSBRs are not requiring the hemodynamic function.
Regional SSBRs. Figure 2 shows that SSBRs were evoked at 0.06-0.065 Hz and 0.1225-0.1275 Hz frequency bands in the left sensorimotor area (SA) and bilateral supplementary motor area (SMA) on both LF and HF conditions. Of note, no other regions had the similar effects with over 20 voxels at the 0.185-0.19 Hz frequency bands on the LF condition. Table 1 summarizes the visual cortex (VC), insula, ventral lateral frontal cortex (VLFC), inferior parietal lobe (IPL), posterior middle temporal gyrus (MTG) and inferior temporal gyrus (ITG) resonated on the HF condition at 0.1225-0.1275 Hz frequency band.
The SSBRs after HRF deconvolution throughout the brain are slightly higher (not statistically significant different) than those before HRF deconvolution. This suggests that the SSBR is independent of HRF convolution, a key hypothesis in the conventional general linear model (GLM). Except for the SA and SMA regions, most areas of the brain showed a little higher power (the ratio of Task/Rest .1 but not arrives significance) on task conditions than on the resting condition ( Figure 3). Although there were no significant differences, to certain extent, the stimuli presented in a constant signal frequency changed the brain state or rhythmic activity. This result, therefore,   supports the perspective of neuroscience that cognitive processing is implemented by widely distributed brain areas acting in collaboration 20 .
The Effect of HRF. Figure 1 suggests that the HRF significantly contributed to the amplitude of BOLD signals. In addition, the HRF deconvolution changed the phase of task signals (1/2 , 3/4 cycle for the HF condition, 1/4 , 3/8 cycle for the LF condition; Figure 4A), because a 4-6 s delay of the neurovascular coupling was eliminated 16 . Due to the phase changes, the task signals showed negative correlations (LF: r 520.353, p ,0.001; HF: r 520.0.783, p ,0.001; Resting: r 50.089, p 50.139) before and after HRF deconvolution. In summary, the HRF changed the amplitude ( Figure 1) and the phase (Figure 4) of BOLD signals, and the patterns of energy distribution along the frequency ranges ( Figure 1), rather than the relative power ( Figure 1, 2) of SSBRs. The results suggest that the effect of SSBRs was clear and independent of neurovascular coupling.

Discussion
With periodic stimuli, we successfully induced the steady-state BOLD responses (SSBRs) in the low frequency ranges. Unlike traditional GLM, the SSBR is independent of the neurovascular coupling and reflects the underlying process of neural oscillations. Therefore, the SSBR may be a meaningful tool to investigate low frequency neural oscillations. We found the SSBRs at the fundamental frequencies of both HF and LF stimuli, the first harmonic of LF stimuli, and a slightly increased SSBR at the second harmonic of LF stimuli. The harmonic phenomenon is thought to stem from the non-linearly coupled neural system 21 . We demonstrated that the non-linear neural oscillation is independent of the neurovascular coupling and differs from the GLM that depends on the linear neurovascular coupling. The between-frequency coupling has been widely observed between low and high EEG frequencies 3 , and between BOLD fluctuations and EEG/LFP oscillations 5,9 . The between-frequency coupling may be an intrinsic feature of the dynamic neural activities. As an integrity system and common basis of EEG and BOLD signals 7 , the neural activities underpin the correlations of the two types of signal in the infra-slow frequency range 7,8 and maybe the harmonic phenomenon they share. However, reliable coupling between 0.0625 Hz and 0.125 Hz cannot be obtained in the current study due to three reasons: first, time points is too few (four or eight) in each cycle to divide a cycle into several bins of phase; second, SSBRs at both frequencies are stimulus-locked and phase-related 10 ; and third, the phase-amplitude coupling is a band-limited rather than a frequency-locked phenomenon 22 , whereas the band width in the infra-slow frequency range is controversial 23,24 . In line with this hypothesis, more studies are required to examine the mechanism of harmonic phenomenon in the low frequency neural oscillations.
In the conditions of LF and HF, the most modulated brain regions are the left SA and bilateral SMA. These regions play roles in tactile sensations, voluntary movements 25 , and action monitoring 26 . The SSBRs in the VC is weaker than those in the SA and SMA; this may be affected by BOLD adaptations, a widely observed phenomenon when stimuli are repeated 27 . While requiring further investigations, it is possible that BOLD adaptations may compete against the SSBRs. Although weaker than in the visual and motor areas, the SSBRs in the MFG, MTG, IPL and insula are considerably evoked at the 0.1225-0.1275 Hz frequency band, especially on the HF condition. These regions are often involved in controlling cognitive activities, such as representing color categories 28 , extracting memories 29 , processing semantic information 30 , and executing disparate cognitive, affective, regulatory functions 31 . Although other regions are regulated to some extent, the SSBRs are strongest in regions that are related to the SRT task and slightly weaker in task control regions, showing a task-related pattern.
Despite the fact that the SSBRs and brain activations are both taskrelated, the underlying mechanisms for them may be different. First, the brain activations require the HRF, whereas the SSBRs may not; Second, periodic visual stimuli reduce the activations 32 but enhance the SSBRs; Third, the effects of mean BOLD signal and BOLD variability are largely irrelative 33 . Furthermore, the BOLD variability is more useful in examining cognitive aging, functional disorders, performance and neuroplasticity 33,34 . Considering the difference between mean BOLD signal and BOLD variability, understanding the mechanism of SSBRs remains an important area for future research. A variaty of tasks, specially those involving in higher-order cognition and deep brain activities, should be performed to uncover the mechanism of the SSBRs. From the point of view of neural oscillations, the SSBRs likely provide us novel insights into how the brain participates in cognitive activities.
The SSBR shows an enlarged variability at the particular frequency when the stimuli are presented at a constant frequency. The increased variability is critical for the neural system to optimally operate. First, recent studies suggested that the neural system works at the ''edge of criticality'' among many possible states or functional network configurations 33,35 . The greater variability of neural oscillations provides a greater dynamic range for the brain to select proper responses from many states or options in an ever-changing environment 33,34 . Therefore, variability is essential to the flexibility, efficiency and adaptability of the neural system 34 . According to the adaptive hypothesis, the SSBRs may provide more insightful evidence than brain activations in studies of life-span development, skill learning, or training. Second, the elevated variability may provide kinetic energy for the neural system to achieve a variety of possible functional states 33 . A low variability can not offer sufficeint capacities for the system to achieve new states (e.g., from the resting state to a task state), resulting in a rigid and maladaptive system 36 . The BOLD variability has been demonstrated to be related to the capacity of swithing from resting state to task state 33 . While requiring further investigation, we suggest that the SSBRs, by increasing neural variability, are also related to the switching capacities. Third, both local information processing and coordination of large-scale networks have their preferred frequency bands, named spectral fingerprints 2 . The spectral fingerprints exist within both the high EEG frequency range 2,37 and the low BOLD frequency range 17,24 for distinct brain regions or functional networks. In the present study, the frequency band bias may be one reason for the distribution of SSBRs in both spatial (e.g., relative strength throughout the brain) and frequency dimensions (e.g., the highest SSBRs at the 0.1225-0.1275 Hz frequency band for both LF and HF conditions). Although high signal-to-noise ratio at 0.125 Hz than at 0.0625 Hz (Figure 1) due to lower intrinsic power within higher frequency band than in lower frequency band may also cause higher and more widespread SSBRs, the spectral fingerprints has been supported by previous studies 17,24 and may be a reasonable hypothesis. It is known that at the natural or dominant frequency, the SSBRs exert the largest effect, and maximize the benefit for neural variability, flexibility and efficiency 37 . Therefore, evoking the SSBRs at the preferred frequency may promote the application of this neural oscillations based method in studies of neuroplasticity.
Although SSBRs and SSEPs share certain characteristics in neural activities, the generative mechanisms of low and high frequency neural oscillations may differ from one another 3,7 . Even within the high frequency ranges, different frequency bands are driven by distinct neural mechanisms 4 . Therefore, more work, such as that using the band-limited method 22 , is necessary to clarify the neural basis of the SSBRs and low frequency neural oscillations.
Finally, we obtained the findings from a sample of young males. In spite of the fact that the steady-state responses were found in both genders and across life span 34,36 , these results should be tested in a wider range of populations. In addition, only two frequencies those are multiple of the TR are tested in the current study. These limited frequencies may reduce power in the analysis and cannot uncover all features of the SSBRs. The characteristics of the SSBRs should be explored with more cognitive tasks and more frequencies. Furthermore, although the HRF deconvolution has been demonstrated to be able to recover neural level signal, simultaneous fullband EEG-fMRI recording is warranted to directly reveal the neural mechanism of SSBRs.

Conclusions
We show that the SSBRs modulate the neural oscillations in a nonlinear fashion. This observation is similar to the effect of the SSEPs, but differs from the linear style of the GLM. The SSBRs are taskrelated and enhance brain efficiency by increasing neural variability. As an important assessment tool in examining the neural activities, the SSBRs can provide us novel and critical insights into neuroplasticity, non-linear brain functions and low frequency neural oscillations during various cognitive tasks.

Methods
Subjects and procedure. Thirty male college students participated in the experiments (From 18 to 20 years of age; 18.30 6 0.53, mean 6 SD). All subjects were right handed confirmed by the Chinese version of Edinburgh Handedness Questionnaire (coefficients . 40). The subjects had normal or correct to normal vision, were free from any medications, neurological and psychiatric disorders. Written informed consent, approved by the research ethical committee of University of Electronic Science and Technology of China, was obtained from each subject before the experiment. The methods were carried out in accordance with the approved guidelines. The data of three subjects were removed from the final analysis due to large head motions (translation .2 mm or rotation .2u). The data of another subject were not used because the MRI data were lost.
In the resting scan of 10 min (the Resting condition; 300 volumes), subjects were required to remain still, focus their eyes on the white crosshair at the center of the gray screen without thinking anything particular during the image acquisition period. In two task runs of 10 min each, subjects were asked to focus their eyes on a green disc (visual angle 3.43u) located in the center of the grey screen. When the disc turned into red (100 ms presentation), the subject pressed a key with the right thumb. Simple visual stimuli were used as in most SSEP studies 10 , while simple reaction was asked to avoid mind wandering or falling asleep. The frequency of the simple reaction time (SRT) task was 0.0625 Hz (once every 16 seconds; the lower frequency condition, LF) or 0.125 Hz (once every 8 seconds; the higher frequency condition, HF) in different runs. These frequencies were located in slow 4 (0.027-0.073 Hz) and slow 3 (0.073-0.198 Hz) frequency bands 23 , respectively. The former has been demonstrated to be suppressed during cognitive processing 38,39 , while the latter gets less attention but is associated with physiological and mental activities 39,40 . The frequency we selected here allowed us to test whether the SRT task would enhance or attenuate BOLD power at both lower and higher frequency bands. The conditions of LF and HF were counterbalanced between subjects with a 10 , 15 min inter-run break. In each task, subjects were asked to remain focused without counting or predicting the stimuli.
Scan acquisition. MRI data were acquired using a 3.0T GE 750 scanner (General Electric, Fairfield, Connecticut, USA) equipped with high-speed gradients. An 8channel prototype quadrature birdcage head coil fitted with foam padding was applied to minimize head movement. Functional images were acquired using a gradient-recalled echo-planar imaging (EPI) sequence. The parameters were as follows: repetition time (TR) 5 2000 ms, echo time (TE) 5 30 ms, 90u flip angle, 39 axial slices (4 mm slice thickness without gap), 64 3 64 matrix, 24 cm field of view.
Data preprocessing. Functional images were preprocessed using the Data Processing Assistant for Resting-state fMRI (DPARSF 2.2, http://restfmri.net/forum/DPARSF). The preprocessing steps included: the 1 st 25 scans (about 3 cycles for LF and 6 cycles for HF) were discarded to allow evoked fluctuations to appear, signal to reach equilibrium and participants to adapt to the scanning noise. The remaining images were slice-time corrected, spatially aligned and then spatially normalized to Montreal Neurological Institute (MNI) EPI template and re-sampled to 3 3 3 3 3 mm 3 isotropic voxels 41 . The head motion was calculated according to the method proposed by Power et al. 42 . There were no differences for the head motion among conditions: F (2, 50) 50.41, p 50.658. The normalized images were spatially smoothed with a 6 mm full width half maximum (FWHM) Gaussian kernel. The linear trend of time courses was removed. Afterwards, six head motion parameters, white matter signal and cerebrospinal fluid signal were regressed out before power analysis. The global signal was not regressed out because global signal regression may remove a global neuronal signal that is caused by a widely distributed ascending input 43 and enhance the neuronal-hemodynamic correspondence 44 .
The Blind HRF Deconvolution. After noise signal regression, the point process analysis was used to detect spontaneous and evoked point events on the resting and task conditions, respectively. BOLD fluctuations of relatively large amplitude (. 1 SD) were collected and the onsets of neural event were saved for HRF reconstruction 19 . The HRF in each voxel was obtained by matching BOLD signal with the canonical HRF and its time derivative. After that, signals at the neural level were recovered by Wiener deconvolution (http://users.ugent.be/,dmarinaz/HRF_ deconvolution.html) 16 .
The Power Analysis. The Power analysis was performed before and after HRF deconvolution. We first calculated the power of the three runs at the whole brain level to test whether SSBRs were evoked by the SRT task at particular frequencies. As an exploratory research, the coarse and holistic analysis was necessary. In addition, brain activities modulated by the cognitive task may be dispersive 20 and may be enhanced 10 or attenuated 38 . The whole brain analysis might facilitate our understanding of these phenomena. The time course of each voxel was converted to the frequency domain without band-pass filtering using the Fast Fourier Transform (FFT). The FFT was adopted because it is the main method to define SSEPs 10 and the same index would benefit the comparison between SSEPs and SSBRs. Frequency resolution was 0.0018 Hz (Sampling rate/sampled data: 0.5 Hz/275). The power of each run for each subject was obtained as the average of all voxels in the gray matter constrained by the Automated Anatomical Labeling (AAL) template without cerebellum 45 . Pairedsample t-tests of LF vs. Resting were carried out at the fundamental frequency (0.06-0.065 Hz, 0.005 Hz bandwidth), the first harmonic (0.1225-0.1275 Hz) and the second harmonic (0.185-0.19 Hz) of the LF stimuli, and of HF vs. Resting was performed at the fundamental frequency (0.1225-0.1275 Hz) of the HF stimuli. The effect of HRF was assessed by comparing the relative power of Task/Resting before and after HRF deconvolution at each of the abovementioned frequency interval because the absolute values were dramatically changed (from several hundreds to about 0.1-0.2; see Figure 1) by the HRF deconvolution.
The frequency bands of 0.06-0.065 Hz, 0.1225-0.1275 Hz and 0.185-0.19 Hz were selected to test the frequency specific power change induced by SRT task at each voxel. The amplitude of low frequency fluctuations (ALFF, the square root of power spectrum) was adopted to examine regional energy change within the three intervals 46 . Regional changes of amplitude for contrasts of LF vs. Resting, HF vs. Resting and the HRF effect were tested by the paired-sample t-test embedded in the SPM8 software (http://www.fil.ion.ucl.ac.uk/spm). All resulting statistic maps were corrected using the family wise error (FWE) method (p ,0.05) for multiple comparisons 47 .
According to the perspective of neuroscience, most brain regions should be modulated during cognitive processing because cognitive processing is implemented by distributed collaboration of widespread brain regions 20 . Considering the rigidity of FWE correction, we further computed the relative power of Task/Resting before and after HRF deconvolution based on the ALFF analysis to examine the manifestation of SSBR throughout the brain without statistical correction. The relative power was used because the absolute values were very different before and after HRF deconvolution.
We further extracted signals from the left sensorimotor area and performed the Pearson correlation analysis to identify the relationship of signals before and after HRF deconvolution. The left sensorimotor area is task-related and modulated by both LF and HF stimuli. Pearson correlation analysis was run on the LF, HF, and Resting conditions, respectively. Before correlation analysis, the signals were normalized by subtracting the mean and dividing by the standard deviation.