Neuronal synchrony and the relation between the BOLD response and the local field potential

The most widespread measures of human brain activity are the blood oxygen level dependent (BOLD) signal and surface field potential. Prior studies report a variety of relationships between these signals. To develop an understanding of how to interpret these signals and the relationship between them, we developed a model of (a) neuronal population responses, and (b) transformations from neuronal responses into the fMRI BOLD signal and electrocorticographic (ECoG) field potential. Rather than seeking a transformation between the two measures directly, this approach interprets each measure with respect to the underlying neuronal population responses. This model accounts for the relationship between BOLD and ECoG data from human visual cortex in V1–V3, with the model predictions and data matching in three ways: Across stimuli, the BOLD amplitude and ECoG broadband power were positively correlated, the BOLD amplitude and alpha power (8–13 Hz) were negatively correlated, and the BOLD amplitude and narrowband gamma power (30–80 Hz) were uncorrelated. The two measures provide complementary information about human brain activity and we infer that features of the field potential that are uncorrelated with BOLD arise largely from changes in synchrony, rather than level, of neuronal activity.


Introduction
Most measurements of activity in the living human brain arise from the responses of large populations of neurons, spanning the millimeter scale of functional magnetic resonance imaging (fMRI) and electrocorticography (ECoG) to the centimeter scale of electro-and magnetoencephalography (EEG and MEG). Integrating results across methods is challenging because the signals measured by these instruments differ in spatial and temporal sensitivity, and in the manner by which they combine the underlying neuronal population activity [1][2][3]. Differences in scale can be partially bridged by bringing the measurements into register. For example, EEG and MEG sensor data can be projected to cortical sources subject to constraints from simultaneously recorded fMRI data [4] or from independent fMRI localizers [5], and ECoG electrodes can be aligned to a high resolution anatomical MRI image [6] and compared to the local fMRI signal.
Yet even when electrophysiological and fMRI data are spatially registered, striking differences in the sensitivity to stimulus and task are often observed, indicating differences in how neuronal responses contribute to the measured physiological signals. For example, the fMRI BOLD signal and EEG evoked potentials differ in which brain areas are most sensitive to visual motion (area MT+ with fMRI [7] versus V1 and V3A with EEG [8]). Within the same visual area, fMRI and sourcelocalized EEG evoked potentials can show different effects of task in similar experimental paradigms, such as the effect of spatial attention on the contrast response function (additive in fMRI [9], multiplicative in EEG [10]). Even when the spatial scale of the two signals is approximately matched at acquisition, such as ECoG electrodes and fMRI voxels (both at ~2 mm), systematically different patterns of responses can be obtained, such as compressive spatial summation in fMRI versus nearly linear summation in ECoG steady state potentials (but not ECoG broadband signals) [11]. Such fundamental functional differences cannot be explained by numerical measurement-tomeasurement transformations.
Rather, these differences must reflect the fact that the measurements are based on different aspects of the neural population response. To explain the differences in measurement modalities requires a computational framework that derives each of these signals from the neuronal responses.
One approach toward developing such a framework has been to measure the BOLD signal and electrophysiological signals simultaneously, or separately but using the same stimulus and task conditions, and to ask how features of the electrophysiological response compare to the BOLD signal. This approach has revealed important patterns, yet after several decades of careful study, some apparent discrepancies remain. A number of studies comparing band-limited power in field potential recordings to the BOLD signal have shown that increases in power between 30 and 100 Hz (gamma band) are more highly correlated with BOLD amplitude than power changes in other bands [12][13][14][15][16][17]. Yet power changes in this band do not fully account for the BOLD signal: very large power changes can occur in the gamma band without a measurable BOLD signal change [18,19], and power changes in lower frequency bands can be correlated with the BOLD signal independently of power changes in the gamma band [20][21][22][23]. It therefore cannot be the case that field potential power in the gamma band is a general predictor of BOLD, even if the two measures are often correlated. Another source of disagreement is that within the gamma band, some reports claim that BOLD is best predicted by synchronous (narrowband) signals (Niessing et al., 2005), and others claim that BOLD is best predicted by asynchronous (broadband) neural signals [11]. Moreover, in some cases, it has been reported that no feature of the local field potential predicts the intrinsic optical imaging signal (closely related to BOLD) as accurately as multiunit spiking activity [24]. Consistent with this claim, a comparison of both motion and contrast response functions measured with single units and with BOLD suggested a tight coupling between BOLD and single unit responses [25][26][27]. To our knowledge, there is currently no single model linking the electrophysiological and BOLD signals that accounts for the wide range of empirical results.
The numerous studies correlating features of electrophysiological signals with BOLD provide constraints in interpreting the relationship between the two types of signals, yet the approach has not led to a general, computational solution. We argue that one reason that correlation studies have not led to computational solutions is that any particular feature of the field potential could be caused by many possible neuronal population responses. For example, a flat field potential (minimal signal) could arise because there is little activity in the local neuronal population, or it could arise from a pair of neuronal sub-populations responding vigorously but in counterphase, resulting in cancellation in the field potential. The same field potential in the two situations would be accompanied by different levels of metabolic demand and presumably different levels of BOLD signal. Similarly, any particular BOLD measurement could be due to many different patterns of neural activity. For example, stimulation of a neuronal population that inhibits local spiking can cause an elevation in the BOLD signal [28], as can stimulation of an excitatory population that increases the local spike rate [29]. In short, there can be no single transfer function that predicts the BOLD signal from the field potential, because the field potential does not cause the BOLD signal; rather, the neuronal activity gives rise to both the field potential and the BOLD signal.
We propose that many of the different claims pertaining to the relationship between BOLD amplitude and features of the field potential can be accounted for by a modeling framework in which BOLD and field potential measurements are predicted from simulated neuronal population activity, rather than by predicting the BOLD signal directly from the field potential (Fig 1a versus  1b). In this paper, we model fMRI and ECoG responses in two stages, one stage in which we simulate activity in a population of neurons, and a second stage in which we model the transformation from the population activity to the instrument measures. We then analyze the simulated instrument measures in precisely the same manner in which we analyze data obtained from fMRI and ECoG measurement of visual cortex. By design, the model employs a minimal set of principles governing how the instruments pool neural activity, rather than a biophysically detailed description of neural and hemodynamic events, enabling us to ask whether this minimal set of principles, together with simulated neuronal population activity, can simultaneously account for the empirical BOLD and ECoG measurements.

Results
We first present an analytic framework to capture basic principles of how the BOLD signal and the field potential pool neuronal signals across a population. Second, we apply this framework to predict ECoG and BOLD responses, and quantify how the simulated ECoG data are related to the simulated BOLD data. Third, we present results from ECoG and fMRI experiments in human V1, V2 and V3 and assess how well the predictions from the simulation can account for the relationship between ECoG and BOLD observed in the data.

Relationship between LFP and BOLD: analytic framework
The fMRI BOLD signal and the local field potential (LFP) measure neuronal population activity in a fundamentally different manner. The goal of this analytic framework is to capture these differences in simple mathematical expressions, and from these expressions derive the relationship between the two instrument measurements. We purposely omit a large number of biophysical details such as cell types, neuronal compartments, the dynamics of blood flow, and so forth, both for tractability and in order to emphasize the basic principles of how different measures integrate neuronal activity. In the sections that follow, we then show that, when coupled to simulated neural responses, the model can account for many important patterns observed in fMRI and ECoG data from human visual cortex.
For this analytic framework, we consider how a population of n neurons responds to a stimulus or task during a brief epoch (time 0 to T), assumed to be on the order of a second. Each neuron will produce a time varying dendritic current, denoted as Ii(t) for the ith neuron, resulting from the trans-membrane potential. We would like to know how these currents, Ii(t), relate to the fMRI BOLD signal and to the LFP signal measured by an ECoG electrode.
We assume that the LFP arises primarily from dendritic membrane currents [2]. We ignore output spikes (although spikes can influence the LFP [30], it is generally thought that the influence is smaller than synaptic and dendritic currents [2]), and including spikes would not change the logic of our arguments). For the ith neuron, the contribution to the LFP is then ! × ! . The constant ! depends on the distance and orientation of the neuron with respect to the electrode, as well as the electrode's impedance. For simplicity, we assume that each neuron is equidistant from the electrode and has the same orientation, like pyramidal neurons perpendicular to the cortical surface, and therefore its contribution to the electrode measurement is scaled by the same constant, . Because currents add, the LFP time series will sum the contribution from each neuron, Field potential recordings are usefully summarized as the power (or band-limited power) in the time series [31]. Here we summarize the LFP response within a short time window as the power in the signal summed over the time window T: (Equation 2) power of sum Importantly, Equation 2 is a linear/ nonlinear (L/N) computation, since the LFP power is computed by first summing the signals (L), and then computing the power (N).
The BOLD signal pools neural activity in a fundamentally different manner because it depends on metabolic demand [e.g., for reviews, see 1,32]. The metabolic demand of each neuron will increase if the cell depolarizes (excitation) or hyperpolarizes (inhibition) [28]. Hence the metabolic demand of a neuron is a nonlinear function of its membrane potential: either a positive or negative change in voltage relative to resting potential causes a current, thereby resulting in a positive metabolic demand. There are many possible nonlinear functions one could assume to summarize the metabolic demand from the dendritic time series, such as the rectified signal (absolute value) or the power (squared signal). For tractability, we assume the metabolic demand of the ith neuron is proportional to the power in the time varying trans-membrane current, integrated over time: , with ! a scaling constant for the ith neuron. (Similar results were obtained if used the absolute value rather than the power). For the entire population of n neurons, we then assume the BOLD signal will sum the metabolic demand of each neuron. For simplicity we use the same scaling constant for each neuron: sum of power Importantly, Equation 3 is a nonlinear / linear (N/L) computation, since the power is computed first (N) and then the signals are summed (L), opposite to the order of operations for the LFP in Equation 3 (Fig 1) (Personal communication from David J Heeger). In other words, we approximate the BOLD signal as the sum of the power, and LFP as the power of the sum, of the separate neuronal time series. The difference in the order of operations can have a profound effect on the predicted signals, as in the simple example with 2 neurons depicted in Fig 1c and 1d. The BOLD signal pooled over the two neurons is the same whether the time series from the two neurons are in phase or out of phase, whereas the LFP power is large when the time series are in phase and small when they are out of phase.

Fig 1. Pooling with different orders of operations can have a large effect on measured brain signals. A)
The approach to directly correlate LFP and BOLD data. B) Current approach to infer the LFP and BOLD from activity in a neuronal population. C) In this illustration the membrane potential of two neurons (x1 and x2) has the shape of a sinusoid with noise, and the sinusoid is in phase between the two neurons. In the simulated electrode measurement the signals are summed and the power is calculated (POWER(SUM) = 2.00). In the simulated measurement of metabolic demand, the power of each of these neurons is first calculated, and then summed across the neurons (SUM(POWER) = 1.01). Here, the POWER(SUM) > SUM(POWER). D) In this illustration the membrane potential of two neurons (x1 and x2) is the same as in panel C except that the two time series are in counterphase. Here, opposite from C), POWER(SUM) < SUM(POWER).
These approximations allow us to make predictions about the relation between LFP and BOLD. By theorem, we know that the power of the sum of several time series is exactly equal to the sum of the power of each time series plus the sum of the cross-power between the different time series (Equation 4): Power of sum Sum of power Sum of cross-power

(Equation 4)
Applying this theorem to Equations 2 and 3 shows the relationship between our model of BOLD and LFP power: We can now see that the LFP power depends on two quantities, one of which is related to the BOLD signal, and one of which is unrelated to the BOLD signal (Equation 5). The first term summarizes the total level of neural activity (summed across neurons), and the second term summarizes the relationship between neural time series (the cross-power, similar to covariance). If and when the second term tends to be large compared to the first, then the LFP power will not be closely related to the BOLD signal.
One cannot deduce from first principles whether the first term in Equation 4 (summed power) or the second term (summed cross-power) will dominate. However, he number of elements contributing to the two terms is quite different: For n neurons, the first term grows with n (the power in each neuron's time series), whereas the second term grows with n 2 (all the pairwise crosspowers). Hence if there is any appreciable covariance, then the LFP power will be dominated by the second term, and the correlation with BOLD will be weak (except in cases where the cross-power and power are highly correlated).

Relationship between LFP and BOLD: simulated experiments
We next used the framework above to compute the BOLD and LFP signals as a function of simulated neuronal population activity.

Simulations of BOLD and LFP responses from neuronal population activity
Cortical neurons receive a large number of inputs from diverse cell types. We simulated three inputs to each neuron (Fig 2a). Input 1 approximated Poisson-like spike arrivals (C 1 , 'Broadband'). Input 2 was a high frequency oscillation, peaked between 40 and 60 Hz, coordinated between neurons (C 2 , 'Gamma'). Input 3 was a low frequency inhibitory signal, 8-12 Hz (C 3 , 'Alpha'). For each simulated neuron i, the total input C on each trial is the sum of these three signals (Fig 2b): We then passed the summed input in each neuron through a leaky integrator to produce the timevarying dendritic current for that neuron (Ii, Fig 2c): The membrane time constant, reflects the time dependence of the trans-membrane current [33,34]. In total, we modeled a population of 200 neurons, each of which produced a one-second time series on each trial. To match the design of our ECoG experiments, a simulated experiment consisted of 240 1-s trials (30 repeats of 8 conditions). Below we explain (1) how the time series was generated for each of the three types of inputs, and (2) how the input levels varied for the 8 conditions. The three inputs were summed in each neuron to produce the total input to the neuron ( ! ). C) The total input was passed through a leaky integrator to produce the dendritic dipole current ( ! ).
Broadband input (C 1 ). Input ! was Gaussian white noise with a small positive bias. The Gaussian white noise approximates Poisson distributed spike arrivals, each of which produces a small positive or negative conductance change, corresponding to excitatory or inhibitory post-synaptic potentials. The small positive bias reflects the assumption of more excitatory than inhibitory synaptic currents, causing a net depolarization. Gaussian white noise was used rather than Poisson distributed synaptic inputs for computational efficiency, but the pattern of results is similar for Poisson or Gaussian distributions. With ! , we simulated high versus low rates of spike arrivals across simulation trials by using large versus small standard deviations of the Gaussian distributions. For all simulations, the coherence was zero between the 200 neurons in the population, meaning the time series in each neuron was independent except for chance correlations. This white noise input, after passing through leaky integration, results in an output whose power spectral density declines with temporal frequency. When the level of this input increases (higher standard deviation of the noise distribution), the result is a broadband elevation in power [34] (Fig 3, left panel), an effect that can be observed in the local field potential [35] as well as intracellular membrane potentials of single neurons in awake macaque [36].
Gamma input ( ! ). Input ! consisted of band pass noise , with fixed amplitude on all trials, and with coherence across neurons that varied between trials. This input approximates the signals giving rise to narrowband gamma oscillations. Across different conditions, we varied the coherence between neurons of ! rather than the amplitude for individual neurons, which was fixed. The motivation for this comes from empirical observations that large gamma oscillations in the LFP tend to reflect increased coherence between neurons [37,38]. This is opposite to the broadband input ( ! ), for which we varied the amplitude in individual neurons across trials, rather than the synchrony between neurons. Narrowband gamma oscillations with a peak between 30 and 80 Hz can be observed in the local field potential [39,40], as well as in the membrane potential of individual pyramidal neurons [41]. When we increase the coherence of ! in our simulations, the result is an increase in the amplitude of the LFP in the gamma band (Fig 3, middle panel), much like narrowband gamma signals observed in microelectrode recordings [42] and human ECoG [43].
Alpha input ( ! ). The alpha input consisted of inhibitory oscillations at approximately 10 Hz, with fixed coherence between neurons, and varying amplitude across trials. The oscillations were inhibitory, i.e. the mean was below the resting potential (compare ! versus ! and ! in Fig 2a).
Because ! was inhibitory, it results in hyperpolarization (or less depolarization), opposite the effect of ! , which resulted in depolarization. This input approximates the signals giving rise to alpha oscillations (Fig 3, right panel). Pyramidal neurons in visual cortex have been hypothesized to receive periodic inhibition, with pulses arriving at approximately 10 Hz [44,45]. Individual neurons in visual cortex can indeed show subthreshold membrane oscillations at frequencies around 10 Hz [46]. From the neuronal population simulations, we computed the LFP and BOLD signals according to the equations above (section "Relationship between LFP and BOLD: analytic framework"). In brief, the LFP was computed by summing the trans-membrane current across neurons (Equation 1), and the BOLD signal was computed by summing the power across neurons (Equation 3). The BOLD signal was averaged across the 30 trials with the same condition. The LFP time series was transformed to a power spectrum, which was averaged across the 30 repeated trials of the same condition. The average power spectrum was summarized by three components: elevation in broadband power, the power of narrowband gamma oscillations , and the power of alpha oscillations (8)(9)(10)(11)(12).
The simulations were structured to approximate the experimental design and the results of our ECoG experiments (Section 2.3). In brief, in the ECoG experiments, there were four grating stimuli of different spatial frequencies, three noise patterns with different power spectra, and one blank stimulus (mean luminance). For each of the 8 stimuli and each of 22 electrodes in V1-V3, we decomposed the measured ECoG responses into the same three components computed from the simulations: broadband, narrowband gamma, and alpha. . The bar plots show the amplitude per neuron and the coherence across neurons for each of the 8 simulated stimuli. The output spectra were computed by summing the time series across the 200 neurons, computing the power spectra on each trial, and averaging the power spectra across 30 repeated trials of the same stimulus. The colors of the traces match the colors of the 8 stimuli in the bar plots. The spectra show broadband changes across all frequencies arising from C 1 , a peak around 50 Hz from C 2 , and a peak around 10 Hz from C 3 . The right panel shows the simulated and measured summary metrics for each of the 8 stimuli: alpha power (green), gamma power (magenta), and broadband (black). The simulated measures are derived from the spectra in the 'Output spectra' panel. The measured metrics are derived from a V1 electrode. The simulated and measured values are close by construction, since the inputs to the simulation were fitted to match the output spectral measurements from this electrode, see a plot of all electrodes in Supplemental Figure S9. B) Same as A, except that the simulation was matched to a V2 electrode rather than a V1 electrode.
The level of the three inputs, C 1 , C 2 , and C 3 were determined so that that the simulated summary metrics (broadband, gamma and alpha) matched the measured ECoG summary metrics for each of the 8 conditions for a particular electrode (Fig 4). Importantly, the determination of the input levels did not take into account the measured BOLD responses. Hence the simulations provided a test: if the inputs were chosen to produce outputs that match the measured ECoG responses, does the simulated BOLD signal accurately predict the measured BOLD signal?
We first review two simulations, one matched to ECoG measurements from a V1 electrode, and one matched to ECoG measurements from a V2 electrode. We then summarize the results of 22 simulations corresponding to the 22 electrodes in V1-V3. As with the ECoG experiments, in the simulation each condition was repeated 30 times for 240 trials.

Broadband and alpha LFP, but not gamma LFP, correlate with BOLD
For each of the 8 stimuli in one simulation, we obtained 4 summary responses, a single BOLD amplitude and three LFP measures (broadband, gamma and alpha, Fig 5A). By design, the values of the LFP components (but not BOLD) were close to the values measured in one of the electrodes.

Fig 5. Simulated correlation between BOLD and LFP. A)
Schematic to show LFP summary metrics derived from spectra: broadband power elevation, narrowband gamma and alpha. Broadband was calculated by the increase in a 1/f n signal, gamma was calculated by fitting a Gaussian on top of the 1/f n line, and alpha was calculated as the difference from baseline in the alpha-frequency-range. B) The levels of C 1 , C 2 , C 3 , and were chosen for each of 8 conditions to approximate the measurements of broadband, gamma, and alpha for 8 different stimuli from a V1 electrode. Red dots indicate that inputs were matched to responses to a grating pattern, blue dots are matched to noise patterns, and black is matched to a blank stimulus. The correlation between broadband and simulated BOLD is high (left), whereas the correlations between gamma and BOLD (middle) and alpha and BOLD (right) are low. Error bars show the 68% confidence interval across bootstraps. C) Same as B except that inputs were chosen so that LFP measures matched ECoG data from a V2 electrode. Here, the correlation is high and positive between broadband and BOLD (left), low between gamma and BOLD (middle), and highly negative between alpha and BOLD (right).
For Simulation 1, we observed three patterns ( Fig 5B). First, the LFP broadband response correlated highly with BOLD (R 2 = 0.88). Second, the power of narrowband gamma oscillations was weakly correlated with BOLD (R 2 = 0.01), and third alpha oscillations were weakly correlated with BOLD (R 2 = 0.06). In this simulation, the LFP broadband response was the only good correlate of the BOLD signal.
In Simulation 2, the broadband LFP was again a good predictor of BOLD, and gamma power was again a poor predictor (Fig 5C). Unlike Simulation 1, the power of alpha oscillations was strongly negatively correlated with BOLD (R 2 = 0.74). The stronger correlation between alpha oscillations and BOLD was not due to a difference in the way the alpha oscillations were simulated, nor to the way BOLD was computed -the simulation algorithm and the analysis code extracting the BOLD and LFP measures were identical. The observations from the two sets of simulations emphasize the fact that the identical mechanism that converts neural activity to BOLD ("neurovascular coupling"), modeled here as power in the time series summed across neurons (Equation 3), can produce very different correlations between BOLD and features of the LFP. The two sets of simulations differed only in the neural population activity, not in the transfer functions relating neural responses to instrument measurements. The higher correlation between alpha power and BOLD in Simulation 2 versus 1 is due to the fact that the alpha power varied more, and the broadband power varied less, across the 8 stimuli in the simulation.

Relationship between ECoG and fMRI: Data
The ECoG responses, like the simulated LFP responses, were separated into three components for each stimulus condition in each electrode: a broadband component, a narrowband gamma component, and an alpha component. The separation into the three components was accomplished using the identical analysis code as in the simulated data. An important feature of this data set is that the 3 components of the ECoG response showed different patterns across stimuli [43]: stimuli comprised of noise patterns caused large broadband increases but little to no measureable narrowband gamma response, whereas grating stimuli elicited both broadband increases and narrowband gamma increases. Gratings and noise stimuli resulted in decreases in alpha power compared to baseline (Supplementary Fig S1). Had all three responses been correlated with each other, it would not be possible to infer how each relates separately to the BOLD signal.
We measured BOLD responses in 4 healthy subjects to the same visual stimuli as used in ECoG (subjects are different from the ECoG subjects), and extracted the signal from regions of interest in visual cortex matched to the previously recorded ECoG electrode locations ( Supplementary Fig S2  and S3). Combining the ECoG and BOLD data sets enabled us to summarize the data with the same four measures for each stimulus that we used to summarize the simulation: a single BOLD response and three ECoG component responses.

Fig 6. Measured correlation between ECoG and BOLD in V1 and V2/V3. A)
The correlation between ECoG and BOLD was calculated for electrodes in V1 and V2/V3. The locations of one sample electrode in V1 and one in V2 are indicated by the enlarged white discs on the cortical surface for subject 1. B) In a foveal V1 site, the broadband ECoG amplitude predicted the BOLD signal (left). Error bars show 68% confidence intervals across bootstraps. Narrowband gamma power (center) and alpha power (right) did not predict BOLD. C) In a V2 site, the broadband ECoG was weakly correlated with BOLD (left). Narrowband gamma did not predict BOLD (middle). Alpha was negatively correlated with BOLD (right). The schematic spectral plots beneath the scatter plots illustrate the three ECoG components. Scatter plots for all other V1 and V2/V3 sites are shown in Supplementary Fig S5. We first consider the V1 site that was used to guide the input levels for Simulation 1 (Fig 6B). We emphasize that the BOLD data for this site was not used to guide the simulation (only the ECoG data were used). Three patterns relating the BOLD signal to the ECoG data for this site match Simulation 1. First, the BOLD amplitude for the different stimuli was highly correlated with the broadband response (R 2 = 0.79). Second, the narrowband gamma response was unrelated to BOLD (R 2 = 0.04). Third, there was little correlation between BOLD and alpha power (R 2 = 0.02). Hence, as in simulation 1, only broadband power was a good predictor of BOLD amplitude.
A different pattern was found in the V2/V3 data, and this pattern was similar to that observed in Simulation 2. For the V2 site guiding Simulation 2, the correlation between BOLD and broadband was smaller than that found for the V1 channel (R 2 = 0.33), and the correlation between alpha power and BOLD was stronger than for V1, and was negative (R 2 = 0.60). Gamma power did not correlate with BOLD, similar to V1 data and to Simulations 1 and 2.
The differences in the BOLD/ECoG responses in V1 compared to V2/V3 could, in principle, arise either because the neurovascular coupling mechanism differs between the two cortical areas, or because the neural responses for these particular experimental conditions differ in the two areas. While we do not have access to the neurovascular coupling mechanisms nor to the neural response at the individual neuron level, we can compare the measured ECoG responses in the two brain areas. We observe one large difference between V1 and V2/V3, which is that the range of broadband responses elicited by the different stimuli was larger for V1 than for V2/V3: the average broadband increase above baseline was 0.43 log units for V1 and 0.26/0.28 for V2/V3 (Supplementary Fig S4). As shown in Simulations 1 and 2, a difference in the neuronal population responses can result in a different pattern of response between features of the ECoG signal and the BOLD signal.

Data model comparison
In the example V1 electrode and in Simulation 1, the BOLD signal was well explained by broadband increases (Fig 5B; Fig 6B, left). In the example V2 electrode and in Simulation 2, the BOLD signal was explained by both broadband increases and alpha decreases (Fig 5C, Fig 6C). By explicitly modeling both the population response and the population-to-instrument transformations, we see that a difference in the relation between instrument measures (BOLD and ECoG) can arise from a difference in the population response, without a difference in neurovascular coupling. We now ask (1) whether these effects are consistent across the measured V1 and V2/V3 electrodes, and (2) how a multiple regression model using broadband, gamma and alpha as predictors fits the BOLD response for both data and simulation.

Predicting the BOLD amplitude from features of the LFP
As we argued in the Introduction, we believe there is no single, general transfer function that can predict the BOLD signal from the LFP. Yet a regression model linking the two measures can be a useful way to summarize the results of a particular experiment or simulation, and to compare results between different experiments or simulations. Here, we fit several regression models to each simulation and to the data. The regression models predicted the simulated or measured BOLD response from either a single LFP component (broadband power, gamma power, or alpha power), or from combinations of LFP components (each of the pairwise combinations, and the 3 components together). These models were fit separately for each of the 9 electrodes in V1 and each of the 13 electrodes in V2/V3 (19 electrodes from subject 1; 3 from subject 2), and for the 22 corresponding simulations. Accuracy of each model was assessed by cross-validation. With a crossvalidation procedure, there is no advantage in accuracy for models with more free parameters.
Simulations. For the simulations, we expect broadband power to positively predict the BOLD signal and alpha power to negatively predict the BOLD signal, because of the construction of the simulations: broadband and alpha power elevations were achieved by increasing the level of inputs per neuron, rather than coherence; the converse was true for gamma. Nonetheless, solving the regression models can be informative because, as seen in Fig 5, simulations with the identical input types and the identical analysis can lead to different patterns, depending on the mixtures of values simulated. Moreover, the regression analyses of the simulated data provide a baseline against which we can compare the regression analyses of the observed data.
The single parameter models across electrodes in V1 and V2/V3 show the same pattern as the single electrode examples in Figs 5 and 6. For simulations fit to V1 ECoG data, broadband alone was a good predictor of BOLD (average cross validated R 2 =0.87 across 9 simulations) while gamma and alpha alone were not (gamma R 2 = 0.17, alpha R 2 = 0.06). For simulations fit to V2/V3 ECoG data, broadband and alpha alone each predicted the BOLD signal with moderate accuracy (broadband R 2 = 0.68, alpha R 2 = 0.37).
For simulations fit to V1 and to V2/V3, the BOLD response was best explained by a regression model combining broadband and alpha (R 2 = 0.95, R 2 = 0.97, Fig 7A and B). The regression coefficients for this model were positive for broadband and negative for alpha. A model that incorporated all three LFP measures -broadband, alpha and gamma -did not explain additional variance for either simulation, confirming our earlier observation that narrowband gamma power was not correlated with BOLD amplitude in simulated data. for each of the 7 model types averaged across the 9 simulations fitted to V1 data. Error bars are standard error of the mean across the 9 simulations. The R 2 was cross-validated (split between even and odd stimulus repetitions), to ensure that the R 2 can be compared between models with different numbers of explanatory variables. Bottom: The regression coefficients show whether the broadband, gamma, and alpha signals were positive or negative predictors of the BOLD signal. B) Same as A, but for the 13 simulations fitted to V2/V3 ECoG data. C) Same as A, but for measured BOLD and ECoG data averaged across the 9 V1 sites. The R 2 was cross-validated (split between subjects for BOLD and stimulus repetitions for ECoG). A red * in the lower plot indicates whether regression coefficients differed significantly from zero by a t-test (p<0.05). D) Same as C, except for the 13 electrodes in V2/V3. (E) Simulated BOLD (x-axis) versus measured BOLD (y-axis) for the 9 V1 sites. Each color corresponds to one site. The R 2 was computed separately for each of the 9 sites, and then averaged. (F) Same as E, but for V2/V3 sites.

Data.
The regression analysis between the measured ECoG and BOLD signal in V1 (Fig 7C) was very similar to the regression analysis in the simulation (Fig 7A). The observed BOLD signal, like the simulated BOLD signal, was best predicted by a combination of broadband and alpha changes (R 2 = 0.73). While narrowband gamma changes were modestly predictive of the BOLD signal across V1 electrodes in a single parameter model (R 2 =0.11), adding gamma as a predictor to either the broadband only model (model 2 versus model 1) or to the combined broadband and alpha model (model 7 versus model 5) led to a slight decrease in accuracy, suggesting that any predictive power of gamma alone was probably due to correlation with the other ECoG predictors, rather than a direct relationship to the BOLD signal.
The model fits of the measured BOLD signal in V2/V3 ( Fig 7D) were qualitatively similar to the model fits for the V2/V3 data (Fig 7B), though with some quantitative differences. The observed BOLD signal, like the simulated BOLD signal, was best predicted by a combination of broadband and alpha changes (R 2 = 0.50) (see also Supplementary Fig S6). As with V1, narrowband gamma in V2/V3 did not explain additional variance in the BOLD signal. Overall, compared to V1, the BOLD signal in V2/V3 was less accurately predicted by the regression models based on the electrophysiological measurements (Fig 7D), and was also less well predicted by the simulated BOLD responses (average R 2 of 35%) ( Fig 7F).
Finally, we directly compared the BOLD responses predicted in the simulations with the measured BOLD data. These were in good agreement, with the BOLD data in the simulations accounting for an average of 80% of the variance in the measured data in V1 and 44% in V2 (R 2 averaged across sites) (Fig 7E-F).
Across simulations and data sets, a general pattern emerges. The broadband signal was significantly and positively predictive of BOLD, and alpha power was significantly and negatively predictive of BOLD. Narrowband gamma had no consistent relation with BOLD. While the relationships between broadband and BOLD and alpha power and BOLD were consistent in terms of sign (the former positive, the latter negative), the level was not always the same. As we noted in the example channels shown in Fig 7, the broadband power was more strongly predictive of BOLD in V1, and alpha power was more strongly predictive of BOLD in V2/V3. An examination of responses to the different stimulus types clarifies the difference between V1 and V2/V3 in these data. Specifically, the BOLD response in V2/V3 to noise patterns was under-predicted by the broadband response alone (Supplementary Fig S5). In V2/V3 alpha decreased more for the noise patterns and this alpha decrease accounted for the BOLD change unexplained by broadband (Supplementary Fig S1). This helps to explain why a model that includes broadband and alpha is much more accurate for V2/V3 than a model that includes only broadband. In contrast, for V1 the BOLD response was well predicted by broadband power in most sites, with little systematic prediction failures, and little room for increased model accuracy when adding predictors such as alpha power.
For each of the 22 simulations, the three inputs C 1 , C 2 , and C 3 defining each of the 8 stimulus conditions were fit to produce the LFP data from the corresponding ECoG electrode. By design, the C 1 (broadband) and C 3 (alpha) inputs were fit to data by varying the level per neuron, whereas C 2 was fit to data by varying the coherence across neurons. In principle, for any of the three inputs, the ECoG data could have been fit by varying either the level per neuron or coherence across neurons. For completeness, we tested the 7 alternative models (Supplementary Fig S7). The most accurate model, quantified as the R 2 between the measured BOLD and the simulated BOLD, averaged across the 22 sites, was the simulation type used in the main text, in which C 1 and C 3 varied in the level per neuron and C 2 varied in the coherence across neurons. The simulation in which broadband coherence rather than level modulated broadband power caused a large drop in R 2 , The models in which the gamma LFP power was modulated by the coherence rather than the level only caused a small drop in R 2 .

Correlation between BOLD and LFP across all frequencies
Simulation. In the previous section, we modeled the BOLD responses as a linear function of three components derived from the LFP. These features -broadband power, narrowband gamma power, and narrowband alpha power -are summary metrics of the power spectrum. We also tested how the power at each frequency in the simulated LFP and in the ECoG data correlated with the BOLD response (simulated and observed). We calculated LFP power for each frequency from 1 to 200 Hz and correlated the power changes from baseline with BOLD changes from baseline (Fig 8). In simulations fit to V1 ECoG data, the LFP correlated well with the BOLD signal across all frequencies except those in the alpha band (8)(9)(10)(11)(12)(13)(14)(15) and below, and those in the gamma band (40-60Hz).
In simulations fit to V2/V3 ECoG data, the pattern was similar, except that the correlation was negative in the alpha band rather than 0, and weaker but still positive in the rest of the spectrum. These patterns match the summary metrics of alpha, gamma, and broadband shown in Fig 7.   Fig 8. The correlation between BOLD and LFP as a function of frequency. A) The correlation between LFP and BOLD for simulations fit to V1 shows that there is a positive correlation across most frequencies, except those including the alpha and gamma changes. B) The correlation between LFP and BOLD for the simulations fit to V2/V3 shows that there is a strong negative correlation around 10 Hz, and a positive correlation across a broad range of frequencies. C) In the V1 data there was a positive correlation between ECoG and BOLD for a broad range of frequencies, except those including the alpha changes. Gray lines represent the 9 individual V1 electrodes, the black line is the average, the red line corresponds to the example electrode shown also in Fig 6B, and to the simulation in panel A. D) In the V2/V3 data there was a strong negative correlation between ECoG and BOLD in the alpha range around 10 Hz and a positive correlation between ECoG and BOLD for a broad range of frequencies. Gray lines represent the 13 individual V2/V3 electrodes, the black line is the average, the red line corresponds to the example electrode shown also in Fig 6C, and to the simulation plotted in panel B. Note that neither the V1 electrodes nor the V2/V3 electrodes show a peak at the gamma frequency.
Data. The correlation between measured BOLD and ECoG power in V1 was qualitatively similar to that found in simulations. In V1, ECoG responses across all frequencies except the alpha band were positively correlated with the BOLD response (Fig 8C), consistent with the regression analyses of the summary metrics, showing that broadband ECoG power was the best predictor of the BOLD signal.
The pattern of correlation between ECoG power and BOLD in V2/V3 was similar to that found in V1, although the overall level of correlation was lower. There were positive correlations between ECoG and BOLD extending across most frequencies that were weaker than in V1, and there was a negative correlation for most sites in the alpha band.
There were some differences between V1 sites. For example, in two sites, the correlation across frequencies dipped in the gamma band (30-80Hz), similar to simulated data. These are also the two sites that showed the largest amplitude gamma responses (sites 8 and 9 in Supplementary Fig S5).
In other words, when cortical sites showed large gamma signals, these signals were uncorrelated with BOLD. The fact that in 7 V1 sites there was a positive correlation between BOLD and LFP power spanning 30-80 Hz might seem inconsistent with our earlier observation that narrowband gamma power was not predictive of the BOLD signal in V1 sites (Fig 6 and 7). However, in this analysis the narrowband and broadband power are not modeled separately and the positive correlation between power at 30-80 Hz in Fig 8 thus likely suggests that broadband power extends into this band, since broadband changes can extend across all frequencies [11,34]. Therefore, if there is little to no narrowband response, we would expect a positive correlation between BOLD and ECoG throughout all frequencies.
There were some site-to-site differences in the correlation between alpha and BOLD. For example, some sites showed a positive correlation with BOLD in the alpha range, and others showed a negative correlation (Fig 8C). These site to site differences depend on the range of responses evoked by stimuli. For example, for electrodes in which stimuli evoked a large of power changes in the alpha band, alpha was more strongly correlated with BOLD. Similarly, for electrodes in which stimuli evoked a large range of broadband responses, broadband was more highly correlated with BOLD ( Supplementary Fig S8). This pattern did not hold for narrowband gamma power changes.

Discussion
This study investigated the relationship between electrophysiological and BOLD measurements in human visual cortex. Our modeling framework decomposed the signals into two stages, a first stage in which we simulated the neuronal population responses (dendritic time series), and a second stage in which we modeled the transfer of the neuronal time series to the BOLD signal and the field potential. This approach differs from the direct comparison of electrophysiological signals and BOLD. The explicit separation into stages clarified both a similarity and a difference between the BOLD amplitude and the field potential power: the two can be approximated as the same operations on the neuronal population activity, but applied in a different order. Specifically, within a brief window, we modeled the BOLD amplitude as the sum of the power in the neuronal time series, and the field potential as the power of the sum of the neuronal time series. Because the order of operations differs, the two signals differ, and each is blind to particular distinctions in the neuronal activity. For example, the BOLD signal (according to our model) does not distinguish between synchronous and asynchronous neural signals with the same total level of activity. In contrast, the field potential does not distinguish counterphase responses from no response. Even if one knew the exact mechanism of neurovascular coupling and the precise antenna function of an electrode, one still could not predict the relationship between the BOLD signal and the field potential without specifying the neuronal population activity that caused both. Hence the relationship between the two types of signals is not fixed, but rather depends on the structure of the underlying responses of the neuronal population.
Although we do not have access to the complete set of individual neuronal responses in any of our experiments in visual cortex, we approximated the responses by specifying the type of signals common to visual cortex. We therefore limited the space of neuronal population responses by modeling the activity as arising from three types of signals, enabling us to compute the complete set of field potentials and BOLD responses to a variety of conditions. Finally, we compared the simulated patterns of BOLD and field potential responses to the actual responses we observed in data from human subjects. These patterns are discussed and interpreted below.

Changes in broadband power predict BOLD
Many studies have reported correlations between BOLD and power in the gamma band LFP (30-130 Hz) (review for human studies: [47]). Yet changes in gamma band power do not reflect a single biological mechanism. For example, several recent studies have emphasized that LFP power changes in the gamma band reflect multiple distinct neural sources, including narrowband oscillations and broadband power shifts, with very different stimulus selectivity and biological origins [43,48,49].
Broadband changes have been proposed to reflect, approximately, the total level of Poisson distributed spiking (or spike arrivals) in a local patch of cortex [34]. In contrast, the narrowband gamma response is caused by neural activity with a high level of cell-to-cell synchrony [50] and likely depends on specialized circuitry [51]. While the two responses are sometimes distinguished as 'high gamma' (referring to broadband signals), and 'low gamma' (referring to oscillatory signals), this distinction is not general. Broadband signals can extend into low frequencies [11,52] so that the two signals can overlap in frequency bands. Hence separating gamma band field potentials into an oscillatory component and a broadband (non-oscillatory) component is not reliably accomplished by binning the signals into two temporal frequency bands, one low and one high, but rather requires a model-based analysis, such as fitting the spectrum as the sum of a baseline power law (to capture the broadband component) and a log-Gaussian bump (to capture the oscillatory component) [43].
There is strong experimental support for the idea that increases in broadband LFP power primarily reflect increases in asynchronous neural activity rather than increases in coherence. First, experiments have shown that broadband power is correlated with multiunit spiking activity [49,53]. Second, unlike the case of narrowband gamma LFP, changes in broadband LFP are not accompanied by changes in broadband spike-field coupling ( [37], their Fig 1A-B). The possibility that neuronal synchrony sometimes affects broadband signals cannot be ruled out, for example as shown in cases of pharmacological manipulations in nonhuman primate [54]. In such cases, there would not be a simple relationship between broadband power and BOLD.
The prior literature has not shown definitively whether broadband LFP, narrowband gamma, or both predict the BOLD signal. The first study that directly compared simultaneously recorded BOLD and electrophysiology showed that both LFP power in the gamma frequency range (40-130 Hz) and multi unit spiking activity (MUA) predicted the BOLD signal [16], and further, that when the LFP power diverged from MUA, the gamma band LFP predicted the BOLD signal more accurately than did spiking. This study however did not separately test whether a narrowband (oscillatory) or a broadband (non-oscillatory) component of the LFP better predicted the BOLD response.
Other studies, too, have shown a variety of patterns when correlating LFP power changes in the gamma band with BOLD. Some reported that BOLD amplitude correlates with narrowband gamma activity [13], while others showed that BOLD correlates with broadband changes [11], and many did not distinguish narrowband from broadband power in the gamma band [55]. Simultaneous recordings of hemodynamic and neuronal activity in macaque V1 showed that BOLD signals from intrinsic optical images can occur in the absence of gamma band LFP changes [56], and that in some circumstances, multiunit activity predicts the BOLD signal more accurately than gamma band LFP [24,57].
Here we separately quantified the broadband power (spanning at least 50-150 Hz) and narrowband gamma power. We found that the amplitude of broadband changes accurately predicted the BOLD signal in V1. The empirical results and the models help resolve the question of why 'high gamma' has been shown to correlate with BOLD, and 'low gamma' sometimes does not [24]. The likely reason is unrelated to the difference in frequency range, nor to the size of the spectral perturbation in the local field potential. In fact, the elevation broadband power is relatively small (2 or 3 fold) compared to the elevation in power often observed in narrowband gamma oscillations (10 x or more) [43]. Instead "High gamma" is predictive of the BOLD signal in many cases not because of the specific frequency range, but because this signal captures the level of asynchronous neuronal response; this signal happens to be most clearly visible in the high frequency range (>100 Hz) where it is not masked by rhythmic lower frequency responses. Hence the distinction in predicting the BOLD response is not about "high" versus "low" gamma, but rather synchronous versus asynchronous responses, and the broadband signal, sometimes labeled high gamma, maps onto the first term on the right hand side of Equation 4, the portion of the field potential measurement which sums the energy demand of each neuron.
Our model fits and data support this view. When we captured the stimulus-related broadband response by simulating a change in broadband coherence across neurons rather than a change in the level of response in each neuron, our predicted BOLD response was highly inaccurate (Supplementary Fig S7).

Changes in narrowband gamma power do not predict BOLD
In contrast, we propose that 'low gamma' often does not predict the BOLD response because 'low gamma' reflects narrowband oscillatory processes, which largely arise from a change in neuronal coherence across the population rather than a change in the level of neuronal population activity. This corresponds to the second term in the right hand side of Equation 4, the portion of the field potential measurement which reflects the cross-power arising from currents in different neurons, and which in our model is independent of the signals giving rise to the BOLD signal.
Our results and model do not argue that narrowband gamma oscillations will never be predictive of the BOLD signal. If in a particular experiment narrowband gamma oscillations were to co-vary with broadband increases, we would expect both signals to correlate with BOLD. This might occur in an experiment with gratings of different contrast; with increasing contrast narrowband gamma responses, broadband responses, and BOLD responses all increase [21,58] and all 3 measures would correlate across stimuli. In such an experiment, if narrowband gamma oscillations had a higher signal to noise ratio than the broadband response, then the oscillatory signal would likely show a higher correlation with BOLD. In contrast, when the choice of stimulus or task can independently modulate broadband power and gamma oscillations so that the two LFP measures are not correlated, as in the experiments presented here and previously [43], then gamma oscillations will not strongly correlate with BOLD.
Our simulation and empirical results are consistent with studies which varied chromatic contrast and spatial frequency, while measuring MEG and BOLD. These studies found that BOLD and narrowband gamma activity did not match in stimulus specificity [18,19]. It is likely that these stimulus manipulations, like ours, independently modulated narrowband gamma power and broadband power, although the studies did not quantify broadband fields, which are more challenging to measure with MEG than with ECoG [59]. We speculate that broadband fields spanning the gamma range would have shown a higher correlation with BOLD.

Neuronal synchrony and the BOLD signal
In our model, the LFP measures are highly sensitive to neuronal synchrony, whereas BOLD is not. In our simulations, increases in neuronal synchrony drove narrowband gamma oscillations in the field potential. There are other cases of population activity with a high degree of neuronal synchrony. One example is the steady state evoked potential associated with a periodic stimulus [first reported by 60,reviewed by 61]. Previous studies have described discrepancies between evoked potentials and the BOLD signal, such as in the case of spatial summation [11], directional motion selectivity [7,8] and spatial attention [9,10]. Our modeling framework suggests that the neural sources generating the steady state potential (synchronous neural activity) are likely to be only weakly related to the BOLD signal (depending largely on asynchronous signals), as these sources will primarily affect the second term on the right hand side of Equation 4 (cross-power). This does not imply that the two measures are always or even usually discrepant; the BOLD signal and steady state potentials are likely to correlate any time that the steady state signals correlate with other measures of neural activity. When measures do dissociate, we do not conclude that one measure is more accurate; instead, the measures offer complementary views of the population activity, emphasizing the degree of synchrony or the average level of the response. An intriguing question is how each of the two signals contributes to perception and behavior.
Neural synchrony can also emerge without being time-locked to the stimulus, often called 'induced synchrony' or 'induced oscillations' [62]. In our simulation, we assumed that narrowband gamma LFP changes were induced by increases in synchrony between neurons, and not by changes in the level of gamma power within the individual neurons. In contrast, we assumed that broadband LFP increases were induced by increased broadband activity in individual neurons, and not by increased broadband coherence between neurons. (In Equation 4, a change in the left hand side, LFP power in the gamma band, can be produced by a change in either the first or second term on the right). This explains why, in our simulation, the broadband power was correlated with BOLD whereas the LFP gamma power was not, findings that were also confirmed by the data. Were our assumptions justified?
In principle, an increase in narrowband gamma power in the LFP could arise because the neurons synchronize in the gamma band, or because ongoing gamma oscillations within each neuron increase in amplitude, independent of coordination between neurons. There is strong experimental support for the former. Experiments which measure both intracellular membrane potential from single neurons and the extracellular LFP show that when there is an increase in narrowband LFP gamma power, the gamma power from individual neurons becomes more coherent with the LFP [41]. Moreover, the coherence between local spiking and the LFP also increases in the gamma band when LFP gamma power increases [37]. These results are consistent with our assumption that a significant part of the increase in gamma LFP power arises from a change in population coherence. To our knowledge, it is not certain whether there is also some increase in the level of gamma signals within individual neurons when the narrowband gamma band LFP power changes. However, since we can attribute a large part of the change in gamma LFP to a change in coherence, we infer that we can only attribute, at most, a small part of the change in gamma LFP to the level of gamma power within neurons.
In our simulation, we made two simple but extreme assumptions. First, we assumed that gamma oscillations occur with no change in the total level of neural activity, and hence no change in metabolic demand or BOLD. Second, we assumed that broadband responses occur with no change in neural synchrony. While these assumptions are likely incorrect at the limit, the simulations nonetheless captured the pattern of ECoG and fMRI results obtained in our datasets. Alternative models in which the broadband response was caused by a change in synchrony were much less accurate (Supplementary Fig S7). Models in which gamma responses were caused by a change in level were only slightly less accurate, and cannot be ruled out entirely (Supplementary Fig S7). However the regression models fit to our data (Fig 7) show that the power of narrowband gamma oscillations does not predict the BOLD response. Hence the most parsimonious explanation is that these responses in the LFP are caused in large part by changes in synchrony.

DC offsets and the BOLD signal
Both our measurements and our simulations showed that broadband electrophysiological responses were related to, but did not fully account for, the BOLD signal. This was especially evident in Simulation 2 and extrastriate data (V2/V3). In these cases, the amplitude of low frequency oscillations (8)(9)(10)(11)(12)(13)(14)(15) was negatively correlated with the BOLD signal, independent of broadband signals. Numerous previous studies have reported that low frequency oscillations are anti-correlated with BOLD, including measurements in motor, visual and language areas [20][21][22][63][64][65]. This result may appear to conflict with the prior discussion, since we argued that oscillations (to the degree that they reflect neuronal synchrony) should have little to no effect on metabolic demand or the BOLD signal. It is therefore important to ask why low frequency oscillations sometimes correlate with the BOLD signal, both in data and in simulation.
One explanation is that alpha oscillations, or a mechanism which generates the oscillations, affect the BOLD signal indirectly, by inhibiting cortical activity. According to this explanation, an increase in alpha power results in a decrease in local spiking activity in turn reducing metabolic demand and the BOLD signal [66]. Alpha oscillations may indeed co-occur with reduced cortical excitation [67]. However, if this coupling between alpha power and spiking were the only explanation for the relationship between alpha power and BOLD, then a more direct measure of neuronal excitation, such as broadband or multiunit activity, would adequately predict the BOLD signal; alpha power would negatively correlate with the BOLD signal, but would provide no additional predictive power. Our data and model do not support this explanation, as we find that for most cortical sites, the most accurate predictor of the BOLD signal is a combined model including both the amplitude of alpha oscillations and broadband power.
We therefore propose that in addition to the indirect effect of modulating cortical excitability, alpha oscillations are also accompanied by a DC shift in membrane potential, making it less depolarized (i.e., closer to the equilibrium potential), and this shift reduces metabolic demand. Indirect evidence for a DC shift comes from MEG and ECoG studies [44,45,68] and can be explained by a simple process: if alpha oscillations reflect periodic inhibitory pulses, then on average they will cause a depolarization. If the neuron was slightly depolarized before the inhibitory alpha pulses, then the pulses would push the neuron toward equilibrium, and hence a lower energy state. In this view, large alpha oscillations reflect larger inhibitory pulses, and hence reducing depolarization. We suggest that this reduced depolarization affects metabolic demand in two ways: by reducing spiking (as discussed above), and by maintaining a less depolarized state, reducing metabolic demand. In our model, the contribution to the BOLD signal from each neuron is the power in the time series (Equation 3), and the mean contributes to power. The idea that a DC shift in the membrane potential affects metabolic demand (in addition to altering excitability) is consistent with the observation that slowly changing currents (<0.5 Hz) correlate with BOLD fluctuations [12,69].

A single modeling framework accounts for patterns of LFP/BOLD correlations across sites
We found that the relationship between the BOLD signal and features of the ECoG data differed across cortical areas. For example, broadband changes in ECoG responses explained more variance in the BOLD data in V1 than in V2/V3. Conversely, low frequency power decreases (alpha, 8-13Hz) explained more variance in the BOLD signal in V2/V3 than in V1. In the absence of a model, we might have interpreted the results as evidence that neurovascular coupling differs across sites. Many previous studies have reported differences in the relation between LFP and BOLD as a function of site or condition, for example between cortical and subcortical locations [70], across cortical regions [71,72], between cortical layers [73], and as a function of medication [74]. Here, we showed that a difference in the relationship between LFP and BOLD need not arise because of a difference in neurovascular coupling. In our results, Simulations 1 and 2, like V1 compared to extrastriate areas, showed differences in the relationship between LFP and BOLD, yet we used the identical model of neurovascular coupling in all simulations. The systematic differences in the two simulations arose because of differences in the neuronal population activity, not because of differences in neurovascular coupling. While our results do not exclude the possibility of differences in neurovascular coupling across locations or states they do caution against interpreting differences in the relationship between field potentials and BOLD as evidence for a difference in neurovascular coupling, since they show that a single model can account for a variety of patterns. More generally, the V1 versus V2/V3 discrepancies bolster the argument that one cannot predict the exact relationship between BOLD and field potentials without also specifying the neuronal population activity.

The role of a simple model in understanding the relation between BOLD and LFP
A complete description of the biophysical processes giving rise to the BOLD signal and the field potential is far beyond the scope of this paper, and is likely premature given the enormous complexity in the nervous system, the vascular system, and the coupling mechanisms between them. Instead, the purpose of our modeling framework was to first begin with a general principle, namely that BOLD and field potentials sum neural activity according to a different sequence of operations; second, to instantiate this principle in simple mathematical rules; third, to combine these rules with a minimal model of neural population activity; and finally, to ask to what extent such a model can account for the patterns in our data. Our model omits many biophysical components, such as compartmentalized neurons, multiple cell types and vessel types, neurotransmitters, the dynamics of blood flow, and so on, and hence it is not a detailed simulation of the nervous system or vascular system. On the other hand, the simplicity of the model facilitates an understanding derived from basic principles, similar to the advantages in building computational, rather than biophysical, models of neural responses [75][76][77][78]. Both types of models and empirical studies are valuable. Here we emphasize that even with a highly simplified model of the BOLD signal, the field potential, and neuronal population activity, we are able to reconcile a wide range of findings in a complicated and technical literature. The model accounts for differences in how broadband field potentials and gamma oscillations relate to the BOLD signal. It can explain differences between cortical areas in the relationship between field potentials and BOLD. The model also provides an explanation for why the amplitude of alpha rhythms is negatively correlated with BOLD, even after accounting for the relationship between broadband signals and BOLD.

Reproducible computations
To test competing computational theories about the relation between the visual input, the LFP and the BOLD response, it is essential to make sample data and code available for others [43,48].
Following standard practices of reproducible research [79][80][81], the Matlab code of the simulation, and sample data and code to reproduce the Figs in this manuscript can be downloaded at https://github.com/WinawerLab/BOLD_LFP.

Conclusions
To understand how the electrophysiology and BOLD responses are related, it is necessary to specify both the manner in which population activity transfers to the two signals, and the neuronal population activity itself. The former shows that the covariance between neuronal time series has a large influence on the field potential and not the BOLD signal. Based on our simulations and empirical results, we made several inferences about the neuronal population responses mediating the BOLD signal and the LFP: that narrowband gamma oscillations in visual cortex likely arise more from synchronization of neural responses than a change in level of the neural response, and hence have a large influence on the field potential and little influence on the BOLD signal; that responses which are asynchronous across neurons manifest in broadband field potentials and an elevated BOLD signal; and that low frequency oscillations observed in field potentials are likely accompanied by a widespread hyperpolarization, which in turn reduces metabolic demand and the BOLD signal.
Our model-based approach brings us a step closer to a general solution to the question of how neural activity relates to the BOLD signal.

Simulated neuronal time series
Simulations were computed for a population of 200 neurons. Each simulation trial was 1 second long with millisecond sampling. The time series for each neuron was derived by summing three inputs, each 1 second long, followed by leaky integration with a time scale of 10 ms to simulate temporal integration in the dendrite (Fig 2). Each simulation was fit to ECoG data from one electrode and consisted of 240 trials, 8 repeats of 30 stimulus conditions. A condition in the simulation was defined by the parameter settings for the 3 inputs (Fig 4): C 1 (broadband), C 2 (gamma) and C 3 (alpha). See Supplementary Methods for details. Variations in these three inputs resulted in power changes in the broadband, gamma, and alpha LFP. The inputs were fit to data such the simulated LFP power changes matched the ECoG data power changes for a particular electrode and stimulus.

C 1 -Broadband input
The first input was a series of random numbers drawn from a normal distribution, with no temporal dependencies and no dependencies between neurons.
Motivation. This input approximates spike arrivals with a Poisson distribution at a fixed rate for a given 1-s trial. A random normal distribution was used rather than a Poisson distribution for computational efficiency. (The pattern of results is the same for either distribution.) The input has a flat (white) power spectrum up to the sampling limit of 500 Hz. When coupled with leaky temporal integration (described in a subsequent section), this input results in a power spectrum that is approximately proportional to 1 ! (brown noise). Several groups have proposed that the approximately 1 ! power spectra observed in field potentials arises from white noise (or Poisson noise) input to individual neurons, coupled to one or more low-pass filters [34,82,83]. Previously proposed sources of filters include an exponentially decaying current response in the synapse following each spike arrival [82], leaky temporal integration in the dendrite [34], and frequency dependent propagation in the extracellular tissue [84], the last of which has since been shown to be unlikely [85]. Regardless of the source of the low-pass filtering, the general proposal makes an interesting prediction, namely that a spectrally broadband increase in field potential power in response to a stimulus is likely to indicate an increase in the rate of spike arrivals following that stimulus [34]. This hypothesis has empirical support, based on correlations between spike rates (single unit and multiunit) and broadband field potentials [49,53], and the fact that a 1 ! baseline spectrum, as well as stimulus-dependent broadband power increases, can be observed in intracellular (single neuron) membrane potentials in awake macaque visual cortex [36]. This hypothesis is the logic behind our choice to model both the baseline 1 ! spectrum and stimulusdependent broadband modulations as arising from spectrally flat inputs followed by low-pass filtering within individual neurons. For computational tractability, we explicitly modeled only one of the low pass filters -leaky integration in the dendrites. We assumed that spectrally broadband signals reflect uncorrelated activity. First, we have shown that the broadband ECoG signal is asynchronous with respect to a visual stimulus, and hence uncorrelated from trial to trial [11].
Here, we extrapolate that within a trial, the contribution to the broadband signal is asynchronous from neuron to neuron. One reason to assume so is based on a physiological model: the broadband signal has been hypothesized to arise from the leaky integration of Poisson distributed spike arrivals [34]. Even if the spike rate is correlated between neurons, the spike timing within a trial is likely have low correlations between neurons.
Parameters. For each simulation the Gaussian distribution defining C 1 always had a mean, = 0.25. The slightly positive mean ensured that in the baseline state, the membrane potential was slightly positive, such that a suppressive signal (described in the section C 3 ) could bring the potential closer to 0, hence reducing the metabolic demand. For the 8 conditions in each simulation, the baseline standard deviation of the distribution was set at = 0.3. A larger results in a larger broadband signal, and can be thought of as reflecting a higher Poisson rate of spike arrivals. The for each condition was calibrated such that the resulting changes in broadband power for each of the 8 stimulus conditions matched the changes in broadband power in the ECoG data (Fig 4,

C 2 -Narrowband oscillations in the gamma band
The second input was band-passed filtered white noise. The white noise was drawn from a distribution with zero mean and fixed standard deviation on all trials and for all neurons, and subsequently band-pass filtered. Unlike C 1 , there were dependencies (coherence) between neurons. The level of coherence varied across the 8 trial types in each simulation.
Motivation. This input approximates a circuit producing narrowband gamma oscillations in the field potential. Parvalbumin positive interneurons project to pyramidal neurons, and can produce fluctuations in the membrane potential of the pyramidal neurons in gamma frequencies from 30-80 Hz [41]. The narrowband rise in gamma power associated with certain stimuli or tasks appears to reflect an increase in synchrony between neurons in this band [86]. Therefore, unlike C 1 , which varied in level but not coherence as a function of condition, C 2 varied in coherence but not level as a function of condition.
Parameters. For all trials and all conditions, the white noise samples were drawn from a normal distribution with = 0 and = 0.2. The covariance of the distributions could range between 0 and 1 (using Matlab's mvnrnd function). The white noise inputs were filtered between 50 Hz and 60 Hz prior to temporal integration in the dendrite: inputs were first zero-padded, then filtered with a 10 th order Butterworth filter in forward and reverse direction. (Fig 4). The covariance for each simulation was calibrated such that the resulting changes in narrowband gamma power for each of the 8 stimulus conditions matched the changes in narrowband gamma power in the ECoG data.

C 3 -Narrowband oscillations in the alpha band
The third input was band-passed filtered white noise, with an added asymmetry such that increased power decreased the mean amplitude. The coherence was the same for all trials and all neurons; the amplitude of the pulses varied by condition.
Motivation. Oscillations in the alpha band (8)(9)(10)(11)(12)(13)(14)(15) are widely observed in visual cortex, with higher amplitudes associated with low sensory stimulation (e.g., eyes closed or zero contrast) or a low level attention. One model of alpha oscillations is that pyramidal neurons in visual cortex receive pulses of inhibition (hyperpolarizing inputs) spaced on the order of 100 ms, generated indirectly by thalamic-cortical loops [44,45]. According to this view, less active states are associated with larger inhibitory pulses, resulting in a time-averaged hyperpolarization, compared to more active states with smaller inhibitory pulses. The inhibition is pulsed rather than continuous, so that the reduced cortical responsiveness is dependent on the phase of the alpha cycle (most reduced following each inhibitory pulse). Individual neurons in visual cortex can indeed show membrane oscillations at frequencies around 10 Hz [46], indicating that it is reasonable to model the alpha oscillation measured in the population as arising from oscillations in individual neurons, rather than arising only from band-limited coherence between neurons.
Parameters. For all trials and all conditions, the white noise samples were drawn from a normal distribution with = 0 and = 1. The covariance of the distributions was fixed at .75 (using Matlab's mvnrnd function). The white noise inputs were filtered between 9 Hz and 12 Hz prior to temporal integration in the dendrite: inputs were first zero-padded, then filtered with a 10 th order Butterworth filter in forward and reverse direction. The envelope was calculated by a Hilbert transform, and added to the signal, and the signal was multiplied by -1, such that increases in power reduced the mean amplitude. The signal was then multiplied by a factor c3, which was calibrated such that the resulting changes in narrowband alpha power for each of the 8 stimulus conditions matched the changes in narrowband alpha power in the ECoG data.

Fitting the LFP power changes to the ECoG power changes
Changing inputs in C 1 , C 2 and C 3 results in a change in LFP power in broadband, gamma and alpha respectively. In order to fit the simulated LFP power changes to the ECoG power changes we quantified the input to LFP output relation, such that a certain change in simulated LFP power could be predicted by change in the input amplitude. Different functions described the relation between the input and LFP. The relation between broadband input and LFP broadband was described as Since gamma and alpha were dependent on broadband amplitude (an increase in broadband noise masks the relative contribution of narrowband oscillations) the following function described the relation between input and gamma and alpha LFP: Parameters a, b, c, d and m were estimated by a separate calibration procedure in which C 1 , C 2 and C 3 were systematically varied and LFP broadband, gamma and alpha were calculated. Supplementary Fig S9 shows that using this procedure the simulated LFP power changes match the ECoG power changes well.

Stimuli and task
Stimuli for ECoG experiments were reported previously [43]. In brief, for one subject, the stimuli came from 8 classes of patterns (30 exemplars per class, 20x20°), including high contrast vertical gratings (0.16, 0.33, 0.65, or 1.3 cycles per degree square wave) noise patterns (spectral power distributions of / ! , / ! , and / ! ), and a blank screen at mean luminance (Supplementary methods and Supplementary Fig S2). For the second ECoG subject, there were the same 8 classes as well as two other stimulus classes -a high contrast white noise pattern and a plaid at 0.65 cpd. The fMRI subjects had the same 10 stimulus classes as the second ECoG subject.

ECoG task
ECoG data were re-analyzed from a previous report [43]. We briefly summarize that experiment here. Subjects viewed static images of gratings and noise patterns for 500 ms each, with 500 ms of zero-contrast (mean luminance) between successive stimuli. Order of presentation was randomized ( Supplementary Fig S2). There were a total of 210 contrast stimuli, shown once each in a single, continuous experiment (and 210 interstimulus blanks). Stimuli were shown on a 15-inch MacBook Pro laptop using Psychtoolbox (http://psychtoolbox.org/). The laptop was placed 60 cm from the subject's eyes at chest level. Screen resolution was 1280x800 pixels (33x21 cm).
Coordinates of the population Receptive Fields (pRF) were obtained from a prior study [11].

fMRI task
The fMRI experiment was a block design, with 12-second stimulus blocks alternating with 12-s blank periods (mean luminance). During the stimulus blocks, images were presented at the same rate as the ECoG experiment: 500 ms duration alternating with 500 ms of zero-contrast (mean luminance) between images (Supplementary Fig S2). All stimuli from each block came from one of the 7 stimulus classes. The exemplars within the block were all different. Subjects participated in 8 scans of 9 blocks each, and block order was randomized using Latin squares. Two subjects (S2 and S3) additionally participated in an identical experiment using lower contrast images, resulting in similar findings. fMRI subjects participated in two pRF runs to identify retinotopy maps. Stimuli for the pRF experiments consisted of a bar (width= 3 deg) that swept across the visual field in 8 directions: the four cardinal directions and the four diagonals. The bar contained a drifting checkerboard with 100% contrast. Images were projected on a screen in the rear of the magnet bore using an LCD projector (LC-XG250, Eiki) with a resolution of 1024x768 (60 Hz refresh rate) and subtending approximately 32x24 visual degrees (32.4x24.3 cm). Subjects viewed the screen with a mirror mounted to the RF coil. The viewing distance was approximately 58 cm.

ECoG procedure
ECoG data were measured from two subjects who were implanted with subdural electrodes (2.3 mm diameter, AdTech Medical Instrument Corp) for clinical purposes at Stanford Hospital. Informed, written consent was obtained from all subjects. ECoG protocols were approved by the Stanford University IRB. In 22 channels in V1 V2 and V3, broadband and narrowband gamma responses were quantified as before [43], and alpha power changes were calculated (see Supplementary Methods).

ECoG recording
ECoG data were recorded at 3052/1528 Hz (ECoG subject 1/ECoG subject 2) from 118/96 electrodes through a 128-channel Tucker Davis Technologies recording system (http://www.tdt.com). Electrodes were localized on a postoperative computer tomography (CT) scan that was co-registered with a pre-operative MRI, and locations were corrected for the brain shift [6]. Electrodes that showed large artifacts or epileptic activity (as determined by the patient's neurologist) were excluded from analysis (7/35 electrodes were excluded in subject 1/subject 2). Off-line, data were re-referenced to the common average, low-pass filtered and resampled at 1000Hz for computational purposes using the Matlab resample function. Line noise was removed at 60, 120 and 180 Hz using a 3 rd order Butterworth filter.

ECoG analyses
Broadband and narrowband gamma responses were quantified as before [43]. We calculated power spectra and separated ECoG responses into broadband and narrowband gamma band spectral power increases. To control for the influence of evoked activity on the spectrum, event related potentials (ERPs) were calculated per condition and the condition specific ERP was regressed from each trial. This procedure makes sure that the broadband increase is not due to a sharp edge in the ERP; the same pattern of results is obtained if this step is omitted. For each condition, the average power spectral density was calculated every 1 Hz by Welch's method [87] from 0 -500 ms after stimulus onset (and 0-500 ms after stimulus offset for the baseline) and a 250 ms Hann window to attenuate edge effects. ECoG power spectra are known to obey a power law and to capture broadband and narrowband gamma increases separately the following function was fitted to the average log spectrum from 35 to 200Hz (leaving out 60Hz line noise and harmonics) from each condition ( Fig 5A): ( ) = ( !"#$%!$&% − x) + !"##$%&"!' (x| , ) In which, The slope of the log-log spectral power function (n) was fixed for each electrode by fitting it based on the average power spectrum of the baseline. For cross-validation, trials were split into even and odd repeats, and broadband and gamma changes were calculated for each. Confidence intervals were calculated by a bootstrap procedure. For each condition C with Nc trials, Nc trials were drawn randomly with replacement and power spectra were averaged. The parameters β were fitted on the average log power spectrum from these bootstrapped trials. This was repeated 100 times, resulting in two sets of distributions of broadband and gamma weights for even and odd trials.
Alpha response amplitude was calculated as follows. Alpha changes are best visible after the initial onset transient and ERP, and we used the power from 250-500 ms to calculate the alpha decreases for each stimulus. Alpha amplitude was calculated by averaging the log-power between 8 and 13 Hz.

Channel selection
Channels for analysis were selected on the basis of three criteria. First, the pRF was located within V1, V2, and V3. Second, the explained variance in a pRF experiment was >15% [11]. Third, the center of the pRF was within the extent of the stimulus (<12 deg) and on the contralateral visual field. Because ECoG subject 2 did not have pRF data, only anatomical estimates of V1, V2, and V3 were used [88]. These criteria yielded 22 channels (19 from ECoG S1, 3 from ECoG S2).

fMRI procedure
Functional MRI data was measured from four subjects (3 female, ages 22-42) with normal or corrected-to-normal vision at the Center for Brain Imaging at NYU. Informed, written consent was obtained from all subjects. The fMRI protocols were approved by the New York University IRB. Functional MRI data were preprocessed and analyzed using custom software (http://vistalab.stanford.edu/software) (see Supplementary Methods). Disc ROIs (radius = 2 mm) were defined in fMRI subjects to match the position of the electrodes in ECoG subjects using a combination of anatomy, pRF centers, and visual field maps. The similarity between the ROI position and electrode position was compared via visual inspection of anatomical images and pRF centers ( Supplementary Fig S3).

fMRI recording
Anatomical and functional MRI data was collected at the Center for Brain Imaging at NYU on a Siemens Allegra 3T head-only scanner with a Nova Medical transmit/receive coil (NMG11) and a Nova Medical phased array, 8-channel receive surface coil (NMSC072).
Two to three T1-weighted whole brain anatomical scans (MPRAGE sequence) were obtained for each subject (voxel size: 1x1x1 mm, TR: 2500 ms; TE: 3.93 ms, flip angle: 8 deg). Functional images were collected using gradient echo, echo-planar imaging (voxel size: 2x2x2 mm, 24 slices, TR: 1500 ms, TE: 30 ms, flip angle: 72 deg). Images were corrected for B0 field inhomogeneity during off-line image reconstruction using a separate field map measurement made half way through the scan session. Slice prescription was set approximately perpendicular to the calcarine sulcus, covering the occipital lobe. In addition, a T1-weighted inplane was collected with the same slice prescription to align functional images to the high-resolution anatomical images.

fMRI analysis
Preprocessing. Anatomical images were co-registered and segmented into gray/white matter voxels using FreeSurfer autosegmentation algorithm (surfer.nmr.mgh.harvard.edu). A 3D mesh of the cortical surface was inflated for ease of visualization. Functional data were preprocessed and analyzed using custom software (http://vistalab.stanford.edu/software). Data were slice-time corrected to adjust for differences in acquisition time among slices in the 1.5-sec frame. Data were motion corrected for both between-and within-scan motion. Finally, data were high-pass filtered for low frequency drift [89] by multiple moving average smoothing (2 iterations, 40 seconds). Data were then converted to percent signal change by dividing each voxel's signal by its mean signal. The first four frames of each run (6 sec) were discarded to allow longitudinal magnetization and the hemodynamic response to reach steady state.
Analysis. Noise was removed from the fMRI data using GLMdenoise, a variant of the standard GLM commonly used in fMRI analysis [90]. In brief, GLMdenoise derives noise regressors for each subject by performing principle components analysis on noise voxels that are unrelated to the task. The optimal number of noise regressors is selected based on improvement in cross-validated R 2 .
The final model is fitted to each voxel's time series and bootstrapped 100 times over 8 runs. Here, the predictors in the GLM were the nine image categories (4 gratings, 4 noise patterns, 1 plaid) and a blank period (a randomly assigned blank block). Voxel bootstraps were averaged across voxels within a region of interest (ROI). The resulting 100 bootstraps per ROI were vector-length normalized and averaged across subjects. The beta estimate for each condition is taken as the median averaged bootstrap and the standard error as one-half the 68% confidence interval.
Population receptive field (pRF) model. The pRF runs were analyzed by fitting a 2D Gaussian to each voxel, modeling its pRF [91]. The pRF is defined by center location (x,y coordinates) and spread (sigma). The resulting maps were used to define retinotopic areas V1, V2, and V3 as in [92].

Predicting fMRI signals directly from ECoG models
The relationship between fMRI and ECoG signals was analyzed using a linear regression model. The cross-validated R2 was used as a metric for model quality and the regression coefficients were used to test whether ECoG predictors (broadband, gamma and alpha) had a positive or negative relation with BOLD (also see Supplementary Methods).
The relationship between fMRI and ECoG signals was analyzed using a linear regression model: Yi = XiBi + intercept, where Y is a vector of fMRI beta estimates for bootstrap i and X is a vector of ECoG broadband, gamma and alpha estimates for bootstrap i. The model was fitted separately for each cortical site (electrode/ROI pair), and for different combinations of predictors -broadband alone, gamma alone, alpha alone, each pairwise combination, and the all three predictors together. To measure the accuracy of the models, we calculated cross-validated coefficient of determination (R 2 ). First, the regression model Yi = XiBi + intercept was fitted for half of the fMRI subjects (1 and 2) and half of the ECoG stimulus repetitions (even repetitions) and then tested on the other half of the fMRI and ECoG data (fMRI subjects 3 and 4 and ECoG odd stimulus repetitions). All R 2 and ΔR 2 values reported in the results are cross-validated in this manner. For model comparison, the same pattern of results were achieved if instead of cross-validation, we solved the models on the complete data sets and computed the R 2 adjusted for the number of regressors.
To test whether different ECoG predictors (broadband, narrowband, alpha) had a positive or negative relation with BOLD, we tested whether the regression coefficient was significantly larger or smaller than zero. The regression coefficient was considered to be significantly different from zero using a two-sided, unpaired t-test across electrodes.