Global signal modulation of single-trial fMRI response variability: effect on positive vs negative BOLD response relationship

Infunctional magnetic resonanceimaging (fMRI),the relationshipbetween positive BOLD responses (PBRs) and negative BOLD responses(NBRs) to stimulation ispotentiallyinformative about the balance of excitatoryand in-hibitorybrainresponses insensorycortex.In this study,weperformedthreeseparateexperiments delivering vi- sual, motor or somatosensory stimulation unilaterally, to one side of the sensory ﬁ eld, to induce PBR and NBR in opposite brain hemispheres. We then assessed the relationship between the evoked amplitudes of contralateral PBR and ipsilateral NBR at the level of both single-trial and average responses. We measure single-trial PBR and NBR peak amplitudes from individual time-courses, and show that they were positively correlated in all experi- ments. In contrast, in the average response across trials the absolute magnitudes of both PBR and NBR increased with increasing stimulus intensity, resulting in a negative correlation between mean response amplitudes. Sub- sequentanalysisshowedthattheamplitudeofsingle-trialPBRwaspositivelycorrelatedwiththeBOLDresponse acrossallgrey-mattervoxelsandwasnotspeci ﬁ callyrelatedtotheipsilateralsensorycorticalresponse.Wedem-onstrate that the global component of this single-trial response modulation could be fully explained by voxel- wise vascular reactivity, the BOLD signal standard deviation measured in a separate resting-state scan (resting state ﬂ uctuation amplitude, RSFA). However, bilateral positive correlation between PBR and NBR regions remained.WefurtherreportthatmodulationsintheglobalbrainfMRIsignalcannotfullyaccountforthispositive PBR – NBRcouplingandconcludethatthelocalsensorynetworkresponsere ﬂ ectsacombinationofsuperimposed vascular and neuronal signals. More detailed quanti ﬁ cation of physiological and noise contributions to the BOLD signal is required to fully understand the trial-by-trial PBR and NBR relationship compared with that of average responses. © 2016 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
BOLD functional magnetic resonance imaging (fMRI) is widely used in human neuroimaging to localise the spatial origin of brain activity in response to experimental tasks or stimuli. The majority of fMRI studies utilise the increase in BOLD signal that occurs following stimulus onset, relative to pre-stimulus or "resting" baseline levels, (termed the positive BOLD response, PBR) to infer that increased neuronal activity occurred in response to the stimulus. This assumption is supported by neurophysiology experiments in both humans and primates which have shown that the fMRI signal is an indirect, vascular correlate of increased neuronal activity in the form of local field potential and multi-unit activity (Heeger et al., 2000;Logothetis et al., 2001;Magri et al., 2012;Mukamel et al., 2005;Viswanathan and Freeman, 2007).
In addition, experimental stimuli often induce a decrease in BOLD signal below the baseline level, termed the negative BOLD response (NBR). For example, when stimuli are delivered unilaterally, such as images presented to one half of the visual field or the movement of one limb, a PBR is induced in the primary sensory cortex contralateral to the stimulation and a NBR is observed in the ipsilateral primary sensory cortex. This lateralisation of PBR and NBR to opposite hemispheres has been reported in primary visual (V1), motor (M1) and somatosensory (S1) cortices (Allison et al., 2000;Bressler et al., 2007;Hlushchuk and Hari, 2006;Kastrup et al., 2008;Newton et al., 2005;Tootell et al., 1998). In cross-modal sensory experiments, visual stimulation induces a NBR in auditory cortex (and vice versa) (Laurienti et al., 2002;Mayhew et al., 2013b), whilst painful stimulation induces a NBR in visual cortex (Derbyshire et al., 1997;Mayhew et al., 2013a). An NBR is not restricted to sensory cortex. Its occurrence has also been reported in posterior cingulate cortex, medial prefrontal cortex and intra-parietal regions (comprising the default mode network, DMN) in response to a wide range of cognitive tasks Northoff et al., 2004;Raichle et al., 2001;Spreng, 2012).
Despite the widespread observation of the NBR, its physiological origin and functional significance remain poorly understood and consequently the NBR is not widely utilised for brain mapping. The BOLD signal arises from a complex neurovascular coupling between cerebral blood flow (CBF), cerebral blood volume (CBV) and cerebral metabolic rate of oxygen consumption (CMRO 2 ) (Buxton et al., 1998) and consequently there are many potential combinations of relative changes in these parameters that can give rise to an NBR (Goense et al., 2012;Pasley et al., 2007;Shmuel et al., 2002).
Investigations into the physiological mechanisms of NBR can be broadly divided into those postulating a purely vascular origin, such as from the 'haemodynamic steal' of blood by an adjacent activated cortical region (Harel et al., 2002;Kannurpatti and Biswal, 2004;Olman et al., 2007;Puckett et al., 2014) and those reporting that NBR originates from local changes in metabolism, such as concurrent decreases in both CBF and CMRO 2 (Devor et al., 2007;Mullinger et al., 2014;Pasley et al., 2007;Schafer et al., 2012;Shmuel et al., 2002) potentially arising from decreases in local field potential neuronal activity in the NBR region (Boorman et al., 2010;Shmuel et al., 2006). A further unknown is whether a combination of these mechanisms may apply in some circumstances. However, in the case of NBR ipsilateral to the stimulus, an origin of blood steal from the contralateral PBR is unlikely to be the sole mechanism due to these regions occupying separate vascular territories (Smith et al., 2004;Tatu et al., 1998). In addition, NBR have further been reported due to local changes in blood volume in large cerebral veins (Bianciardi et al., 2011) and cerebrospinal fluid in the ventricles (Bright et al., 2014) of humans, as well as due to increases in neuronal activity and metabolism without a compensatory increase in CBF during hippocampal seizures in rats (Schridde et al., 2008).
Given the wide variety of contexts in which NBRs have been observed an important open question is whether comparable or different mechanisms underlie their generation and whether NBRs represent different physiological processes in different scenarios. Of key interest is ascertaining the circumstances in which NBR provides a useful neuroimaging marker of cortical inhibition, whether this reflects increases in local inhibitory neuron activity or decreases in excitatory input (Cauli et al., 2004;Ferbert et al., 1992;Lauritzen et al., 2012) given that increases in inhibition have been shown to lead to both increases (Enager et al., 2009;Pelled et al., 2009) and decreases (Devor et al., 2007) in BOLD signal.
Evidence for an association between NBR and measures of cortical inhibition comes from reports that individuals with higher baseline concentrations of the gamma-aminobutyric acid (GABA) inhibitory neurotransmitter in anterior cingulate cortex have been shown to display the largest magnitude (absolute value) of NBR in the same region (Northoff et al., 2007). Increased NBR magnitude ipsilateral to median nerve stimulation has been linked to increases in the perception threshold for electrical stimuli delivered to fingers of the contralateral hand (Kastrup et al., 2008;Schafer et al., 2012), which is thought to form a behavioural manifestation of ipsilateral cortical inhibition. Additionally, single-trial NBR amplitudes have been shown to correlate with the power of simultaneously recorded 8-13 Hz EEG oscillations in the somatosensory cortex (Mullinger et al., 2014), providing further evidence of a link between NBR and inhibitory neuronal processes (Jensen and Mazaheri, 2010;Mathewson et al., 2011).
The NBR displays many of the stimulus-response properties that characterise the PBR. The average magnitude of the NBR increases with increasing stimulus intensity and duration (Klingner et al., 2010;Shmuel et al., 2002), suggesting that NBR reflects neuronal inhibition required to optimise task performance, by reducing sensitivity and allocation of processing resources to the unattended or irrelevant part of the sensory field.
In addition to the average response, single-trial responses can be measured from the peak amplitude of each trial's fMRI timecourse. Trial-by-trial amplitude variability is commonly treated as "noise" by conventional general linear modelling analyses, which assume a consistent amplitude response across trials, but has been widely reported to contain information which is behaviourally relevant to the dynamics of network processing Scheibe et al., 2010).
Taking into account all the caveats stated above, here we hypothesize that important functional information may be contained in the relationship between contralateral PBR and ipsilateral NBR amplitudes at the single-trial level, such as regarding the balance of excitation and inhibition within a cortical network. Therefore we will investigate how the single-trial amplitudes of PBR and NBR, evoked concurrently by the same individual stimulus, relate to each other. The degree of single-trial PBR and NBR amplitude variability, the frequency with which they display signal polarity which is opposite to the average response polarity, and the relationship between PBR and NBR are currently uncharacterised.
Here, we use unilateral visual, motor and somatosensory stimulation to induce contralateral PBR and ipsilateral NBR in primary visual (V1), motor (M1) and somatosensory (S1) cortices respectively, to allow an assessment of the generalisability of findings across these three sensory modalities. First, we investigate the unknown relationship between natural single-trial variability in the PBR and NBR amplitude and compare this to the relationship between the mean PBR and NBR response amplitudes with increasing stimulus intensity. Second, we investigate the unclear role that global fMRI signals and resting-state haemodynamic signal properties play in modulating single-trial fMRI signals and the PBR-NBR relationship.

Materials and methods
Three fMRI experiments were performed in different subject cohorts. visual: 14 subjects (4 female, 27.8 ± 5.4 years); motor: 17 right-handed subjects (7 female, 26 ± 4 years); and somatosensory: 18 right-handed subjects (8 female, 27 ± 3 years). All data were collected with approval from the local ethics committee and informed consent was obtained from all subjects. These data were initially collected for other purposes, however all data contained the primary experimental condition (unilateral stimulation of primary sensory cortex) required to answer the scientific questions which we pose in this work. Table 1 summarises the key parameters for the three fMRI experiments.

Stimulus paradigms
Visual Black/white checkerboard stimuli were presented to the left visual hemi-field at either high (100%, HC) or low (25%, LC) contrast levels in a pseudo-randomised order. Individual experimental trials consisted of a one-second duration checkerboard presentation, with phase reversal at 500 ms. Trials were separated by an inter-stimulus interval (ISI) of either 16.5, 19 or 21 s Porcaro et al., 2010). Five, 11-min experimental runs were acquired in each subject, each consisting of 17 trials per contrast, providing 85 trials per contrast in total. Subjects were instructed to continually fixate on a central cross that was displayed throughout. Within the same scanning session a three-minute resting-state scan was also acquired, during which subjects were instructed to lie still, keep their eyes open and think of nothing in particular.

Motor
Subjects performed an isometric contraction of a compliant rubber bulb, by opposing the thumb with the first two fingers of their righthand. An increase in the contraction force applied to the rubber bulb increased the pneumatic pressure inside a rubber tube, which was translated into an analogue electrical signal by in-house-built electronics and recorded by a Ni-DAQ (National Instruments) (van Wijk et al., 2009). Subjects were instructed to perform the isometric contraction for 5-s at one of two force levels, 10% or 30% of maximum voluntary contraction (MVC). Real-time visual feedback of contraction force was provided via a projector screen to inform the subject of task performance. Subjects were required to match the vertical position of a centrally displayed force indicator to a target force level of either 10% or 30% of MVC. Individual's MVC was measured prior to the experiment using a hand dynamometer. Two 12-min experimental runs were used to acquire 60 trials at each force level in a pseudo-random order, with an ISI of 5, 7 or 9 s. A six-minute resting-state scan was also acquired, during which subjects were instructed to lie still, keep their eyes open and think of nothing in particular.

Somatosensory
Median-nerve stimulation (MNS) was applied via two electrodes placed on the right wrist using square wave pulses of 0.5 ms duration (Digitimer DS7A). The stimulation current amplitude was set just above the subject's motor threshold (range 2.6-7 mA, mean 4.6 ± 1 mA) to cause a small thumb distension. In a single experimental run MNS was applied at 2 Hz in 40 trials, each comprising a 10 s stimulation period followed by a 20.5-21 s duration ISI (Mullinger et al., 2013;Mullinger et al., 2014).

fMRI data acquisition
All data was acquired on a Philips 3T Achieva MRI scanner using an eight channel receiver-array head coil. The three experiments were acquired with the following imaging parameters.
To facilitate co-registration of functional data, whole-head, T 1weighted anatomical images with 1 mm isotropic resolution were acquired on each subject. The subject's cardiac and respiratory cycles were continuously recorded throughout the motor and somatosensory experiments using the scanner's inbuilt pulse-oximeter and pneumatic bellows which were attached to the subject's index finger of the left hand and positioned around their ribcage respectively, but not during the visual experiment.

fMRI preprocessing
Physiological data were temporally aligned with the BOLD data acquisition. RETROICOR (Glover et al., 2000) was used to reduce physiological noise in the motor BOLD data. The BOLD data from the motor and visual experiments were then motion corrected using MCFLIRT (Jenkinson et al., 2002).
Somatosensory DABS data were separated into BOLD and ASL data sets for subsequent analysis. The BOLD data were then physiologically corrected using RETROICOR (Glover et al., 2000). The use of background suppression, which nulls the contribution of the background tissue signal to the ASL images and removes cardiac and respiratory sources, meant that this physiological noise correction was not required for the ASL data (Garcia et al., 2005). All data were then motion corrected using FLIRT (Jenkinson et al., 2002) and linearly interpolated to an effective TR of 2.6 s. Tag-control ASL image pairs were then subtracted to create cerebral blood flow (CBF) data. The BOLD-weighted image pairs were averaged to produce mean BOLD-weighted data.
Prior to statistical analysis, BOLD and CBF data from all experiments were spatially smoothed (5 mm FWHM Gaussian kernel), high-pass temporally filtered (100 s cutoff) and registered to high-resolution anatomical and MNI standard brain images using FSL 4.1.8 (www.fmrib.ox. ac.uk/fsl). At this point three subjects were excluded from the motor and five from the MNS datasets due to excess motion (N3 mm).

fMRI data analysis
The sequential fMRI GLM analyses performed in this study is illustrated schematically in Fig. 1 and described in detail in the following section.
A series of three separate GLM analyses were performed on each dataset for the respective purposes of: 1) PBR and NBR ROI identification, followed by response amplitude measurement and PBR-NBR correlations; 2) investigation of the contribution of trial-by-trial PBR modulations and global signal modulations to BOLD response variability across the whole-brain; and 3) investigation of how voxel-wise variability in the BOLD response amplitude can be explained by intrinsic BOLD signal parameters such as temporal fluctuation and global signal contributions.

GLM analysis 1: PBR and NBR ROI identification
A GLM analysis (Fig. 1, GLM 1) was performed to identify those brain regions where stimulation induced, on average, a significant increase or decrease in BOLD signal compared to the passive baseline periods between stimuli. For each experiment, first-level design matrices were constructed for each subject using regressors of stimulus timings to model the main effect of stimulation. Visual: 1) 100% (HC) contrast trials and 2) 25% (LC) contrast trials; motor: 1) 10% contraction trials and 2) 30% contraction trials; and somatosensory: all MNS trials. All stimulus regressors were convolved with the canonical double-gamma haemodynamic response function.
To control for potential differences in heart-rate and depth of respiration between trials and between experimental conditions unrelated to the neuronal response to the task, the respiration-per-volume-time (RVT) (Birn et al., 2008) and the variation in the heart-rate interval (HRI) (Chang et al., 2009;de Munck et al., 2008) were computed from the physiological data, for all experimental runs of the motor and somatosensory data. These data were down-sampled to form continuous time-courses with one sample point per TR period and convolved with the respiration-response function (Birn et al., 2008) and cardiacresponse function (Chang et al., 2009) respectively. The RVT and HRI regressors and the six motion parameters characterising head translation and rotation were incorporated into the first-level design matrix as confounds of no interest.
First-level statistical analyses were performed using FILM modelling (Woolrich et al., 2001) in FEAT 6.00 (www.fsl.ox.ac.uk). Positive and negative contrasts were computed on all regressors. For the somatosensory data, first-level results were combined across all subjects using fixed-effects to identify group PBR and NBR regions. BOLD and CBF somatosensory datasets were analysed separately. For the visual and motor datasets, for each subject first-level results were combined across all experimental runs using a second-level, fixed effects analysis to calculate an average response per condition and per subject. These results were then combined across all subjects at the third, group-level using a fixed effects analysis. Separately for the visual and motor experiments, third-level contrasts combined statistics across both conditions (visual HC and LC; motor 10% and 30% MVC) to identify group PBR and NBR regions which were common across stimulus intensities as well as subjects. All BOLD Z-statistic images were threshold using a Z N 2.3 and cluster corrected significance threshold of p b 0.05. The CBF Z-statistic images were less stringently threshold at Z N 1.6, uncorrected due to the inherently lower contrast-to-noise ratio of CBF data than BOLD data.

ROI timecourse extraction and single-trial response measurement
For the motor and somatosensory data, respiratory (RVT) and cardiac (HRI) trends were removed from each voxel's BOLD timecourse by scaling the RVT and HRI regressor by their respective GLM regression coefficient and subtracting them from the data. For each dataset, we defined group-level regions of interest (ROI) from the cortical areas that exhibited a significant contralateral (c) PBR or ipsilateral (i) NBR to the stimulus, from which response timecourses were extracted for each experimental run. ROIs were defined from a 3 × 3 × 3 voxel cube centered on the peak-statistic BOLD voxel in PBR and NBR regions of V1, M1, S1 and primary auditory (A1) cortex.
Subsequently for each subject and each ROI, single-trial BOLD response timecourses were extracted based upon stimulus timings. To define the baseline signal level for each experimental run, the average BOLD response timecourse was calculated for each run and the mean of the final three time-points (two in the motor experiment due to shorter ISI) was used as the baseline measure. The visual/motor/ somatosensory data contained 17/30/40 trials per run respectively, meaning that the baseline correction was based upon 41/60/120 time points respectively to ensure a robust correction.
Single-trial BOLD response timecourses were converted to percent signal change relative to this baseline value. For each subject, the timeto-peak of the PBR and NBR signal change was then found from the mean response timecourses. Single-trial PBR and NBR amplitudes were measured as the peak signal change within the time window of: mean time-to-peak ± (2*TR). These measures were used to assess the linear correlation of single-trial PBR-NBR amplitudes for each subject as follows: visual: cV1-iV1, cV1-iA1 for both HC and LC visual trials; motor: cM1-iM1 for both 10% and 30% MVC trials; and somatosensory: cS1-iS1. An equivalent analysis was also performed on the somatosensory CBF data, using the BOLD ROIs to also extract CBF responses, allowing direct comparison between CBF and BOLD responses from the same spatial location. For each subject, single-trial positive (PCBF) and negative (NCBF) CBF response amplitudes were measured from the cS1 and iS1 ROIs respectively and their linear correlation computed.
For each experimental run the standard deviation of the BOLD signal time series was calculated on a voxel-wise basis, a measure we term 'task fluctuation amplitude' (TFA), analogous to the 'resting state fluctuation amplitude' (RSFA, (Kannurpatti and Biswal, 2008;Kannurpatti et al., 2012)). TFA was averaged across runs to obtain a TFA value per voxel, per subject. The voxel-wise matrices of TFA were then concatenated across subjects and used as a covariate in subsequent group-level GLM analysis. Here for each voxel, the vector of TFA values forms a regressor containing between-subjects variability in the effect size. This regressor therefore models how the response at each voxel varies between individuals depending on the TFA, effectively removing the component of the group response (in a voxel wise manner) that correlates with TFA. Such group-level voxel-wise covariates have been previously used to correct for differences in vascular reactivity (Murphy et al., 2010) and also for structural confound modelling across subjects (Oakes et al., 2007;Spisak et al., 2014).
Additionally, for each experimental run the BOLD time series were averaged across all brain voxels giving a single 'task global signal' timecourse (TGS) for use in a subsequent GLM.

Analysis of resting-state data acquired during visual and motor experiments
Resting-state BOLD data were preprocessed equivalently to the task data (physiological noise correction (for motor experiment only), motion correction, temporal filtering (high pass N 0.01 Hz) and 5 mm spatial smoothing) and physiological trends of RVT and HRI regressors were removed by linear regression. Subsequently, the standard deviation of the resting-state BOLD signal time series in each voxel (RSFA, (Kannurpatti et al., 2012)) was calculated. RSFA is thought to provide a measure of a voxel's vascular reactivity as it has been shown to positively correlate with the amplitude of the voxel's BOLD response to both motor and breath-hold tasks (Kannurpatti et al., 2012). Additionally, the BOLD time series were averaged across all brain voxels to give a single 'resting-state global signal' timecourse (RSGS). The linear correlation between the RSGS and each voxel's resting-state BOLD timecourse was calculated using the voxel-wise Pearson correlation coefficients as a measure of the similarity between each voxel's BOLD signal and the RSGS, which can be thought of as approximating the degree to which each voxel contributes to the resting global signal. Voxel-wise brain volume matrices of this correlation with RSGS and RSFA were separately concatenated across subjects and used as a covariate in subsequent group-level GLM analyses.
GLM analysis 2: Investigating trial-by-trial response modulation across the whole brain GLM analyses were performed to identify those brain regions where the trial-by-trial variability in the BOLD response to the stimulus was correlated with the contralateral PBR amplitude (Fig. 1, GLM2). Separately for the visual, motor and somatosensory data-sets, the singletrial measurements of cPBR amplitude were mean subtracted and used to form regressors to model the parametric modulation of the BOLD response to the stimulus. This regressor had an amplitude that was either positive or negative for each trial depending on whether each single-trial BOLD response was larger or smaller, respectively, than the mean response. Two separate GLM analyses were then performed ( Fig. 1), incorporating the following additional regressors into the design matrices used in the initial analyses: GLM2A single-trial cPBR modulations; GLM2B single-trial cPBR modulations and the TGS timecourse (which was not convolved with any HRF).
The purpose of including these two GLMs was to dissociate between the contribution of sensory network PBR modulations and the contribution of global signal modulations to BOLD response variability in primary sensory cortex. We aim to discover whether including the TGS regressor in the model accounts for all of the response modulation in sensory cortex, as well as throughout the rest of the brain. The effect of interest is the total variance explained by each model, not the specific loading on its constituent factors. Therefore in GLM2B the single-trial cPBR modulations and the TGS regressors were not orthogonalised. By adopting this approach, the GLM can either allow the TGS to explain all the cPBR variability or separate the response variability into sensory network and global components. Group-level, fixed effects statistical maps were calculated using the same methods stated above.
For the somatosensory data, an additional GLM analysis was performed to investigate correlations between trial-by-trial PBR fluctuations and the CBF response, to test whether the CBF response was modulated in a comparable manner to the BOLD response. This GLM therefore analysed the CBF data using the same first-level design matrix that was used to model the somatosensory BOLD data.
GLM analysis 3: Group-level analyses to investigate how task response modulations are explained by vascular reactivity and global signal fluctuations Finally, we investigated the extent to which single-trial variability in the amplitude of the BOLD response to stimulation can be explained at the voxel-wise level by intrinsic, predominantly vascular parameters of global signal or signal fluctuation amplitude (Fig. 1, GLM3). For each dataset, additional third-level, fixed-effects analyses were performed on the subject's second-level contrast of the single-trial PBR modulation regressor (result of GLM 2). In separate analyses we incorporated (where available): 1) RSGS, 2) RSFA, or 3) TFA as a voxel-wise 4-dimensional (4th dimension = subject) covariate-of-no-interest. These covariates were included to investigate how voxel-wise variability in the GLM fit to the cPBR modulations was explained by fMRI signal properties derived from the task data and also from an independent resting-state dataset. Using the TFA GLM (Fig. 1, GLM 3C) enabled us to estimate the effect of intrinsic voxel signal properties on response modulation in the somatosensory dataset where no resting-state data was available, and to compare results between the TFA GLM and the RSFA GLM (Fig. 1, GLM 3B) in the visual and motor experiments.

Group average fMRI responses
Significant cPBR and iNBR were observed in the relevant primary sensory cortex in each experiment. The main effect fMRI responses to visual, motor and somatosensory stimulation in V1, M1 and S1 are illustrated in Fig. 2A, B and C. For the somatosensory data, the group conjunction ROIs show that significant contralateral PCBF and ipsilateral NCBF responses were observed in primary somatosensory cortex and that these were highly comparable to the regions exhibiting PBR and NBR (Fig. 2C). In addition, MNS also induced significant fMRI responses in M1 but fMRI timecourses were extracted from ROIs centered on S1, where the peak response voxel was located. During somatosensory stimulation, additional NBRs were observed in small regions of contralateral S1, adjacent to the PBR. Additionally, visual stimulation evoked a NBR bilaterally in A1 ( Fig. 2A) as well as medial prefrontal cortex and posterior cingulate cortex (Mayhew et al., 2013b). Here we focus our investigation on the relationship between contralateral PBR and ipsilateral NBR response regions as a marker of the balance of functional activity between directly stimulated and unstimulated sensory network regions, as there is a reasonable body of evidence to suggest that the ipsilateral NBR doesn't arise from purely vascular origins and may reflect cortical inhibition. We include the A1 NBR to provide a comparison to a region outside of the directly stimulated sensory network, but to constrain our analysis to a reasonable number of regions, we omit the NBRs from the contralateral S1 and midline regions that displayed the lowest signal-to-noise ratio. Fig. 2D shows that in both V1 and M1 the average magnitude (absolute value) of both cPBR and iNBR increased with increasing stimulus intensity for the motor and visual experiments. This negative correlation, such that more positive mean PBR occurred concurrently with a more negative mean NBR, is in agreement with previous studies (Klingner et al., 2010;Shmuel et al., 2002), suggesting that the magnitude of ipsilateral inhibition increases with the magnitude of increasing contralateral excitation.
Modelling variations in the depth of subject's breathing and heartrate in motor and somatosensory data showed that the BOLD signal was significantly correlated with these physiological fluctuations in widespread areas of grey matter, including primary sensory cortex. Fig. 3 shows that the RVT regressor explained BOLD variance bilaterally in parietal, visual and motor cortices, in agreement with previous work (Birn et al., 2008). The HRI regressor explained BOLD variance in the dorsal surface of the brain, bilateral parietal areas, the precuneus region adjacent to the sagittal sinus and in the junction between the rostral visual cortex and the cerebellum as previously reported (Chang et al., 2009).

Relationship between single-trial PBR-NBR amplitudes
We observed substantial trial-by-trial variability in both PBR and NBR amplitudes in each experiment. The mean ranges of response amplitudes pooled across conditions and subjects were: visual: cV1 PBR −0.5%-2.67%, iV1 NBR −1.58%-0.89%; motor: cM1 PBR −0.6%-1.81%, iM1 NBR −1.52%-0.91%; and somatosensory: cS1 PBR −0.27%-1.83%, iS1 NBR −1.13%-0.49%. The BOLD response in ipsilateral primary sensory cortex was negative on average, but not consistently negative on a single-trial basis. On average across all subjects, a bilateral PBR was observed in 24 ± 5% of visual, 32 ± 5% of motor and 16 ± 6% of somatosensory trials. These bilateral PBR occurred in the trials that exhibited  Fig. 3. Group-level correlations between BOLD signal and RVT (i) and HRI (ii) for the motor (A) and the somatosensory (B) experiments. All maps are shown cluster corrected at a threshold of p b 0.05 on the MNI standard brain. Note that the somatosensory DABS data was not acquired over the whole head, but from an imaging slab centred on the somatosensory cortex. the largest contralateral PBR amplitudes, as outlined in greater detail in subsequent analyses.

B) Somatosensory A) Motor
We observed a significant (p b 0.05) positive, linear correlation between single-trial cPBR and iNBR amplitudes in the majority of individual subjects for each of the data sets. Specifically, this correlation was seen in the visual HC data for 12 out of 14 subjects; visual LC = 11/14 subjects; motor 10% MVC = 11/13 subjects; motor 30% MVC = 11/13 subjects; and somatosensory = 10/11 subjects. We also observed that contralateral visual PBR was significantly (p b 0.05) correlated with auditory NBR at the single subject level (in the visual HC data for 11/14 subjects; and in the visual LC data for 10/14). Fig. 4 displays a boxplot for each condition and each experiment that shows the distribution and median value of the correlations between single-trial cPBR and iNBR amplitudes for the individual subject's data. No instance of a significant negative correlation was observed in any subject. Therefore at the single-trial level, we observed that larger magnitude (more positive) cPBR occurred concurrently with smaller magnitude (less negative) iNBR. Fig. 5 displays example data taken from the subject with the median correlation strength in each experiment. In each case a significant positive linear relationship between cPBR and iNBR was observed. For the visual and motor data (Fig. 5A, B, C) the low (LC, 10% MVC) and high (HC, 30% MVC) intensity stimulus conditions are plotted in blue and red respectively. We observe that the negative correlation between the average cPBR and iNBR amplitudes across stimulus intensities is preserved at the individual subject level, as illustrated by the black lines connecting the mean values in Fig. 5A, B and C.

Relationship between single-trial cPCBF-iNCBF and cPBR-cPCBF amplitudes
In the somatosensory data, the relationship between contralateral and ipsilateral CBF responses was comparable to that observed in the BOLD data. A significant positive correlation (p b 0.05) between single-trial cPCBF and iNCBF amplitudes was found in four subjects (R range = −0.12-0.82, median = 0.18). No significant (p N 0.05) negative cPCBF-iNCBF correlations were observed. A significant positive correlation between single-trial PBR and PCBF amplitudes was also found in seven subjects (R range = − 0.04-0.83, median = 0.26). No significant negative correlations between cPBR and cPCBF were observed. Additionally we observed that the correlation between a subject's single trial cPBR and iNBR amplitudes was itself significantly linearly correlated (p b 0.001) with the strength of the relationship between that subject's cPBR and cPCBF response amplitudes (Fig. 6). The strength of the cPBR-iNBR relationship was not significantly correlated with the mean or standard deviation of a subject's cPBR or cPCBF amplitude suggesting this effect does not arise from differences in the data signal-to-noise ratio between subjects. This implies that the degree of inter-hemispheric correlation between PBR and NBR is strongly related to the local relationship between the BOLD and the perfusion signals in contralateral S1. This result aids the interpretation of global signal effects discussed below.
Spatial pattern of single-trial response modulations across the whole brain GLM2A (Fig. 1) showed that single-trial contralateral PBR amplitude was significantly positively correlated with the BOLD response to the stimulus across widespread brain regions in all three experiments (Fig. 7). These results indicate that the positive correlation with cPBR amplitude (Figs. 4 & 5) was not spatially specific to the iNBR region, but occurred globally in the grey matter. However, the statistical maps show that in each experiment the most significant single-trial modulations occurred bilaterally in the directly stimulated primary sensory cortex, V1, M1 and S1.
This finding is further illustrated in Fig. 8 which displays the results from the GLM including regressors formed from both single-trial cPBR amplitude and TGS (GLM 2B). Here, significant correlations with single-trial cPBR modulations were observed only in bilateral sensory    Fig. 5. Single-trial correlation between cPBR and iNBR amplitude in subjects with the median R value (from Fig. 4) for: A) visual, cV1-iV1; B) visual, cV1-iA1; motor, cM1-iM1; and C) somatosensory, cS1-iS1 data. Blue and red colours (A-C) represent singletrial data points for low and high stimulus intensities respectively. The black lines drawn between the centres of the linear fit lines in A-C indicate the negative direction of the correlation between mean response amplitudes and the low and high stimulus intensities.  6. Positive relationship between inter-hemispheric BOLD coupling and contralateral PBR-PCBF coupling in the somatosensory experiment. A significant positive linear fit between the single-trial cPBR-iNBR correlation and the single-trial cPBR-cPCBF correlation across the group is seen. cortex (Fig. 8A), as the global signal regressor accounted for response modulations over the rest of the brain (Fig. 8B). The remaining correlation between cPBR amplitude and the BOLD signal in bilateral primary sensory cortex suggests that the respiration, heart-rate and global brain signals, cannot account for all of the response variability in the ipsilateral sensory cortex. The differences between the maps shown in Fig. 8A and B illustrate the distinction between sensory network specific signal fluctuations and global fluctuations over the rest of the brain.
We further observed that variability in the CBF response to somatosensory stimulation correlated with the single-trial amplitude of the cPBR over widespread cortical grey matter regions (Fig. 7D). These correlations occurred with lower statistical significance than in the BOLD data, as was also observed for the single-trial correlations, consistent with the lower contrast-to-noise ratio of CBF data (Aguirre et al., 2002;Detre and Wang, 2002).

Effect of group-level covariates upon single-trial cPBR response modulation
In the visual and motor experiments, where resting-state data was available, we found that including group-level voxel-wise covariates of RSGS or RSFA (Fig. 1, GLM 3A & B) reduced both the statistical significance and the spatial extent of single-trial correlations with the cPBR amplitude. This reduction is evident by comparing the whole brain correlations shown in Fig. 7 with the reduced area of correlation shown in Fig. 9. Including RSGS (Fig. 9A) reduced the significance of the correlation with cPBR (note the lower Z-stats in Fig. 9 compared with Fig. 7), but response modulations across widespread brain regions were still observed. In comparison, after voxel-wise removal of RSFA (Fig. 9B), the significant single-trial response correlation was reduced to a number of focal cortical areas, with the most significant modulations remaining in bilateral primary sensory cortex. In the visual data, significant modulations remained in inferior lateral visual areas and anterior auditory cortex. In the motor data, modulations remained in anterior cingulate cortex and precuneus. Highly comparable results were observed when using TFA (GLM 3C) as a voxel-wise covariate (Fig. 9C) with significant correlations with cPBR amplitude remaining in bilateral primary sensory cortex in all three experiments. Due to the strong similarity of the TFA results in all three experiments, we suggest that the RSGS and RSFA results obtained for the visual and motor experiments could also be extended to the somatosensory experiment.

Discussion
Improving understanding of the origin and functional significance of BOLD signal fluctuations during both task execution and rest is of vital importance to neuroimaging research. Here, we addressed this issue using unilateral stimulation in three sensory modalities. We demonstrate that the average amplitudes of contralateral PBR (cPBR) and ipsilateral NBR (iNBR) were negatively correlated as previously reported and thought to result from increased ipsilateral inhibition occurring concurrently with increased contralateral excitation (Klingner et al., 2010;Shmuel et al., 2002). However we do not observe the same relationship between single-trials where the natural fluctuations in cPBR-iNBR amplitudes were positively correlated. Further investigation revealed that the single-trial variability in the cPBR was positively correlated with the amplitude of the BOLD response to the stimulus across the grey matter in the whole brain, not just in iNBR regions. Voxelwise contribution of the global brain signal (TGS) accounted for the global grey matter component of single-trial response variability, however positive correlations between bilateral primary sensory cortex remained. The voxel-wise temporal standard deviation of the BOLD signal, whether measured at rest (RSFA) or during the task (TFA), was found to account for a large proportion of the global response modulation, suggesting that much of the single-trial variability arises from intrinsic vascular effects. However, the bilateral modulation between primary sensory cortices remained when these additional confounds were modelled. Therefore we find that the positive correlation of single-trial BOLD responses between cPBR and iNBR regions of primary sensory cortex could not be completely accounted for by heart-rate or respiratory variability, global signal modulations, RSFA or TFA. Other potential sources of this relationship could perhaps be attributed to common neuronal responses to stimulation between the sensory cortices, as discussed below.

Potential sources of response variability
Here, we present results from three experiments each using a different subject cohort, stimulus modality, stimulus paradigm and fMRI acquisition. These factors all contribute to the between-subject and between-experiment variability in our measurements. However, we consider that these differences strengthen the findings of this work, as they demonstrate the generality of our results. We focussed on robust fMRI responses in primary sensory cortex and found that the number of subjects exhibiting the positive, linear relationship between singletrial cPBR and iNBR amplitudes, and the strength of this relationship, was highly comparable across different stimulus conditions and modalities (Figs. 4, 5 & 7). This somewhat surprising finding, given the documented negative correlation of mean cPBR-iNBR amplitudes with stimulus intensity (Fig. 2D and (Klingner et al., 2010;Shmuel et al., 2002)), led us to investigate the contribution of local and global intrinsic vascular signals to BOLD response variability measured from the bilateral sensory network across all experiments. We now discuss the possible sources of variability which may explain the observed results.

Physiological and "noise" signals
The resting-state BOLD signal, measured from any given voxel in the brain, possesses a degree of positive correlation with every other brain voxel (Aguirre et al., 1998;Desjardins et al., 2001). For example, right M1 will exhibit strong positive correlation with all other regions of the motor network and weaker positive correlation with all other brain regions. This arises because of shared BOLD signal variance between all brain voxels, which is not restricted to those sharing functional or anatomical connections, a phenomenon called the 'global brain signal' (Aguirre et al., 1998;Desjardins et al., 2001;Murphy et al., 2009). This signal is usually estimated by averaging BOLD signal across all voxels in the brain. It has been shown that physiological variability in heart-rate, depth and rate of respiration, and the level of expired carbon dioxide induce changes in CBF that modulate the BOLD signal independent from changes induced by neuronal activity. These factors, along with head motion and equipment noise, are all components of the conglomerate "noise" which contributes to the variance in BOLD signal fluctuations and the global brain signal (Birn, 2012;Birn et al., 2006;Chang and Glover, 2009b;Power et al., 2012;Van Dijk et al., 2012). Therefore the global signal can be simply conceptualised as a summation of local-network and long-range vascular and neuronal signals, plus noise. In the current study, we observe that global CBF fluctuations display comparable behaviour to BOLD signals by correlating with the contralateral response to stimulation (Fig. 7). We also observe that the strength of the cPBR-iNBR correlation is related to the coupling between cPCBF and cPBR responses (Fig. 6). Under the assumption that global signal modulations contribute substantially to the positive correlation between cPBR and iNBR, this relationship further suggests that modulations in contralateral CBF responses also bear a relationship to the global signal, i.e. modulations in contralateral CBF correlate with iNBR.
The global brain signal is conventionally regressed out of restingstate data in order to remove common variance between brain regions and improve the spatial specificity of functional connectivity measurements . However recent work has cast doubt on the validity of this procedure of blindly removing an unknown mixture of BOLD signal components which are of uncertain origin (Gotts et al., 2013;Murphy et al., 2009;Saad et al., 2012). The exact contribution of the physiological, vascular and neuronal components to the BOLD signal is dependent upon both within-and between-subject factors, e.g. voxel proximity to major blood vessels, individual physiology and the experimental paradigm. Therefore evidence combined across different experiments may help to disentangle the various contributions.
Our paradigms employed different stimulus durations (visual 1 s, motor 5 s, somatosensory 10 s), intensities (visual contrast, motor contraction force), mean ISIs (visual 18 s, motor 7 s, somatosensory 20 s) and task-demands (passive for visual and somatosensory, active response for motor) ( Table 1). As different TR were used in each experiment the sampling of physiological noise and therefore its removal by RETROICOR would also vary across datasets. The consistency of our findings, despite this between-experiment variability, suggests the single-trial global modulations are not dependent upon experimental conditions, but represent a general response characteristic. Additional analysis showed that cPBR and iNBR amplitudes were independent of the trial position within the experimental session as well as the pre-and post-stimulus BOLD signal amplitude in each experiment (data not shown). Whilst it could be suggested that the significance and spatial extent of the correlation between the BOLD response and single-trial cPBR amplitudes increased with task demand (lowest and least widespread Z-statistics for brief, passive visual stimulation compared with strongest and most widespread Z-statistics in the motor experiment, Fig. 7), we propose that the most important factor was the duration of the ISI which determined the sampling of the global signal. The single-trial cPBR amplitudes reflect a combination of both local and more widespread neuronal, and vascular signals, plus noise. We therefore hypothesize that the strength of the correlation between the cPBR and the BOLD response across the whole brain arises due to the cPBR amplitudes effectively sub-sampling the modulation of the global-brain signal from the stimulated region, at the frequency of the stimulation. Therefore for datasets with the shortest ISI, more samples of the global signal (relative to the total number of sample points) are modelled resulting in better fits of the cPBR amplitude with the global BOLD signal over the whole brain. This is reflected by our observation that the statistical significance of the response modulation increases with decreasing ISI such that the stronger effect is observed in the motor (peak Z-stat = 23.4; ISI = 7-9 s) than the visual (peak Z-stat = 13.9; ISI = 16.5-21 s) or somatosensory (peak Z-stat = 16.8; ISI = 20 s) experiments. The contribution of the global signal to the local BOLD response variability and the correlation between cPBR and iNBR can be judged by comparing Fig. 8A with Fig. 7. This shows a large reduction in the extent and strength (Z-statistic) of the correlation within each sensory network (Fig. 8A) when the global signal is modelled. However, we still found that the single-trial cPBR-iNBR responses were more strongly correlated with each other than with the global signal, as demonstrated by the remaining bilateral sensory network correlation when the global signal was included in GLM2 (Fig. 8A). Taken together these results suggest that a number of signal modulations, of different neurophysiological origins acting over different spatial scales (global, whole network or local), are linearly superimposed within the stimulated primary sensory cortex, contributing to the responses observed.

Intrinsic vascular reactivity
We also found that the global, voxel-wise variability in the strength of the BOLD signal correlation with the cPBR can be strongly predicted by RSFA, a measure of vascular reactivity derived from an independent resting-state experiment. The RSFA of voxels in the motor cortex has been previously shown to predict the amplitude of the BOLD response in those voxels to a finger tapping task and a hypercapnic breathhold (Kannurpatti and Biswal, 2008;Kannurpatti et al., 2012). Here, we extend this work outside of primary sensory cortex by demonstrating that the amplitude of the BOLD response to stimulation in any brain voxel is strongly predicted by its RSFA. We also show that the BOLD signal standard deviation calculated from the task data (TFA) provides a highly comparable calibration of the voxel-wise BOLD response (Fig. 9). In both cases, after accounting for RSFA, global BOLD signal correlations with cPBR are minimal, but bilateral correlations with the amplitude of the iNBR regions remain.
In comparison, the level of similarity between a voxel's resting-state BOLD signal and the resting-state global brain signal (RSGS) did not predict the voxel-wise BOLD response to stimulation (Fig. 9A). It is of further interest to note that, for the data shown in this work, removal of the global brain signal from the task time series data does not transform the single-trial cPBR-iNBR correlation from positive into negative. The recomputed cPBR-iNBR single trial correlations after voxel-wise removal of the global brain signal timecourse using linear regression (Macey et al., 2004) are displayed in Fig. S1. Global signal removal has the general effect of decreasing correlation values, in agreement with previous functional connectivity work (Murphy et al., 2009), and introducing some negative correlations. However, as is clear from Fig. S1, the median value remains positive, except between V1 and A1. The observation of a negative correlation only arising between areas which belong to different brain networks (as opposed to cPBR vs iNBR within the same network) raises the possibility that in this case global signal removal has accounted for confounding effects in the BOLD data and reveals the functional coupling between two regions. However, inverse correlations after global signal removal have to be interpreted with caution as the preprocessing has been reported not only to introduce a negative bias but also to bias correlations differently in different areas depending on the true, but unknown, underlying relationship structure (Gotts et al., 2013;Murphy et al., 2009;Saad et al., 2012). A better understanding of the global brain signal's constituents would enable full separation of specific network-related BOLD signal fluctuations from the conglomerate noise, allowing negative single-trial PBR-NBR coupling to emerge.
Taken together, these results demonstrate that global signal modulations do not account for all of the positive single-trial cPBR-iNBR relationship. TGS would have been expected to account for voxel-wise variability in the strength of the BOLD signal correlation with the cPBR in the event that this response variability was driven only by modulations of the global brain signal that were independent of task state (i.e. resting or stimulation). We do not observe such behaviour which suggests that not all of the positive cPBR-iNBR single-trial correlation can be attributed to non-sensory specific, global vascular effects. As we currently lack an accurate, informed method of separating the global signal from local effects and neuronal from vascular "noise" effects, it is not possible to discern the exact underlying single-trial PBR-NBR relationship. Recent work using a dual-echo BOLD recording offers a promising way of addressing the issues of removing global signal confounds by measuring this "noise" at short echo times (Bright and Murphy, 2013;Kundu et al., 2012). Further insight into neurovascular coupling and separation between different haemodynamic signals, such as capillary signals from venous draining effects could be provided by acquisition of fMRI data at high spatial-resolution (Iranpour et al., 2015) and analysis omitting spatial smoothing (Mikl et al., 2008;Scouten et al., 2006).

Neuronal contributions
We have shown that fluctuations in the global BOLD signal contribute to positive single trial cPBR-iNBR correlations, but the mean cPBR-iNBR amplitudes are negatively correlated. We suggest that when comparing average BOLD responses between brain regions the effect of the global signal, which is not time-locked to the stimulus presentation, is averaged out leaving only correlations driven by the stimulus. The negative correlation of mean cPBR-iNBR is driven by increasing stimulus intensity (Fig. 2D). It therefore reflects local activity and is thought to result from increased ipsilateral inhibition occurring concurrently with increased contralateral excitation (Klingner et al., 2010;Shmuel et al., 2002). The variability in response amplitude between trials arises from neuronal, as well as vascular sources. For example variability may be related to differences in subject's arousal, attention and spontaneous brain network activity which contribute to variations in bottomup afferent input and top-down modulation of stimulus processing (Debener et al., 2006;Mayhew et al., 2012Mayhew et al., , 2013aScheibe et al., 2010). Whilst we have shown that non-neuronal physiological signals, largely indexed by the global signal modulations, form a large contribution to the positive single trial cPBR-iNBR correlations (Figs. 4, 5 & 7), these vascular effects cannot fully explain the response variability and the cPBR-iNBR relationship observed (Figs. 8 & 9). The question of which neurophysiological factors explain the remaining BOLD signal variance is an important topic of ongoing research. Therefore the findings of this study have consequences for understanding the spatial specificity of BOLD signal fluctuations, the interpretation of differences and similarities in BOLD responses between networks and of the composition of the BOLD signal.

Spontaneous neuronal fluctuations
In addition to the global fluctuations in fMRI signal which have a non-neuronal vascular or noise-related origin, spontaneous fluctuations in neuronal activity are widely observed. During the passive state of resting in the scanner, BOLD signals exhibit spontaneous fluctuations which are highly correlated between brain regions that share a common functional specialisation (Biswal et al., 1995;Lowe et al., 1998). These fluctuations correspond to the activity of resting state networks and electrophysiological studies in primates suggest that a neuronal component underlies this BOLD signal variance that is exploited for functional connectivity measurements (Leopold and Logothetis, 2003;Miller et al., 2009;Vincent et al., 2007), complementing magnetoencephalography studies in humans (Brookes et al., 2011;de Pasquale et al., 2010). Further work has shown that coherence in neuronal activity is not restricted to functionally specialised networks, but occurs over widespread brain areas (Scholvinck et al., 2010). This, combined with reported relationships between the global signal and EEG vigilance measures (Wong et al., 2013), provides evidence that a component of the global BOLD signal is driven by neuronal activity.
In the light of these previous findings, it is not possible for the present study to distinguish between haemodynamic and neuronal BOLD signal components which underlie the positive correlation between single-trial cPBR and iNBR amplitude. However, by accounting for head-motion and implementing RETROICOR, RVT and HRI correction we have taken considerable steps, in line with current best practice in the literature, to discount the possibility that physiological noise confounds can explain the response modulations observed. Previous work has shown how physiological correction improves the estimation of both task fMRI activation and resting-state functional connectivity Glover, 2009a, 2009b;Khalili-Mahani et al., 2013). Cardiac and respiratory signals were not recorded in the visual experiment, but the observation of similar results between the visual experiment and the motor/somatosensory experiments, for which physiological correction was performed, suggests that these factors do not dominate the effects we observe. The observation of similar modulations between BOLD and CBF signals suggests a haemodynamic, rather than artefactual, origin, but whether those CBF changes are induced by underlying fluctuations in neuronal activity or reflect purely vascular mechanisms cannot be conclusively determined. Given that we have accounted for several vascular factors and yet correlations, although weaker, still persist in sensory areas (Figs. 8,9 & S1), the presence of neuronal mechanisms must be considered.
It has been suggested that coherent fluctuations in neuronal activity are not specific to the resting state and persist during task performance, therefore forming an intrinsic feature of brain activity (Smith et al., 2009). Several reports have shown that task responses are superimposed upon this so-called "background" activity (Arieli et al., 1996;Fox et al., 2007). Previous work by Fox et al., found that the trial-by-trial variability in the amplitude of the BOLD signal evoked in left M1 by a right hand button press could be accounted for by variability in the BOLD signal from the right M1 Fox et al., 2006). Furthermore these authors link this bilateral BOLD signal variability to behavioural variability in the force of the button press. This finding is cited as evidence for superposition of task-evoked and spontaneous activity within the motor network. Interestingly, despite employing a unilateral motor task, Fox et al. do not report an ipsilateral NBR. Instead they state that ipsilateral M1 showed signs of small amplitude PBR on average. Their method of calibrating the cPBR M1 response, by subtraction of a scaled version of the ipsilateral M1 BOLD signal, yields a similar effect to that demonstrated in the present study. We report that the largest (smallest) amplitude contralateral BOLD responses are observed concurrent with the largest (smallest) amplitude ipsilateral BOLD responses respectively. The Fox et al. subtraction shifts all contralateral response trials towards the mean amplitude. However, Fox et al. only investigated the response modulation within the motor network and did not consider either global or physiological contributions. Here we demonstrate that the trial-by-trial amplitude variation in the cPBR is correlated with the BOLD signal fluctuation across widespread brain regions, and furthermore even after accounting for global and physiological signal contributions, a positive relationship remains between the responses of contralateral and ipsilateral sensory network regions which exhibit opposing average BOLD responses (cPBR vs iNBR), providing complementary results to those presented by Fox et al. In conclusion, contralateral PBR and ipsilateral NBR regions that exhibit negative correlation between their mean response amplitude show positive correlations between their trial-by-trial amplitudes. We observe that this positive relationship largely reflects a global modulation of response amplitude and complements recent findings of widespread activity during basic visual stimulation (Gonzalez-Castillo et al., 2012), thus providing new insight into the spatial extent and temporally dynamic complexity of brain responses during simple tasks. Whilst we have shown that the global BOLD response amplitude variability can be explained by physiological (e.g. RVT and HRI) and vascular effects (e.g. RSFA, TFA), substantial unexplained trial-by-trial variance remains in the stimulated sensory cortical areas which preserves the positive cPBR-iNBR correlation. This is a surprising result, given the known negative correlation between mean cPBR and iNBR amplitudes, and further investigation is required to distinguish the neurophysiological sources of inter-trial response variability. Promising approaches to this include use of multi-echo MRI data to fully separate spatially homogenous sources of BOLD signal noise from network specific activity (Bright and Murphy, 2013;Kundu et al., 2012) and characterisation of haemodynamic and neuronal markers of changes in brain state (Saka et al., 2010).