Spatially Extended fMRI Signal Response to Stimulus in Non-Functionally Relevant Regions of the Human Brain: Preliminary Results

1. INTRODUCTION

Analysis of data sets in functional neuroimaging via magnetic resonance (fMRI) is based on the haemodynamic response function (HDR) and the phenomenon of blood-oxygenation level dependency (BOLD). In this paper we demonstrate that for one dataset, a spatially varying HDR model can be fit across the whole brain, leading to a coherent global spatial structure of positive and negative responses in the HDR parameters. We present the results from this dataset as interesting preliminary findings, providing evidence for the existence of a neurally driven blood redistribution mechanism.

Haemoglobin is paramagnetic (has positive magnetic susceptibility in an external magnetic field) when de-oxygenated, but is diamagnetic (has negative susceptibility) when oxygenated. Therefore, the magnetic susceptibility of blood is altered depending on the relative concentration of deoxy-haemoglobin [1]. Brain activity requires the local consumption (metabolisation) of oxygen, forcing an increase in the flow rate and volume of cerebral blood that is much greater than the rate at which oxygen is metabolised, thereby leading to an increase in venous blood oxygenation. An increase in concentration of oxy-haemoglobin causes an increase in BOLD signal (reduced attenuation due to de-oxyhaemoglobin) and the change in this concentration over time defines the BOLD HDR (Fig. 1).

In fMRI brain activation studies, the (BOLD) HDR is the core phenomenon being measured with the assumption that it represents underlying brain activity. The initial dip in oxygenation levels, detectable at high field strengths [2-5] is thought to be due to highly local micro-vascular depletion, though there is some debate as to whether it is robustly detectable with fMRI [6-8]. This is followed by the large over-compensation in oxygenated blood delivery to the site of activation detectable at lower field strengths. The HDR then falls a little more slowly than during the rise and subsequently dips below the baseline for a while (negative overshoot) before gradually returning to baseline. The source of the negative overshoot is not fully understood, but physiological mechanisms that have been proposed [9] include: 1) increased oxygen extraction rate due to sustained post-stimulus elevation in de-oxyhaemoglobin; 2) cerebral blood flow overshoot after stimulus ends; or 3) delayed recovery of venous cerebral blood volume.

Negative BOLD signals not attributable to the initial dip or negative overshoot have been observed in numerous activation studies [10, 11], and understanding the properties of such negative signals has been the object of considerable research [12-14]. Furthermore, the physiological interpretation of these negative BOLD signals is a debated topic, and suggested mechanisms include: 1) local vascular compensation or blood stealing [12]; 2) local decrease in neural activity [13, 14]; or 3) neurally controlled constriction of vessels as described in the discussion of Smith et al. (2004) [13]. As noted therein, if such a neurally controlled phenomenon genuinely exists, then it might be described as blood sharing rather than blood stealing since it reflects a global redistribution of blood through neuronally constricted and dilated vessels (haemodynamic balancing). Considerable human and animal research provides evidence that key physiological ingredients for a system capable of performing haemodynamic balancing exist [15, 16], with a range of possible neuronal triggering mechanisms [17-19].

When a single global delay is used for fitting the HDR, the magnitudes at voxels are inefficiently estimated by use of delays either too short or long, depending upon the difference between the particular delay imposed and the voxel’s ‘true’ delay. We have shown [20] that within ‘conventionally activated’ areas (i.e. regions of significant stimulus-responding activation) that there is up to a 4 s difference in delay of the HDR between voxels. In addition we found that magnitude estimates across conventionally activated areas are highly correlated with delay estimates. However, there is a large potential gain in statistical sensitivity for detecting stimulus-responding regions when fitting variable delay models to voxels with small response magnitudes, i.e. voxels outside the conventionally activated areas. In regions with smaller response magnitudes, the time to reach the peak of the response is relatively short. The ensuing discrepancy between the assumed fixed delay model and the true delay leads to systematic under-estimation of magnitude (which is already small). Fig. (2) illustrates the problem by showing the fit of a 10 s delay Poisson curve (defined later in the methods) to a 5 s delay Poisson curve with time to repetition (TR) of 1 s; the fitted curve’s magnitude is only 45% of that for the true response.

Given the underestimation from inappropriate delay assumptions, and the typically lower signal-to-noise ratios (SNRs) beyond these brain regions predicted to be activated by the stimulus, the signal from small magnitude voxels may not emerge from the noise to be detected as significantly active using conventional voxel-by-voxel thresholding-based tests, e.g. using Gaussian random fields [21]. (For conventionally restricted purposes of detecting direct stimulus-responding activation this may be a good thing.) Furthermore, the voxels that we show later to have a predominantly negative response (with delays of the order of that for the negative overshoot) would be missed by HDR models focused on matching the earlier positive response; the negative overshoot can continue up to 25 s beyond the end of the stimulus period - see e.g. Buxton (2002) [22] p. 418.

Previous attempts to detect spatial variation of delay across voxels [11, 23-27] were unable to detect the global response described here, probably because either 1) the form of delay variation they were capable of detecting was restricted by convolution assumptions and the need to choose canonical basis functions; 2) they were not looking for the possibility of extended regions of negatively responding voxels; 3) they did not directly estimate baseline levels of signal that were uncontaminated by negative overshoots from the response of a previous stimulus; 4) the effects would not appear significant using voxel-by-voxel thresholding-based tests. It is therefore informative to attempt to directly fit and assess a variable delay HDR model across all voxels in the brain to determine whether additional (possibly smaller magnitude) patterns of response emerge in both predicted and unpredicted brain areas; and if such patterns do emerge, to examine the possible implications on the underlying biology as well as on inferences derived from conventional activation and brain connectivity studies.

2. MATERIALS AND METHODOLOGY

Parameter Estimation

We fitted a HDR model directly to the average response curve at each voxel. This method contrasts with the usual approach of fitting a version of the stimulus function convolved with a canonical HDR or set of HDR basis functions to the voxel time courses. Avoiding convolution is preferable here in order to avoid the imposed restrictions on the shape of the HDR [20]. The HDR was modelled by a reference function h_s(τ), where τ was the time from the onset of activity and s was the voxel location index. The fitting procedure has four steps:

1. Estimate the baseline signal at each voxel by least squares using the two minute resting condition occurring immediately before and after the listening experiment and subtract the baseline estimate from the voxel data.
2. Average the data across cycles to produce a ‘mean response’ cycle for each voxel.
3. Fit a curve h_s(τ) to each voxel's mean response cycle using least squares estimation.
4. Create summary maps for the parameter estimates of the fitted curves.

In Kornak et al. (2009) [20] we compared a range of reference functions (gamma, constrained gamma, polynomials, and double gamma) for this direct fitting procedure. Here we forewent the options therein and instead fitted a simpler Poisson curve using the optimisation function nlregb in Splus (http://www.insightful.com/products/splus/default.asp).

The Poisson curve is defined as:

(1)  $h_s\left(\right)=\left(\begin{array}{cc}{k}_s\frac{{}_s^1{e}^{{}_s}}{\left(1\right)!}& =1,2,...,T1\\ 0& =0\end{array}\right)$

where T is the number of timepoints in the stimulus cycle and λ_s is the mean of the underlying distributional form for the Poisson curve at voxel s , providing for a measure of delay (in terms of scan intervals) at each voxel.

The Poisson curve, as defined here, shifts the functional form of the standard Poisson probability distribution one unit to the right; the first timepoint's response value is set to zero and each subsequent timepoint τ corresponds to the value of the distributional form at τ - 1. This shifted function is then scaled by a constant k_s defined to be the response magnitude at voxel s . An implicit assumption of this model is that at time 0 of each cycle there is no response due to the stimulus.

The Poisson curve imposes an inherent constraint that the mean and variance are equal; the Poisson distribution is defined by the single parameter λ that equals both the mean and variance. Consequently, the Poisson curve imposes a specific non-linear relationship between delay and width (when defined in terms of FWHM or SD), and is therefore not invariant to changes in scale of the time axis. However, we found that the Poisson curve adequately modelled the main rise and fall of the HDR, while providing reliable fits for the design and dataset considered here. Examples of Poisson curve fits (Fig. 3). Summary parameter estimates from the Poisson curve fits θ_s =(k_s,λ_s) (magnitude and delay, respectively) of h_s(τ) were used to represent the nature of activation in each voxel.

It should be noted that non-stimulus locked activity and noise will be reduced relative to the signal in the empirically averaged cycles. Any systematic effects will be reduced because they are not synchronised to the stimulus and therefore do not behave additively when averaging, but rather tend to cancel each other out. The Poisson curve least squares fit implicitly assumes that noise exists in the average cycles, where the noise effectively includes all non-stimulus-locked forms of signal because non-stimulus-locked biological sources and signal noise cannot be separated. In addition, there will also be a component of signal due to the scanner noise that occurs every 2.33 s. However, since this signal will cycle at the frequency the data are acquired it should manifest as a constant additive effect. This constant additive effect should have been subtracted prior to the fitting of the Poisson because of the subtraction of the average of the two minutes of baseline at the beginning and end of the experiment.

The flexibility of allowing a spatially varying delay within a HDR model allows the Poisson curve to fit different types of responses. This flexibility of fit led to improved sensitivity for detecting spatial pattern in the magnitude and delay parameter maps, thereby enabling the findings reported here.

3. RESULTS

The baseline signal acquired in the 2 minute of ‘scanner noise only’ at the beginning and end of the experiments was very stable and proved important for obtaining stable estimates of haemodynamic response across all voxels. When comparing mean baseline estimates for the two minutes before and after experiments across all voxels, the mean absolute difference of ‘before’ and ‘after’ was 0.08 with a standard deviation of 0.06. This compares with an average absolute empirical signal change of 1.2 across all voxels (based on the maximal signal change of the 15 timepoints in the averaged cycles).

Fig. (4) shows parameter estimates for one coronal slice containing primary auditory cortex. For those voxels within core sites of the auditory cortex, a positive correlation between the Poisson curve magnitude parameter estimates (magnitude map) in panel (a) and the the corresponding map of estimates of delay (delay map) in panel (b) was seen over the first 10 s. These ‘conventionally active’ voxels can be seen in panel (a) as the set of high response magnitude voxels (dark region) on the left side of the map. However, beyond the principle regions of activity it was noticeable that the fitted magnitude parameters were spatially clustered into positive and negative magnitude voxels and that the delay map was inversely correlated with the pattern of the magnitude map; this pattern would be inconsistent with fitting to random noise. (The negative voxels, i.e., negative relative to baseline, would be detected as negative activation if found significant in a conventional fMRI analysis.) Since the Poisson curve fitting procedure from which the delay was obtained is identical for both positive and negative magnitudes, there was no asymmetry in the delay fitting procedure about zero magnitude, and therefore we would not expect the negative correlation in Fig. (4) between the magnitude (a) and delay (b) maps.

Additional maps were generated from the original parameter estimates to illustrate this phenomenon more clearly. Panel (c) of (Fig. 4) displays a version of the magnitude map thresholded at zero with black voxels corresponding to positive magnitudes and white to negative. Panel (d) shows the delay map thresholded at its median value (black voxels consist of delays above the median), and panel (e) gives the sum of the two thresholded maps (c) and (d). The map in panel (e) of (Fig. 4) tends to homogeneous mid-grey (value =1). This is in contrast to what would be expected in the absence of a global haemodynamic effect. Under a null-hypothesis of no correlation between the length of the delay and whether a voxel is positive or negative responding you would expect to see an even split of positive and negative responding voxels above and below the median delay. Therefore, you would expect to see approximately 25% black ( = 2 ), 25% white ( =0 ) and 50% mid-grey (=1) in panel (e) of (Fig. 4). This result provides some insight into the reason for thresholding the delay at the median. If there is a 50-50 split of positive and negative responding voxels (which is approximately what we see), thresholding at the median of the delay is the optimal way to try and detect an association in a 2 × 2 table, i.e., we are looking for whether one 50-50 split can be assumed independent of another 50-50 split.

Fig. (3) shows typical empirical and fitted HDRs for the four main classes of voxel responses that we observed: high-magnitude positive- and negative-responding voxels, and low-magnitude positive- and negative-responding voxels. The first thing to note is that the empirical HDRs appear to exhibit both a positive and negative component. The Poisson curve is either completely above or below zero depending on the value of k_s and is therefore only capable of fitting either the positive or negative component of the observed HDR.

The positive-responding voxels had a short delay compared with the negative ones. The fitted mean delays of positive responding voxels generally lay in a range matching the positive rise and fall of the HDR for conventionally responding voxels, whereas the fitted mean delays for the negative-responding voxels more closely matched typical delays for the negative overshoot; up to 28 s post stimulus e.g., Buxton (2002) [22] p. 418. Hence, over the time range of the delay estimates, there was a strong negative correlation between magnitude and delay. This contrasts with the smaller positive correlation seen when only the positive magnitude fits are considered. Two further observations are worthy of note. First, the observed pattern of effects occurred across the whole slice and was stimulus-locked, i.e. was temporally correlated with the stimulus paradigm. Second, the effect extended across brain regions that would not have been predicted a priori and in a conventional analysis would not have been detected as activity.

It is straightforward to confirm statistically that this global BOLD effect is not simply spatially correlated random variation. The test was performed based only on the voxels in the right half of the slice in Fig. (4) to avoid the left side that contained the peak activity. This region contained 643 voxels, whereas the whole slice contained 1236 voxels. If the null hypothesis were true, i.e. if the region only contained time courses formed from random noise (possibly temporally correlated) with zero mean, then the probability of obtaining a negative magnitude would be 0.5 for every voxel regardless of delay estimate.

Table 1 displays a 2 × 2 table of zero-thresholded magnitude against median-thresholded delay. A two-sided Fisher's exact test emphatically rejected the null hypothesis of no association between the sign of the magnitude and the size of the delay ( p<10^-15 ). The cross-ratio measure of association for 2 × 2 contingency tables [30] was 499, implying a very high degree of association.

Fig. (5) displays a scatter plot of delay against magnitude for all voxels, i.e., the complete set of slices covering the whole brain. The plot further indicates the strong asymmetry of the delay fits about the zero magnitude line. Those voxels designated as ‘conventionally activated’ comprised only a small proportion of voxels, and correspond to the dark region of high-magnitude positively responding voxels on the left side of (Fig. 4) (a). This set of voxel responses contrasts with the pattern expected if the remaining parameter estimates were based on white or temporally correlated noise, which would be approximately symmetric about the zero magnitude line. A further interesting feature of the plot is the group of low-magnitude negative responding voxels with a delay below 2 s. Spatially the largest group of these voxels lies to the bottom right of the main auditory area of activation -- see the light region below the main activated area in the delay map (b) of (Fig. 4). These voxels could potentially be interpreted as reflecting voxels with a dominant initial oxygenation dip. We defer comment on this to the discussion.

There are a few voxels with delays greater than 35 s because the Poisson curve fitted delay is not constrained to be less than 35 s. In addition, the fitted curve is only evaluated up to 35 s and therefore there is no signal beyond that to constrain the fitted curve to return quickly to zero. If the observed data across the cycle looks like it is a steady rise or decline it could be fit by a Poisson curve with a delay longer than the stimulus cycle.

Note that most of the voxels would not reach a level of conventional activation based on a “massively univariate” modelling approach, such as that used in SPM (i.e. voxel-by-voxel modelling followed by thresholding). The evidence that the voxels in this dataset do not consist of noise is primarily driven by the observation that most positive responding voxels have short delay and most negative ones have long delay (there is a consistent pattern across voxels that could not occur by chance). The highly significant Fisher's exact test result (which is still conclusive after conventionally activated voxels are removed) indicates that most of the brain is responding in a synchronised way to the stimulus. Because the Fisher's exact test as used here is essentially a global test of activation, we cannot say equivocally that any specific voxel is activated, but the very strong relationship in the 2 x 2 table indicates a very large proportion of voxels are “responding” in synchronicity with the stimulus.

Some alternative approaches within the Poisson curve fitting framework were examined, in order to relax the forced dichotomisation into all-positive or all-negative responses. First, the data were modelled using the sum of two Poisson curves, such that one was constrained to have a positive magnitude and the other a negative one. However, this formulation contained some non-identifiability, in that it was possible for the two curves to fully or partially cancel each other out. This led to the algorithm not converging in a high proportion of voxels, with the magnitude parameters tending to infinity for both the negative and positive directions. Even when the positive and negative delays were constrained to be within different ranges (according to expected positive delay and expected negative overshoot delay ranges), convergence problems persisted. The second approach was to force the fits for all voxels to be positive. Surprisingly, problems of non-convergence again occurred in many voxels, particularly those which had negative magnitudes in the original unrestricted fits. This suggests that although there appeared to be components of both positive and negative responses in most voxels, they could still be characterised in some sense as predominantly positive or negative responding.

4. DISCUSSION

5. CONCLUSION

The results in this paper for a single-subject dataset at high SNR show a global structure of HDR parameter estimates that is inconsistent with fitting to noise (even if the noise is temporally correlated). These results suggest the possibility of a haemodynamic balancing effect occurring outside of conventionally activated areas to compensate for highly localised oxygenated blood delivery. We are unable to determine any possible confounding effect which may produce these results artefactually, and assume that the effect becomes discernible through the extra flexibility obtained when employing a suitable spatially varying model for the HDR with adequate estimates of baseline BOLD signal. We acknowledge the limitation that the study is based on a single subject. It is a primary goal of this paper to encourage independent replication and further testing in terms of alternative sensory stimuli (visual, motor etc.) and different subject populations. If properties of this spatio-temporal effect are found to be generalisable, it raises far-reaching issues for techniques aiming to optimally localise activation and quantify brain connectivity in fMRI.

REFERENCES