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 functional magnetic resonance imaging (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, V2, and 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.
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 magneto-encephalography (EEG and MEG). Integrating results across methods is challenging because the signals measured by these instruments differ in spatial and temporal sensitivity, as well as in the manner by which they combine the underlying neuronal population activity [1–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 blood-oxygen-level dependent (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 source-localized 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 approximately 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-to-measurement 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 Hz and 100 Hz (gamma band) are more highly correlated with BOLD amplitude than power changes in other bands [12–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–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 [13], 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 (LFP) 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–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 subpopulations 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. 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. By design, the model employs a minimal set of principles governing how the instruments pool neuronal activity, rather than a biophysically detailed description of neuronal and hemodynamic events. This approach enables us to ask whether this minimal set of principles is sufficient to guide simulations of neuronal population activity, such that the parameters of the simulations are fit to ECoG measurements from human visual cortex, and the output of the simulations predicts fMRI BOLD responses in the same regions for the same stimuli.
The fMRI BOLD signal and the 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, I (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 αi × Ii (t). The constant αi 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, α. These neurons act together like a single, equivalent circuit, and hence the LFP time series will sum the contribution from each neuron, LFP (t) =α∙∑inIi (t) (Eq 1) Field potential recordings are often 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: LFPpower=α∙∫0T ([∑inIi (t) ]2) dt]powerofsum (Eq 2) Importantly, Eq 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]. (Recent work on the neurobiology of neurovascular coupling indicates that much of the BOLD signal is not caused directly by changes in the level of metabolic products such as glucose, but rather by signaling molecules that tend to correlate with metabolic demand [33]). For simplicity, we discuss the BOLD signal throughout in terms of metabolic demand and return to this issue explicitly in the Discussion, (3. 6 The role of a simple model in understanding the relation between BOLD and LFP). The metabolic demand of each neuron will increase if the cell depolarizes (excitation) or hyperpolarizes (inhibition) [28,34,35]. 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: βi × (POWER (Ii (t) ) or βi×∫0T (Ii (t) 2) dt, with βi a scaling constant for the ith neuron. (Similar results were obtained if we 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: BOLD=β∙∑in (∫0T (Ii (t) 2) dt) sumofpower (Eq 3) Importantly, Eq 3 is an 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 Eq 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 two 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. 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 (Eq 4): ∫0T ([∑inXi (t) ]2) dt=∑in (∫0T (Xi (t) 2) dt) +∑i≠jn (∫0T (Xi (t) ∙Xj (t) ) dt) PowerofsumSumofpowerSumofcross-power (Eq 4) Applying this theorem to Eqs 2 and 3 shows the relationship between our models of BOLD and LFP power: LFPpower=αβ∙BOLD+α∙∑i≠jn (∫0T (Ii (t) ∙Ij (t) ) dt) (Eq 5) 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 (Eq 5). The first quantity summarizes the total level of neural activity (summed across neurons), and the second quantity 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 Eq 4 (summed power) or the second term (summed cross-power) will dominate. However, the number of elements contributing to the two terms is quite different: For n neurons, the first term has n numbers (the power in each neuron’s time series), whereas the second term has nearly n2 numbers (all the pairwise cross-powers). 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). To see how these equations translate to quantitative measures of BOLD and LFP, we consider a small neuronal population whose time series conform to a multivariate Gaussian distribution. We assume that each neuron’s time series has the same mean, m; the same variance, σ2; and all of the pairwise correlations have the same value, ρ: X∼N (μ, Σ) μ= (m⋮m) Σ=[σ2⋯σ2ρ⋮⋱⋮σ2ρ⋯σ2] (Eq 6) X is the population time series, μ is the mean of each time series, and ∑ is the covariance matrix. We can now rewrite the simulated BOLD signal (the sum of the power) and the LFP (power of the sum) in terms of the parameters of the multivariate Gaussian (and arbitrary scaling factors α, β), BOLD=β∙[n∙ (m2+σ2) ]LFPpower=α∙[n∙ (m2+σ2) + (n2−n) (m2+σ2ρ) ] (Eq 7) where n is the number of neurons. This enables us to visualize how the BOLD signal and the LFP power depend on just three values: the variance, correlation, and mean in the neural time series, rather than on all the individual time series (Fig 2). For these neuronal time series, the LFP, modeled as the power of the sum of neuronal time series (panel A), is dominated by the neuronal cross-power (panel C). The BOLD signal, modeled as the sum of the power in the neuronal time series (panel B), makes little contribution to the LFP, except when the correlation between neurons is low (ρ is close to 0); in this case, there is no cross-power, and BOLD and LFP power are correlated. In section 2. 1, we proposed formulae to derive instrument measures from neuronal population activity. Here, we ask how we might simulate neuronal activity with a small number of parameters. A low-dimensional characterization of the population activity is useful since we normally do not have access to the time series of an entire population of neurons. Moreover, a low-dimensional representation can lead to better understanding and generalization even when high-dimensional data are available [36,37]. After simulating the population activity, we then use the analytic framework from section 2. 1 to compute the BOLD and LFP signals. The parameters for the simulations were fit to ECoG recordings from human V1, V2, and V3 [38]. Because there were recordings from multiple electrodes and multiple stimuli, we ran multiple simulations fit to the different ECoG responses. We then used these simulations to predict the BOLD signal and compared these predictions to the measured BOLD signal for the same stimuli and same cortical locations (but in different observers). The steps for simulating the neuronal population data and the derived LFP and BOLD, and for comparing the simulations to empirical data, are summarized in Table 1. 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, V2 and V3, we decomposed the measured ECoG responses into three spectral components: broadband, narrowband gamma, and alpha (Fig 3). An important feature of this dataset is that the three components of the ECoG responses showed different patterns across stimuli [38]: 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 both resulted in decreases in alpha power compared to baseline (also see S1 Fig). Had all three responses been tightly correlated with each other, it would not be possible to infer how each relates separately to the BOLD signal. The simulations were structured to approximate the experimental design and the results of our ECoG experiments. To match the design of our ECoG experiments, a simulated experiment consisted of 240 trials, each of which were 1 second long (30 repeats of 8 conditions). The LFP time series were transformed to power spectra, which were averaged across the 30 repeated trials of the same condition. The simulation parameters—i. e. , the level of the three inputs, C1, C2, and C3—were fit to the measured ECoG summary metrics (broadband, gamma, and alpha) for each of the 8 conditions for a particular electrode (Fig 6). To verify the validity of this procedure, we asked whether the simulations using the fitted parameters produce simulated spectra, which, when analyzed like the ECoG spectra, reproduce the original values of broadband, gamma, and alpha. In other words, do we close the loop from measured spectral components (broadband, gamma, and alpha) to inferred input parameters (C1, C2, C3) to simulated population activity, to simulated spectral components (broadband, gamma and alpha)? The original values are not reproduced exactly because the simulations are stochastic, but overall, the original broadband, gamma, and alpha values are recovered with high accuracy (S9 Fig). As described above, the fitting of the parameters for C1, C2, and C3 was constrained by the assumptions that for C1, the correlation between neurons was 0 (and the amplitude was varied for fitting); for C2, the amplitude was fixed at a nonzero value (and the correlation was varied for fitting); and for C3, the correlation was fixed at a nonzero value (and the amplitude was varied for fitting). Results from alternative models with different constraints show poorer fits and are described briefly below and more extensively in S7 Fig. Importantly, the parameter fits did not take into account the measured BOLD responses. Hence the simulations provided a test: if the input parameters were chosen to produce outputs that match the measured ECoG responses (training data), does the simulated BOLD signal accurately predict the measured BOLD signal (test data)? We measured BOLD responses in four 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 (S2 Fig and S3 Fig). For an example V1 site, the predicted BOLD signal accurately matched the measured BOLD signal, with 89% of the variance in the measured BOLD signal explained by the prediction, as quantified by R2—the coefficient of determination (Fig 7A). Across V1 sites, the predicted BOLD signal from the simulations accounted for a median of 80% of the variance in the measured data (Fig 7C). For an example V2 site, the predicted BOLD signal also matched the measured BOLD signal (R2 = 0. 74, Fig 7B). Across V2/V3 sites the simulations explained a median of 40% of the variance in the data. The explained variance in V2/V3 is substantial but lower than in V1. One likely reason for the higher variance explained in V1 is that for the particular stimuli used in these experiments (gratings and noise patterns), the BOLD response reliability was higher in V1. For example, the median R2 computed by using half the BOLD data as a predictor for the other half (split half by subjects) was 86% for V1 and 63% for V2/V3. Similarly, the stimulus-evoked BOLD responses in V1 were larger than in V2 and V3, with more stimulus-related variance to explain: a mean of 1. 8% signal change in V1 versus 1. 2% in V2 and 0. 8% in V3 (S4 Fig). It is possible that a stimulus set more tailored to extrastriate areas, such as textures or more naturalistic scenes, would have evoked more reliable responses in extrastriate cortex. For each of the 22 simulations, the three input parameters C1, C2, and C3 defining each of the 8 stimulus conditions were fit to produce the LFP data from the corresponding ECoG electrode. By design, the C1 (broadband) and C3 (alpha) inputs were fit to ECoG data by varying the level per neuron, whereas C2 was fit to data by varying the correlation 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 correlation across neurons. For completeness, we tested all 8 combinations of models (S7 Fig). The most accurate model, quantified as the R2 between the measured BOLD and the simulated BOLD (median across the sites in V1 or in V2/V3), was the simulation type used in the main text, in which C1 and C3 varied in the level per neuron and C2 varied in the correlation across neurons. Models in which the broadband correlation rather than level was used to fit the ECoG broadband power were much less accurate. The models in which the gamma LFP power was fit by modulating the level rather than the correlation in the simulated population caused a small drop in R2. The previous analysis showed that when simulations were fit to ECoG data, the simulated BOLD response predicted the measured BOLD response. Here, we used regression analyses to assess how the simulated LFP predicted the simulated BOLD response and how the measured LFP predicted the measured BOLD response.
Many studies have reported correlations between BOLD and power in the gamma band LFP (30–130 Hz) (review for human studies: [52]). 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 [38,53,54]. Broadband changes have been proposed to reflect, approximately, the total level of Poisson-distributed spiking (or spike arrivals) in a local patch of cortex [40]. In contrast, the narrowband gamma response is caused by neural activity with a high level of cell-to-cell synchrony [55] and likely depends on specialized circuitry [56]. 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,57] so that the two signals can overlap in frequency bands. Hence, separating gamma band field potentials into an oscillatory component and a broadband (nonoscillatory) 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) [38]. 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 [54,58]. Second, unlike the case of narrowband gamma LFP, changes in broadband LFP are not accompanied by changes in broadband spike-field coupling (Figure 1A-B in [43]). The possibility that neuronal synchrony sometimes affects broadband signals cannot be ruled out, for example, as shown in cases of pharmacological manipulations in nonhuman primates [59]. 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 multiunit 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 (nonoscillatory) 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 [60]. 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 [61] and that, in some circumstances, multiunit activity predicts the BOLD signal more accurately than gamma band LFP [24,62]. 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 LFP. In fact, the elevation in broadband power is relatively small (2- or 3-fold) compared to the elevation in power often observed in narrowband gamma oscillations (10x or more) [38]. 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) in which 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 Eq 4, the portion of the field potential measurement that 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 (S7 Fig). 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 synchrony across the population rather than a change in the response level per neuron. This corresponds to the second term in the right-hand side of Eq 4, the portion of the field potential measurement that reflects the cross-power arising from currents in different neurons, 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 covary 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,63], and all three 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 [38], then gamma oscillations will not strongly correlate with BOLD. Our simulation and empirical results are consistent with studies that 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 [64]. We speculate that broadband fields spanning the gamma range would have shown a higher correlation with BOLD. 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 [65,66]. 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 Eq 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 [67]. Neural synchrony can also emerge without being time-locked to the stimulus, often called “induced synchrony” or “induced oscillations” [68]. 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 Eq 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 that 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 [47]. Moreover, the coherence between local spiking and the LFP also increases in the gamma band when LFP gamma power increases [43]. 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 (S7 Fig). Models in which gamma responses were caused by a change in level were only slightly less accurate and cannot be ruled out entirely (S7 Fig). However, the regression models fit to our data (Fig 9) 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. 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–15 Hz) was negatively correlated with the BOLD signal, independent of broadband signals. Numerous previous studies have reported that low-frequency oscillations are anticorrelated with BOLD, including measurements in motor, visual, and language areas [20–22,69–71]. 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 that 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 [72]. Alpha oscillations may indeed co-occur with reduced cortical excitation [73]. However, if this coupling between alpha power and spiking was 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 mean 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 mean shift comes from MEG and ECoG studies [49,50,74], which refer to alpha oscillations as being asymmetrical (i. e. , they are not centered at 0—there is a mean shift). This can be explained by a simple process: if alpha oscillations reflect periodic inhibitory pulses, then on average, they will cause a hyperpolarization (or less 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, 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 (Eq 3), and the mean contributes to power. The idea that a mean 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,75]. Moreover, if alpha oscillations are associated with a mean shift in membrane potential, this would explain why cortical excitability depends on the phase of the alpha cycle: at one phase, the membrane potential is more depolarized, and hence cortex is more excitable, and in the opposite phase, cortex is more hyperpolarized and hence less excitable. This is consistent with the observations that the threshold for eliciting a phosphene with TMS changes with alpha phase [76,77] and that the alpha phase at the time of stimulus presentation influences the size of the BOLD response in visual cortex [78]. Inhibition takes two neurons—one to inhibit and one to be inhibited. In our simulations, the alpha oscillations (C3) were associated with inhibitory fluctuations in the membrane potential (mean below 0), which in turn was associated with decreases in BOLD. It is important to note that these fluctuations are meant to capture the effect of local inhibition on the postsynaptic neurons (the neurons being inhibited). The inhibitory neurons themselves are presynaptic, and the action of inhibiting other neurons is presumably an active process that consumes energy. Therefore, inhibition is expected to increase energy demand in some neurons (the presynaptic neurons) and decrease energy demand in other neurons (postsynaptic neurons). We did not model the inhibitory neurons explicitly; however, the neural activity associated with active inhibition would be expected to contribute to the measured broadband signal in the ECoG data and is implicitly included in the broadband inputs in our simulations (C1). More complex models (see paragraph 3. 6) in which the circuitry of excitatory and inhibitory neurons is explicitly represented (such as [63,79,80]) may provide insight into how the balance between excitation and inhibition influences the field potential and the BOLD signal. 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–13 Hz) 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 [81], across cortical regions [82,83], between cortical layers [84], and as a function of medication [85]. 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. 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; hence, it is not a detailed simulation of the nervous system or vascular system. We modeled the BOLD signal as a function of dendritic currents summed across neurons within an imaging region. The logic motivating this is straightforward. Neuronal activity consumes a large amount of energy, and this energy demand is dominated by the cost of restoring the membrane potential following ionic flows from synaptic potentials and action potentials [86,87]. As a result, the increased energy demand from neuronal responses is related to the dendritic currents. Neurovascular coupling is the process of increasing blood flow to meet this energetic demand; a failure of the hemodynamic response such as stroke can cause neuronal and glial cell death, highlighting the importance of the relationship between blood supply regulation and neuronal activity [33]. We modeled the hemodynamic response as being proportional to the energy demand from dendritic currents. This model was proposed as a computational summary of the approximate relationship between the BOLD signal and neuronal activity, not as a hypothesis about a causal mechanism. Recent work suggests that energy consumption per se (e. g. , the change in the cerebral metabolic rate of oxygen consumption) is not the triggering mechanism for the increased blood flow, rather neurotransmitters and other molecules associated with synaptic events are part of a complex cascade that causes vasodilation and changes in blood flow [33,88,89]. The exact biological mechanism responsible for neurovascular coupling is an area of highly active, ongoing research [35]. The key assumptions in the model—that the BOLD signal is correlated with changes in membrane potential and that the order of operations differs for the BOLD signal and the LFP—makes accurate predictions for our dataset. A separate and important research question is how closely the biophysical mechanisms match this computational-level description, and what these mechanisms are. The simplicity of the model has benefits. It facilitates an understanding derived from basic principles, similar to the advantages in building computational, rather than biophysical, models of neural responses [90–93]. 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. We note that drastic simplifications are the norm in many fields of neuroscience, such as receptive field modeling of visual neurons; most such models omit fixational eye movements, optical properties of the eye, retinal and cortical circuitry, etc. , instead modeling responses as a few simple mathematical computations of the stimulus (filtering, thresholding, and normalization) [94]. These highly simplified models will certainly fail under some conditions [95], yet they have proven to be of immense value to the field [93], in part due to their simplicity and in part because the alternative in which the responses of visual neurons are computed from a complete, neurobiologically realistic model of the nervous system simply does not exist. 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 [38,53]. Following standard practices of reproducible research [96–98], 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/dorahermes/Paper_Hermes_2017_PLOSBiology. 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 that 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.
Informed, written consent was obtained from all subjects. The fMRI protocols were approved by the New York University IRB and ECoG protocols were approved by the Stanford University IRB, according to the principles expressed in the Declaration of Helsinki. 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 milliseconds to simulate temporal integration in the dendrite (Fig 4). 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 three inputs (Fig 4): C1 (broadband), C2 (gamma), and C3 (alpha). Variations in these three inputs resulted in power changes in the broadband, gamma, and alpha LFP. The inputs were fit to data such that the simulated LFP power changes matched the ECoG data power changes for a particular electrode and stimulus. Stimuli for ECoG experiments were reported previously [38]. 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 k/f4, k/f2, and k/f0), and a blank screen at mean luminance (S2 Fig). 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 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 electrodes in V1 V2 and V3, broadband and narrowband gamma responses were quantified as before [38], and alpha power changes were calculated. fMRI data was measured from four subjects (three 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. fMRI data were preprocessed and analyzed using custom software (http: //vistalab. stanford. edu/software). Disc regions of interest (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 (S3 Fig). The relationship between fMRI and ECoG signals was analyzed using a linear regression model. The cross-validated coefficient of determination (R2) was used as a metric for model accuracy, and the regression coefficients were used to test whether ECoG predictors (broadband, gamma, and alpha) had a positive or negative relation with BOLD. The relationship between fMRI and ECoG signals was analyzed using a linear regression model: y=Xb+c+ε where y is a vector of fMRI amplitudes (beta estimates), with n entries for the n different stimuli; X is a matrix of ECoG responses, n by 1,2, or 3, where the columns correspond to one or more of broadband, gamma, and alpha estimates; b are the 1,2, or 3 beta weights for the broadband, gamma, and alpha estimates; c is a constant (the y-intercept); and ε is the residual error term. 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 all three predictors together. The n stimuli in the regressions included the contrast patterns and the blank stimulus. Inclusion of the blank stimulus is important for capturing the sign of the mean response. For example, if all contrast patterns induced a BOLD response of a particular level (say, +1) and induced ECoG responses of a particular level (say, –1) and we did not include the blank stimulus in the regression, then after subtracting the mean from each measure, all beta estimates would be approximately 0. This would mask a systematic relationship between ECoG and BOLD measures (in this example, an anticorrelation) arising from viewing stimuli with contrast compared to viewing a blank screen. Models were evaluated by split-half cross-validation. First, the regression model yi = Xibi + ci + ε was fit using half of the fMRI subjects (1 and 2) and half of the ECoG stimulus repetitions (even repetitions). To cross-validate this model, the beta values (bi) were then applied to the left out half of the ECoG data (odd stimulus repetitions) to predict the left out half of the fMRI data (fMRI subjects 3 and 4). The same procedure was applied by reversing the training and testing data. This resulted in two testing datasets with BOLD responses predicted from ECoG for each stimulus condition (Xibi + ci) and an actual measured BOLD value. For each cortical site, the coefficient of determination (see below) was calculated between the concatenated predictions and BOLD data values of the two test sets. All R2 values reported in the results are cross-validated in this manner. The same pattern of results was achieved if instead of cross-validation, we solved the models on the complete datasets and computed the R2 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 0. The regression coefficient was considered to be significantly different from 0 using a bootstrap statistic: for each model, the median of the beta values across sites was calculated after resampling 10,000 times. If <2. 5% of the resampled statistics were smaller than zero, the beta values were considered significantly positive, and similarly, if <2. 5% of the resampled statistics were greater than 0, the beta values were considered significantly negative. All model predictions were quantified using the coefficient of determination on cross-validated predictions. For predicting BOLD data from simulations of population neuronal activity (Fig 7, S7 Fig), the predicted BOLD has arbitrary units. In these cases, the observed BOLD and the predicted BOLD were both normalized by subtracting the mean and then dividing by the vector length. When predicting BOLD responses from features of the LFP data (broadband, gamma, and alpha) by regression, the predicted BOLD data were in the same units as the measured BOLD, and no normalizing or rescaling was done. To quantify the accuracy of the models, we calculated the cross-validated coefficient of determination, R2: R2=1–SSresidualsSSdata SSresiduals=∑i (yi−fi) 2 SSdata=∑i (yi−y¯) 2 where y are the data values and f are the prediction values. Because the model fits are cross-validated, it is possible for the model errors (residuals) to be larger than the data values, hence R2 can be lower than 0, and spans (−∞, 1]. In the case in which the model predictions and the data are unrelated and each are normally distributed with equal variance, R2 will tend to –1.
There are several methods for measuring activity in the living human brain. Here, we studied functional magnetic resonance imaging (fMRI), which depends on the vascular response to neuronal activity, and surface field potentials, which measure electrical activity from many neurons. These two widely used measurements of human brain activity often provide different and potentially conflicting results. We propose a quantitative model for how these two measurements integrate activity from neuronal populations. The fMRI signal is highly sensitive to the average level of local neuronal activity but not the degree of synchrony between neurons. In contrast, the field potential is most sensitive to synchronous neuronal signals. Our model accounts for several observations seen in fMRI and field potential data: some very large features of field potential recordings, such as gamma oscillations, can occur with little to no associated fMRI signal. The model predicts this because the gamma oscillations result more from increased neuronal synchrony than increased neuronal activity. Other field potential signals, such as broadband changes, which are likely driven by the level of neuronal activity rather than a change in synchrony, are highly correlated with fMRI. The two measures thus provide complementary information about human brain activity.
medicine and health sciences
action potentials
diagnostic radiology
functional magnetic resonance imaging
membrane potential
brain
electrophysiology
neuroscience
magnetic resonance imaging
simulation and modeling
brain mapping
electrocorticography
neuronal dendrites
neuroimaging
research and analysis methods
imaging techniques
animal cells
visual cortex
cellular neuroscience
radiology and imaging
diagnostic medicine
cell biology
anatomy
physiology
neurons
biology and life sciences
cellular types
neurophysiology
Neuronal synchrony and the relation between the blood-oxygen-level dependent response and the local field potential