PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of ploscompComputational BiologyView this ArticleSubmit to PLoSGet E-mail AlertsContact UsPublic Library of Science (PLoS)
 
PLoS Comput Biol. 2010 September; 6(9): e1000919.
Published online 2010 September 9. doi:  10.1371/journal.pcbi.1000919
PMCID: PMC2936513

A Dynamic Neural Field Model of Mesoscopic Cortical Activity Captured with Voltage-Sensitive Dye Imaging

Valentin Markounikau, 1 , 2 , * Christian Igel, 1 , 2 Amiram Grinvald, 3 and Dirk Jancke 1 , 2 , 3 , 4
Olaf Sporns, Editor

Abstract

A neural field model is presented that captures the essential non-linear characteristics of activity dynamics across several millimeters of visual cortex in response to local flashed and moving stimuli. We account for physiological data obtained by voltage-sensitive dye (VSD) imaging which reports mesoscopic population activity at high spatio-temporal resolution. Stimulation included a single flashed square, a single flashed bar, the line-motion paradigm – for which psychophysical studies showed that flashing a square briefly before a bar produces sensation of illusory motion within the bar – and moving squares controls. We consider a two-layer neural field (NF) model describing an excitatory and an inhibitory layer of neurons as a coupled system of non-linear integro-differential equations. Under the assumption that the aggregated activity of both layers is reflected by VSD imaging, our phenomenological model quantitatively accounts for the observed spatio-temporal activity patterns. Moreover, the model generalizes to novel similar stimuli as it matches activity evoked by moving squares of different speeds. Our results indicate that feedback from higher brain areas is not required to produce motion patterns in the case of the illusory line-motion paradigm. Physiological interpretation of the model suggests that a considerable fraction of the VSD signal may be due to inhibitory activity, supporting the notion that balanced intra-layer cortical interactions between inhibitory and excitatory populations play a major role in shaping dynamic stimulus representations in the early visual cortex.

Author Summary

Understanding the functioning of the primary visual cortex requires characterization of the non-linear dynamics that underlie visual perception and of how the cortical architecture gives rise to these dynamics. Recent advances in real-time voltage-sensitive dye (VSD) imaging permit recording of cortical population activity with high spatial and temporal resolution. This wealth of data can be related to cortical function, dynamics, and architecture by computational modeling. Here we used a mesoscopic neural field model to describe brain dynamics at the population level as measured by VSD imaging. Introduced in 1972 by Wilson and Cowan, these models are derived from statistical mechanics to analyze the collective properties of large numbers of neurons. For simplicity, the cortical planar tissue is assumed to contain only two types of homogeneously distributed neurons (excitatory and inhibitory) that interact via recurrent lateral connections. This study shows 1) how a concise neural field model can simulate VSD data quantitatively in space and time by identifying the underlying non-linear dynamics, 2) how such a model can support hypotheses about visual information processing, and 3) how the model can be linked to the neuronal architecture.

Introduction

Visual cortical activity does not exclusively mirror visual input but rather reflects the contribution of additional recurrent processes involving lateral and local feedback couplings. Understanding cortical processing requires a theoretical understanding of the underlying activity dynamics, which can be attained by modeling at various levels of abstraction. Naturally, the chosen level should match the level at which neuronal recordings are made [1]. The activity patterns observed using voltage-sensitive dye (VSD) imaging reflect population activity at the mesoscopic (intra-areal) level [2], [3]. This suggests the application of mean-field models in which large numbers of neurons are averaged. Moreover, we are interested in the relation of neuronal dynamics to the spatial dimensions of the cortical sheet (and more generally to metric embeddings spanned by more abstract parameters, see [4], [5]). Neural field (NF) models [1], [6][11], in which the efficacy of synaptic coupling depends on the notion of distance between neurons or ensembles of neurons, are therefore our preferred choice. Here, we show that a minimalistic multiple-layer NF models can simulate mean VSD data in space and time with high accuracy. The model is an abstract functional description of VSD-recorded dynamics. Thus, it is in the first place phenomenological. However, its interpretation in biological terms allows to link its structure and parameters to the neuronal functional architecture.

The imaging data that we model showed: i) Two stationary stimuli (a square and an elongated bar) presented in rapid succession produce a pattern that signals propagation of activity across the bar's retinotopic representation in early visual cortex. ii) The obtained pattern was different from activity when the bar was flashed alone, and did not match the simple superposition of activities evoked by individually-presented square and bar stimuli. iii) Rather, we observed propagation of a wave front of activity that was also found when a square stimulus moved physically in visual space [12].

Based on the VSD imaging data [12], we hypothesized that a two-layer neural field [7][9] model can account for the findings i–iii. If so, this would imply that the feedback from higher brain areas is not a principal requirement to produce motion patterns across primary visual cortex upon presentation of a square and a bar flashed in rapid succession (as debated, e.g., in [12][17]).

Voltage sensitive dye imaging measures relative fluorescence changes that are linearly correlated to changes in membrane potentials [18],[19]). This technique currently allows recording of in vivo cortical activity at sub- as well as suprathreshold level with at least 10 ms temporal resolution and a spatial resolution of 50 An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e001.jpgm across several millimeters of cortex. Hence, it is well suited to capture the real-time dynamics of millions of neurons at once. However, the signal does not distinguish between excitatory or inhibitory contributions to the overall activity. Therefore, our model explicitly assumes that the VSD signal reflects a mixture of activity from excitatory and inhibitory neuronal populations, and thus contrasts with studies that interpret VSD data as mainly reflecting excitatory activity (e.g., [20][22]). We expect to gain insights about the relative contributions of excitatory and inhibitory activities in the VSD signal, because the model allows separate inspection of its inhibitory and excitatory layers.

In the following, we first describe the underlying data and its preprocessing, our model structure, and our parameter identification procedure. Then we present the results including the model fit in comparison to further simplified models, the model prediction regarding similar yet novel stimuli, and the results of a standard linear stability analysis of the homogeneous solution of the model. We then follow with a discussion of the findings in relation to alternative modeling approaches, the physiological interpretation of our model, and the role of excitation and inhibition in the model. Finally, we consider our results in the context of hypotheses concerning the cortical representation of motion and the origin of the line-motion illusion.

Methods

Data

The data underlying this study were recorded by Jancke et al. at the Department of Neurobiology at the Weizmann Institute of Science in Israel using VSD imaging of cat visual cortex [12], [18], [23]. Animals were initially anesthetized with a mixture of ketamine (15 mg An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e002.jpg) and xylazine (1mg An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e003.jpg). After tracheotomy, animals were respirated and anesthetized with 1.5% halothane (0.8% during recordings) in a 1[ratio]1 mixture of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e004.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e005.jpgO. The animals were paralyzed with pancuronium bromide (0.2 mg An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e006.jpg, intravenously). Area 18 was stained for 2.5–3 h with the voltage-sensitive dye RH-1691. Its molecules bind to the external surface of excitable membranes and transform changes in membrane potential into changes in fluorescence intensity, which is correlated linearly with membrane potentials of layer 2 and 3 cortical neurons [19], [24]. Using a high-speed camera, the VSD signals were recorded with a temporal resolution of 9.6 ms. Stimuli were presented binocularly. The projection of the area centralis to the monitor screen was determined using a fundus camera. If necessary, projection of the eyes was converged using a prism in front of the eye that was ipsilateral to the recorded hemisphere. To control for possible eye drift during the experiment, the position of the area centralis and receptive field positions were measured repeatedly. The stationary stimuli were a square of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e007.jpg and a bar of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e008.jpg. Both stimuli were aligned with their upper edges. The line-motion (LM) stimulus consisted of a square briefly presented before the bar. Additionally, squares were moved vertically with 4, 8, 16, and 32 deg/s. Stimuli were placed at An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e009.jpg eccentricity. In the LM setting the square was presented for 50 ms followed by an inter-stimulus interval of 10 ms, after which the bar was presented for 130 ms. In single presentations the square and the bar were displayed with the same individual timings as in the LM condition.

The imaged signals (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e010.jpg) reflect relative changes in fluorescence compared to baseline: Imaging data were normalized by its DC level during pre-stimulus period (200 ms) for each pixel; heart-beat and respiration-related artifacts were removed by dividing by the average of “blank” signals recorded in absence of stimulation. For more details about the imaging methods we refer to [12].

We restrict our analysis to the first 250 ms of the VSD recordings during which propagation of activity was observed. The VSD data are sequences of frames that image 3 mmAn external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e011.jpg7 mm of cortex. This area was discretized into 24An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e012.jpg50 pixels (px). We denote the fluorescence change represented by a pixel at position An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e013.jpg at time An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e014.jpg by An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e015.jpg. As the stimulus representations in this study vary only along the posterior-anterior cortical axis we reduced the dimensionality of the data by averaging An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e016.jpg along the medial-lateral cortical axis (as done, e.g., in [25]), see Figure 1. For each vertical position An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e017.jpg and time step An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e018.jpg, we define An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e019.jpg, where An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e020.jpg (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e021.jpg) and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e022.jpg (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e023.jpg). To eliminate the low signal-to-noise at the border of the representations, activity was averaged across the central pixels 4–20. These values were defined as mean An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e024.jpg standard deviation of a Gaussian function fitted to the distribution of activity values across x-axis in the flashed bar condition (see Figure S1 in Text S1). However, the exact choice of these positions did not affect the results. Given this dimension reduction, the variable An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e025.jpg is from hereon referred to as An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e026.jpg for clarity. Finally, the VSD data were additionally normalized using the mean level of activity when no stimulus was presented as reference. That is, for each stimulus condition the activity was averaged over the first 20 ms and all An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e027.jpg and the result was subtracted from the data.

Figure 1
Deriving a space-time diagram of the VSD signals.

Model

We aimed for a model that quantitatively captures the spatio-temporal cortical dynamics observed by VSD imaging in response to the stimuli described above. The model should have as few parameters as possible, and these parameters should allow for functional interpretations (e.g., in terms of lateral interactions and time constants). The model should be at the same level of abstraction as the neuronal data. As our data reflect the dynamics of neuronal populations and describe the spread of activity across the cortical sheet, NF models are an appropriate choice [1], [4], [5], [21], [26].

The model is one-dimensional as the dynamics of interest evolve in one dimension along the (apparent) movement direction as described above. Our two-layer NF is governed by the following system of integro-differential equations:

equation image
(1)

equation image
(2)

where An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e033.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e034.jpg denote the mean membrane potentials of model neurons at cortical position An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e035.jpg and time An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e036.jpg in the excitatory and inhibitory layer, respectively, and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e037.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e038.jpg are the corresponding time constants. The resting potentials are determined by An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e039.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e040.jpg. The transfer functions (or axonal response functions) An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e041.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e042.jpg are represented by sigmoidal functions of the membrane potential, which relates depolarization to firing rate. For instance, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e043.jpg is given by

equation image
(3)

The parameter An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e045.jpg describes the slope of the response function and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e046.jpg its threshold value. The transfer function An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e047.jpg is parameterized analogously. The efficacy of the synaptic connectivity between two positions An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e048.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e049.jpg is assumed to be translation invariant and isotropic. It is given by gain factors An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e050.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e051.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e052.jpg and synaptic strength (weight) functions An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e053.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e054.jpg of Gaussian shape. For example, the coupling An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e055.jpg from the excitatory to the inhibitory layer is modeled by

equation image
(4)

The parameter An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e057.jpg controls the width and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e058.jpg controls the strength of the interactions. Eccentricity in the imaged region of area 18 (see legend Figure 1) is approximately linear along the vertical meridian [27], [28] mapping one degree of the visual field to An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e059.jpg mm of cortex. The external visual stimulus mapped accordingly to cortical coordinates is denoted by An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e060.jpg, and the afferent input is computed by smoothing the stimulus by convolution with a Gaussian function An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e061.jpg. The afferent input couples only into the excitatory layer, and the inhibition is coupled into the excitatory layer locally.

There are several reasons for the local coupling from the inhibitory to the excitatory layer. From a technical point of view, we want our model to have as few parameters as possible. Equations 1 and 2 can be viewed as the minimal extension of the original Amari neural field model [9] to distinct dynamics for excitatory and inhibitory neurons (see also Eqn 10), allowing for different time constants and making it consistent with Dale's law (we added the response function An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e062.jpg such that all communication between neurons can be interpreted as being encoded by mean firing rates). The resulting model is similar to the dynamics considered in [4], and it can be argued that the local coupling resembles physiology because the majority of interneurons projects locally [29].

The presented NF model is parameterized by 15 values. However, because all model neurons communicate only through the transfer function An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e063.jpg (see Eqns 1 and 2), shifting of resting potentials An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e064.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e065.jpg can be compensated by changing the transfer function thresholds An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e066.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e067.jpg in Eqn 3 such that the dynamics are not affected. Thus the parameters An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e068.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e069.jpg are redundant.

The VSD data reflect the activity of both excitatory and inhibitory populations of neurons. Accordingly, we define the modeled VSD signal An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e070.jpg to be a combination of the field's excitatory and inhibitory activities. We assume an affine linear mixture

equation image
(5)

with the non-negative coefficients An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e072.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e073.jpg controlling how strongly excitatory and inhibitory activity is reflected in the dye signal, and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e074.jpg is some offset [30], [31]. We do not explicitly consider units of measurement to keep the notation uncluttered. Formally, the dye signal An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e075.jpg is measured in change of fluorescence intensity An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e076.jpg and the potentials in mV. Therefore, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e077.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e078.jpg have units An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e079.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e080.jpg is measured in An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e081.jpg.

The time delay between the stimulus presentation and the response onset in the VSD data is a sum of the time the neuronal signal needs to travel from the retina to the primary visual cortex and the time the neuronal populations in the primary visual cortex need to build up a detectable activity. A fixed retino-cortical time delay of two time frames (19.2 ms) was used in our model to align the response onsets of the model and the VSD data.

In our numerical experiments, we iterated the dynamical systems defined by Eqns 1 and 2 starting from the initial conditions An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e082.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e083.jpg for all An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e084.jpg. Then we let the system relax for some time period in the absence of afferent input. When we present our results, the time An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e085.jpg is after this relaxation phase. In our model, we had to discretize the spatial dimension. We simulated 150 spatial positions. The center 50 positions were mapped to the 50 pixels of the VSD image data. The other positions were added to avoid boundary effects and are not shown in the results.

Model parameters identification procedure

In order to explore the system parameter space we adopted a grid search procedure on a reduced parameter set. The reduced set contained 10 parameters An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e086.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e087.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e088.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e089.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e090.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e091.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e092.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e093.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e094.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e095.jpg, where the excitatory and inhibitory resting potentials were assumed to be equal (i.e., An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e096.jpg). We defined a grid with 3 points (chosen by educated guess, see Table 1, third column) for each parameter, which results in total to An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e097.jpg parameter configurations. The remaining parameters An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e098.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e099.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e100.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e101.jpg were fixed to physiologically plausible values. The transfer function thresholds An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e102.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e103.jpg were both set to −40 mV (see discussion of parameter redundancy in the previous section). The feed-forward gain An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e104.jpg was set to 70 to yield maximal amplitude of the input signal of about 60 mV. According to the literature, one retinal position is represented in our region of interest by a population of cortical neurons distributed around the mean position with a standard deviation of approximately 0.6 mm [28]. The parameter An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e105.jpg, determining the feed-forward smoothing of the input (as parameter of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e106.jpg, see Eqns 1 and 4), was therefore fixed to 4 pixels in the discretized model, which corresponds to a value of 0.51 mm.

Table 1
Summary of the model parameters.

After we have simulated our NF model with a given set of model parameters and obtained the spatio-temporal patterns of the excitatory and inhibitory layer in response to the stimuli used for system identification, we can compute the values An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e123.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e124.jpg, and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e125.jpg, see Eqn 5, using the ordinary least-squares (OLS) method under the constraints An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e126.jpg. That is, we solve

equation image
(6)

where (with a slight abuse of notation) An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e128.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e129.jpg, and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e130.jpg refer to the concatenated signals of all stimulus configurations considered in the optimization procedure. This yields the optimal values for An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e131.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e132.jpg, and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e133.jpg in terms of the mean squared error between aggregated model signal An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e134.jpg and observed dye patterns An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e135.jpg for the given spatio-temporal patterns of the excitatory An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e136.jpg and inhibitory An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e137.jpg layer. From the solution, we get the aggregated signal

equation image
(7)

and obtain the mixture coefficient

equation image
(8)

of excitation and inhibition. The An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e140.jpg coefficient indicates that the simulated signal is comprised of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e141.jpg excitation and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e142.jpg inhibition (under the assumption that the values of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e143.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e144.jpg vary in the same interval). For optimization of the system, only four of the seven available stimulus configurations were used, namely the flashed square, the flashed bar, the LM stimulus, and a square moving at 32 deg/s. The other three stimulus conditions, squares moving at different speeds, were used to test whether the model generalizes to unseen, but related stimuli.

We selected parameter configurations from the sets produced by the grid search procedure according to the following criteria: (1) the system with given parameters is stable (see section “Analysis of Stability”); (2) after stimulus presentation, the simulated activity eventually decays to zero; (3) the correlation coefficient between the dye and the simulated data is larger than 0.8; (4) the rate of activity change during the onset period (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e145.jpg20–70 ms for the flashed square, LM and moving square stimuli; An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e146.jpg80–130 ms for the flashed bar stimulus) is similar in the modeled and measured responses.

In order to show that the model fit can be considerably improved by fine-tuning, we adjusted some model parameters using a randomized direct optimization algorithm. We used the selected parameter set from the grid search procedure (see Table 1) as starting point. We defined an objective function that quantifies the correlation coefficients between dye and simulated data as well as inhibition/excitation ratio:

equation image
(9)

The correlation coefficients An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e148.jpg were computed between the dye and the simulated data for each stimulus configuration considered in the optimization and the trade-off parameter is set to An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e149.jpg.

This objective function is non-convex and multi-modal (i.e., there are undesired local optima). Our method of choice for such problems is the covariance matrix adaptation evolution strategy (CMA-ES) [32], [33]. The CMA-ES is an iterative, direct, stochastic optimization method and one of the most efficient biologically inspired search heuristics for real-valued optimization [34]. It has been successfully used to adapt neural field models [35][37] and is explained in detail in Text S1, section A.

Results

VSD imaging of cortical responses to the line-motion paradigm

The imaged visual cortical area included the complete retinotopic representation of the applied stimuli. Figure 1 shows the observed evolution of activity in response to the LM paradigm within single time frames as originally recorded (upper row). After the square was presented (frame zero), activity emerged and reached suprathreshold amplitudes around the thalamic input location (see region colored yellow/red in frames 48–58 ms; verified through spike recordings at various electrode positions [12]). After 60 ms the bar was presented, giving rise to activity that was gradually drawn-out (68–116 ms) along the retinotopic bar representation in the anterior cortical direction (see upper arrows and legend for more details). In contrast, the lower part of the responses showed no propagation (see lower row of arrows). The anteriorly propagating activity was interpreted as a motion signal that correlates to the perceived illusory line-motion as found in psychophysical studies [13]. Thus, instead of representing the bar at once, activity propagated across the retinotopic map. As shown in the original data, the same characteristics were obtained in multiple other experiments [12]. Since this effect occurred along the posterior-anterior cortical y-axis, we averaged each camera frame along the x-axis (see Figure 1) to enable a one-dimensional model approach (see 2nd row of Figure 1). Using these 1D frames, we depict activity in space-time diagrams that allow inspection of entire time courses (3rd row).

Aggregated activity dynamics of both excitatory and inhibitory layers – emergence of propagating activity

The NF parameters found by grid search are summarized in Table 1. Figure 2 shows both the data and the NF responses. Activity in the model started with a delay of 19.2 ms (two time frames, see section “Model parameters identification procedure”), following low amplitude activity (light blue-green colors) that spread rapidly across several millimeters of cortex. With increasing amplitudes (yellow, red), the speed of the spread gradually decreased.

Figure 2
Visual inputs, modeled, and VSD responses.

For the flashed square (Figure 2B), activity spread symmetrically around the stimulus input location. With either the flashed bar (Figure 2C) or LM stimulus (Figure 2D), nearly identical responses were found in the lower part (0–3 mm) of the space-time diagram, as there were no differences in stimulus inputs in this direction. In the upper part, however, as outlined by rectangles, there were crucial differences between these conditions around the elongated bar input location. In this region, we observed a gradual propagation of high-amplitude activity – from the lower left to the upper right for the LM stimulus but not for the single flashed bar (Figure 2C, a more detailed analysis is provided in Figure S2 in Text S1). Note that such propagation of high-amplitude activity was also observed both in the VSD and modeled responses to a moving square (Figure 2A). Thus, similarly to the LM stimulus, real moving input gave rise to a propagating wave front of cortical activity reporting physical motion across the retinotopic map (see Figure S2 in Text S1). In the model the LM effect results from the following mechanism: The activity in response to the flashed square is sustained after the stimulus offset and serves as a starting point for the response to the flashed bar. The bar-evoked excitation then propagates from this region producing a gradual spread of activity, drawn-out from the highest activity amplitudes. Hence, the LM effect persists in the model as long as the inter-stimulus interval is in the range of activity decay times after the square stimulus presentation.

Figure 3 shows the evolution of activity at the center of the cortex image (the model responses are given by dashed lines). Although the response to the square stimulus was sustained longer than the measured dye signal and the maximum amplitude was not fully reached in the NF, the model in general captured the time courses of the cortical responses to the different stimuli used. Moreover, a small adjustment of the parameters removed these discrepancies: Using evolutionary optimization (see “Methods” and Text S1, section A) the duration of the model response to the flashed square was greatly reduced (Figure 3, red lines), and the correlation coefficient for this single configuration increased to 0.81. The evolutionary optimization changed the parameters to An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e151.jpg mV, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e152.jpg mV, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e153.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e154.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e155.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e156.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e157.jpg (see Text S1, Table S1).

Figure 3
Time courses of aggregated activity of a single model neuron.

Nonlinear space-time interactions

We next analyzed whether nonlinear interactions are required to explain our findings to explain the space-time responses of the LM condition. The superposition of the responses to the single square and the bar alone differed from the response to their combined presentation in the LM stimulus (Figure 4). Note that the superposition also gave rise to propagation of activity as can be seen by the temporal offsets between the time courses at different locations (compare stippled red and green curves, Figure 4). However, there were marked deviations from linearity. The VSD LM response (Figure 4A) exhibited facilitation (around An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e159.jpg ms) and suppressive effects (around An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e160.jpg ms) as compared to the superposition of the single square and the bar alone. The significance of these effects was tested as follows. We defined three consecutive time intervals. The pre-stimulus interval, frames 1–4, the “facilitatory” interval, frames 5–12, and the “suppressive” interval, frames 13–20. First, we found no significant difference between all pixel values of the LM and the superposition in the pre-stimulus interval (Mann-Whitney U = 19808, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e161.jpg = An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e162.jpg = 200, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e163.jpg, two-tailed). Second, we showed that the difference between the LM and the superposition is significant in the “facilitatory” interval (Mann-Whitney U = 64195, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e164.jpg = An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e165.jpg = 400, PAn external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e166.jpg, two-tailed) as well as in the “suppressive” interval (Mann-Whitney U = 64369, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e167.jpg = An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e168.jpg = 400, PAn external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e169.jpg, two-tailed). These effects are emphasized by thin red lines showing the difference of the LM response and the superposition (positive values indicate facilitation, negative values suppression). The model (Figure 4B) showed both of these effects, although the facilitatory effect (blue ellipse) was considerably smaller than in the measured data and was not observed at the center of the pattern (green line). We conclude that non-linear dynamics are necessary to explain the deviation from a simple superposition of the single square and bar representations.

Figure 4
Comparison of the LM condition to the superposition of responses to the square and the bar alone.

Comparison to reduced models

In order to verify that the different components of our model are necessary to reproduce the VSD-recorded dynamics, we tested several simplified models. We asked whether the spatio-temporal patterns could be simulated by a simple feed-forward model without lateral interactions, in which Gaussian smoothing of the input stimulus and low pass filtering produce the spread of activity (similar to the “G-waves” of Baloch and Grossberg [15]). It turned out that switching off lateral couplings and optimizing the parameters of the reduced model did not give satisfactory results.

Next, to investigate the role of inhibition and the necessity of different time constants for inhibition and excitation, we considered a single layer Amari-type NF [9],

equation image
(10)

where the kernel function An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e178.jpg was either purely excitatory (Gaussian) or had a Mexican-hat shape. The latter models excitation and inhibition as the difference between two Gaussian kernels. Applying the grid search procedure to the reduced Mexican-hat model produced fairly good fits to the data with correlation values of up to 0.75 (compared to 0.85 for the full two-layer model), the individual correlation coefficients computed between the simulated and measured responses to particular stimulus conditions were: 0.55 for the square moving at 32 deg/s, 0.51 for the flashed square, 0.81 for the flashed bar, and 0.86 for the LM stimulus condition. The amplitudes of activity in response to smaller sized stimuli (i.e., flashed and moving squares) were too low. Testing this model on the moving stimuli produced responses that were too weak and too prolonged. The overall correlation coefficient between simulated and measured data for the four conditions with moving squares was 0.70. The individual correlation coefficients computed between the simulated and measured responses to particular stimulus conditions were 0.78 for the square moving at 4 deg/s, 0.84 for 8 deg/s, 0.69 for 16 deg/s, and 0.55 for 32 deg/s. As indicators for the goodness-of-fit, we computed the AIC (Akaike information criterion [38]) for the two-layer and the single layer NF. The AIC is given by An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e179.jpg, where An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e180.jpg denotes the residual sum of squares, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e181.jpg the number of data points (4 stimulus configurations, 50 spatial positions, and 25 time steps), and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e182.jpg the number of parameters [38]. Note that the absolute values of the AIC cannot be interpreted. Although the two-layer model has more parameters, its AIC is smaller (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e183.jpg compared to An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e184.jpg) indicating a better fit. These results were obtained using the Mexican-hat kernel, the results for the simple Gaussian kernel were even worse.

In addition to the coarse grid search, we performed global optimization of the Amari-type field model using the CMA-ES (see Text S1, section A) without coming close to the fit quality of our two-layer model. Clearly, the fact that the system identification procedure did not find suitable parameters for a different model class does not prove that no suitable parameters for these alternative models exist (as such, non-existence proofs can be regarded as a general challenge in computational neuroscience). However, we regard the failure to fit the reduced models as a strong indication that lateral “intracortical” interactions and inhibition evolving with independent time constants are indeed necessary for the best fit of our data.

Decomposing the model aggregated activity: Separate inspection of layers reveals their high correlation

In Figure 5 the decomposed excitatory and inhibitory NF responses showed the observed gradual spread in both components.

Figure 5
Decomposing NF activity into excitatory and inhibitory components.

The parameter An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e187.jpg, which reflects the mixing ratio of excitatory to inhibitory signals, was optimal for a value of 0.54 as revealed by our fitting procedure. This suggests that the aggregated signal was caused to a large extent by inhibitory processes. If we used for this parameter configuration either only An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e188.jpg (or only An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e189.jpg) in the regression Eqn 6 (i.e., set An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e190.jpg or An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e191.jpg, respectively), the quality of the fit significantly declined compared to the use of both layers' aggregated activities An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e192.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e193.jpg (verified by Wald-test [39], [40]. This indicated that elimination of the variable An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e194.jpg had a significant impact on the goodness-of-fit, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e195.jpg). We computed the residuals An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e196.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e197.jpg, see section “Model identification procedure”. Let An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e198.jpg be the number of data points used for fitting the models (4 stimulus configurations, 50 spatial positions, and 25 time steps). We compared An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e199.jpg to the An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e200.jpg distribution [39], [40]. When interpreting this result one has to keep in mind that the residuals of our models are correlated because of the spatio-temporal structure of the data and their variances are not homogeneous.

However, the choice of the mixing ratio An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e201.jpg could vary over a wide range without a considerable decrease in fit quality as indicated by the high correlation coefficients for different mixing ratios (see Figure 6). This result was on the one hand surprising as we expected to reveal a critical value for this parameter. Consequently, we cannot estimate directly the relative fraction of excitatory and inhibitory signals in the VSD data. On the other hand, our results strongly suggest that both inhibitory and excitatory activity contribute significantly to the measured VSD responses as demonstrated by our statistical goodness-of-fit analysis. Moreover, our control experiments with simpler single layer models (see “Comparison to reduced models”) indicated that balanced interactions of the strongly coupled excitatory and inhibitory layer are required to produce the high correlation between the model aggregated spatio-temporal patterns and the VSD dynamics.

Figure 6
Different mixing ratios of excitatory and inhibitory activity.

Generalization to moving squares with different speeds

We finally tested the derived model using different stimuli that were similar to those used in the grid search procedure. Novel square stimuli moving at different speeds (4, 8, 16 deg/sec; the 32 deg/s stimulus was used for optimization) were fed into the model. The results, summarized in Figure 7, revealed that our model also reproduced these new spatio-temporal patterns. The correlation coefficient between simulated and measured data for the four moving square conditions was on average 0.79.

Figure 7
Model and VSD responses to moving squares with different speeds.

Analysis of stability

For a better understanding of the dynamics of our model, we performed a standard linear stability analysis in the absence of afferent input (e.g., see [41]). We linearized the system governed by Eqns 1 and 2 near its homogeneous solution, which is described by An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e209.jpg with scalars An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e210.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e211.jpg.

The spatial dimension is considered in frequency space. The analysis yields two constraints that have to be met for the linearized system to be asymptotically stable:

equation image
(11)

equation image
(12)

for all spatial frequencies An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e214.jpg, where An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e215.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e216.jpg are the spatial Fourier transforms of the weight functions An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e217.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e218.jpg (e.g., An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e219.jpg) and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e220.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e221.jpg are the derivatives of the transfer functions. Substituting the parameters listed in Table 1 into the Inequalities 11 and 12 reveals that the homogeneous solution of our model is indeed asymptotically stable, which is a minimal requirement (for detailed derivation of Inequalities 11 and 12 see Text S1, section B).

Discussion

This is the first time in vivo VSD patterns have been modeled quantitatively in space and time on a mesoscopic population level. Our proposed dynamical system is an abstract functional description of in vivo recorded VSD dynamics that correlate with changes in potentials across neuronal membranes and reflect the mass activity of a large pool of neurons.

Despite its rather few degrees of freedom, the NF model matches the VSD activity patterns evoked by briefly flashed visual stimuli including the LM paradigm. Our statistical analysis indicate that the two-layer structure of our model as well as the assumption that both excitatory and inhibitory activity contribute to the dye signal are necessary to obtain a fit with the highest correlation to the VSD imaging data.

The question arises whether our model is realistic in the sense that it could be implemented by the cortical network. We argue that the relevant model parameters are indeed within physiologically plausible ranges (see “Physiological interpretation of the model”). In this context, we suggest a significant contribution of inhibitory activity to the VSD dynamics.

Finally, the here presented modeling of V1 implies that feedback from higher brain areas is not necessary to produce activity patterns resembling the percept of illusory motion.

Relation to alternative large-scale model

First, we compare our model to other in silico simulations that addressed the same experimental findings. Rangan et al. [20] simulated the VSD data from Jancke et al. [12] using a large-scale integrate-and-fire model, which consists of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e222.jpg neurons modeling three types of channel conductances (NMDA, GABA, AMPA) and two types of connections (isotropic short-range and orientation-specific long-range connections) [42]. Importantly, both their model and ours simulate the LM effect without additional modeling of higher brain areas [20]. While our NF model was specifically designed to capture the LM data, the large-scale network developed by Rangan et al. [20] incorporated additional aspects of cortical processing on a different level of abstraction, as their large-scale model accounts for further experimental observations (for instance correlations between spontaneous activity and stimulus orientation [42], [43]) that cannot be addressed by our model, which is only as complex as needed to match our data. Therefore, it is not surprising that our model allowed a better fit per se.

Comparing the simulations in Rangan et al. [20] with our study, it becomes apparent that both approaches have a spatial and temporal integration problem: activity in response to the square stimulus is much less integrated over space and time than for the longer bar stimulus. The large-scale model [20] therefore included additional pre-processing in the lateral geniculate nucleus (LGN). Using such an LGN model that implements normalization, as proposed by several authors [25], [30], [44][47], may indeed help to enhance our fits but was not explicitly tested here.

In the large-scale model, the response following the bar arises about 20 ms earlier than the VSD-recorded responses. In contrast, the response to the square stimulus was An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e223.jpg ms delayed. Using the constant retino-cortical time delay in our model, the simulated and measured signal onsets were aligned in all stimulus conditions. Our model thus has the advantage of capturing the timing of VSD signal onsets more accurately (see Figure 3).

Possible extensions of our NF model

The model responses to both the flashed bar and LM condition fitted the observed VSD measurements. In contrast, for the single flashed and the moving squares the model revealed a discrepancy to the VSD data in the extent of lateral spread. One reason for this effect is our simple Gaussian smoothing that we used as a model for the retino-cortical processing. Increasing the kernel width An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e224.jpg resulted indeed in a wider activity spread, but the tested widths were inappropriately large to match the common experimental findings. As another straightforward solution, we increased the widths of the coupling kernels An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e225.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e226.jpg, however, the grid search did not find models with an accurate fit for such wider kernels.

For the flashed square, prolonged activity was observed compared to the data. Importantly, tuning the gains An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e227.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e228.jpg, and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e229.jpg, the resting potentials An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e230.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e231.jpg, and the steepness of the transfer function An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e232.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e233.jpg using evolutionary optimization eliminated the discrepancy (see Figure 3).

Finally, the flashed and moving square stimuli evoked model responses that were lower in amplitude than that measured. As stated in “Relation to alternative large-scale model”, using a normalization method [25], [30], [44][47] in the retino-thalamic processing step could be a suitable solution to adjust the amplitudes of the responses to the small- and large-sized stimuli. As we aimed at capturing essential nonlinear activity dynamics with the simplest form of our NF model, we also neglected further specific mechanisms of cortical processing like short-term synaptic plasticity and did not parameterize axonal propagation speed (see, e.g., [22], [31]). It should be noted that increasing the model complexity by introducing more parameters would certainly enhance the fit to the data. For instance, in a comparable version of the model, to account for spiking population data [48], we added a “shunting inhibition” term that allowed to produce a rapid spread of activity ahead of a moving square without pushing the model into the active mode [12]. However, the ability of the present model with its minimal complexity to fit and generalize suggests that the principles of two-layered architecture with lateral connectivity are sufficient to account for VSD-recorded dynamics.

Physiological interpretation of the model

Our model is functional in the sense that it characterizes the dynamics of the VSD signal. Still, the model must also be plausible from a physiological point of view in the sense that it can be implemented by the underlying brain structures. In the following, we therefore discuss the ability of the cortical architecture to give rise to the model dynamics.

The NF model is a graded response mean-field approximation. As shown in the work of Eggert and van Hemmen [49], graded response models are able to describe the evolution of a population of spiking neurons in the case of slow dynamics. The VSD data used in this study lack these very fast ms-dynamics due to the frame duration and averaging and therefore, fulfills this requirement. The model accounts for the neuronal activity measured in cortical layers 2/3 and ignores other cortical structures. Furthermore, as it matches averages over several trials and numerous cell types with different connectivities and different neuronal response profiles, the neglected details affect the obtained model parameters. In this sense the model dynamics implicitly simulate all the factors that influence VSD-recorded dynamics.

The parameters An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e234.jpg mV and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e235.jpg mV in the model can be associated with the mean resting membrane potentials of excitatory and inhibitory neurons, respectively. The physiological values reported in the literature span from An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e236.jpg mV to An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e237.jpg mV with a standard deviation of about 10 mV [50][54] for excitatory and inhibitory neurons in the cat primary visual cortex. Our modeled membrane potentials are in this range. However, it should be noted that since all model neurons communicate only through the transfer function An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e238.jpg (see Eqns 1 and 2), shifting of resting potentials An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e239.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e240.jpg can be compensated by simply changing the transfer function thresholds An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e241.jpg such that the system dynamics are not affected, see section “Model”. In the model, the mean membrane time constants An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e242.jpg ms and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e243.jpg ms are in agreement with experimental data [51], [55].

The lateral connectivity is determined in the model by the Gaussian interaction kernels An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e244.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e245.jpg (see Eqn 4). As shown in Table 1, the values An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e246.jpg mm and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e247.jpg mm were used. The absolute values of the axonal extension of excitatory cells found in the cat primary visual cortex have been reported to reach up to 3.5 mm [56] or even up to 6–8 mm [57]. However, if our An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e248.jpg-values are interpreted as measures of lateral extent of axonal-dendritic connections, they can be compared with the results of a recent quantitative study that measured and modeled the spatial and orientation preference distribution of labeled axonal bouton density of excitatory neurons in area 18 of the cat [58]. Modeling of the spatial distribution by a single Gaussian function revealed a An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e249.jpg-value of 0.6 mm, which is smaller than in our NF. However, their additional approach designed to differentiate between oriented and non-oriented components suggested values up to 1.1 mm for the non-oriented component, which is close to our results.

In contrast, inhibitory interneurons in the primary visual cortex act mostly locally [59], [60]. These interneurons are activated by widely spreading lateral excitatory connections (see [61] and [62], [63] for results for macaque and cat, respectively). This has been implemented in the NF by the local coupling of the inhibitory term into the excitatory field equation (see Eqns 1 and 2). Thus, although the coupling from the inhibitory to the excitatory layer is only local in our model, inhibition can act over a wide range.

Excitation and inhibition in the VSD signal

We considered a mixing ratio of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e250.jpg, which stresses that the dye signal reflects both excitation and inhibition [30], [31]. This is in contrast to studies in which the signal is interpreted as caused by excitation only, but also to studies that presume only 25% inhibition [20]. The latter is based on the fact that about 25% of the neurons in the cat primary visual cortex are inhibitory (GABA-immunoreactive) [64]. We cannot exclude the possibility that the dye binds more strongly to inhibitory neurons, the number of active cells however, does in any case not necessarily reflect their functional impact. The NF dynamics are the result of processing across the closely coupled excitatory and inhibitory layers. We showed that the resulting dynamics were only achievable if both layers interact.

For example, as summarized in their review article, Ferster and Miller [65] pointed to the discrepancy between the observed contrast invariance of orientation tuning curves of simple cells and the notion of only weak cortical inhibition. In fact, recent intracellular in vivo recordings in cat have demonstrated strong inhibitory input [54] that may counteract cortical excitatory inputs in a push-pull manner [66].

Relation to hypotheses of cortical representation of motion

Different amplitudes of the dye signal correspond to different levels of the degree of depolarization across the observed neural populations. Thus, with higher amplitudes, the probability rises that the signal reflects supra-threshold activity [12], [19]. Information about the stimulus trajectory should essentially be encoded at high amplitudes of activation, most likely as a propagating wave of spiking activity [12], [48]. In contrast, low amplitude activity may reflect initial passive spread [18] without a close coupling to the input speed of the stimulus. Therefore, we were particularly interested in the relationship between the speed of lateral propagation and the level of activity.

We found that the model did not reproduce cortical axonal conduction speeds of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000919.e251.jpg m/s [12], [67], [68]. Part of the problem to measure speeds of lateral spread is due to the fact that in vivo VSD imaging (as well as in vivo intracellular recordings) must measure weak deviations from baseline levels to capture the earliest input and are therefore confronted with low signal-to-noise-ratios. These initially very low activity levels observed in the VSD experiments (see [12] and also [68]) were indeed not high enough in the data to significantly contribute to our simple grid-based search algorithms if not explicitly parameterized. As we tried to keep our model as simple as possible we therefore did not further optimize these lowest activity levels. However, it is crucial to note that the measured propagation speeds are of postsynaptic origin. Thus, propagation of activity signals integrative properties of the neurons and the surrounding network, rather than reporting true axonal conduction speeds [25]. These network properties are reflected in the continuously rising VSD signal which we indeed captured in our NF model. Hence, in the model the speed of the initial low-amplitude activity spread (coded by the light blue color in Figure 7) demonstrates initially emerging activity within the network and shows only weak dependency on the stimulus speed as observed experimentally [12].

In conclusion, the model captures the main signatures of the spread of activity observed in the new data set including real motion. This generalization to novel stimuli similar to those used to determine the model parameters supports our choice of the neural field model architecture.

It remains an open question as to which neuronal mechanisms lead to the line-motion illusion. The dichoptic experiments of Hikosaka et al. [13] demonstrated that the retina and the LGN cannot be the processing stages where the LM effect arises. Instead, Hikosaka et al. argued that attentional effects are responsible for the line-motion sensation [13], [14]. Another proposal was that since area MT plays a major role in motion integration in humans and monkeys, the LM effect should arise along the dorsal pathway. For instance, Baloch and Grossberg [15] discussed a number of processing steps in their model of the V1–V2–MT–MST pathway that could give rise to the LM illusion.

A recent human fMRI study by Larsen et al. [16] demonstrated that MT+ activation in response to true motion was similar to activation following the presentation of a corresponding illusory moving stimulus. The three-stage theory proposed by these authors suggests that in the first stage higher areas solve the “correspondence problem” (i.e., identify two images as two successive views of the same object), while in the second stage MT+ computes the motion trajectory between these two object representations, and finally the computed trajectory is back-projected to V1 and filled in by a sequence of visual representations of the object [16]. Ahmed et al. [17] reported such feedback activation in area 17/18 from area 19/21 in ferrets using VSD imaging. The activation was motion-dependent and locked to the offset of the first stimulus in their apparent-motion paradigm. This feedback activated the path between subsequent retinotopic stimulus representations in area 17/18 and was interpreted as to play an important role in the computation of continuous motion.

In contrast, our model does not explicitly account for an impact of back projections from higher cortical areas. Nevertheless, these back propagating waves are slower than the initial LM effect characterized here by the immediate drawing-out of activity representing rapid motion [17]. However, we cannot exclude that particularly later parts of the VSD responses may involve feedback from higher visual areas.

The study by Jancke et al. [12] suggested that bottom-up processes are the main source of the initial line-motion activity. Our model confirms that lateral interactions in primary visual cortex are sufficient to generate responses to the illusory LM stimulus that are nearly indistinguishable from the responses to true motion.

Supporting Information

Text S1

A) Covariance matrix adaptation evolution strategy; B) Linear stability analysis of the homogeneous solution.

(0.97 MB PDF)

Acknowledgments

The authors thank Jennifer Meyer for the implementation of a preliminary version of the model. We thank John Lipinski and Benedict Ng for helpful comments on the manuscript. We thank Gregor Schöner, Sebastian Schneegans, and Eray Basar for fruitful discussions. The authors are grateful to the three anonymous reviewers for useful remarks on the manuscript.

Footnotes

The authors have declared that no competing interests exist.

The authors acknowledge support from the German Federal Ministry of Education and Research (BMBF) within the Bernstein group “The grounding of higher brain function in dynamic neural fields”, project number 01GQ0704. The funder had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

References

1. Deco G, Jirsa VK, Robinson PA, Breakspear M, Friston K. The dynamic brain: From spiking neurons to neural masses and cortical fields. PLoS Comput Biol. 2008;4:e1000092. [PMC free article] [PubMed]
2. Freeman WJ. Neurodynamics: an exploration in mesoscopic brain dynamics. Springer; 2000.
3. Dinse HR, Jancke D. Time-variant processing in V1: From microscopic (single cell) to mesoscopic (population) levels. Trends Neurosci. 2001;24:203–205. [PubMed]
4. Jancke D, Erlhagen W, Dinse HR, Akhavan AC, Giese M, et al. Parametric population representation of retinal location: Neuronal interaction dynamics in cat primary visual cortex. J Neurosci. 1999;19(20):9016–9028. [PubMed]
5. Erlhagen W, Bastian A, Jancke D, Riehle A, Schöner G. The distribution of neuronal population activation (DPA) as a tool to study interaction and integration in cortical representations. J Neurosci Methods. 1999;94:53–66. [PubMed]
6. Beurle RL. Properties of a mass of cells capable of regenerating pulses. Philos Trans R Soc Lond B Biol Sci. 1956;240:55–94.
7. Wilson R, Cowan D. Excitatory and inhibitory interactions in localized populations of model neurons. Biophys J. 1972;12:1–24. [PubMed]
8. Wilson R, Cowan D. A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue. Kybernetik. 1973;13:55–80. [PubMed]
9. Amari SI. Dynamics of pattern formation in lateral-inhibition type neural fields. Biol Cybern. 1977;27:77–87. [PubMed]
10. Jancke D. Orientation formed by a spot's trajectory: A two-dimensional population approach in primary visual cortex. J Neurosci. 2000;20: RC86:1–6. [PubMed]
11. Erlhagen W, Jancke D. The role of action plans and other cognitive factors in motion extrapolation: A modelling study. Vis Cogn. 2004;11:315–340.
12. Jancke D, Chavane F, Na'aman S, Grinvald A. Imaging cortical correlates of illusion in early visual cortex. Nature. 2004;428:423–426. [PubMed]
13. Hikosaka O, Miyauchi S, Shimojo S. Focal visual attention produces illusory temporal order and motion sensation. Vision Res. 1993;33:1219–1240. [PubMed]
14. Hikosaka O, Miyauchi S, Shimojo S. Voluntary and stimulus-induced attention detected as motion sensation. Perception. 1993;22:517–526. [PubMed]
15. Baloch A, Grossberg S. A neural model of high-level motion processing: Line motion and formotion dynamics. Vision Res. 1997;37:3037–3059. [PubMed]
16. Larsen A, Madsen KH, Lund TE, Bundesen C. Images of illusory motion in primary visual cortex. J Cogn Neurosci. 2006;18:1174–1180. [PubMed]
17. Ahmed B, Hanazawa A, Undeman C, Eriksson D, Valentiniene S, et al. Cortical dynamics subserving visual apparent motion. Cereb Cortex. 2008;18:2796–2810. [PubMed]
18. Grinvald A, Lieke E, Frostig RD, Hildesheim R. Cortical point-spread function and long-range lateral interactions revealed by real-time optical imaging of macaque monkey primary visual cortex. J Neurosci. 1994;14:2545–2568. [PubMed]
19. Petersen CH, Grinvald A, Sakmann B. Spatiotemporal dynamics of sensory responses in layer 2/3 of rat barrel cortex measured in vivo by voltage-sensitive dye imaging combined with whole-cell voltage recordings and neuron reconstructions. J Neurosci. 2003;23:1298–1309. [PubMed]
20. Rangan AV, Cai D, McLaughlin DW. Modeling the spatiotemporal cortical activity associated with the line-motion illusion in primary visual cortex. Proc Natl Acad Sci U S A. 2005;102:18793–18800. [PubMed]
21. Blumenfeld B, Bibitchkov D, Tsodyks M. Neural network model of the primary visual cortex: from functional architecture to lateral connectivity and back. J Comput Neurosci. 2006;20:219–241. [PMC free article] [PubMed]
22. Grimbert F, Chavane F. Neural field model of VSD optical imaging signals. 2007. Technical Report 6398, INRIA.
23. Sharon D, Jancke D, Chavane F, Na'aman S, Grinvald A. Cortical response field dynamics in cat visual cortex. Cereb Cortex. 2007;17:2866–2877. [PubMed]
24. Sterkin A, Lampl I, Ferster D, Grinvald A, Arieli A. Realtime optical imaging in cat visual cortex exhibits high similarity to intracellular activity. Neurosci Lett. 1998;51:41.
25. Sit Y, Chen Y, Geisler W, Miikkulainen R, Seidemann E. Complex dynamics of V1 population responses explained by a simple gain-control model. Neuron. 2009;64:943–956. [PMC free article] [PubMed]
26. Trappenberg T. Tracking population densities using dynamic neural fields with moderately strong inhibition. Cogn Neurodyn. 2008;2:1476–1492. [PMC free article] [PubMed]
27. Tusa R, Rosenquist A, Palmer L. Retinotopic organization of areas 18 and 19 in the cat. J Comp Neurol. 1979;185:657–678. [PubMed]
28. Albus K. A quantitative study of the projection area of the central and the paracentral visual field in area 17 of the cat. Exp Brain Res. 1975;24:159–179. [PubMed]
29. Somogyi P, Freund TF, Cowey A. The axo-axonic interneuron in the cerebral cortex of the rat, cat and monkey. Neuroscience. 1982;7:2577–2607. [PubMed]
30. Meyer J, Igel C, Jancke D. Modelling dynamic activity patterns in early visual cortex based on voltage sensitive dye experiments. In: Hossmann KA, editor. Symposium “Neuro-Visionen 3”, Perspektiven in Nordrhein-Westfalen. Verlag Ferdinand Schöningh; 2006. pp. 193–195.
31. Symes A, Wennekers T. Spatiotemporal dynamics in the cortical microcircuit: A modelling study of primary visual cortex layer 2/3. Neural Netw. 2009;22:1079–1092. [PubMed]
32. Hansen N, Ostermeier A. Completely derandomized self-adaptation in evolution strategies. Evol Comput. 2001;9:159–195. [PubMed]
33. Suttorp T, Hansen N, Igel C. Efficient covariance matrix update for variable metric evolution strategies. Mach Learn. 2009;75:167–197.
34. Beyer HG. Evolution strategies. Scholarpedia J. 2007;2:1965.
35. Igel C, Erlhagen W, Jancke D. Optimization of dynamic neural fields. Neurocomputing. 2001;36:225–233.
36. Igel C, von Seelen W, Erlhagen W, Jancke D. Evolving field models for inhibition effects in early vision. Neurocomputing. 2002;44–46:467–472.
37. Igel C, Glasmachers T, Heidrich-Meisner V. Shark. J Mach Learn Res. 2008;9:993–996.
38. Akaike H. A new look at the statistical model identification. IEEE Trans Automat Contr. 1974;19:716–723.
39. Wilks SS. The large-sample distribution of the likelihood ratio for testing composite hypotheses. Ann Math Stat. 1938;9:60–62.
40. Fan J, Zhang C, Zhang J. Generalized likelihood ratio statistics and Wilks phenomenon. Ann Stat. 2001;29:153–193.
41. von der Malsburg C, Cowan JD. Theory of ontogenesis of orientation domains – Intracortical dynamics part. 1981. Technical Report 81–1, Department of Neurobiology, Max-Planck-Institut for Biological Chemistry, Göttingen.
42. Cai D, Rangan AV, McLaughlin DW. Architectural and synaptic mechanisms underlying coherent spontaneous activity in V1. Proc Natl Acad Sci U S A. 2005;102:5868–5873. [PubMed]
43. Kenet T, Bibitchkov D, Tsodyks M, Grinvald A, Arieli A. Spontaneously emerging cortical representations of visual attributes. Nature. 2004;425:954–956. [PubMed]
44. Heeger DJ. Normalization of cell responses in cat striate cortex. Vis Neurosci. 1992;9:181–97. [PubMed]
45. Carandini M, Heeger DJ, Movshon JA. Linearity and normalization in simple cells of the macaque primary visual cortex. J Neurosci. 1997;17:8621–8644. [PubMed]
46. Graham N, Sutter A. Normalization: contrast-gain control in simple (fourier) and complex (non-fourier) pathways of pattern vision. Vision Res. 2000;40:2737–2761. [PubMed]
47. Freeman T, Durand S, Kiper D, Carandini M. Suppression without inhibition in visual cortex. Neuron. 2002;35:759–771. [PubMed]
48. Jancke D, Erlhagen W, Schoner G, Dinse H. Shorter latencies for motion trajectories than for flashes in population responses of cat primary visual cortex. J Physiol. 2004;556:971–982. [PubMed]
49. Eggert J, van Hemmen JL. Unifying framework for neuronal assembly dynamics. Phys Rev E Stat Nonlin Soft Matter Phys. 2000;61:1855–1874. [PubMed]
50. Hirsch JA, Gilbert CD. Synaptic physiology of horizontal connections in the cat's visual cortex. J Neurosci. 1991;11:1800–1809. [PubMed]
51. Borg-Graham LJ, Monier C, Fregnac Y. Visual input evokes transient and strong shunting inhibition in visual cortical neurons. Nature. 1998;393:369–373. [PubMed]
52. Carandini M, Ferster D. Membrane potential and firing rate in cat primary visual cortex. J Neurosci. 2000;20:470–484. [PubMed]
53. Tamás G, Somogyi P, Buhl EH. Differentially interconnected networks of GABAergic interneurons in the visual cortex of the cat. J Neurosci. 1998;18:4255–4270. [PubMed]
54. Monier C, Chavane F, Baudot P, Graham LJ, Frégnac Y. Orientation and direction selectivity of synaptic inputs in visual cortical neurons: A diversity of combinations produces spike tuning. Neuron. 2003;37:683–680. [PubMed]
55. Hirsch JA, Alonso JM, Reid RC, Martinez LM. Synaptic Integration in Striate Cortical Simple Cells. J Neurosci. 1998;18:9517–9528. [PubMed]
56. Kisvárday ZF, Toth E, Rausch M, Eysel UT. Orientation-specific relationship between populations of excitatory and inhibitory lateral connections in the visual cortex of the cat. Cereb Cortex. 1997;7:605–618. [PubMed]
57. Gilbert CD, Wiesel TN. Columnar specificity of intrinsic horizontal and corticocortical connections in cat visual cortex. J Neurosci. 1989;9:2432–2442. [PubMed]
58. Buzás P, Kovács K, Ferecskó AS, Budd JM, Eysel UT, et al. Model-based analysis of excitatory lateral connections in the visual cortex. J Comp Neurol. 2006;499:861–881. [PubMed]
59. LeVay S. Patchy intrinsic projections in visual cortex, area 18, of the cat: Morphological and immunocytochemical evidence for an excitatory function. J Comp Neurol. 1988;269:265–274. [PubMed]
60. Matsubara J. Local, horizontal connections within area 18 of the cat. Prog Brain Res. 1988;75:163–72. [PubMed]
61. McGuire BA, Gilbert CD, Rivlin PK, Wiesel TN. Targets of horizontal connections in macaque primary visual cortex. J Comp Neurol. 1991;305:370–392. [PubMed]
62. Kisvárday ZF, Martin KAC, Freund TF, Maglóczky Z, Whitteridge D, et al. Targets and quantitative distribution of GABAergic synapses in the visual cortex of the cat. Exp Brain Res. 1986;64:541–552. [PubMed]
63. Beaulieu C, Somogyi P. Targets and quantitative distribution of GABAergic synapses in the visual cortex of the cat. J Comp Neurol. 1990;2:296–303. [PubMed]
64. Gabbott PLA, Somogyi P. Quantitative distribution of GABA-immunoreactive neurons in the visual cortex (area 17) of the cat. Exp Brain Res. 1986;61:323–331. [PubMed]
65. Ferster D, Miller K. Neural mechanisms of orientation selectivity in the visual cortex. Annu Rev Neurosci. 2000;23:441–71. [PubMed]
66. Troyer T, Krukowski A, Priebe N, Miller K. Contrast-invariant orientation tuning in cat visual cortex: thalamocortical input tuning and correlation-based intracortical connectivity. J Neurosci. 1998;18:5908–5927. [PubMed]
67. Bringuier V, Chavane F, Glaeser L, Frégnac Y. Horizontal propagation of visual activity in the synaptic integration field of area 17 neurons. Science. 1999;29:695–699. [PubMed]
68. Benucci A, Frazor R, Carandini M. Standing waves and traveling waves distinguish two circuits in visual cortex. Neuron. 2007;55:103–117. [PMC free article] [PubMed]

Articles from PLoS Computational Biology are provided here courtesy of Public Library of Science