Extracting novel information from neuroimaging data using neural fields

Abstract

We showcase three case studies that illustrate how neural fields can be useful in the analysis of neuroimaging data. In particular, we argue that neural fields allow one to: (i) compare evidences for alternative hypotheses regarding neurobiological determinants of stimulus-specific response variability; (ii) make inferences about between subject variability in cortical function and microstructure using non-invasive data and (iii) estimate spatial parameters describing cortical sources, even without spatially resolved data.

Introduction

This paper reviews some recent applications of neural field models in the analysis of neuroimaging data. Neural fields model current fluxes as continuous processes on the cortical sheet, using partial differential equations (PDEs). The key advance that neural field models offer, over other population models (like neural masses), is that they embody spatial parameters (like the density and extent of lateral connections). This allows one to model responses not just in time but also over space. Conversely, these models are particularly useful for explaining observed cortical responses over different spatial scales; for example, with high-density recordings, at the epidural or intracortical level. However, the impact of spatially extensive dynamics is not restricted to expression over space but can also have profound effects on temporal (e.g., spectral) responses at one point (or averaged locally over the cortical surface) [1]. This means that neural field models may also play a key role in the modelling of non-invasive electrophysiological data that does not resolve spatial activity directly. In what follows, we attempt to shed light on these different uses of neural fields and put forward three reasons why these models can be useful in the analysis of neuroimaging data. Each of these motivations is demonstrated by analysing a particular dataset obtained using three different modalities: electrocorticography (ECoG), magnetoencephalography (MEG) and local field potential recordings (LFPs). We argue that neural fields allow one to: (i) compare evidences for alternative hypotheses regarding the important neurobiological determinants of stimulus-specific response variability; (ii) make inferences about between subject variability in cortical function and microstructure using non-invasive data and (iii) obtain estimates of spatial parameters describing cortical sources in the absence of spatially resolved data. Our analyses exploit dynamic causal modelling [2] and include model space explorations that embody different hypotheses about the generation of observed responses in relation to model evidence - obtained using Variational Bayes [3]. This model comparison uses a variational free-energy bound to furnish optimized models in a manner similar to fitting empirical spectra with AR and ARMA models, see e.g. [4]. The advantage this approach has over other optimization criteria is that it provides an optimal balance between model fit and complexity; yielding models that are both parsimonious and accurate. The analyses presented here showcase particular instances where neural field models serve as a mathematical microscope, allowing us to extract information that is hidden in electrophysiological data.

Review

In the following, we present analyses of three different datasets obtained with ECoG, MEG and LFP recordings respectively; these illustrate how neural fields combined with dynamic causal modelling can disclose important neurobiological properties of cortical microcircuitry to which we do not have direct access through conventional analysis techniques.

Comparing alternative hypotheses about the determinants of stimulus-specific gamma peak variability

In our first study [5], we exploited high density spatially resolved data from multi-electrode electrocorticographic (ECoG) arrays to study the effect of varying stimulus contrast on cortical excitability and gamma peak frequency. These data were obtained at the lab of Pascal Fries at the Ernst Strungmann Institute for Neuroscience in collaboration with the Max Planck Society in Frankfurt. Details of the task and the surgical procedures are described elsewhere [6, 7]. Neuronal signals were recorded during a task where a monkey first fixated its gaze (in a window with 1° radius) and later released a lever when it detected a colour change at fixation (see Figure 1). During the fixation period, a stimulus was presented at 4° eccentricity. Crucially, the contrast of the stimulus varied between trials.

Figure 1

Experimental Task . The monkey released a lever when it detected a colour change at fixation. The colour change could occur at any moment in time following fixation onset, starting from 200 ms post-fixation onset to 3 seconds post fixation onset. While the monkey fixated, a stimulus located at 4° of eccentricity was displayed (a physically isoluminant sinusoidal grating; diameter: 1.2°; spatial frequency: 0.4-0.8 cycles/deg; drift velocity: 0.6 deg/s drifting within a circular aperture). Stimulus contrast varied between zero and an 82% contrast value.

We analysed the activity before the fixation colour change and the subsequent behavioural response. We obtained model predictions using a canonical microcircuit field model introduced in [8], see also [9] and Figure 2. Table 1 describes the priors over model parameters. During model optimisation, these prior means are multiplied by scale factors with lognormal priors.

Figure 2

The Canonical Microcircuit (CMC) neural mass model . This figure shows the evolution equations that specify a CMC field model of a single source. This model contains four populations occupying different cortical layers: the pyramidal cell population of the Jansen and Rit (JR) model is here split into two subpopulations allowing a separation of the sources of forward and backward connections in cortical hierarchies. As with the JR model, second-order differential equations mediate a linear convolution of presynaptic activity to produce postsynaptic depolarization. This depolarization gives rise to firing rates within each sub-population that provide inputs to other populations.

Table 1 Prior expectations of model parameters

The canonical microcircuit model we used conforms to known neuroanatomy and produces observed differences in gamma responses that reflect varying stimulus conditions. It provides a probabilistic mapping from hidden neuronal states to observed spectra, which are modelled as a mixture of induced activity, channel noise and sample noise. The schematic in Figure 2 summarizes the equations of motion or state equations that specify a canonical microcircuit neural mass model of a single source. This model contains four populations, each loosely associated with a specific cortical layer. The second-order differential equations describe changes in hidden states (e.g., voltage) that subtend observed local field potentials or EEG signals. These differential equations effectively mediate a linear convolution of presynaptic activity to produce postsynaptic depolarization. Average firing rates within each sub-population are then transformed through a nonlinear (sigmoid) voltage-firing rate function σ(⋅) to provide inputs to other populations. These inputs are weighted by connection strengths. In Figure 3, electrophysiological responses and model predictions are shown in dashed and solid lines. The real and imaginary parts are shown in the left and right panels respectively. In this figure, lines of different colours correspond to different contrast conditions, where the peak frequency ranges between 48 and 60 Hz and contrast levels vary between 0 and 82% of some maximum value. There are 9 contrast conditions; where in the first two there is no prominent gamma peak. It is known that gamma peak frequency is highly contrast-dependent, see e.g. [10, 11]. This is in accordance to what we observe here, namely a monotonic relation between peak frequency and stimulus contrast.

Figure 3

Model predictions and data for all conditions simultaneously. Empirical responses (cross spectral densities) and model fits are shown in dashed and solid lines. The real and imaginary parts of model fits and observed spectral responses are shown in the left and right panels respectively. The two sets of curves in the left panel correspond to auto- and cross-spectral densities from bipolar sites. Subsequent model comparisons rely on model evidence in data from all conditions.

Bipolar differences were extracted from ECoG sensors covering a large part of the primate brain. From the modeller’s vantage point, the important issue here is that a single source in the visual cortex is sampled by several sensors. In our approach we have allowed for the optimization of both the parameters governing topographic aspects of cortical activity and the deployment of the lead fields (their dispersion and location). In brief, we endowed our observation model with sufficient degrees of freedom to characterize unknown aspects of sensor sensitivity.

In this study, we wanted to delineate important mechanisms underlying contrast effects. We tested these effects by allowing for trial-specific changes in model parameters. In particular, we asked which model can best explain spectral responses for different contrast conditions, see Figure 4. Our candidate DCMs allowed for contrast dependent changes in the strength of recurrent connections, the strength of intrinsic connections or the extent of intrinsic connections (or combinations thereof). The ability of each model to explain induced responses was evaluated in terms of their Bayesian model evidence; which provides a principled way to evaluate competing or complementary models or hypotheses.

Figure 4

Candidate models. We consider three putative mechanisms of gain control that lead to models or hypotheses about trial-specific modulations of (G:) recurrent connections of neuronal populations, (L:) horizontal connections between distinct neuronal subpopulations and (E:) spatial dispersion of horizontal connections.

The first set of parameters comprised the gains of neural populations that are thought to encode precision errors. These gain parameters correspond to the precision (inverse variance or reliability) of prediction errors in predictive coding. This fits comfortably with neurobiological theories of attention under predictive coding and the hypothesis that contrast manipulation amounts to changes in the precision of sensory input, as reviewed in [12]. These changes affect hierarchical processing in the brain and control the interaction between bottom up sensory information and top down modulation: here, we focus on the sensory level of the visual hierarchy.

The second set of parameters comprised intrinsic connection strengths among pyramidal and spiny stellate cells and inhibitory interneurons. This speaks to variations in cortical excitability, which modulates the contributions of different neuronal populations under varying contrast conditions – and a fixed dispersion of lateral connections. This hypothesis fits comfortably well with studies focusing on the activation of reciprocally connected networks of excitatory and inhibitory neurons [13–15].

The last set included the spatial extent of excitatory and inhibitory connections. From single cell recordings, it is known that as the boundary between classical and non-classical receptive field depends crucially upon contrast [16]. As stimulus contrast is decreased, the excitatory region becomes larger and the inhibitory flanks become smaller. The hypothesis here rests on a modulation of the effective spatial extent of lateral connections; effective extent refers to the extent of neuronal populations that subserve stimulus integration - that are differentially engaged depending on stimulus properties (as opposed to the anatomical extent of connections).

In Figure 5, we see the results of Bayesian model comparison [17] between seven candidate models. These are the three models depicted in Figure 3 and their combinations. Model comparison pools the evidence for models in data from all conditions. By construction, the free energy ℱ(G, q) ≈ ln p(G|m) approximates the log-evidence (marginal likelihood) of the data. In other words, the free energy reports the probability of obtaining the data under any given model. Variational free energy was introduced by Richard Feynman in the context of path integral formulations of quantum mechanics and has been used extensively in machine learning to finesse the difficult problem of exact Bayesian inference. The models we compared allowed only a subset of parameters to vary with contrast level, where each model corresponds to a hypothesis about contrast-specific effects on cortical responses. We found that the model involving modulations of all candidate parameters has the highest evidence (with a relative log evidence difference of 16 with respect to the model that allows for modulations of all but the extent parameters – model 6). This simply means that we can be confident that all three candidate mechanisms contribute to the modulation of spectral responses, which is consistent with a plethora of studies that have considered visual contrast effects in terms of adaptive (extraclassical) gain control mechanisms and predictive coding in the visual cortex.

Figure 5

Results of Bayesian Model Comparison . The model involving modulations of all parameters (model 7) has the highest evidence with a relative log evidence of 17 with respect to the model that allows for modulations in all but the extent parameters (model 6). The first three models correspond to hypotheses (i), (ii) and (iii) while models 4 and 5 to combinations (i) + (ii) and (ii) + (iii). In the bottom panel we include model posterior probabilities. These results suggest that we can be almost certain that all three synaptic mechanisms contribute to the formation of cross spectral data features.

Explaining inter-subject variability of peak gamma frequency of visually induced oscillations

Our second case study also involved spectral responses obtained from the visual cortex during a perception experiment [8]. This study was however based on non-invasive MEG data. These data were obtained during a task whose focus was on understanding how activity in the gamma band is related to visual perception [18]. In earlier work, the surface of the primary visual cortex was estimated via retinotopic mapping and was found to correlate with individual gamma peak frequencies. A similar visual MEG experiment found a correlation between gamma peak frequency and resting GABA concentration, as measured with MR spectroscopy [19].

Our focus here was on obtaining mechanistic explanations for this intersubject variability in peak gamma frequency, observed during visually induced oscillations. In particular, we wanted to disclose the origins of individual gamma peak variability and understand whether this variability can be attributed to cortical structure or function (the level of cortical inhibition as expressed by resting GABA concentration). The generative model we used in this study was similar to that of the study described above (but allowing for patchy horizontal connectivity, see [1, 20]). This neural field DCM was particularly useful to discern the roles of cortical anatomy and function as it parameterises both structure and functional excitation-inhibition balance. Our model included parameters describing the dispersion of lateral or horizontal connections which we associate with columnar width and kinetic parameters describing the synaptic drive that various populations are exposed to. We focussed on estimates of the parameters describing columnar width and the excitatory drive to inhibitory interneurons. We then looked at the three-way relation between these estimates, the size of V1 and peak gamma frequency.

Our results show that both structure and function contribute to intersubject variability in gamma peaks. We found correlations between columnar width and gamma peak; this confirmed the crucial role of cortical anatomy (i.e., macroscopic V1 size) in explaining individual differences, see Figure 6 (top panel). We also found correlations with the parameter describing the drive to inhibitory interneurons, see Figure 7. This underlined the importance of cortical inhibition. This correlation remained significant when controlling for V1 size and columnar width, suggesting that this correlation cannot by accounted for by spatial effects at the microscopic or macroscopic level. The important thing to note here is that we have estimated microscopic properties of cortical sources (microanatomy and local inhibition) non-invasively. We also found a significant correlation between columnar width and V1 surface which implies that subject with a larger V1 might have bigger macrocolumns, see Figure 6 (bottom panel). This result is consistent with the previous empirical results and establishes the construct validity of DCM in this neural field setting.

Figure 6

Correlations with columnar width. (Top panel): We found a correlation between the log scaling of the posterior estimate of columnar width and peak gamma frequency (Pearson r = 0.271, p = 0.06, 30 d.f., one-tailed test). This suggests that increases in peak gamma frequency across subjects can be attributed to a greater columnar width. (Bottom panel) Correlation between the log scaling of the posterior estimate of columnar width and V1 surface (Pearson r = 0.36, p = 0.02, 30 d.f., one-tailed test). This result suggests that a larger V1 is constituted by bigger macrocolumns. This finding, together with that reported in the top panel, confirms our earlier empirical finding that gamma peak and V1 size are correlated. Furthermore, our use of an underlying generative (mechanistic) model, offers insights into local cortical microstructure that would be hard (or impossible) to disclose otherwise. Confidence intervals are plotted as dotted lines.

Figure 7

Correlation between the log scaling of the connectivity estimate between pyramidal cells and inhibitory interneurons and peak gamma frequency (Pearson r = − 0.37, p = 0.02, 30 d.f.,two-tailed test). Confidence intervals are plotted as dotted lines.

Estimating the spatial properties of sources when there is no explicit spatial information

A common theme in the examples discussed so far has been the modelling of lateral interactions in sensory cortices. Both studies considered above exploited spatially resolved data (either invasive or non-invasive). However, it is possible to obtain estimates of parameters that describe the topographic properties of cortical sources like the extent of lateral connections, even when using spatially unresolved data, like data from a single LFP electrode. This is precisely what we did in this final case study [21]. This study considered LFP data obtained from the rat auditory cortex under anaesthesia. Local field potentials were recorded from primary (A1) auditory cortex in the Lister hooded rat, following the application of the anaesthetic agent Isoflurane; 1.4 mg (see [22] for details). Ten minutes of recordings were used to obtain frequency domain data-features g _Y (ω) using a vector autoregression model [23].

In this analysis we first established the identifiability of our model using simulated data to ensure that Bayesian inversion scheme was able to recover veridical estimates. Almost universally, the fitting or inversion of Dynamic Causal Models optimizes a free energy bound on the log-evidence for a model m. This bound is optimized with respect to a variational density q(θ) over unknown model parameters. By construction, the free energy bound ensures that when the variational density maximizes free energy, it approximates the true posterior density over parameters, q(θ) ≈ p(θ|G, m). The (approximate) conditional density and (approximate) log-evidence are used for inference on parameters and models respectively (see also Figure 8).

Figure 8

Summary of the approach. We assume the measured signal is a mixture of predicted spectra, channel noise (a mixture of white and pink noise with weights alpha and beta) and Gaussian observation noise epsilon. Bayesian inference uses a Laplace approximation to the posterior density over unknown model parameters to obtain approximate conditional density and log-evidence – that are then used for inference on parameters and models respectively.

Here, we examined the conditional mean of the transit time when the simulated log-scaling deviates from its prior expectation (of zero). This revealed how our model inversion operates in regions of parameter space that are remote from prior assumptions. Inferring axonal conduction speeds using single channel data illustrates our point that suitably informed models enable one to access the spatial attributes of neuronal infrastructures, even in the absence of spatially resolved data. We synthesized observed spectra by setting the values of the parameters equal to their prior expectations, with the exception of transit time, which was varied over a range of 36% to 272%. The synthetic data were inverted and the resulting conditional estimates are shown in terms of the conditional mean and 90% confidence intervals in Figure 9 (left panel). One can see that there is a remarkable agreement between the conditional mean and true values and that the (relatively precise) confidence intervals include the true value (with the exception of the most extreme deviations). This is an interesting result that establishes the identifiability of the model, at least in relation to conduction speed. Clearly this only establishes face validity (the inversion does what it is meant to) but does not speak to the physiological validity of the model, which is addressed using empirical data below. Having said this, this sort of validation speaks to the possibility of estimating small differences in conduction velocity that might be elicited experimentally through pharmacological or other manipulations.

Figure 9

(Left panel) Sensitivity analyses of synthetic data. This figure shows the conditional estimates of the inverse conduction speed parameter, υ. These estimates were obtained from simulated data sets generated using υ as per Table 1, with a log-scaling from −1 to 1. We illustrate the agreement between the true parameter value (dotted line) and its conditional estimate. The bars represent 90% conditional confidence intervals. (Right panel) Parameter estimates obtained using LFP data. We show here the conditional estimates (90% confidence intervals) of the neural field model parameters using empirical data from a single trial recording of auditory responses in the rat. The estimates lie close to the prior values (a log-scaling of zero), which were motivated by animal physiological studies [24].

Finally, we used the empirical LFP data described above to obtain conditional estimates of the (log-scaling of the) parameters. Figure 9 (right panel), presents these estimates in terms of their conditional means and 95% confidence intervals. Most confidence intervals include a log – scaling of zero, with the exception of inverse intrinsic conduction speed v and neuronal gain g, which increased by 160 and 120% respectively with respect to their prior values. This is not surprising as the priors were selected to reproduce spectra that are typical of this experimental setup. There key thing to notice here is that that we were able to obtain fairly precise estimates of topographic parameters like the extent of intrinsic connections c despite the fact that we are only using data from a single electrode.

Conclusions

This paper reviews some recent results regarding the modelling of spatially extended neuronal dynamics using neural field models. We have focused on fitting these models to empirical data and entertained questions about cortical structure and function. Addressing these questions relies on exploiting neuroimaging data to invert neural field DCMs. Optimizing neural field models can be a hard task due to their intrinsic nonlinearities. We have here shown that using dynamic causal modelling it is possible to validate neural field models in relation to empirical data. Our approach comprises three important components: first, a careful selection of priors that conform to known neurobiology; where parameters about which we have little prior knowledge are given flat priors – and the parameters of the observation model are optimized concurrently. Second, an appropriate cost function: here a variational free energy bound on log model evidence. Third, we call on linearity and ergodicity assumptions that allow one to summarize cortical activity efficiently, in terms of frequency domain responses.

In our approach, neural fields serve as generative models that allow us to define a probabilistic mapping from free parameters to observed cross spectra. This assumes that the measured signal is a mixture of predicted spectra, channel and observation noise and can furnish predictions for conventional measures of linear systems; like coherence, phase delay or cross correlation functions, as detailed in [25]. In brief, there is a mapping between model parameters (effective connectivity) and spectral characterizations (functional connectivity) that provides a useful link between the generative modelling of biophysical time series and dynamical systems theory. Neural field models prescribe a likelihood function; this function taken together with the priors specify a dynamic causal model that can be inverted using standard variational procedures [3]. This Variational Laplace scheme approximates model evidence with a variational free energy. The (approximate) posterior density and (approximate) log-evidence are used for inference on parameters and models respectively. In other words, one can compare different models (e.g., neural field and mass models) using their log-evidence and also make inferences on parameters, under the model selected.

In brief, our approach is based on a combination of neural fields with appropriate observation models that are optimised in relation to observed data and are – crucially – compared in terms of their evidence. This provides a principled way to adjudicate among different models or hypotheses about functional brain architectures and the physiological correlates of neuronal computations.

Software note

The procedures described in this review can be accessed as part of the SPM academic freeware (in the DCM toolbox: http://www.fil.ion.ucl.ac.uk/spm/). This code has been written in a modular way that allows people to select from a suite of neural mass and field models to analyse their data – or indeed analyse simulated data that can be generated by the routines. It has also been written in a way that allows people to specify their own models – and connectivity architectures – in terms of Matlab routines (using an equation of motion and a nonlinear mapping from hidden states to observed measurements).