Introduction

Since its discovery in the early 1990s, functional magnetic resonance imaging (fMRI)1 with blood oxygenation level dependent (BOLD) contrast2,3,4,5 has become the mainstay of human brain imaging. BOLD fMRI reflects changes in cerebral blood flow, volume and oxygenation. Flow increases contribute to positive BOLD signals, diluting paramagnetic, signal-reducing deoxyhæmoglobin6. Control of blood flow involves modulations of vessel diameters, thus changing the relative sizes of intravascular (IV) and extravascular (EV) compartments. When the IV compartment expands, the EV compartment shrinks, meaning that the local fMRI signal also depends on the volume of the EV tissue. Potentially, conformational changes other than vasodilation, such as swelling of activated brain cells7,8,9, can influence the local volume fractions of tissue and cerebrospinal fluid (CSF) during neuronal activation10,11.

The commonplaceness of fMRI has given the impression that this imaging method can be used in an “automatic” mode, without worrying about the biophysical basis of fMRI signal changes. Some recent studies do not report even the basic acquisition parameters (APs), such as echo time (TE), repetition time (TR), or flip angle (FA), as if these were well-established defaults with no appreciable influence on the data. However, even though an informed choice of APs can maximize BOLD contrast, it remains unclear which APs are the most suitable for inferring the actual neuronal function from the signal and whether some APs may generate sham activations or eliminate neuronal activation-related signals. Although the EV signal that varies under the direct influence of the IV compartment must be regarded as an integral part of the functional contrast mechanism, the diversity of APs employed in fMRI warrant an inspection of BOLD fMRI signals to evaluate effects confounding their interpretation in terms of the neuronal origins.

Short TE (e.g. ~15 ms at 3 T) yields little BOLD contrast compared with longer (≥30 ms) TE, meaning that data acquired at multiple TEs could differentiate between BOLD and non-BOLD signals12,13. On the other hand, flip angle (FA) in general alters tissue contrast. In a recent single-TE study14, FA seemed irrelevant to the functional contrast of the fMRI signal (expressed as percent signal change, %ΔS). A low FA has many beneficial effects (reviewed in Ref. 14), including reduced artifacts due to blood inflow, subject movement and physiological noise (e.g. respiration and cardiac pulsation), reduced radiofrequency energy deposition and improved tissue contrast, with apparently no detrimental effects on BOLD contrast. Thus, the traditional choice of FA (at the Ernst angle) was challenged and vastly reduced FAs were suggested14. Notwithstanding, since stimulus-related conformational changes remain potential confounds, we varied both TE and FA to assess the effects of tissue contrast weightings on fMRI signals elicited during visual stimulation and to disambiguate BOLD from non-BOLD signal changes.

Results

Thirteen healthy subjects were imaged at 3 T with four FAs (12.6°, 22.5°, 50.0° and 90.0°, subsequently labeled as FA12, FA22, FA50 and FA90) and a standard TE (30 ms) while they were viewing a reversing circular checkerboard. To identify BOLD changes, scans were repeated with a shorter TE (16.8 ms) for the highest and lowest FAs (TE values later referred to as TE16 and TE30). The TR was 0.8 s except for two subjects whose data were acquired at TR = 0.5 s and 1.0 s. The data of five subjects had to be rejected for motion- and physiological-related confounds.

Figure 1 presents the wide-spread aggregate effects of APs on fMRI responses. The yellow–blue Venn diagrams, with green overlap areas, indicate that for the two different FAs presented, the sets of activated volume image elements (voxels) are largely disparate for three of the seven anatomical atlas-based regions of interest (ROIs) in both the left and the right hemispheres [posterior cingulate cortex (PCC), precuneus and posterior superior temporal gyrus (pSTG)] but similar in the four visual areas [primary visual cortex (V1), cuneus, lingual gyrys (LingG) and middle occipital gyrus (MOG); these atlas-derived areas have some overlap]. Fig. 1 further shows the linear regression of the β-coefficients of the two FAs being evaluated in statistical comparison. We applied three different masks to restrict the analysis to the voxels that responded at both or either FA (see the legend of Fig. 1). The responses from voxels that were responding on both or either FA were different in the majority of cases (slope ≠ 1) indicating that, on the whole, the two FAs produced quantitatively different results.

Figure 1
figure 1

Mean slopes m (across eight subjects, ± SEM) of regression lines βFA90 = m βFA12 + c for each region of interest in the left (top) and right (bottom) hemispheres.

Slopes were calculated for significantly-activated voxels common to both FA12 and FA90 (circles), significantly-activated voxels extracted using the FA12 activation map (squares) and significantly-activated voxels extracted using the FA90 activation map (triangles). Slopes statistically significantly (asterisks; two-tailed t-tests, p < 0.05) different from 1 (dashed horizontal line) indicate FA-dependent activation strengths. The colored Venn diagrams indicate the proportions of significantly activated voxels for FA12 (yellow), FA90 (blue) and voxels significantly activated at both FAs (green). The wide error range in left PCC for the FA90 mask is due to one subject for whom only 5 data points were available, with m = −4; without that subject m would be 1.14 ± 0.24.

Full size image

Figure 2 shows the differences of the β-coefficients for all subjects in several ROIs separately as well as the totals across subjects. These voxel-by-voxel differences, where the β-values at FA90 were subtracted from those at FA12, show that the shapes of the distributions for the visual areas (cuneus, middle occipital gyrus, lingual gyrus; V1 not shown separately) are regular and quite different from those for PCC and precuneus. For example, the PCC distribution is skewed toward negative difference (in all subjects except for subject S2 scanned at TR = 0.5 s). This result is also illustrated in the Supplementary Fig. S1, where many PCC voxels are in the negative region and more negative at FA12 than at FA90.

Figure 2
figure 2

Histograms portraying the FA12–FA90 difference between task-related beta coefficients in 5 ROIs; results are given separately for the individual subjects S1–S8 (the applied TR indicated below each subject label) and for the data pooled across all subjects (“All”; top row).

Voxels with significant activations [p < 0.05, false discovery rate (FDR)-corrected] at either or both of the FAs were extracted from each ROI, shown in 21 non-overlapping bins centered at −0.0050, −0.0045, … 0.0050. A voxel with equal beta coefficients for the two FAs would land in the centermost bin (indicated by a white vertical line). Values are provided for the number of voxels used in the analysis (N), the median of the difference of the beta values (, also illustrated by the blue-pink border), the scale of the y-axis (y) and the p-value for equal medians for the two FAs obtained by the Wilcoxon rank sum test; significantly different medians (p < 0.05) are emphasized.

Full size image

Of the responding voxels, less than 50% reacted at both FA12 and FA90 in PCC, whereas ~80% reacted at both FAs in cuneus, in part due to the higher average response amplitude in cuneus. Interestingly, however, the proportion of “negative” and “positive” voxels responding only at one FA differred according to FAs: in PCC, the proportion of negatives per positives was 5:1 at FA12 and 1:1 at FA90. In cuneus, the proportions were 4:5 and 1:5, respectively. In both areas, these FA dependencies were statistically significant (p < 0.0007, χ2 test; see Supplementary Tables S1–S2 for voxel counts and test results).

Functional connectivity analysis performed using left V1 as a seed also provided different results for different FAs, with the number of significantly responding voxels varying by a factor of two or more in several ROIs, with unlike trends for positive and negative responses (see Supplementary Table S3).

While the above analyses show crude differences between FAs at the group and ROI level, we next focus on data from individual subjects and voxels. To start off, the stimulation-related activation maps in Fig. 3 derived from data acquired from one subject in a single fMRI run, with two FAs interleaved in 76-s long blocks (see Fig. 7 in Methods), reveal striking differences in the precuneus and PCC. For example, the large negative cluster found at FA12 (white dashed oval) is absent at FA90—an effect that is not explained by the temporal signal-to-noise ratio (tSNR), which was found to slightly favor FA90 (Fig. S2). The tSNR assesses data quality: an equal tSNR warrants the use of a single statistical threshold, but, if anything, the slightly poorer tSNR at FA12 should yield fewer statistically significant voxels. This one-subject result is congruous with the group results and once again underlines that fMRI results are strongly dependent on APs: not only is one set of APs more or less sensitive to detect changes, but different sets reveal different attributes. Interestingly, the deviant voxels correspond closely to the areas of dissimilar tissue contrast at the FA pairs we compared.

Figure 3
figure 3

Activation maps.

Task-related activation maps showing statistically significantly responding voxels (activated in red, deactivated in blue) in two slices acquired with FA12 and FA90 at TE30 for the subject (S4) shown in Fig. 4. Activations are thresholded at p < 0.05 (FDR-corrected), corresponding to t-values of 2.38 and 2.32 for FA12 and FA90, respectively. White dashed line encircles an area negatively responding at FA12 but with no significant signal change at FA90.

Full size image
Figure 7
figure 7

Flip angle switching, data sorting and stimulation timing.

Full size image

Finally, we inspect the voxel-by-voxel signal properties. Figure 4A shows, for one subject, how the FA/TE settings affect tissue contrast. Quite saliently, for instance between FA12 and FA90, the appearance of CSF gradually changes from bright to dark. In the enlarged slice of Fig. 4B (FA12/TE30), the corresponding %ΔS traces for each voxel of the slice are plotted, color-coded by the frames of the slices at the top of the Figure, in panel A. The four traces in Fig. 4B correspond to FA12 and FA90 at both TE16 and TE30; TE30 with FA22 and FA50 are included in insets E and F.

Figure 4
figure 4

Anatomical and time-series results in subject S4.

(A) The anatomical tissue contrast and signal level max(S) of an oblique slice (D) imaged with six different FA/TE combinations as indicated. The signal (%ΔS) time courses (B) correspond to the different APs (color-coded as in (A) and overlaid on the 1343 voxels of the FA12/TE30 baseline image (red square in (A)). The stimulation period and the amplitude scaling (C) are shown for each voxel separately. The insets (E) and (F) enclose two areas whose voxels exhibit a wide range of responses (with cartoon schematics in (G–K)), ranging from classical BOLD to non-BOLD. The bands in (E) and (F) represent mean ± SEM values.

Full size image

In different voxels, the response behavior ranged from “Classical BOLD” to “non-BOLD”, illustrated schematically in Figs. 4G–K (the traces in this row are not data but cartoons). Classical BOLD traces (G, f(TE)) remain stable for all FAs at TE30 (blue, green, yellow and red) but are smaller for TE16 (dark violet and brown traces). Non-BOLD responses deviate from classical BOLD in terms of their FA/TE sensitivities, being either FA- and TE-sensitive (H, f(TE, FA)) to the extent of disappearing at some FA (I, f(TE, FA)), FA- and TE-insensitive (J), or displaying inverted polarity (K, f(FA)), etc. To guide the eye, FA-sensitive responses are emphasized with cyan shading in B: the larger the %ΔS discrepancy between FA conditions, the larger the shaded area.

These response patterns are exemplified in two enlarged, numbered voxel grids enclosing an area in the visual cortex (inset E) and a region of variable tissue contrast around the parieto-occipital sulcus (inset F). Additionally, some classical BOLD responses (type G) are indicated by the dashed boxes in B. Traces with similar %ΔS across APs (type J) were found in voxel 8 and f(TE,FA) traces with no abolished responses (type H) were found in voxels 6–7, 11–14 and 16–17. The non-BOLD responses in inset F include those of inverted polarities (type K) in voxels 43–44 and abolished signals in voxels 28–29 and 46–47. Responses of types H and I are especially disturbing, as they would evade the TE-based detection of non-BOLD responses.

Since activation-dependent displacement of unequally-weighted tissues obtained at a given set of scan parameters could produce sham BOLD changes, we further examined the differential effect of TE and FA on %ΔS, curious whether it could be related to changes in voxel composition and thus AP-dependent variation in image intensity (appearance of different tissues is influenced by both FA and TE; see Fig. 4A). Overall tissue contrast was poorest at FA50 and most enhanced at the lowest FA (FA12); correspondingly, we found that voxels located anterior and posterior to the sulcus (bright area in inset F) were highly responsive in the high-contrast setting (red, voxels 42–44) but the traces were flattened in the low-contrast setting (green). Interestingly, the high-contrast traces were of inverted polarities on opposing sides of the sulcus (e.g. voxel 33 vs. 43, FA12), suggesting that, in these locations, differentially bright tissue or CSF was displaced back and forth at the stimulus periodicity. These responses, similar to FA-dependent signals with opposite polarities at the borders of ventricles (Fig. S3), deviate from classical BOLD predictions and would be mistaken for BOLD changes in standard fMRI.

Fig. 5 shows, for the same slice as in Fig. 4, different components of signal change reflecting distinct tissue properties. The area of non-BOLD responses in Fig. 4 is clearly influenced by both M0 and T1 effects instead of showing purely T2*-driven responses. Fig. 5 also provides a semi-quantitative T1-map (bottom-right) as a way to locate CSF-containing voxels. Additionally, to illustrate the potential signal changes due to dynamic differences in tissue composition, Fig. 5 (bottom-right) provides the %ΔS expected for a hypothetical shift of tissue boundary by 50 μm from gray matter into voxels containing either CSF or white matter.

Figure 5
figure 5

Contributions of different tissue parameters on %ΔS signal change due to a simulated 50-μm shift in tissue boundary and putative location of CSF based on T1.

Top: The black squares show the locations of the two ROIs of Fig. 4. Bottom-left: The computed signal change due to the boundary of gray matter (GM) shifting toward CSF and white matter (WM) is shown for different FAs (with color convention as in Fig. 4). The initial concentration of gray matter [GM]0 is shown, the remainder being filled with the other tissue. The %ΔS was simulated with the signal equation (Methods) using parameters modeled from data obtained from voxels containing relatively pure samples of each tissue class. Bottom-right: T1 map with ranges of values colored to pinpoint possible CSF partial volume, dark purple indicating a low fraction of CSF in the voxel.

Full size image

Discussion

While a large share of fMRI studies is conducted with “conventional” APs (including TR: 2–3 s, TE: 25–40 ms, FA: Ernst angle of gray matter, full or nearly full brain coverage and about 3 × 3 × 3 mm3 voxel size at 3 T), some of the most rigorous work connecting fMRI signal and electrophysiology employ markedly different APs, e.g. with one or a few imaging slices, a TR of 0.1–1 s, a much higher spatial resolution and an FA chosen appropriately for the purpose, e.g. < Ernst angle to reduce flow effects. Thus, the fMRI signals of these two types of experiments derive from unequal biophysics, with potential ramifications for the interpretation of conventional fMRI.

If fMRI signal changes reflected similar physiological and biophysical phenomena at different FAs, the proportions of negatively vs. positively responding voxels should stay approximately equal across the FAs. In our data, this seems not to be the case as FA-dependent effects prevail. The responses appeared largely BOLD in the visual cortices, with non-BOLD features emerging in the other ROIs. Effects of M0 and T2* components remained largely unaffected as FA changed (TE theoretically only modulating the T2* component). However, at higher FAs, the T1 component influenced the total %ΔS. The largest observed deviations from BOLD responses coincided with the CSF-containing voxels (see T1-map in Fig. 5); another consideration is that for non-BOLD signals to arise, responses due to parameters other than T2* must be of at least comparable amplitudes. Detailed examination of single voxels indicated that all of the parameters had accumulative effects in some brain areas, whereas in other areas the total responses resulted from almost equally strong but counteracting contributions, especially at FA90.

PCC and precuneus are nodes of the “default-mode” resting-state network, which is deactivated during goal-directed behavior as has been shown by positron emission tomography15, electrophysiology16,17 and fMRI (negative BOLD signals)18. In our data, PCC and precuneus showed prominent non-BOLD effects (activations of inverted polarity in Fig. 4F are from those areas), meaning that the default-mode signal decrease in fMRI might include non-BOLD effects. Thus, the task-related dampening of fMRI signals in PCC might be a less specific marker for neuronal suppression than is often implied. Besides, these cortical areas have previously been observed to be prone to physiological noise derived from CSF fluctuations19; on the other hand, large draining veins running through CSF-containing voxels have been shown to bring about signals that anti-correlate with visual-cortex activations at zero or very short latency, insensitive to TE and not related to physiological noise20. In our data, many partially CSF-filled voxels contained physiological-like noise, but with non-zero mean. Physiological noise, comprising BOLD-like and non-BOLD-like components, is known to reduce the gain of SNR due to increased FA, with BOLD-like noise at maximum near the T2*21. The origin of the BOLD noise is not fully understood, but it is generally associated with hemodynamic and metabolic fluctuations in gray matter22. Here the noise did not differ between stimulation and control conditions in responding voxels (see SI and Supplementary Fig. S4) and was thus not driving the sustained signal changes.

Aside from physiological noise, inflow-related signals depend on APs. The inflow signals and their relation to APs are modulated differently depending whether the inflowing blood enters the imaging partition from another imaging partition or from outside the imaging volume: On one hand, if blood is displaced by an equal amount of blood that has been targeted by full-amplitude RF excitations within the latest TR, the MRI signal decreases. On the other hand, if the newly entering blood comes from outside the imaging volume and has not been subjected to RF pulses, signal increases (omitting any previous spin history effects). Our data are prone to both kinds of inflow effects; however, since our imaging volume was thin, the effect from the slow flow from neighboring partitions was assumed to be small and was not considered further. To set a reference with respect to inflow effects, FA12 gives the least corrupted measurements because, with small FA, the magnetization of spins in the imaging volume is more similar to that of the blood inflowing from outside the imaging volume, while for FA90, inflow effects increase ΔS for the opposite reason—the difference in magnetization between inflowing and present spins is larger. Thus, the additional signal increase for FA90 with respect to the lower FAs seen in a number of positively-responding voxels in our data can in principle be due to increased inflow. However, the ΔS should be amplified by the inflow effect for negative signal changes too, which was generally not the case and inflow effects cannot explain the ΔS in voxels where the low FAs result in the greatest ΔS.

Additionally, bulk head movement could induce, conceivably with a stimulus-related temporal pattern, tissue-contrast-based signal changes. While our data remain partly ambivalent to the nature of this effect, head movement parameters captured by motion correction do not support the bulk movement hypothesis and the measured voxel timecourses are not consistent with stimulus-locked head motion, e.g. high-contrast borders did not systematically fluctuate in a time-locked fashion, as they should if they were motion-induced.

Practically, scan parameters that result in as poor tissue contrast as possible will minimize functional contrast arising from conformational, partial volume-related non-BOLD changes. However, fMRI images with anatomical detail are more convenient for analysis. Unfortunately, without obtaining data at different FAs, activation maps are prone to misinterpretation. Task-related non-BOLD signals may themselves provide some information about neuronal activation, despite their putative non-specificity. In future, non-BOLD effects may be further discriminated by, for example, the use of flow-combating bipolar gradients in the direction of the assumed flow and especially by acquiring data at a higher spatial resolution.

The prominence of non-BOLD signals in our results suggests that their identification, e.g. by modeling different signal components as done here, would crucially improve the extent to which (i) positive or negative fMRI responses are related to neuronal activity and to electrophysiological data, (ii) various fMRI studies with different scan parameters are commensurate, (iii) functional networks are reliably detected and (iv) fMRI signal changes portray BOLD responses. However, our results are likely less central for traditional group-level analyses of activation locale (brain mapping) because the discrepancies tend to average out due to anatomical variability and even localized non-BOLD changes still represent functional signal. Also, the main positive BOLD response is bound to survive due to its strength with respect to e.g. T1-related effects and only some AP-dependent inconsistency will remain. Still, averaging out the details of the signal prevents gaining a more accurate picture of the complex spatio-temporal patterns of brain function, e.g. whole-brain stimulus-locked activations at individual level23. These types of analyses are essential for furthering the utility of fMRI in brain research.

The central implications of our results are that the BOLD vs. non-BOLD contributions to fMRI signals can be disambiguated by applying multiple APs in a single experiment and that without disambiguation, the BOLD assumption underlying interpretation of fMRI results—whether portraying brain dynamics, individual activation locale, or connectivity—is at stake.

Methods

Subjects

Thirteen healthy young adults (5 females, 8 males) reporting no history of neurological disorders participated in the fMRI experiments. The study was approved by the Ethics Committee of Helsinki and Uusimaa Hospital District and undertaken after the written informed consent of the subjects. Five subjects were excluded from analysis: head motion was considered too large (≥1.3 mm) in three subjects and the time series of two subjects were contaminated by large fluctuations throughout the brain that could not be explained by cardiac and respiratory cycles. The strict inclusion criteria were maintained as our purpose was to examine the method rather than to generalize brain reactivity beyond the cohort. The remaining eight subjects (4 females, 4 males) were all right-handed and within 20–35 years of age (mean 25).

Stimuli

During the control condition, the subjects fixated on a small cross centered on a gray (63 cd m−2) circle of 30° diameter. During the stimulation condition, the subjects continued to fixate on the cross but the gray background was replaced by a circular black- (0.6 cd m−2) and-white (126 cd m−2) “checkerboard” (Fig. 6) reversing at 8 Hz.

Figure 6
figure 6

Stimulation.

The visual stimulus consisted of a circular checkerboard (left) with the checks switching between black and white at a rate of 8 Hz. During the control condition, the subjects viewed a grey circle. The total eccentricity of the stimulus was 15°, but the scanner setup limited the view to a smaller portion, as illustrated by the red-shaded boundaries. Due to differences in head size and shape, each subject viewed a slightly different portion of the scene.

Full size image

The stimuli were delivered with Presentation software (Neurobehavioral Systems Inc., Albany, CA, USA), projected on a semitransparent back-projection screen with a data projector (Christie X3, Christie Digital Systems Ltd., Mönchengladbach, Germany) and viewed via a mirror placed outside the measurement coil.

Fig. 7 illustrates the stimulus timing: each stimulation period lasted 20 s, alternating with 56-s control periods; the runs started and ended in 28-s control periods. To ease the tedium of the experiment, subjects were instructed to silently cycle through the alphabet at their own leisurely pace, thinking of words beginning with each letter.

Data acquisition

FMRI was performed using a 3-T whole-body MRI scanner (Signa VH/i 3.0 T HDxt, GE Healthcare Ltd., Chalfont St Giles, UK) with a 16-channel head coil (MR Instruments Inc., MN) and gradient-echo echo-planar acquisitions (field-of-view = 19.2 cm2, matrix = 64 × 64). Fourteen (ten for Subject 1, nine for Subject 2) oblique axial slices were defined in-plane with the sinus rectus, crossing the subjects' calcarine sulcus and posterior cingulate cortex (slice thickness = 3 mm, spacing = 0 mm). The data were acquired with ASSET (SENSE parallel imaging technique) R-factor at 2 and second-order shimming was used. Two experimental runs were acquired (TE = 30 ms) using alternating, but fixed, FA pairs (90.0° and 12.6°, or 50.0° and 22.5°). The order of and within the fixed pairs was randomized across the 13 subjects. Additionally, a third set of dual-FA (90.0° and 12.6°) fMRI data were acquired at TE = 16.8 ms starting with Subject 3. The FA was switched exactly in the middle of the control condition following each stimulus block (see Fig. 7). Four trials (consisting of a pre-stimulus control period, the stimulus block and a post-stimulus control period) were acquired for each FA in the pair for a total of eight trials per run. Each experimental run for Subject 3 onwards contained 760 time points acquired at a repetition time TR = 800 ms; for Subject 1, TR = 1000 ms (608 time points) and for Subject 2, TR = 500 ms (1216 time points).

T1-weighted images were acquired for anatomical reference with a 3D SPGR sequence at a spatial resolution of 1.0 × 1.0 × 1.0 mm3.

Data analysis

The data for the six FA/TE conditions (90/30, 50/30, 22/30, 12/30, 90/16 and 12/16°/ms, decimals omitted) were extracted from the three functional runs with alternating FAs between trials (see Fig. 7 for a schematic of the imaging parameter sequencing, data sorting and stimulation timing). In the beginning of the image acquisition and after each FA change, a portion of data was discarded to allow the longitudinal magnetization to stabilize and the actual data used began with an abbreviated control block, followed by 20 s of stimulation and 28 s of the control stimulus. Thus, each of the dual-FA 608-s runs were split and sorted into two concatenated single-FA pseudo runs. Between-AP variation was controlled by within-run FA switching and the limited number of trials per condition was balanced by having two adjacent parameter combinations for each TE30 condition on the FA/TE plane, enabling examination of gradual changes in observables.

Statistical region-of-interest preprocessing (Figures 1,2,3, S1–S2 and Tables S1–S3)

The data were analyzed using the Analysis of Functional Neuroimages (AFNI)24 software and Matlab scripts. Concatenated runs were 320 points long (256 points for Subject 1 and 512 points for Subject 2). The concatenated runs were corrected for head motion, spatially registered to the first volume of the FA90 functional run and de-obliqued. A brain mask was applied to exclude any non-brain signals. Each trial within the concatenated runs was then individually intensity-normalized by dividing the time series in each trial by its mean. The time series in each processed, concatenated run was inspected and, where necessary, short segments containing head motion spikes too large to be corrected by the 3D-registration were excluded from further analysis.

Anatomical and EPI datasets were aligned; for one subject, the AFNI program align_epi_anat.py was used to improve the alignment. The anatomical datasets were then transformed into Talairach space. This transformation was subsequently applied to all processed EPI datasets and their activation maps.

Analysis of beta coefficients (Figures 1 and 2)

A quantitative measure of the relationship between activation amplitude expressed as general linear model coefficients for the (below) stimulus response function (β) and FA was determined by statistical analysis across individual subject datasets. Anatomically-defined region of interest (ROI) masks based on the Talairach atlas installed in AFNI (whereami program) were used to extract all voxels in seven ROIs in the left and right hemispheres. We plotted the voxel-wise beta coefficients for FA90 vs. FA12 and calculated the linear regressions (βFA90 = m βFA12 + c, where m = slope, c = intercept) of each of the resulting scatter plots using the polyfit function in Matlab for each region; only significantly active voxels were used, as described in the figure captions. Slope values for each fit in each ROI were collected for all subjects and two-tailed t-tests were computed to determine whether the mean slopes of the subject group in each ROI differed statistically significantly from 1. Under ideal BOLD conditions, m = 1; that is, FA has no effect on the beta coefficients.

Task-related activation analysis (Figure 3)

Individual subjects' activation maps for FA12 and FA90 (TE30, TR = 800 ms) conditions were computed using the 3dDeconvolve AFNI program. Head motion parameters and a reference waveform computed by convolving the stimulus time course with an impulse response function in the AFNI waver program (default parameters: delay time 2 s, rise time 4 s, fall time 6 s, post-stimulus undershoot 0.2 times the peak amplitude) were used as regressors. The 3dREML AFNI program was then used to account for temporal autocorrelations. Statistically-significant maps were thresholded to p < 0.05 (false discovery rate (FDR)-corrected).

Time course analysis (Figures 4 and S3)

The data presented are from a subject that was scanned at TR = 800 ms; the rest of the subjects were analyzed in the same way and the results shown are representative. The concatenated runs were corrected for motion, realigned to the first data point of the FA12/TE30 scan, first-order intensity trends were removed, %ΔS values from four trials were averaged and time courses temporally filtered. See the Supplementary Information for details about data processing and presentation. The resulting time courses of the six FA/TE combinations were plotted for every brain voxel in an image plane (plane 11 of 14 in Fig. 4, chosen to be far from the edge of the image stack, sinus rectus and ventricles to avoid artifacts), overlaid on the gray value of the FA12/TE30 scan for anatomical reference.

Tissue parameter modeling (Figure 5)

Tissue parameters M0, T1 and T2* were modeled from the TE30 data using Matlab least-squares optimization for the signal equation for spoiled gradient-echo acquisition25:

where M0, T1 and T2* are the longitudinal equilibrium magnetization (in arbitrary units), longitudinal relaxation time and effective transverse relaxation time. T2* was first fitted for , A being a constant, using the FA90 data. The parameters were extracted from the average signals across 8 s preceding the start (control) and end (stimulation) of each stimulation period, i.e. altogether 2 × 40 data points per AP condition. The parameters thus acquired were then used (without changing other conditions) to produce the %ΔS contribution of each parameter for each FA. Further, the putative location of CSF was mapped using voxel-vise T1 as a guide, assuming that a T1 of 1.5–2.0 s likely indicates the presence of some CSF in the voxel and a T1 > 2.0 s reflects a substantial partial volume of CSF.