The Influence of Mexican Hat Recurrent Connectivity on Noise Correlations and Stimulus Encoding

Noise correlations are a common feature of neural responses and have been observed in many cortical areas across different species. These correlations can influence information processing by enhancing or diminishing the quality of the neural code, but the origin of these correlations is still a matter of controversy. In this computational study we explore the hypothesis that noise correlations are the result of local recurrent excitatory and inhibitory connections. We simulated two-dimensional networks of adaptive spiking neurons with local connection patterns following Gaussian kernels. Noise correlations decay with distance between neurons but are only observed if the range of excitatory connections is smaller than the range of inhibitory connections (“Mexican hat” connectivity) and if the connection strengths are sufficiently strong. These correlations arise from a moving blob-like structure of evoked activity, which is absent if inhibitory interactions have a smaller range (“inverse Mexican hat” connectivity). Spatially structured external inputs fixate these blobs to certain locations and thus effectively reduce noise correlations. We further investigated the influence of these network configurations on stimulus encoding. On the one hand, the observed correlations diminish information about a stimulus encoded by a network. On the other hand, correlated activity allows for more precise encoding of stimulus information if the decoder has only access to a limited amount of neurons.


MODEL DETAILS
We used adaptive exponential integrate and fire neurons (AEIF) (Brette and Gerstner, 2005). Each neuron of population X ∈ {E, I} indexed by i comprised two differential equations describing the development of the membrane potential V X,i and an adaptation current w X,i over time: with −g L,X (V X,i − E L ) defining the leak term and g L,X ∆ T exp being an approximation to exponential rise of the sodium (Na + ) current of an action potential, assuming that the activation of Na +channels is instantaneous and inactivation can be neglected (Fourcaud-Trocme et al., 2003). Moreover, C denotes the capacitance, g L,X is the leak conductance of population X, E L the leak reversal potential, V T the threshold parameter, ∆ T the slope factor, a X the sub-threshold adaptation of population X, τ w the adaptation time constant, and I XY,i denotes the afferent and recurrent currents. Moreover, each neuron exhibited a threshold condition: if V X,i ≥ V cut : V X,i ← E L (with clamping for τ refr,X ), (4) w X,i ← w X,i + b X , i.e. each neuron was reset to the leakage potential in case the membrane potential rose beyond the cut off point V cut . In addition, this crossing was considered to be the particular point in time a neuron spiked. Then w X,i was increased by a constant term b X , modeling spike-triggered adaptation. Moreover, the dynamic of V X,i was clamped for a period of τ refr,X to E L after the crossing to simulate an absolute refractory period.
Note that we used different refractory periods τ refr,X for excitatory and inhibitory neurons. We assumed a shorter period for fast spiking inhibitory neurons of 2 ms in comparison to the excitatory cells with 3 ms. Similarly, we used different leak conductances g L,X among the two populations to account for the experimental observation that the conductance is about twice as large in inhibitory cells (Karayannis et al., 2007).
Furthermore, the currents I XY,i with Y ∈ {A, E, I} were modeled as conductance based: Afferent currents were mediated by fast exponentially decaying AMPA-like synapses. Recurrent input consisted of excitatory as well as an inhibitory part. The excitatory current depended on a mixture of fast AMPA and slow bi-exponential NMDA-like synapses, whereas the inhibitory current was based on fast GABA A -like synapses: where α ∈ [0, 1] determines the ratio between AMPA and NMDA receptors and t k j is the time point of the kth spike of pre-synaptic neuron j. In addition,ḡ XY denotes the maximum conductance, i.e. the synaptic coupling strength, τ ... are time constants, E Y is the synaptic reversal potential, and ∆ ij is the synaptic delay between neuron j and i as given below. Figure 1: Maximum absolute deflection of membrane voltage at resting potential (-65 mV) for single presynaptic spikes as function of connection strength. Excitatory post-synaptic potential (EPSP) size on the left and inhibitory post-synaptic potential (IPSP) on the right; in blue for excitatory neurons and inhibitory neurons in dotted green.

SYNAPTIC DELAYS
Synaptic delays ∆ ij were distance dependent and computed as with a basic time delay ∆ ij that was sampled uniformly from the interval [0.2, 1.2] ms. Moreover, v denotes the signal velocity and d ij is the distance between cells i and j. We assumed a signal velocity of v = 200 kpx/s for the one-dimensional networks (see below) and v = 13.3 kpx/s in the two-dimensional setting which corresponds to 0.2 m/s in our cat cortex scaling.

NUMERICAL SIMULATION
We integrated all equations using a simple Euler scheme with a particular step size dt. The runtime was optimized by pre-compiling Python code into faster C code using the Numba package (The Numba Development Team, 2015).
To allow fast computation of the exponential upstroke as in g L, , we precomputed the function. We used 10 6 sampling points and linearly sampled membrane voltages in a range between 92.5 mV and V cut . We did not interpolate in between points. For membrane voltages falling in between two sampling points, we used the result of the closest point. For values smaller than the lower bound, we simply took a voltage of 92.5 mV. Randomly sampled membrane potentials V X,i from a uniform distribution over the interval [E L , V T ] were used as initial conditions. All other dynamic parameters were set to 0. Moreover, before any measures or statistics were computed a simulation was run for at least 1 second to allow the network to settle into a stationary state.

PARAMETER SETTINGS
In here we provide all parameter settings, see table 1. The relation between potential sizes and recurrent connections strength is shown in figure 1. The pinwheel map of the heterogeneous input for the 2D networks is also depicted in figure 2, showing the preferred orientations s P O i of 100 × 100 neurons.  Figure 3: One-dimensional network realization with an excitatory as well as inhibitory population. Recurrent connectivity are sampled from two Gaussian with widths σ E and σ I with periodic boundary conditions.

NOISE CORRELATIONS IN ONE-DIMENSIONAL NETWORKS
We turned to simulations of ring networks similar to the model by Hansen et al. (2012) to test our hypothesis in a simpler topology as well. Such a one-dimensional model is depicted in figure 3. Our networks consisted of N E = 10,000 excitatory and N I = 2500 inhibitory neurons placed on a ring with a length scale of 10 kpx. Due to periodic boundary conditions the maximum distance between two cells is therefore 5 kpx. Each neuron received K E = 200 excitatory and K I = 100 inhibitory recurrent inputs. In this one-dimensional setting we only considered homogeneous input of 15 Hz with K A = 100 afferent connections per cell.
In figure 4A the spike count correlation coefficient as a function of distance between cell pairs is shown for a Mexican hat network close to self-sustained activity (σ E = 125 px < σ I = 250 px,ḡ EE = 0.35 nS, g IE = 0.6 nS). Each blue dot corresponds to the r SC value for an individual cell pair calculated over 50 stimulus presentations of 1 second length. For comparison the noise correlations of an inverse Mexican hat (σ E = 250 px > σ I = 125 px) are shown in the small inset. No distance dependency of correlations was observed in this case.
The distance dependency of spike count correlations in one-dimensional Mexican hat networks were fitted using an exponentially decaying sine wave with where r SC (d) is expectation across neuron pairs with the same distance, C r SC is the correlation fit for d = 0, λ r SC is the mean decay lifetime, and ω r SC denotes the spatial frequency. The fits were calculated using curve fitting from SciPy's optimization library (Jones et al., 2001). This optimization was repeated 10 times with random initial values and the parameters of the best fit were chosen.
The average value for bins of 20 px is plotted in green in figure 4A. For the Mexican hat network this average correlation coefficient can be neatly fitted with a damped sinusoidal wave (see equation 14); shown on top of the green curve in red. For close by cells correlation is relatively large and we fitted an initial value C r SC ≈ 0.64. Furthermore, the exponential decay was quite slow with λ r SC ≈ 3.4 kpx. Lastly, the spatial frequency was fast with ω r SC ≈ 1.4 per thousand pixels.
Moreover, we investigated whether the spatial scales of the noise correlations correspond to the spatial scales of the moving bumps. The average spatial Fourier spectrum alongside the average spatial autocorrelation is shown in figures 4C and D, respectively. Indeed, the highest power was measured for 1.4 cycles per thousand pixels, which nicely fits the spatial scale of the noise correlations. Moreover, the spatial power spectrum does not show a spiky peak at only 1.4 cycles, but a slight broadening around this value. This matches the observations one makes looking directly at the network activity. Usually we counted 14 bumps within the entire network of 10,000 neurons, but occasionally due to splits and merges of bumps one more or less could be found.
Furthermore, we observed an increase of correlated variability with increase in the size of the integration window. As seen in figure 5A, longer integration time windows yielded larger amplitudes of the damped sine wave. Moreover, increasing the integration window beyond 1 second yielded a saturation of average correlation of close by cells as depicted in figure 5B. We recall that the noise correlation coefficient is defined as the ratio between the covariance and the product of individual neurons' standard deviations of spike counts. In this line of thought, we took a look at the relation between the mean covariance of nearby neurons and the geometric mean spike count variance of individual neurons. The results are depicted in figure 5D. We observed that for smaller window sizes the covariance was much smaller in comparison to the variance. Hence, much of the cells' covariation was hidden by the spiking threshold.
Moreover, as shown in the temporal power spectrum (figure 5C), the movement of the bumps was rather slow. This is indicated by the pronounced power in the low frequency band. The peak of the temporal power spectrum was found for 0.2 Hz.

DEPENDENCY OF NOISE CORRELATIONS ON NETWORK PARAMETERS
We explored the dependency of the fitted parameters from equation 14 as a function of different connection strengthsḡ EE as well as different Mexican hat sizes in terms of inhibitory connection spread σ I . The fitted parameters as well as some other basic statistics such as the average firing rate, the coefficient of variation (CV), and the average noise correlation coefficient (r SC ) for cells at most 100 px apart are shown in figure 6.
As one can see in figure 6A, the average firing rate increased slowly with increasing excitatory to excitatory strength until self-sustained activity is reached as depicted by the orange border. Beyond the border the average firing rate quickly diverged. Moreover, figure 6B shows an increase of the coefficient of variation (CV) averaged across all excitatory cells with an increase ofḡ EE . As figure 6B illustrates with increasing strength the average CV exhibited a rather sudden jump beyond 1 after which a seemingly exponential growth happened. Of note, the value of 1 is a standard reference point because it is also the CV for any homogeneous Poisson process. Hence, the spiking activity of the cells was more random and irregular than Poisson spiking. CVs larger than 1 are usually indicative of bursting activity (Christodoulou and Clarkson, 1995). Bursting neurons show periods of no activity followed by periods with strong repetitive firing. We could identify that such bursting behavior was indeed present for each individual cell by looking at the spike and voltage traces (data not shown). Yet, if one takes the joint neural activity into account, one can understand that the bursts are due to the moving activity bumps. Thus, cells may be participating in a bump, i.e. bursting, or are silent if they are part of an inactive region. Figure 6C shows the noise correlation coefficient r SC over 50 stimulus presentations using a time window of 1 second and averaged over cell pairs being at most 100 px apart. Accordingly, the coefficient is averaged among close by cells to get an estimate of the magnitude or initial amplitude of the distance dependent noise correlations. For lower connection strengthḡ EE , the average noise correlation was essentially zero. For valuesḡ EE ≥ 0.2 nS there was a rapid increase followed by a plateau around 0.4 close to the bifurcation to self-sustained activity and a further increase beyond this point. When inspecting the raster plots for small values ofḡ EE < 0.2 nS, we observed no bump activity. Rosenbaum and Doiron (2014) showed that the homogeneous fixpoint with asynchronous activity can still be stable even in case of Mexican hat coupling due to a finite network size. Indeed, this was the case here. We tested different network sizes, i.e. we increased the neuron density, and scaled the number of afferent and recurrent connections and the connection spread proportionally while keeping the connection strength fixed. Using low values ofḡ EE = 0.15 nS and 0.2 nS, we observed an increase in noise correlations between neurons that are close by for larger network sizes, as depicted in figure 7.
The middle row (D-F) of figure 6 shows the fitted parameters for equation 14 for different connection strengths. The parameters were only fitted forḡ EE ≥ 0.2 nS.
With increase in connection strength, the frequency of the noise correlations ω r SC decreased slightly in figure 6D. Likewise, by looking at raster plots of the corresponding runs, we could observe a decrease in bump frequency, too (data not shown).
Furthermore, the relation betweenḡ EE and λ r SC was non-monotonic. First, in figure 6E we identified an increase in λ r SC which means a slower decay of average correlations, followed by a decrease towards instability. Overall the exponential decay of correlations with distance was rather slow, similar to what was observed in figure 4A. Hence, there were even considerable positive and negative correlations between neurons with wide distances of 4000 or 5000 pixels apart.
As expected, the shape of the curve for the fitted initial correlations C r SC in figure 6F strongly resembles the one above ( figure 6C) showing the empirical average noise correlation coefficient r SC . Since the empirical one is based on cell pairs at most 100 px apart, its magnitude is a bit lower than compared to C r SC . As a reminder, C r SC fits the hypothetical average noise correlation value for cell pairs with no distance in between them.
The bottom row (G-I) of figure 6 displays the parameter fits as a function of inhibitory spread σ I for a fixed connection strength ofḡ EE = 0.4 nS. Similar to increasing connection strength, widening the inhibitory connection spread yielded a decrease in the spatial frequency of noise correlations (G). Surprisingly, there was no effect of inhibitory connection spread on the decay λ r SC (H) and only a minor one on the amplitude parameter C r SC (I). We observed a weak decrease in C r SC for wider spreads.

BOUNDARY CONDITIONS
We further investigated the question: Do Mexican hat networks yield spatial activity patterns for other boundary conditions as well?
We ran simulations with absorbing as well as reflecting boundary conditions in one-dimensional networks. In the former case, neurons close to one of the boundaries received less inputs from their neighbors because outgoing connections beyond the boundary were simply cut off. In the latter case connectivity at the boundaries was denser than in the center of the network. Connections that were supposed to reach over the boundary were mirrored and projected back into the network. Figure 8 shows network activity for these two different boundary conditions. As for circular boundary conditions, we observed the emergence of spatial patterns in Mexican hat networks close to self-sustained activity. In turn, these spatial patterns yielded noise correlations (data not shown). Likewise, no patterns and correlations were observed for balanced and inverse spreads (data not shown).

NOISE CORRELATIONS AND ADAPTATION
Clearly, adaptation plays an important role in the emergence of noise correlations. Adaptation yields movement of the bump solutions. Accordingly, removing adaptation from the spiking neurons reduced the magnitude of noise correlations for nearby neurons considerably (figure 9).
As shown by Hansel and Sompolinsky (1998)   The smallest window we use is 20 ms because for smaller windows only a tiny minority of cells exhibits at least a single spike within 50 stimulus presentations, which strongly biases our calculation of the spike count correlation coefficient. Increasing window size increases the amplitude of the damped wave. Top right (B): Average noise correlation for cells at most 100 px apart as a function of integration window size. This corresponds to less than the first fourth of a cycle of the damped sine waves one the top left. The envelope shows the standard deviation among all cell pairs. Increasing the window size beyond 1 second yields a saturation effect. Bottom left (C): Average temporal power spectrum. Movement of bumps happens on slow time scales as the peak and pronounced power on very low frequencies suggest. Additionally, there is another shallow peak in the gamma frequency band, as shown in the inset. Bottom right (D): Contribution of spike count variance and covariance to the average r SC shown above. For smaller window sizes the covariance is much smaller than the neuron pairs' gemoetric mean variance of spike counts. For increasing window size the difference is reduced and both lines become parallel. The expected number of spikes per neuron is shown as a black dotted line. show parameter explorations of different excitatory to excitatory connection strengths.ḡ EE is explored using a Mexican hat topology with σ E = 125 px > σ I = 250 px, andḡ IE = 0.6 nS. The orange line marks the bifurcation to self-sustained activity. All values are computed from 50 stimulus presentations and averaged across 10 network realizations. Values are computed from excitatory neurons only. Error bars mark the standard deviation across network realizations. The top row shows from left to right: (A) the average firing rate as a function of connection strength, followed by (B) the average coefficient of variation (CV), and (C) the average noise correlation coefficient (r SC ) for cells with at most 100 px distance between them. The r SC is computed over an integration window of 1 second and based on 50 stimulus presentations. Middle row shows the curve parameters from fitting distance dependend noise correlation with equation 14. From left to right: (D) Average revolutions per thousand pixels (ω r SC ), (E) decay constant (λ r SC ), and (F) maximum average noise correlation (C r SC ). Values were computed only forḡ EE ≥ 0.2 nS because for smaller values no distance dependent noise correlations were observed. Bottom row (G-I) same as in the middle but for exploring inhibitory connection spread σ I usinḡ g EE = 0.4 nS and σ E = 125 px.   Figure 8: Spiking activity of a Mexican hat network (σ E = 125 px < σ E = 250 px,ḡ EE = 0.3 nS, g IE = 0.6 nS) for reflecting (top) as well as absorbing (bottom) boundary conditions. A small black dot represents a spike of a particular neuron at a particular point in time. Neurons are ordered according to their position on the ring. Small plots show the activity histograms over time (horizontal) as well as over space (vertical). As for circular boundary conditions spatial patterns are observed. Figure 9: Average noise correlation within a Mexican hat network close to self-sustained activity (σ E = 125 px < σ E = 250 px,ḡ EE = 0.4 nS,ḡ IE = 0.6 nS). Values are averaged across pairs at most 100 px apart over 50 stimulus presentations of 1 second each with homogeneous input. Error bars mark standard deviations across 10 network realizations. On the left without adaptation a X = b X = 0 and on the right with adaptation a E = 10a I = 2 nS and b E = 10b I = 50 pA. Removing adaptation significantly reduces correlations (Wilcoxon rank-sum test, p < 0.001, 10 network samples).