Origin of Synchronized Low-Frequency Blood Oxygen Level–Dependent Fluctuations in the Primary Visual Cortex

BACKGROUND AND PURPOSE: Low-frequency (<0.08 Hz) fluctuations in spontaneous blood oxygen level–dependent (BOLD) signal intensity show synchronization across anatomically interconnected and functionally specific brain regions, suggesting a neural origin of fluctuations. To determine the mechanism by which high-frequency neural activity results in low-frequency BOLD fluctuations, I obtained measurements of the effects of neurovascular coupling on the frequency content of BOLD fluctuations. MATERIALS AND METHODS: 3T recordings of BOLD signal intensity in the primary visual cortex were obtained in response to visual stimuli presented at varying temporal frequencies to determine which stimulus frequencies were successfully transmitted to BOLD signal intensity. Additional BOLD time series recordings were performed in a resting state and during natural visual stimulation, and frequencies comprising BOLD fluctuations were measured. Magnetoencephalographic (MEG) time series recordings were obtained in a resting state to measure which components of MEG signal intensity best correlated in frequency distribution to observed BOLD fluctuations. RESULTS: Visually driven oscillations in BOLD signal intensity were observed up to 0.2 Hz, representing a mismatch between low-pass filter properties of neurovascular coupling and observed frequencies of spontaneous BOLD fluctuations, which are <0.05 Hz in the primary visual cortex. Visual stimulation frequencies of >0.2 Hz resulted in frequency-dependent increases in mean BOLD response. Amplitude modulation of high-frequency neural activity was measured in MEG time series data, which demonstrated 1/frequency distribution with the greatest power comprising frequencies <0.05 Hz, consistent with the distribution of observed BOLD fluctuations. CONCLUSION: Synchronized low-frequency BOLD fluctuations likely arise from a combination of vascular low-pass filtering and low-frequency amplitude modulation of neural activity.

S ince the discovery of synchronized low-frequency (Ͻ0.08 Hz) fluctuations in blood oxygen level-dependent (BOLD) signal intensity among functionally related brain regions, 1 the application of functional MR imaging (fMRI) to studying functional connectivity has led to increasing interest in identifying the physiologic source of these fluctuations. Early speculation about the technique questioned the possibility that fluctuations might be generated by artifacts such as aliasing of cardiac or respiratory pulsations in the low-frequency range, vasomotor oscillations, or effects of attention. 2 Subsequent investigations now provide strong evidence that synchronous low-frequency fluctuations have a neural source. First, time series measurements of BOLD signal intensity at a high temporal frequency (TR ϭ 400 ms) showed only minimal components of the BOLD signal intensity at cardiac and respiratory frequencies in recordings from several cortical regions. 3 Subsequent analysis of respiration artifact showed significant artifact only in regions of prominent vessels such as the circle of Willis. 4 Second, a growing literature has demonstrated high functional specificity of networks exhibiting synchronized BOLD fluctuations. Resting (no task) recordings have convincingly delineated synchronized networks that are intrinsically active in the absence of a task, [5][6][7] with anatomic boundaries corroborated from studies using a wide range of tasks by means of fMRI, positron-emission tomography, and electroencephalographic (EEG) techniques. [8][9][10] Other reports have shown synchronized networks involving precisely those regions known to be involved in visual, auditory, and language function. 11,12 Subcortical networks involved in memory 13 and thalamo-hippocampal 14 circuits also showed loci of synchrony that corresponded to known neuroanatomic pathways. Large-scale analyses have demonstrated the functionally organized topology of synchronized networks. 15 Third, networks with synchronized BOLD fluctuations have shown coherent modulation with behavior and alteration in numerous pathologic states. Direct recording from the primary auditory cortex in neurosurgical patients has demonstrated tight correlation between BOLD fluctuations and separately acquired local field potentials, phase-locked to a common stimulus. 16 Brain regions involved in episodic memory function have shown increased synchrony during performance of a memory-encoding task. 17 In a button-pressing task, a large fraction of the trial-to-trial variability in BOLD response was accounted for by synchronized fluctuations. 18 Differences in amplitude of low-frequency oscillations have been observed between eyes-open and eyes-closed conditions. 19 Alteration in synchrony of spontaneous BOLD fluctuations have been observed in autism, 20,21 Alzheimer disease, 22 mild cognitive impairment, 23 multiple sclerosis, 24 blindness, 25 schizophrenia, 26 and agenesis of the corpus callosum. 27 Yet, the mechanism by which neural activity results in BOLD fluctuations that are orders of magnitude slower than the underlying neural activity is not understood. Measurement of frequencies contributing to synchronized BOLD activity have shown distribution in which nearly all of the power in the fluctuations arises from frequencies of Ͻ0.08 Hz. 3,28 To evaluate what mechanisms might contribute to lowfrequency BOLD fluctuations, I recorded fMRI time series data at a high temporal frequency in the primary visual cortex in response to visual stimuli presented at prescribed frequencies to measure the temporal structure of neurovascular coupling. Subsequently, magnetoencephalographic (MEG) recordings of the occipital lobes were obtained to evaluate the temporal relationship between spontaneous MEG activity and BOLD fluctuations.

Subject Characteristics
Experimental procedures were approved by the University of Utah institutional review board. Informed consent was obtained for all subjects. 3T fMRI recordings were obtained from 19 healthy volunteer subjects (16 men, 3 women; average age, 29.5 Ϯ 3.6). All fMRI subjects were fitted with MR imaging-compatible lenses to allow comfortable reading of the 8-point text within the scanner. An additional 6 volunteers (5 men, 1 woman; average age, 35.3 Ϯ 6.2) were then studied with MEG, described below.

fMRI Acquisition
Images were acquired on a 3T Trio scanner (Siemens, Erlangen, Germany). Subject monitoring was performed throughout the examination by real-time eye tracking by using an infrared camera mounted on a 12-channel head coil. Scanning protocol consisted of an initial 2-mm isotropic magnetization-prepared rapid acquisition of gradient echo (MPRAGE) acquisition for the anatomic template. A field map sequence was acquired for off-line distortion and magnetic inhomogeneity correction. BOLD echo-planar images (TR ϭ 1 second, TE ϭ 28 ms, generalized autocalibrating partially parallel acquisition [GRAPPA] with an acceleration factor of 2, 11 sections at 3-mm section thickness) were obtained for 4 minutes during presentation of a flashing checkerboard stimulus shown to alternating left and right hemifields in 20-second intervals with a block design. Motion correction was performed during BOLD imaging with a prospective acquisition-correction technique sequence.
On-line thresholding was performed in real time by using a general linear model with Neuro3D software (Siemens) to obtain activation maps. Three-or 5-section regions in an oblique axial plane including the pericalcarine primary visual cortex and either the lateral geniculate nuclei (LGNs) or the frontal eye field foci of activation were selected for all subsequent BOLD series, and additional 4-minute acquisitions with the same visual stimuli were performed to use as a template for visual regions of interest in off-line postprocessing. Subsequent BOLD series used the same sections and protocol (3 or 5 sections, 3-mm section thickness, 64 ϫ 64 matrix with in-plane resolution of 3 ϫ 3 mm, TR ϭ 200 ms for 3 sections in 8 subjects, TR ϭ 300 ms for 5 sections in 11 subjects, TE ϭ 28 ms to central k-space, GRAPPA ϭ 2).
Nineteen subjects then were assigned 1 of 3 protocols. In 6 subjects, 14-minute BOLD acquisitions were performed in a resting state condition (eyes closed, subjects instructed to "remain awake and let thoughts pass through their mind without focusing on any one thought") and then during natural visual stimulation (Bugs Bunny cartoon, Looney Tunes).
In 7 subjects (data from 1 subject were discarded because of an inability to remain awake through the examination), visual stimuli were flashed for a single screen refresh on an LCD projector at several prescribed frequencies. Each run consisted of a single frequency condition with a 7-minute duration. Between 4 and 6 runs were performed per subject.
In 6 subjects, visual stimuli were flashed at 6 Hz during ON periods and an isoluminant gray screen was displayed with a central fixation spot during OFF periods. ON and OFF periods alternated at multiple frequencies (from the start of the ON period to the start of next ON period). Each frequency condition lasted for 1 minute, and 7 conditions were presented in pseudorandom order during 1 run. Each subject was shown the same frequency conditions in varying order between 4 and 6 runs. A total of 9 frequencies were shown to the subjects (not all frequency conditions were shown to all subjects.)

Stimulus Presentation
Visual stimuli were generated with custom-designed software by using the Psychophysics Toolbox 29 platform for Matlab (MathWorks, Natick, Mass). A stimulus-presentation computer was synchronized with the scanner by acquisition of an optical pulse emitted by the scanner at the start of each BOLD sequence via an optic-electric converter and a digital acquisition board.
Stimuli were displayed on an LCD projector with a luminous RE-6011 shutter (Avotec, Stuart, Fla), and submillisecond accuracy was confirmed with use of a photometer and oscilloscope. All presentation data were logged with a frame-by-frame audit at the level of each screen refresh.
Stimuli were viewed on an 8-inch screen mounted in the bore of the scanner approximately 6 inches from the patient's eyes and viewed by a mirror mounted on the top of the head coil. All subjects were fitted with magnetically inert goggles with appropriate vision corrective lenses.

fMRI Postprocessing
Off-line postprocessing was performed in Matlab (MathWorks) by using SPM5 toolbox with MarsBaR toolbox extension. 30 A field map sequence was used for distortion correction, and all images were motion corrected by using a realign and unwarp procedure. BOLD images were then coregistered to the MPRAGE anatomic image sequence. No image normalization or smoothing was performed.
Activation maps were generated from dedicated 4-minute visual stimulation (Fig 1) BOLD data by using a general linear model (false discovery rate, P ϭ .005). Regions of interest corresponding to the primary (pericalcarine) visual cortex, extrastriate visual cortex, LGN, and frontal eye fields were obtained for each subject by selecting the dominant cluster in each region. This technique allowed voxel selection independent of frequency conditions in subsequent BOLD acquisitions.
Time series were extracted with the MarsBaR toolbox for each region of interest by using average data across the region of interest for each BOLD sequence obtained in that subject. Time series data were processed by linear detrend to correct for scanner drift. Averaging for FUNCTIONAL ORIGINAL RESEARCH each frequency condition and power spectral analysis was performed with custom-designed software for the Matlab platform.
The resting fMRI time series data were extracted from the visual cortical, LGN, and frontal eye field clusters described previously by using MarsBaR toolbox. A linear detrend function was used to correct for scanner drift over each 14-minute acquisition. Power spectral analysis was performed with Time Series Tool software in Matlab.

MEG Analysis
MEG data were collected by using a Neuromag (Elekta Instruments, Stockholm, Sweden) whole-head system with 306 MEG sensors comprising 102 magnetometers and 204 gradiometers. The MEG sensor elements are in 102 locations. There are 2 orthogonally placed planar gradiometers and 1 magnetometer at each sensor location. The Elekta Neuromag whole-head MEG system has a single dewar, and the sensors are immersed in liquid helium to keep the SQUIDS in cryogenic temperature. The MEG system is located in a 2-layered magnetically shielded room. The data were collected with a sampling rate of 600 Hz and a high-pass filter of 0.01Hz, with a 32-bit resolution.
Six subjects were studied with MEG, consisting of a 10-minute recording with their eyes closed. Subjects were instructed to remain awake and to try to clear their minds, letting thoughts pass through rather than focusing on any 1 mental activity. A gradiometer overly-ing the left occipital pole was selected for each patient, and raw time series data were extracted.
Raw data were then bandpass-filtered in Matlab with the Time Series Tool by using the Ideal Pass Filter in ␦, ␣, ␤, and ␥ ranges. Data were then full-wave rectified, and the envelope function was computed by replacing each point in the time series by the maximal value of the time series within 1 period corresponding to the lower bound of the bandpass frequency range. Thus, for the ␣ range filter (8)(9)(10)(11)(12)(13), the sliding maximal function with a 125-ms window (corresponding to 8 Hz) was used to calculate the envelope function. Power spectral analysis was performed with the Time Series Tool in Matlab.

Synchrony between Functionally Related Brain Regions
In 6 subjects, BOLD time series data in a resting state (eyes closed, no task) and during natural visual stimulation were acquired for 14 consecutive minutes each. From these subjects, statistically significant clusters were seen in postprocessed images in the LGN in 3 subjects, in the frontal eye fields in 2 subjects, and in the primary visual cortex and extrastriate visual cortex in all subjects.
The contralateral primary visual cortices showed highly synchronous BOLD fluctuations (Fig 2) in all subjects. Synchrony was also observed between the visual cortex and the LGNs (Fig 2C) in each of the subjects measured. Synchrony between the bilateral frontal eye fields was seen with a slower time course, indicated by a broad slow peak on the cross-correlation plot in 1 subject and a smaller narrower peak in the other subject ( Fig 2C). The very slow crosscorrelation peak in frontal eye fields in 1 subject was distinctly different from the correlation observed in other visual areas, and it is possible that this represents an artifactual physiologic signal intensity such as from eye movements during the acquisition.

Effect of Neurovascular Coupling on the BOLD Response
In 13 subjects, BOLD responses in the primary visual cortex were recorded to flashing visual stimuli presented at varying frequencies. Two separate stimulation paradigms were used. In the first paradigm (6 subjects, Fig 3A right), radial checkerboard stimuli were flashed to the bilateral visual field for a single screen refresh (13 ms) at 6 frequencies ranging from 0.1 Hz to 4 Hz.
In the remaining 6 subjects (Fig 3A left), a 6-Hz flashing checkerboard stimulus was alternated ON and OFF at several of 9 frequencies ranging from 0.017 to 0.5 Hz. In this paradigm, the same total number of flashes of the stimulus was constant for every frequency condition.
For both stimulation paradigms, the F1(peak-to-peak amplitude) component of the driving stimulus was preserved in the BOLD responses for frequencies up to approximately 0.15 Hz, beyond which there was a gradual decline in peak-to-peak amplitude between 0.15 and 0.3 Hz (Fig 3B). Although 2 separate paradigms were used, 1 with weak stimuli and 1 with strong stimuli, virtually identical measurements of filter properties of neurovascular coupling were observed in the frequency range over which both paradigms overlapped.
Also for both stimulation paradigms, there was a uniform increase in mean BOLD signal intensity with increasing stim- ulus frequency (Fig 3C), observed even in the paradigm for which the total number of stimulus presentations was constant for each frequency condition. This suggests that periods of higher frequency neural activity than can be reproduced in BOLD fluctuations (ie, Ͼ0.2 Hz) may be transmitted as an increased mean BOLD response. If such high-frequency activity was modulated in amplitude or frequency at relatively slow rates, these data would imply that the resulting BOLD response may show a resulting slow fluctuation.

Frequency Components of BOLD Fluctuations
Power spectra were obtained to analyze frequencies present in the spontaneous (eyes closed) BOLD fluctuations from the initial 6 subjects, normalized and averaged across the 6 subjects ( Fig 4A). None of the subjects exhibited any significant contribution in visual cortical BOLD time series data from frequencies attributable to cardiac (ϳ1 Hz) or respiratory (ϳ0.2 Hz) pulsations, similar to observations by Cordes et al. 3 Rather, virtually the entire power distribution of the BOLD fluctuations was concentrated at frequencies Ͻ0.05 Hz. When these frequencies are compared with low-pass filter properties observed in Fig 3 from which stimulus frequencies (and presumably frequencies of neural activity) are successfully transmitted to BOLD signal intensity, there is a mismatch. Frequencies between 0.05 and 0.15 Hz are virtually absent from BOLD fluctuations (Fig 4B). Comparing power spectra from the BOLD time series obtained in the eyes-closed condition with the natural visual stimulation condition shows higher frequencies in the eyesclosed condition (Fig 5). Average power was calculated between 0 and 0.02 Hz (P Ͻ .02) and between 0.02 and 0.05 Hz (P Ͻ .03) for each subject. A pair-wise 2-tailed t test (n ϭ 6) demonstrated statistical significance between the 2 conditions in each frequency range (Fig 5). The difference between conditions would be consistent with the hypothesis that frequency distribution of BOLD fluctuations is influenced by underlying neural activity and that fluctuations do not merely represent correlated physiologic noise.

Amplitude Modulation of Neural Activity
To evaluate whether underlying neural activity may show slow modulation at frequencies likely to contribute to the observed low-frequency BOLD fluctuations, I studied an additional 6 subjects with MEG recordings in the eyes-closed condition. Processing of MEG data is illustrated in Fig 6. Raw data from the occipital gradiometer were first bandpass-filtered into each of the ␦, ␣, ␤, and ␥ bands, and the mean value of the data was subtracted. Each trace was then full-wave rectified, 31 and the envelope function was calculated. The resulting amplitude modulation curve of ␣ activity for 1 subject is illustrated with accompanying power spectrum in Fig 6. MEG data were obtained in the eyes-closed condition, and the full-range power spectrum of the data was dominated by ␣ activity, as might be expected for an occipital recording (Fig 7A). For this subject, the power spectrum of the envelope function for ␣ filtered activity matches precisely a 1/frequency curve (Fig 7B).
Summary data for each of the frequency bands measured are presented in Fig 7C-F. All bands show very high power at low frequencies, with curves approximating scaled versions of 1/frequency, though frequency bands show variable rates of decay in power with increasing frequency. For all frequency bands, there is a striking increase in power at frequencies below ϳ0.05 Hz.

Discussion
An accurate understanding of the mechanism whereby spontaneous neural activity translates to much lower frequency  BOLD fluctuations is increasingly important as techniques using these synchronized fluctuations are being applied to a wide range of pathologic conditions. [17][18][19][20][21][22][23][24] In particular, it would be useful to know whether synchronized fluctuations reflect correlated neural noise or whether there is neurophysiologic information accessible within the temporal structure of the fluctuations. Data presented here indicate that a large part of the mechanism responsible for translation of neural activity to BOLD signal intensity is mediated by filter properties of neurovascular coupling. Low-pass filtering is demonstrated to restrict transmission of oscillations at frequencies of Ͼ0.3 Hz and results in increasing damping of frequencies between 0.1 and 0.3 Hz.
Yet neurovascular coupling alone is unlikely to explain the frequency distribution of observed BOLD fluctuations, in which frequencies Յ0.05 Hz are actively represented, but frequencies between 0.05 and 0.2 Hz are not represented to the extent expected by filtering alone. From these data, it is expected that some other mechanism contributes to the observed frequency range of spontaneous BOLD fluctuations. Although low-frequency fluctuations have traditionally been described as Ͻ0.08 Hz, 27 the visual cortical fluctuations observed in this report were almost exclusively Ͻ0.05 Hz.
A contribution of high-frequency neural activity to syn-chrony in BOLD signal intensity, and not merely of lowfrequency filtered components or harmonics of the neural activity, is also suggested from several other reports. Examination of local field potential (LFP) in the monkey visual cortex demonstrated that though synchrony in the raw LFP diminished rapidly with cortical distance, band-limited power in the ␥ frequency range showed synchrony over much larger cortical distances. 29 In a simultaneous measurement of LFP and optical hemodynamic response, there was tight correlation of power of local field potential oscillations in the ␥ range with hemodynamic response. 32 Gamma-range LFP and BOLD responses showed high temporal correlation and coincided with periods of synchrony among interneurons, in a study measuring individual unit activity, LFP, and BOLD response. 33 Other reports have also shown correlation between lowfrequency synchronized BOLD fluctuations and other ranges of higher frequency neural activity. An investigation of simultaneous EEG and fMRI revealed a correlation between BOLD fluctuations and power in ␦, , ␣, ␤, and ␥ ranges, with variable contributions in different spatially separated distributed networks. 34 An analysis of resting-state networks by using epidural EEG and fMRI in anesthetized rats showed the highest correlation between the ␦ range power in the contralateral somatosensory cortex. 35 Higher frequency ranges have also been shown to be dependent, with coupling between ␥ and è oscillations in 1 study, 36 supporting the idea of distributed synchrony among temporal frequency ranges.
One possibility for generating slow fluctuations in hemodynamic responses is infrequent high-amplitude epochs of spontaneous neural activity, or neural avalanches. 37 Temporal inhomogeneity of spontaneous neural activity is well described. Many neocortical neurons exhibit, for example, slow 2-state fluctuations in membrane potential and spiking rates. 38 Such fluctuations participate in sensory-information encoding 39 and are synchronized over wide cortical networks. 40 Discrete states at slow time scales among more distributed neural networks have also been proposed as a mechanism for neural computation. 41 Yet analysis of time series data in BOLD and MEG responses does not show discrete periods of activity followed by relative inactivity but shows continuous slow fluctuations with somewhat variable frequency and amplitude within a narrow range. Thus, temporal inhomogeneity may contribute to low-frequency fluctuations, but data appear more consistent with a slow modulation of frequency and amplitude of neural activity than infrequent, salient, and discrete periods of activity.
The data previously mentioned demonstrate that at frequencies of Ͼ0.3 Hz, increasing frequencies of presented stimuli (and neural activity) are transmitted as increased mean BOLD signal intensity, suggesting that slow modulation in the time of the BOLD signal intensity could be produced by a similar modulation in the time of frequency or amplitude of neural activity. MEG data confirm that such slow modulation in the amplitude of given frequencies is indeed present. Moreover, across all frequency bands, envelope functions of amplitude modulation demonstrate 1/frequency distribution, with the greatest power below 0.05 Hz.
Such modulations would be an excellent candidate for amplification of very low frequencies seen in BOLD fluctuations. They would preserve synchrony of high-frequency power that other studies have shown to be correlated to BOLD lowfrequency synchrony and could explain the mismatch between low-pass-filter properties of neurovascular coupling shown previously and observed BOLD frequency distribution.
The 1/frequency distribution seen in amplitude modulation of neural activity is a pervasive property of numerous physical systems, ranging from light emission by quasars to flow through rivers and current flow through resistors, 42 and is seen extensively in biologic systems including fMRI time series data. 43 The presence of 1/frequency distribution in the modulation of the band-limited amplitude of neural activity suggests that this modulation is generated by a stochastic process, though this does not exclude the possibility that parameters of the noise may reflect static or dynamic properties of the neural system.
The data in this report are drawn from different populations of volunteers, in particular the subjects from whom fMRI and MEG data were obtained. It might strengthen the agreement in frequency distribution between the MEG amplitude modulation data and the fMRI time series data presented previously if the results had been drawn from the same population. Yet the frequency ranges observed in visual cortical fMRI time series are remarkably consistent with those reported in other studies, 3 and observing similar results for frequency ranges of resting state data across multiple modalities, subject populations, and experimental designs strengthens the argument that amplitude modulation of neural activity produces frequency distribution that is conserved across subjects.

Conclusion
Synchronized BOLD fluctuations in the primary visual cortex comprised frequencies of Ͻ0.05 Hz, but visual stimulation with frequencies Յ0.2 Hz can drive oscillations in the BOLD signal intensity. Therefore, although low-pass filtering from neurovascular coupling likely makes a large contribution to the frequency components observed in BOLD fluctuations, it is insufficient to explain the observed frequencies by itself. Amplitude modulation of neural activity appears to occur at predominantly low frequencies, with the greatest power Ͻ0.05 Hz and 1/frequency distribution. Such changes in amplitude of neural activity are likely reflected in the BOLD response by a slow fluctuation and would be an ideal component for amplification of very low frequencies in the BOLD response. Amplitude modulation, together with vascular low-pass filtering, could provide a mechanism for the observed frequency distribution of spontaneous BOLD fluctuations.
If synchronized BOLD fluctuations, which have been shown to correlate with neural responses, brain pathology, and functional connectivity, are generated at least in part from amplitude modulation of epochs of synchronized high-frequency neural activity, then the implication that information may be encoded in such low-frequency modulation suggests a target to explore in future studies. Moreover, because the distribution of low frequencies of BOLD fluctuations was altered with visual stimulation, the timing of amplitude modulation of neural activity may represent a useful metric of regional connectivity that could be used in neurophysiologic and pathologic investigations.