Functional magnetic resonance imaging (fMRI), with blood oxygenation level-dependent (BOLD) contrast, is a widely used technique for studying the human brain. However, it is an indirect measure of underlying neuronal activity and the processes that link this activity to BOLD signals are still a topic of much debate. In order to relate findings from fMRI research to other measures of neuronal activity it is vital to understand the underlying neurovascular coupling mechanism. Currently, there is no consensus on the relative roles of synaptic and spiking activity in the generation of the BOLD response. Here we designed a modelling framework to investigate different neurovascular coupling mechanisms. We use Electroencephalographic (EEG) and fMRI data from a visual stimulation task together with biophysically informed mathematical models describing how neuronal activity generates the BOLD signals. These models allow us to non-invasively infer the degree of local synaptic and spiking activity in the healthy human brain. In addition, we use Bayesian model comparison to decide between neurovascular coupling mechanisms. We show that the BOLD signal is dependent upon both the synaptic and spiking activity but that the relative contributions of these two inputs are dependent upon the underlying neuronal firing rate. When the underlying neuronal firing is low then the BOLD response is best explained by synaptic activity. However, when the neuronal firing rate is high then both synaptic and spiking activity are required to explain the BOLD signal.
Functional magnetic resonance imaging (fMRI) is an extensively employed neuroimaging technique that allows the non-invasive recordings from human brain of neuronal activity with relatively high spatial resolution. However, the blood oxygenation level-dependent (BOLD) contrast on which fMRI is based is only an indirect measure of this activity. The processes that link the underlying neuronal activity to the BOLD signals are still a topic of much debate. In particular, there is no consensus on the relative roles of synaptic and spiking activity in the generation of BOLD signals. In order to relate findings from fMRI research to other measures of neuronal activity it is important to understand the underlying neurovascular coupling mechanism [1]. Most of our present knowledge about neurovascular coupling comes from animal experiments. These studies have combined hemodynamic measures such as cerebral blood flow (CBF), with electrical measurements such as local field potentials (LFPs) and single/multi-unit activity (S/MUA). LFPs correspond primarily to weighted averages of synchronised dendro-somatic components of synaptic signals in a neuronal population, whilst S/MUA measures the action potentials of a single cell or population of cells, respectively [2]. In a pioneering study [3], found, in monkey visual cortex, that although both LFPs and MUA correlate with the BOLD response, this response could be predicted more accurately from LFPs. This result has been confirmed in awake animals [4]. On the other hand, [5], [6] and [7] found strong positive correlations between blood flow and spiking activity. More recently, [8], [9] and [10] have shown that when synaptic and spiking activity is uncoupled (by drug injection in [8], [10] and using a stimulus that elicits only synaptic activity in [9]), changes in CBF do not reflect underlying spiking activity and relate closer to LFPs. This growing body of evidence (Table 1) therefore supports the hypothesis that BOLD signals are more closely coupled to synaptic input and processing activity than to the output spikes of a population of neurons. In addition, this work (Table 1) provides support to a growing consensus in which the BOLD signal is thought to result from pre-synaptic activity and the release of neurotransmitters, in particular glutamate [11], as well as vasodilatory substances, such as nitric oxide [12], [13] and [14]. An increase in pre-synaptic activity and concomitant release of glutamate induces fluctuations in transmembrane potential at the post-synaptic neuron, and these fluctuations are measured with LFPs. This activity is also thought to be responsible for triggering the release of vasodilatory agents to the extracellular medium, which induce changes in blood flow and consequently the BOLD response [11]. However, when it comes to the human brain the number of studies directly addressing the question of how BOLD relates to synaptic versus spiking activity is relatively smaller (Table 1 and 2), and the data in these studies comes exclusively from neurosurgical patients, whose physiology may be compromised (Table 2). Of the few such studies [15], observe significant correlations between BOLD signals and both synaptic and spiking signals in auditory cortex, whilst [16] found no correlation between BOLD signals and neuronal firing in the hippocampal area. The link between neuronal activity and the BOLD response has not only been investigated at a microscopic level, using invasive co-localised recordings, but also at a macroscopic scale using fMRI and Electroencephalography (EEG). EEG (and Magnetoencephalography (MEG) ), are well established non-invasive techniques that are well suited to studying neuronal activity since they provide direct (not confounded by the hemodynamic response) measurement of post-synaptic potentials (magnetic fields) in cortical pyramidal cell populations with high temporal resolution [17]. Studies using both EEG and fMRI in humans have focused on correlations between BOLD signals and oscillatory EEG power measured in different frequency bands. For example, [18], [19] and [20] have shown that reductions in ongoing-scalp EEG alpha power (8–13 Hz) correlate with increases in BOLD activity in human occipital cortex. Using intra-cranial recordings in epileptic patients [21], have found a close spatial correspondence between regions of fMRI activation and sites showing EEG energy variation in the gamma band (). The main conclusion of this body of work is that increases in EEG frequency are associated with increases in BOLD signal. Even though these studies do not address our question (input versus output) directly, they seem to point in the direction of the biological hypothesis constructed from animal evidence (see above): increases in pre-synaptic activity, decrease effective membrane time-constants and result in faster oscillatory dynamics; at the same time more neurotransmitters are released (e. g. glutamate), which lead to increases in BOLD signal [14]. Here we design a powerful and efficient modelling framework to explicitly investigate competing hypotheses for the relationship between neuronal activity and the BOLD response in the healthy human brain. We use this framework to explore the relative contribution of synaptic and spiking activity to the generation of fMRI signals in visual cortex. The participation of healthy subjects prohibits the use of invasive electrophysiological measures. Therefore we use a mathematical modelling framework that allows us to non-invasively infer the degree of local synaptic and spiking activity, together with EEG-fMRI data, in which subjects were exposed to a reversing checkerboard of varying frequencies. This is similar in spirit to the use of ‘virtual electrodes’ in EEG analysis [22], but provides more specific biophysical information. This framework consists of a biophysically informed forward model from neuronal activity to the observed EEG and fMRI signals. Models linking neuronal activity to EEG/MEG signals have been proposed by [23], [24] and [25], to mention a few. These models usually use one or two state variables to represent the mean electrical activity of neuronal populations at the macro-column level, and are referred to as neural mass models [26]. Models linking ‘neuronal activity’ to BOLD signals include the metabolic models proposed by [27], [28] and the Balloon model, proposed by [29]. The Balloon model describes how evoked changes in blood flow are transformed into the BOLD response and has been extended by [30], who introduced a blood flow-inducing signal relating ‘neuronal activity’ and CBF and by [31], where different metabolic pathways have been proposed for supporting excitatory and inhibitory synaptic activity. In the above work ‘neuronal activity’ is usually not explicitly modelled and often corresponds to the stimulus input functions. Models linking a common underlying neuronal substratum to both EEG and fMRI signals have also been developed [32]. Some models are phenomenologically motivated, such as the ‘Heuristic’ proposed by [33]. This model aims to explain empirical results which relate frequency-specific power changes in EEG with fMRI signals and predicts that increases in the BOLD contrast reflect increases in the Root Mean Squared (RMS) frequency of EEG. We have validated these predictions in previous work [34] using simultaneous EEG-fMRI data in humans with a visual flicker stimulation task. As predicted by [33], the RMS frequency significantly explained more BOLD activity than the total time-varying spectral power or any linear combination of frequency-band amplitude modulations (e. g. alpha or gamma power). Biophysically motivated models include [35]–[37]. Most of these theoretical frameworks combine the neural mass model approach for EEG with the Balloon model for fMRI, but the coupling between neuronal activity and blood flow differs from model to model. For instance [35], propose that the squared post-synaptic membrane potential from both excitatory and inhibitory cells from a cortical area drives increases in cerebral blood flow, whilst [37] consider all the incoming action potentials from populations within and outside the voxel to be the input to the BOLD response. In [36] this input is proportional to the total concentration of nitric oxide (NO) synthesised by neurons in the cortical unit. The parameters of this model have been estimated using EEG-fMRI data from the visual cortex of one subject exposed to a reversing checkerboard with varying frequency [38]. Despite these theoretical efforts, the existing modelling frameworks have not yet been used in conjunction with real electrophysiological and hemodynamic data to compare different neurovascular coupling mechanisms, although important steps in this direction have been taken by [36], [39]. In [39], the authors have compared different models to investigate the role of excitatory and inhibitory activity in the generation of BOLD signals, using fMRI data from one subject. They found BOLD signals to be best explained by excitatory activity alone. Here we use the forward model proposed by [36] and embed it within a Bayesian framework. Using EEG and fMRI data in combination with Bayesian inference allows us to estimate the underlying synaptic and spiking activity, along with other biophysical model parameters. These quantities are computed using the variational Laplace method described in [40]. This optimisation scheme has been successfully applied to other input-state-output systems, such as [41], [42]. However, inverting generative models using multi-modality datasets, can be a technically demanding task, if the temporal characteristics of the datasets are very different, which is the case for EEG-fMRI data. Here we develop a computationally efficient scheme for model inversion. Instead of inverting the model in a single (computationally demanding) step we adopt a ‘multi-step inversion’ approach. This approach is based on partitioning model inversion into multiple, independent and computationally efficient steps that are motivated by the time-scales of data involved. This is a general procedure that can be used with other datasets and in other multimodal studies, such as with MEG-fMRI or LFP-fMRI data. Finally, once equiped with this mathematical and computational framework we posit models embodying different hypotheses about neurovascular coupling and adjudicate between them using Bayesian model evidence [43]. We compare three models. The first assumes that blood flow depends on the amount of vasodilatory substances (e. g. nitric oxide) released as a result of synaptic activity (synaptic input model), as proposed by [36]. The second assumes blood flow is driven by the firing rate of pyramidal cells from the same unit (spiking output model). These hypotheses are then compared against a third model where both these quantities contribute to the BOLD response (mixture model). In the long term, we anticipate that this modelling framework will be used to test neurovascular coupling hypotheses in a variety of experimental contexts with a range of subject cohorts.
We use a realistic biophysical model, proposed by [36], of how electrical and vascular dynamics are generated within a cortical unit. The unit comprises three subpopulations of cells: two layer IV GABAergic interneuron populations (the transmission and feedback interneurons (INs) ) and a layer V pyramidal cell (PC) population (Figure 1 (a) ). Interneurons are modelled as single compartment neurons, whilst the pyramidal cell has three compartments (soma, basal and apical tuft dendrites). Here we briefly describe the forward model. A summary of all the equations and parameters of the model can be found in Text S1. For a more detailed description please consult the original work [36]. A neural mass model (NMM) characterises the population dynamics of electrical states such as the membrane potentials in the somas of the neurons and electric currents flowing in the neuropil. This modelling framework is appropriate for data that reflect the behaviour of neuronal populations, such as EEG and fMRI data. The neural mass model can be viewed as a special case of ensemble density models, where the ensemble density is summarised with a single number representing mean activity [44]. Assuming that the equilibrium density of the neuronal states has a point mass (i. e. , a delta function), we can reduce the density dynamics to the location of that mass. What we are left with is a set of non-linear differential equations describing the evolution of this mode. The time variations of membrane potential in the individual compartments of the pyramidal cell and single compartment interneurons, , are determined by the differential equation for a simple voltage source circuit: (1) where is the effective membrane resistance of the compartment, and is cell-type and compartment specific. is the membrane time constant (same for all cells and compartments). The current, , that flows through the membrane of the cell depends on the connections between different elements of the cortical unit and its external inputs (Figure 1 (a) ). The cortical unit receives external excitatory input in different subpopulations, whilst its sole output is the firing rate of the pyramidal cells, . The excitatory inputs to the transmission interneuron, , and basal dendrites of the pyramidal cell, , correspond to thalamo-cortical afferent projections. The input to the apical tuft dendrites, , mediates cortico-cortical interactions. These currents can be found in Figure 1 (a). In terms of synaptic connections within the cortical unit, the total inhibitory synaptic effect on the pyramidal cell is given by: , where is the transmission inhibitory current and the feedback inhibitory current. The inhibitory synaptic currents depend nonlinearly on the membrane potential of the GABAergic cells through a threshold function: . The excitatory synaptic current generated by the pyramidal cell has the same form: : (2) The parameters are set to and to ensure that the output stays between and. The and parameters determine the voltage sensitivity by setting the membrane potential maximum growth and growth rate, respectively. These parameters are estimated from the data. determines the membrane potential near the asymptote where maximum growth occurs. The threshold function, , is also used to construct the firing rate coupling model (see below). The equations for the membrane potential at the soma of the three-compartment pyramidal cell, as well as the extracellular potential along its apical dendrites can be determined from the potentials and currents at the individual compartments (given by Eq. 1). These equations can be found in Text S1. The apical dendrites of the layer V pyramidal cells are arranged in parallel to each other and perpendicularly oriented to the surface of the cortex. This geometry facilitates the summation of electric currents in the neuropil. The mesoscopic effect resulting from the spatial average of these extracellular currents corresponds to the electrical signal measured with EEG. The state variables, , and parameters, , of the neural mass model described above are summarised in Table 3 of the main text and Tables 1 and 2 in Text S1. The coupling between local neuronal activity, described by the neural mass model, and subsequent changes in vascular dynamics is our question of interest. These changes are expressed in the BOLD signal and have previously been modelled in an extended Balloon approach [30], in which a set of four ordinary differential equations comprise the hemodynamic forward model from ‘neuronal activity’ to hemodynamic responses. The full derivation of these equations can be found in [29] and [30]. In brief, for a particular region, neuronal activity, , causes an increase in a vasodilatory signal, , that is subject to auto-regulatory feedback. Inflow, responds in proportion to this signal with concomitant changes in blood volume and deoxyhemoglobin content. These equations are summarised in Text S1. The hemodynamic parameters, , comprise the rate constant of the vasodilatory signal decay, the rate constant for autoregulatory feedback by blood flow, transit time, Grubb' s vessel stiffness exponent, and the resting oxygen extraction fraction, respectively. The whole dynamic system is driven by the input. Different inputs, , correspond to different aspects of neuronal activity and consequently different coupling hypotheses between neuronal activity and the BOLD response. A summary of the hemodynamic model' s state variables, , and parameters, , can be found in Table 3 of the main text and Tables 1 and 2 in Text S1. In the next section we specify the neurovascular coupling mechanisms we are interested in comparing. The original electro-vascular model proposed by [36] is represented by a set of stochastic differential equations describing the dynamics of the neuronal and vascular states, . In [36] the stochastic aspect of the model is instantiated by incorporating an additive multidimensional Wiener process to model physiological noise. In this paper, however, we use a deterministic version of the model. This means that the dynamics are completely determined by the state of the system and stochastic effects enter only at the observation level (Eq. 3). This deterministic approach resulted in very similar frequency-response curves to those in [36] (see Results: synthetic data below) and allows us to use standard Bayesian estimation routines, widely used with deterministic forward models for EEG (e. g. [42]) and fMRI (e. g. [40]). The observation equations for EEG, , and fMRI, , data are then given by: (3) where the errors are assumed to be i. i. d. , . The temporal variations of the EEG signal are well approximated by the extracellular electric current in the neuropil, , obtained from the NMM multiplied by the lead field matrix, . This matrix contains information about the geometry and conductivity of the head, and is therefore employed to map the distributed electric sources within the brain to scalp EEG recordings [45]: (4) The observation function for fMRI is a static nonlinear function of the cerebral blood volume and the concentration of deoxyhemoglobin directly [30]: (5) The factors, and are dimensionless but depend on the characteristics of the fMRI recording system. For 1. 5 T and TE of 40 msec, . is the resting blood volume fraction. To link the two main components of the biophysical model, the neural mass model and the Balloon model, we specified three different biologically plausible neurovascular coupling mechanisms based on previous empirical results. These mechanisms are described below: Using EEG-fMRI data in combination with Bayesian inference allows us to estimate the underlying synaptic and spiking activities, along with other parameters of the biophysical framework. Additionally, we can compare the different neurovascular coupling hypotheses using Bayesian model evidence. In Bayesian inference, prior beliefs about parameters, , of model are quantified by the prior density, . Inference on the parameters, , after observing data, , is based on the posterior density. These densities are related through Bayes' rule: (12) where is the probability of the data (likelihood) conditioned upon the model and its parameters. The normalisation factor, , is called the model evidence and plays a central role in model comparison (see below). The posterior density is an optimal combination of prior knowledge and new observations, weighted by their relative precision (i. e. , inverse variance), and provides a complete description of uncertainty about the parameters. Generally, the choice of priors reflects either empirical knowledge (e. g. , previous measurements) or formal considerations (e. g. , biological or physical constraints). Here we use empirical knowledge for both the neural mass model parameters and the coupling/hemodynamic parameters, based on estimates obtained by [36]. Under Gaussian assumptions, also known as a fixed-form Laplace approximation [55], the problem of estimating the posterior density reduces to finding its first two moments, the conditional mean and conditional covariance. The prior density is also assumed to be Gaussian with mean and covariance (see Table 3 for a list of prior mean values). A non-linear model, such as the local electro-vascular (LEV) model used here, Eq. (3), can be linearised by expanding the observation equation about a working estimate of the conditional mean: (13) such that, and. In this paper, the error covariance is assumed isotropic over the EEG and fMRI predictions. The linearised model, Eq. (13), can be used in a Variational Laplace (VL) optimisation scheme that iteratively updates the moments of the conditional density, . VL is a generic approach to estimate the posterior density, and can be formulated by analogy with statistical physics as a gradient ascent on the ‘negative Free Energy’, , of the system. The full derivation of the algorithm is described in [56]. The maximisation of with respect to in effect maximises a lower bound on the log model evidence, , [57]: (14) The model evidence is the probability of obtaining observed data, , given model, , and is at the heart of Bayesian Model Selection (BMS). The last term in Eq. (14) is the Kullback-Leibler (KL) divergence between the approximate posterior density, , and the true posterior, . This quantity is always positive, or zero when the densities are identical, and therefore is bounded below by. Through the iterative optimisation described above, the KL divergence is implicitly minimised and becomes an increasingly tighter lower bound on the log-model evidence. Model comparison can then proceed using as a surrogate for the log-model evidence. This approximation to the posterior density has been evaluated using Markov Chain Monte Carlo (MCMC) [58]. These schemes are more computationally intensive but allow one to estimate the posterior density without assuming it has a fixed form. Comparison between the model evidence obtained by MCMC methods and by variational approaches showed similar estimates, confirming that the approximations entailed by the variational approach lead to accurate model selection [55]. Again through Bayes' rule we can relate the model evidence to the model posterior probability, : (15) where is the prior distribution over models. Selecting the optimal model corresponds to choosing the model that maximises the posterior. If no model is favoured a priori then is a uniform distribution, and the model with the highest posterior probability is also the model with the highest evidence, . Given two models, and, we can compare these models using Bayes Factors, [59], which are defined as the ratio of the corresponding model evidences, or equivalently the difference in their log-evidences: (16) Bayes factors have been stratified into different ranges deemed to correspond to different strengths of evidence. ‘Strong’ evidence, for example, corresponds to a BF of over 20 (log-BF over 3) [59] in favour of model when compared to model. The equivalent posterior model probability is greater than 0. 95 [43]. Here we use Bayes factors to compare the neurovascular coupling models defined in the previous section.
In this section, simulations are used to explore the behaviour of the model and its ability to reproduce EEG and BOLD data under the experimental conditions described in the previous section. The response of the three neurovascular coupling models to changes in stimulus frequency is also shown. These synthetic signals are used to test the model inversion routines and to verify that Bayesian model comparison can be used to infer the correct coupling model. The LEV model was numerically integrated using the multi-step Adams-Bashforth-Moulton predictor-corrector algorithm implemented in the MATLAB (The MathWorks, Inc.) function ode113. The integration step used was 1 msec (1000 Hz) for the electrical and vascular states. The integrated signals were then downsampled to 100 Hz in the EEG case and to 0. 3 Hz for the BOLD signal. The input to the model is described below. Finally we fit the electro-vascular model with the three different coupling mechanisms to the EEG and fMRI data. We used the same ‘multi-step’ inversion procedure described in the previous section. Figure 9 shows the model predictions for EEG, as well as predictions of the coupling models and the BOLD response.
In this paper we used EEG-fMRI data and a biophysically informed mathematical model to investigate the relationship between neuronal activity and the BOLD signal in human visual cortex. In particular, we explored the contributions of synaptic input and spiking output activities to the generation of the BOLD response. We have provided preliminary evidence that the BOLD signal is dependent upon both synaptic and spiking activity but that the relative contribution of these two factors are dependent upon the underlying neuronal firing rate. When the underlying neuronal firing is low then BOLD signals are best explained by synaptic input, in agreement with previous animal studies, such as [3]. This result is also in line with more recent studies, such as [9] and [10], which show that the BOLD response is only affected by changes in synaptic-related activity (measured with LFPs) and not by changes in spiking activity (measured with MUA) when these two signals can be dissociated. However, when the neuronal firing rate is high then both synaptic and spiking activity are required to explain the BOLD signal, as observed in, for example [6], and [15]. We were particularly encouraged to find that a combination of synaptic input and spiking output frequency response curves (Figure 6) can explain the doubly-peaked BOLD response observed by [6] and replicated in our own data. One possible explanation for the increased performance of the output model with higher frequencies comes from neuroenergetic studies such as e. g. [7] and [63]. In these studies brain metabolism was found to depend strongly on neuronal spiking, with increases in oxygen consumption reflecting higher firing rates. More recently [64], have found that differences in the BOLD response between different brain areas (motor cortex and thalamus) could be explained by underlying differences in the firing rates of the corresponding neuronal populations. Our results also support the conclusion that the relationship between synaptic activity, spikes and BOLD signals depends on the specific neuronal circuitry engaged in task processing. Moreover, one can speculate that different coupling mechanisms involving different types of cells and molecules could come into play depending on the task in question. Despite our initial concern about the small number of fMRI samples compared to EEG, our initial results with synthetic data showed that it is possible to make inferences on different hypotheses for the neurovascular coupling using a generative modelling framework and Bayesian model comparison. The issue of different time-scales was addressed by partitioning the estimation of electrical and vascular states into a multi-step approach. In this approach we first estimated the electrical states and parameters from the EEG data and then integrated the full electro-vascular model using these estimates. From the integrated model we extracted the input time-series to the Balloon model, which we then inverted using BOLD data. The last two steps were repeated for each coupling model. This method significantly increases the computational efficiency of the model inversion. However, this multi-step approach is only possible with a deterministic model. In this work we used a deterministic version of the stochastic electro-vascular model proposed by [36]. Under different experimental conditions, which do not induce a large sensory response, the introduction of stochastic effects might be essential to reproduce the empirical data. In this case, other Bayesian inversion frameworks can be employed to estimate the model parameters, such as [65] and [66]. It is also worth noting that despite the fact that the mixture model had more parameters than the input and output models, this extra complexity did not provide a significantly better fit to the data in the low-frequency analysis than the input model. This complexity is correctly penalised using Bayesian methods, such as the one used here. One concern about the coupling models defined here regards the definition of NO concentration. As mentioned in the Methods section, NO is thought to have a pre-synaptic synthesis [12], [13]. However, here and in [36] the concentration of NO is modelled through post-synaptic quantities such as the transmembrane capacitive currents. Although in principle these two phenomena are directly related (increases in pre-synaptic activity mean larger post-synaptic effects) this is not always the case. Changes in transmembrane currents at the post-synaptic level can be caused by different processes such as chemical-gated channels, electric-gated channels, and passive leakage, not all of them being related to pre-synaptic activity. Therefore the transmembrane currents are an indirect way of quantifying the amount of NO released during synaptic activity. However, this issue is also encountered in experimental measures of synaptic activity, such as local field potentials. This signal is a surrogate post-synaptic signal, which is also affected by other slow potentials occurring at the cellular level that do not have a purely pre-synaptic origin. A natural extension to this work is the inclusion of multiple cortical units in the model representing multiple brain areas. For instance, sub-cortical areas such as the thalamus and other cortical areas activated by the experimental task could be included. Having more than one area would facilitate the differentiation between input and local processing synaptic activity, such as in [37]. In a recent study [67], have decomposed the effect of these two types of synaptic activity on hemodynamic signals by reducing the thalamic input to a rodent' s cortex. The authors found that although both input and local neuronal processing contribute to BOLD signals, as previously found, this contribution is larger from local processing. Another extension would be to probe the contribution of excitatory and inhibitory neuronal populations to the generation of BOLD signals, such as in [39]. This model-driven approach could, for instance, be used to study the findings of [68], where a negative BOLD response in deeper cortical layers, adjacent to positive-BOLD areas, was found to be associated with a reduction in local neuronal firing. Very recently [64], have optically driven genetically modified inhibitory cells and measured a negative BOLD signal in response to this stimulation, in the rat cortex. This result can inform the development of new generative models of neurovascular coupling. To our knowledge this paper presents the first quantitative model comparison of different biologically plausible mechanisms for neurovascular coupling in human cortex using EEG-fMRI data and a realistic biophysical model. However, even though our results were consistent across the three subjects and the majority of sessions, the case study approach adopted here has its limitations. Namely, it does not quantitatively address the issue of inter-subject variability and it therefore precludes inferences at the population level. With a larger sample of subjects, inter-subject variability can be accommodated using the Random-Effects (RFX) model selection approach developed by [69]. This approach fits a Bayesian hierarchical model to group model evidence data to obtain the frequencies with which each model is used in the population. This approach can be combined with the methodology developed in this paper. We hope that future studies with other datasets and different experimental conditions will employ our modeling approach so that a balance of evidence can be reached that clearly disambiguates between different hypotheses concerning neurovascular coupling. Understanding the underlying biophysical mechanisms behind the coupling between neuronal activity and the BOLD response is vital not only for improving the interpretability of the BOLD response, but also for relating findings from fMRI research with results from other neuroscientific disciplines.
Functional magnetic resonance imaging (fMRI), with blood oxygenation level-dependent (BOLD) contrast, is a widely used technique for studying the human brain. However, the relationship between neuronal activity and blood flow, the basis of fMRI, is still under much debate. A growing body of evidence from animal studies suggests that fMRI signals are more closely coupled to synaptic input activity than to the spiking output of a neuronal population. However, data from neurosurgical patients does not seem to support this view and this hypothesis hasn' t yet been tested in the healthy human brain. Here we design a powerful and efficient modelling framework that can be used to non-invasively compare different biologically plausible hypotheses of neurovascular coupling. We use this framework to explore the contribution of these two aspects of neuronal activity (synaptic and spiking) to the generation of hemodynamic signals in human visual cortex, with Electroencephalographic (EEG) -fMRI data. Our results provide preliminary evidence that depending on the frequency of the visual stimulus and underlying firing rate, fMRI relates closer to synaptic activity (low-frequencies) or to both synaptic and spiking activities (high-frequencies).
mathematics
applied mathematics
biology
computational biology
nonlinear dynamics
biophysics
neuroscience
Bayesian Comparison of Neurovascular Coupling Models Using EEG-fMRI