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

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.


Text S1
A Covariance matrix adaptation evolution strategy Identifying suitable parameters of a neural field model is a non-linear optimization problem. In this study, we used a semi-automatic approach to solve it. We split the optimization problem into a linear and non-linear part: 1. After we have simulated our neural field model with some model parameters x ∈ R 13 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 λ u , λ v , and c using the ordinary least-squares (OLS) method under the constraints λ u , λ v ≥ 0, see Eq 6 in the main paper. This yields the optimal values for λ u , λ v , and c in terms of the mean squared error between aggregated model signal and observed dye patterns for the given the activities of the excitatory and inhibitory layer u and v, respectively.
2. For additional fine-tuning of the model parameters g uu , g uv , g vu , β u , β v , h u , and h v we use a randomized direct optimization algorithm described below.
The objective function in the non-linear part of our model identification procedure is nonconvex, non-differentiable, and multi-modal(i.e., there are undesired local optima). Our method of choice for such problems is the covariance matrix adaptation ES (CMA-ES) explained in the following.
The covariance matrix adaptation evolution strategy. Evolution strategies are randomized direct optimization algorithms [1][2][3][4]. They are one of the major branches of evolutionary algorithms, a class of algorithms drawing inspiration from principles of neo-Darwinian evolution theory. Evolution strategies are iterative optimization methods. In every iteration (also referred to as generation), they sample a set of candidate solutions (the offspring) from a probability distribution over the search space, evaluate these points using the objective function f , and construct a new probability distribution over the search space. In typical ESs, this search distribution is parameterized by a set of candidate solutions, the parent population, and by parameters of the variation operators that are used to create the offspring from the parent population.
Evolution strategies are most frequently applied to real-valued optimization. Arguably the most elaborate ES for real-valued optimization is the covariance matrix adaptation ES (CMA-ES, [5][6][7][8][9]). It relies on on Gaussian mutations, that is, solutions are modified by adding random vectors drawn according to multi-variate zero-mean normal distributions. The CMA-ES is a variable metric algorithm adapting the shape and strength of its Gaussian search distribution. "[The] CMA-ESs represent the state-of-the-art in evolutionary optimization in real-valued R n search spaces" [4]. This claim is backed up by many performance comparisons across different suites of benchmark problems [5,7,[10][11][12][13].
The CMA-ES has been successfully applied for adapting neural dynamics in previous studies. In particular, the algorithm has been compared with the BFGS (Broyden - Fletcher -Goldfarb -Shanno) method based on analytically derived gradients for learning the parameters of neural fields, and the CMA-ES clearly outperformed the quasi-Newton method [14,15]. An application of the CMA-ES to modelling of human car driving behavior using neural attractor dynamics is presented in [16]. We can only expect very weak convergence results for optimization of non-differentiable and multi-modal objective functions in real-valued search spaces. Loosely speaking, if we store the best solution found so far in the optimization process, impose a lower bound on the step size 1 , and assume that we just need to determine the optimal solution with an accuracy of > 0, then the CMA-ES converges if the number of iterations goes to infinity, see [17] for details.
The CMA-ES in detail. For completeness, we briefly describe the CMA-ES as used in this article closely following [9]. For more detailed information, we refer to the literature [5][6][7][8][9].
In each generation k of the CMA-ES, which is shown in Algorithm 1, the lth offspring x (k+1) l ∈ R n , l ∈ {1, . . . , λ}, is generated by additive multi-variate Gaussian mutation and weighted global intermediate recombination: l:λ denotes the lth best individual of the λ offspring ranked according to the values y ) of the objective functions f : R n → R. Considering the best µ of the offspring in the recombination implements non-elitist, rank-based selection. We use the standard choice for the recombination weights w l ∝ ln(µ + 1) − ln(l), The CMA-ES adapts both the n×n-dimensional covariance matrix C (k) of the normal mutation distribution as well as the global step size σ (k) ∈ R + . The covariance matrix update has two parts, the rank-1 update considering the change of the population mean over time and the rank-µ update considering the successful variations in the last generation. The rank-1 update is based on a low-pass filtered evolution path p (k) of successful (i.e., selected) steps and aims at changing C (k) to make steps in the promising direction p (k+1) more likely by morphing the covariance towards p is a normalization constant. The rank-µ update aims at making the single steps that were selected in the last iteration more likely by morphing Putting both updates together, we have The constants c cov and µ cov are fixed learning rates. The learning rate of the covariance matrix update c cov = 2 (n+ √ 2) 2 is roughly inversely proportional to the degrees of freedom of the covariance matrix. The parameter µ cov mediates between the rank-µ update (µ cov → ∞) and the rank-one update (µ cov = 1). The default value is µ cov = µ eff .
An appropriate adaptation of the mutation strength (step size adaptation) is of utmost importance to balance exploration and exploitation during search and to reach an optimum with high accuracy [1,2,18]. The global step size σ (k) is adapted on a faster timescale. It is increased if the selected steps are larger and/or more correlated than expected and decreased if they are smaller and/or more anticorrelated than expected: where E[ N (0, I) ] =χ n is the expected length of an n-dimensional random vector drawn from a zero mean Gaussian distribution with covariance matrix equal to the unit matrix and the (conjugate) evolution path is with learning rate c σ = µ eff +2 n+µ eff +3 and damping factor The matrix C − 1 2 is defined as BD −1 B T , where BD 2 B T is an eigendecomposition of C (B is an orthogonal matrix with the eigenvectors of C and D a diagonal matrix with the corresponding eigenvalues) and sampling N (0, C) is done by sampling BDN (0, I).
The CMA-ES is robust in the sense that it does not rely on tweaking of hyperparameters. The values of the learning rates and the damping factor are well considered and have been validated by experiments on many basic test functions [6]. They need not be adjusted dependent on the problem and are therefore no hyperparameters of the algorithm. Also the population sizes can be set to default values, which are λ = max(4 + 3 ln n , 5) and µ = λ 2 for offspring and parent population, respectively. If we fix C (0) = I, the only hyperparameter to be chosen problem dependent is the initial global step size σ (0) .

B Linear stability analysis of the homogeneous solution
We perform a standard linear stability analysis of the homogeneous solution of the DNF equations 1 and 2 of the main paper in the absence of afferent input. We follow the classic approach of von der Malsburg [19] (see also [20,21]). The field equations without input are written as Here * denotes convolution andu(x, t) , t)). We assume the weight functions to be translation invariant and isotropic.
We restrict our considerations to the homogeneous solution of the system with u(x, t) = u 0 and v(x, t) =ṽ 0 for all x and t: with W uu = w uu (x )dx and W vu = w vu (x )dx . We consider small perturbations (x, t) = u(x, t) −ũ 0 and η(x, t) = v(x, t) −ṽ 0 at time step t: By replacing the transfer function f by its fist-order Taylor approximation the system is linearized: where the dependencies of x and t are removed for simplicity. We have By subtracting the homogeneous solution 4 and 5 from 10 and 11 we get Now we apply the Fourier transform in the spatial dimension. Letĝ denote the Fourier transform of some function g. Eqns 12 and 13 become whereŵ andŵ vu is computed analogously. Writing Eqns 14 and 15 in vector form gives with If the matrix A(k) has two distinct eigenvalues, the homogeneous linear differential equation 18 can be solved as with the initial values (ˆ 0 (k),η 0 (k)) T . The columns of the matrix T (k) are the eigenvectors, where λ(k) + and λ(k) − are the corresponding eigenvalues. These can be compute by setting the characteristic polynomial λ(k) 2 − tr(A(k))λ(k) + det(A(k)) to zero: Solving this quadratic equation gives For convenience, we write If the real part of the largest eigenvalue is negative, then the system is asymptotically stable. In our system, we have to consider a Hopf bifurcation, in which a pair of complex conjugate eigenvalues (λ ± ) of the linearization crosses the imaginary axis of the complex plane, and a saddle point bifurcation, where two fixed points, one stable and one unstable, collide and annihilate each other.
For the first case, we assume that the imaginary part of the eigenvalues is non zero. This means B 2 + C < 0, which implies C < 0. The bifurcation occurs when the real part of λ is zero (i.e., λ ± = 0 ± iω, with ω ∈ R, ω = 0). In this case, we have and therefore and applying this to the Gaussian kernel gives Thus, under the assumption that the discriminant B 2 + C is negative, the system is asymptotically stable if B < 0 (the real part of the largest eigenvalue has to be negative, see above). That is, the following inequality has to be true Now we look at the second case where B 2 + C ≥ 0. Then the system is obviously not asymptotically stable for B ≥ 0, therefore we consider B < 0. The bifurcation point is λ ± = λ = 0 + 0i. From Eqn 20 we see that this equality holds if C = 0: For non negative discriminant and B < 0 this leads to the stability condition To summarize, the analysis reveals that the system is asymptotically stable if both, inequality 26 and inequality 28, are fulfilled.  Tables   Table S1. Summary of the optimized model parameters. The second column gives the description of the parameters. The third column tells whether the parameter was optimized. The parameter values are shown in the last column.