|Home | About | Journals | Submit | Contact Us | Français|
Shared neural variability is ubiquitous in cortical populations. While this variability is presumed to arise from overlapping synaptic input, its precise relationship to local circuit architecture remains unclear. We combine computational models and in vivo recordings to study the relationship between the spatial structure of connectivity and correlated variability in neural circuits. Extending the theory of networks with balanced excitation and inhibition we find that spatially localized lateral projections promote weakly correlated spiking, but broader lateral projections produce a distinctive spatial correlation structure: Nearby neuron pairs are positively correlated, pairs at intermediate distances are negatively correlated and distant pairs are weakly correlated. This non-monotonic dependence of correlation on distance is revealed in a new analysis of recordings from superficial layers of macaque primary visual cortex. Our findings show that incorporating distance-dependent connectivity improves the extent to which balanced network theory can explain correlated neural variability.
The spiking activity of cortical neurons is often characterized by their average response over a large number of trials, prompting a wealth of theoretical studies relating the structure of neuronal networks to their trial averaged firing rate dynamics1. However, trial averages do not capture the stochastic and irregular dynamics characteristic of cortical populations and the nervous system in general2. Indeed, trial-to-trial fluctuations are central to contemporary theories of cortical computation3,4. A deep mechanistic understanding of neuronal variability remains an open challenge.
Early theoretical studies deduced that variable spiking activity could arise through a balancing of strong, yet opposing, excitatory and inhibitory synaptic inputs5,6. Expanding on this conjecture, van Vreeswijk and Sompolinsky7 showed that networks of recurrently coupled model neurons robustly create a state where strong excitation is approximately balanced by inhibition, creating a push-pull dynamic that generates irregular spiking activity. More recently, balanced networks have been implicated in theories of optimal coding8, working memory9, and stimulus tuning10. Numerous experimental studies have established that excitation is often approximately balanced by inhibition in cortical circuits11,12,13,14,15,16,17. In sum, balanced networks provide a parsimonious model of the irregular spiking activity observed in cortical circuits.
Early balanced network models produced asynchronous activity through sparse connectivity7,18. However, several experimental studies reveal that local cortical networks are densely connected, with connection probabilities between nearby neurons sometimes exceeding 40 percent19,20,21,22. These data imply substantial overlap between local synaptic inputs, which could, in principle, synchronize cortical networks. However – counter to intuition – balanced networks with dense connectivity show weak spike train correlations23. This “asynchronous state” results from the correlated excitatory (e) or inhibitory (i) afferents to neuron pairs being actively cancelled by a strong negative e-i correlation establishing weak correlations even when connectivity is not sparse23.
Consistent with the predicted asynchronous state, some multi-unit extracellular recordings show noise correlations that are nearly zero on average24. However, a majority of population recordings in cortex reveal comparatively large correlations25,26. Several studies suggest that the magnitude of noise correlations is dependent on many factors27, including arousal28, attention29, anesthetic state23,24,30,31 and cortical layer32,33. Finally, while in vivo whole cell recordings reveal strong positive e-e and i-i correlations coexisting with strong e-i correlations13, these correlation sources do not always perfectly cancel as predicted by some theoretical models28. Taken together, these studies show that cortical circuits can exhibit both weak and moderate noise correlations, at odds with predictions from the current theory of balanced networks23.
In this study, we generalize the theory of correlations in densely connected balanced networks to include the widely observed dependence of synaptic connection probability on distance34,21. We show that spatially broad recurrent projections disrupt the asynchronous state, producing a signature spatial correlation structure: Nearby pairs of neurons are positively correlated on average, pairs at intermediate distances are negatively correlated and distant pairs are weakly correlated. These positive and negative correlations cancel so that the average correlation between pairs sampled randomly over a large range of distances is nearly zero. We uncover this non-monotonic dependence of correlation on distance in recordings from superficial layers of macaque primary visual cortex, but only after correcting for a latent source of shared fluctuations. Our findings decouple balanced excitation and inhibition from asynchronous network activity, greatly extending the applicability of balanced network theory to explain cortical dynamics.
We consider a network of excitatory and inhibitory exponential integrate-and-fire model neurons. Neurons provide recurrent, lateral synaptic input to one another and receive feedforward synaptic input from a non-local presynaptic population. A detailed mathematical analysis of correlations in the limit of large network size is provided in Supplementary Note S.1 and Supplementary Figures 1–3. Below, we provide an outline of these theoretical results and confirm their predictions using computer simulations. We first use a simplified network model to demonstrate how the asynchronous state considered in previous theoretical work23 is broken by heterogeneous input correlations35. We then consider a more realistic model where neurons belong to a continuous spatial domain and connection probability depends on distance.
To demonstrate the mechanisms affecting correlations in recurrent networks, we first simulated a simplified network of N=20,000 neurons (half excitatory and half inhibitory) that all receive the same fluctuating feedforward input and are each connected with probability 0.25 (Figure 1a). Despite the fact that neuron pairs share all of their feedforward input and 25% of their recurrent synaptic input on average, spiking activity was asynchronous, with an average pairwise spike count correlation of 6.8×10−5 (Figure 1b,c).
This small average correlation is a defining characteristic of the asynchronous state. Mathematically, this state is realized when spike count covariances in the network satisfy23
where CSS denotes the average spike count covariance between pairs of neurons in the recurrent network, N is the number of neurons in the network and O(1/N) denotes asymptotic proportionality to 1/N for large N. Note that covariance and correlation scale identically with network size in balanced networks, so we discuss them interchangeably23.
Our theoretical analysis proceeds by noting that spike count covariance is inherited from synaptic input covariance26 and therefore the two scale similarly with N in the asynchronous state,
where CII denotes the average covariance between neurons’ synaptic inputs.
Synaptic inputs can be decomposed into their feedforward and recurrent sources, I=F+R, so that neurons’ input covariances decompose as
where CFF is the average covariance between neurons’ feedforward input currents, CRR between their recurrent inputs and CRF between one neuron’s recurrent and the other neuron’s feedforward synaptic input. Recurrent synaptic input, R, is composed of positive contributions from lateral excitatory synaptic inputs and negative contributions from inhibitory (R=e−i).
Shared input fluctuations are visualized by averaging the inputs to several neurons, so that the unshared contributions average out (Figure 1d). Overlapping inputs cause CFF and CRR to be positive (Figure 1d – the blue and red curves have large fluctuations). If feedforward input correlation is moderate, CFF~O(1), then recurrent input tracks the feedforward input so that CRF is negative and nearly perfectly cancels the positive sources of correlations (i.e., 2CRF=−(CFF+CRR)+O(1/N); see Supplementary Note S.1). As a result, the covariance between the total synaptic inputs is weak, CII~O(1/N), (Figure 1d – the black curve has small fluctuations). This cancellation arises naturally in the balanced state and does not require a precise tuning of model parameters23. Since spiking correlations are inherited from synaptic input covariance26, this cancellation of input covariances leads to small, O(1/N), spike count correlations.
To study the impact of heterogeneity on correlations in balanced networks, we modified the above model by dividing the neurons into two populations. Each population received a separate feedforward input (Figure 1e). The two feedforward input sources were statistically identical but uncorrelated. Recurrent connectivity was not changed: Neurons were randomly connected without respect to population membership (identical to Figure 1a). This input heterogeneity dramatically changed the structure of correlations in the network. Pairs of neurons in the same population had strongly positive spike count correlations on average (0.34), while neuron pairs from opposite populations were negatively correlated with a nearly identical correlation magnitude (−0.34), and the average correlation between all pairs was nearly zero (4.2×10−4; Figure 1f,g).
The mechanism responsible for this change in correlations can be understood by again separating the synaptic input covariance into its recurrent and feedforward sources, but generalizing the decomposition to account for neuron “distance” to obtain
Here, CII(d) is the average covariance between input currents to pairs of neurons separated by distance d where d=0 for neurons in the same population and d=1 for opposite population pairs, and similarly for the other terms. Feedforward input is only correlated between neurons in the same population, so CFF(0)>0, but CFF(1)=0.
In contrast, recurrent connections do not respect population membership and thus neither do the statistics of recurrent input,
Since covariances CRR(d) and CRF(d) do not depend on d, but CFF(d) does, cancellation cannot be achieved in Eq. (1) for both d=1 and d=0 simultaneously. In other words, the one “copy” of shared recurrent synaptic input cannot cancel both versions of the feedforward synaptic inputs. The loss of cancellation causes the total synaptic current shared by neurons in the same population to inherit shared fluctuations from their feedforward inputs (Figure 1h; compare blue and black), giving rise to positive O(1) correlations between same-population pairs. A competitive dynamic introduces negative correlations between neurons in opposite populations (Figure 1f,g). A similar mechanism was considered in a recent theoretical study35.
For illustrative purposes, we considered a simplified network model with discrete subpopulations, but correlations and connectivity in many cortical circuits depend on continuous quantities like physical distance or tuning similarity36,19,21. Next, we generalize these findings to more biologically realistic networks with connection probabilities that depend on neuron distance.
We next considered a network of Ne=40,000 excitatory and Ni=10,000 inhibitory model neurons arranged on a square-shaped domain modeling a portion of a cortical layer. The neurons receive feedforward synaptic input from a separate layer of Poisson-spiking neurons and are connected with a probability that decays with distance (Figure 2). Specifically, the probability of a connection between two neurons in the recurrent network obeys
where d is the distance between the neurons measured along the two-dimensional network, ~ denotes proportionality, g(d;σ2)~exp(−d2/(2σ2)) is a Gaussian-shaped function and αrec approximately represents the average length of a recurrent synaptic projection. Similarly, the probability of a synaptic projection from a neuron in the feedforward layer to a neuron in the recurrent layer decays with distance like a Gaussian with width parameter αffwd, where distance is measured parallel to the cortical surface (Figure 2). We next show that the asynchronous state is realized when αrec<αffwd, then show that the asynchronous state cannot be realized when αrec>αffwd.
As above, asynchronous spiking requires a cancellation between input covariances, cf. Eq. (1), except that d now represents continuous instead of binary distance. Therefore, conditions on asynchrony require first an understanding of how input covariances depend on pairwise neuron distance.
Overlapping feedforward synaptic projections introduce O(1) correlations between the feedforward inputs to neuron pairs. Since nearby pairs share more feedforward inputs, these correlations are distance-dependent. Specifically, synaptic divergence causes feedforward input correlations to be twice as broad as synaptic projection widths (Figure 2b),
The fact that CFF(d)~O(1) at first seems to preclude the possibility of an asynchronous state since CFF(d) is one component of CII(d) in Eq. (1) and the asynchronous state requires CII(d)~O(1/N). However, the asynchronous state is realized under a cancellation between positive (CFF and CRR) and negative (CRF) sources of correlations in Eq. (1). Cancellation at all distances requires that all correlation sources have the same shape (Supplementary Note S.1), meaning that
The implications of this requirement on the spatial profile of spiking correlations are clarified by noting that recurrent synaptic input is generated by spike trains in the. Synaptic divergence causes the correlations between neurons’ recurrent synaptic inputs to be broader in space than the correlations between spike trains according to (Figure 2c)
Here, σRR is the width of correlations between neurons’ recurrent synaptic input currents and σSS is the width of spike train correlations in the recurrent network. In general, we use α’s to denote the widths of synaptic projections and σ’s to denote the widths of correlations.
Correlations between recurrent inputs are constrained by the cancellation required in the asynchronous state. Specifically, Eq. (2) requires that the width of correlations between recurrent synaptic inputs satisfies
Combining the two expressions for σRR above yields
The existence of a real solution to Eq. (3) requires that αffwd>αrec; in other words, the spatial width of the recurrent projections must be narrower than the width of feedforward projections for the asynchronous state to exist. Further, Eq. (3) implies that σ2SS>σ2FF−α2rec so that spike train correlations are spatially narrower than correlations between feedforward input currents. Thus, recurrent dynamics actively sharpen the spatial profile of correlations in the asynchronous state (compare to the sharpening of tuning curves in previous work37).
To test these theoretical findings, we performed network simulations with feedforward synaptic projections broader than recurrent projections (Figure 3a). The simulations confirmed that CRR(d) and CRF(d) decayed similarly with distance to CFF(d) (Figure 3b red, purple and blue). This allowed a cancellation between positive and negative sources of correlations, so that correlations between neurons’ total synaptic currents and between their spike trains were weak over all distances (Figure 3b–e). Despite their small average, spike count correlations had a larger standard deviation (stdev=0.11; Figure 3d), consistent with results for non-spatial networks23 (Figure 1b). Neurons in the network receive strong excitation that is canceled by strong inhibition on average (Figure 3f and Supplementary Figure 4a,b), confirming that the network maintains a balanced state. Correlations computed from simulations agreed with closed-form mathematical predictions (Figure 3e, compare black and dashed red; see Supplementary Note S.1 for equations). Additional simulations confirmed that mean correlations decay toward zero at increasing network size (Supplementary Note S.2 and Supplementary Figure 4c,d).
As noted above, the cancellation between positive and negative correlations necessary for the asynchronous state cannot be realized when recurrent projections are broader than feedforward (αrec>αffwd) since Eq. (3) cannot be solved in this case. Instead, neuron pairs inherit correlations from overlapping feedforward inputs so that CSS(d)~O(1).
This prediction was confirmed by numerical simulations identical to those discussed above, but with recurrent projections broader than feedforward (Figure 4a). As predicted, recurrent input correlations (Figure 4b, red and purple) were too spatially broad to cancel with the more sharply decaying feedforward correlations (Figure 4b, blue) so that the total input correlation between nearby neurons was large (Figure 4b, black; compare to Figure 3b). This effect introduced moderately strong correlations between nearby spike trains (Figure 4c,d,e) that do not decay to zero at increasing network size (Supplementary Note S.2 and Supplementary Figure 4e,f). Nevertheless, the network maintained excitatory-inhibitory balance (Figure 4f and Supplementary Figure 4a,b).
Since recurrent inputs must cancel feedforward inputs in balanced networks, CRF(d) is negative (Figures 1, ,33 and and4;4; Supplementary Note S.1). Moreover, broad recurrent projections cause CRF(d) to decay slowly with distance (Figure 4b, purple). Through Eq. (1), this imparts a non-monotonicity in the dependence of CII(d) on d (Figure 4b, black) and spike count correlations inherit this non-monotonic shape (Figure 4e).
Following the same argument made for the homogeneous network, the spike count correlations averaged over neuron pairs at all distances is O(1/N) (Figure 4e, dashed gray and Supplementary Note S.1). However, as noted above, the average correlation over each distance cannot be O(1/N). Hence, there must be a cancellation between positive and negative correlations at different distances. As in Figure 1e–h, a competitive dynamic causes nearby neurons to be positively correlated and more distant neurons to be negatively correlated. This competitive dynamic does not extend beyond the reach of recurrent projections, so sufficiently distant neurons are weakly correlated. Hence, correlation decreases then increases with distance. This non-monotonicity can be explained more precisely using a mathematical theory of correlation transfer, (Figure 4e red dashed curve and Supplementary Note S.1). The heterogeneity of positive and negative correlations at different distances increases the standard deviation of pairwise correlations, but only modestly (stdev=0.16; Figure 4d, compare to Figure 3d).
In summary, when recurrent projections are spatially narrower than feedforward projections (αrec<αffwd as in Figure 3), correlations are weak between pairs of neurons at all distances. When recurrent projections are broader than feedforward (αrec>αffwd as in Figure 4), nearby neurons are positively correlated, neurons at moderate distances are negatively correlated and distant neurons are weakly correlated. Moreover, the average correlation between pairs of neurons sampled randomly at all distances is small. The non-monotonic dependence of correlation on distance is a distinct signature of correlations arising from broad recurrent projections. We next investigate whether this correlation structure predicted by our theory is present in cortical recordings.
We next asked whether our theoretical characterization of correlations in spatially extended networks can explain correlations in a cortical circuit. Layers 2/3 and layer 4C of macaque primary visual cortex (L2/3 and L4C) provide an ideal circuit for testing our predictions. Pairs of neurons in L2/3 exhibit moderately large noise correlations that decay with distance, but neurons in L4C, which are a primary source of inter-laminar input to L2/3, exhibit extremely weak pairwise noise correlations32,33 (Figure 5 with data from previous studies36,33).
Neurons in L4C receive much of their feedforward input from thalamic projections, which form spatially broad synaptic fields, around 1 mm in diameter, but lateral projections within macaque L4C form narrower, sub-millimeter, synaptic fields34. Thus, our theoretical prediction that correlations are weak when αffwd>αrec is consistent with the weak pairwise correlations observed between L4C neurons in vivo (as in Figure 3).
Inter-laminar projections from L4C to L2/3 have a similar sub-millimeter width as excitatory intra-laminar projections within L4C and lateral projections from inhibitory basket cells in L2/3 form sub-millimeter synaptic fields, similar to those in L4C34. Excitatory neurons in L2/3, however, form long-range lateral synaptic projections with synaptic fields spanning several millimeters in diameter34. Our theoretical results can be generalized to this setting where inhibitory and excitatory projections have different spatial profiles (Supplementary Note S.1). This extension predicts the same correlation structure reported in Figure 4. However, correlations measured in L2/3 are positive on average over a broad range of distances36 (Figure 5b), in disagreement with this prediction.
We hypothesized that this inconsistency is explained by recent studies showing that much of the correlated variability measured in L2/3 arises from a low-dimensional shared source of latent variability38,31,30,39,40. We conjectured that this shared variability increases pairwise correlations in L2/3 at all distances, thereby “washing out” the negative correlations predicted by our theory.
To search for low-dimensional variability in our data, we used Gaussian process factor analysis41,31 (GPFA), a statistical algorithm that extracts shared fluctuations from a population of spike trains (see Materials and Methods). Applying this algorithm to our L2/3 recordings reveals a source of one-dimensional covariability that decays with distance (Figure 5c). This distance-dependence implies that nearby neurons are affected similarly by the latent variable.
To test whether one-dimensional latent variability explains the discrepancy between our theoretical predictions and data, we built a two-layer network model representing a 10 mm by 10 mm square of cortex (Figure 6a). The first layer, representing L4C, was similar to the model in Figure 3 with the profile of feedforward and recurrent projections chosen to match experimentally constrained thalamic and lateral projection widths34. The second layer, representing L2/3, was similar to Figure 4, with feedforward synaptic input from excitatory neurons in the L4C model and recurrent projection widths also chosen to match anatomical measurements. To capture latent variability in L2/3, the feedforward synaptic input to each neuron in the second layer was modulated by a time-varying, multiplicative gain modulation. We chose a multiplicative source of variability to be consistent with the properties of low-dimensional variability previously reported in macaque V131, but an additive source of latent variability would produce similar overall results. The gain modulation contributes an O(√N) source of covariance to the feedforward inputs that the recurrent network cannot cancel23. To capture the distance-dependence of latent variability (Figure 5c), the magnitude of the gain modulation was heterogeneous across the network in such a way that nearby neurons received similar modulations and more distant neurons received less similar modulations.
Simulations of this two-layer model revealed that correlations between neurons in L4C were extremely small on average (Figure 6b,c gray), consistent with our theoretical predictions (Figure 3) and consistent with in vivo recordings (Figure 5). Correlations in the model L2/3 layer were moderately large and positive over all distances (Figure 6b,c black), comparable to in vivo recordings (Figure 5).
Thus, our model recovers the coarse structure of correlations in L4C and L2/3. However, our explanation of positive correlations in L2/3 is unsatisfying because the addition of globally shared variability destroyed the distinct non-monotonic relationship between correlation and distance predicted by our theory (compare Figure 4e to Figure 6b). We next asked whether this structure could be recovered by filtering out globally shared variability. To accomplish this, we computed the “residual correlation” matrix estimated by GPFA. Residual correlations approximate the spike count correlations with the contribution from low-dimensional variability removed31.
Residual correlations computed between the simulated L2/3 spike trains exhibit the predicted non-monotonic dependence on distance, corroborating the ability of the GPFA algorithm to extract low-dimensional variability and leave the structure of residual correlations intact. We next computed the mean residual correlation in macaque L2/3 as a function of electrode distance. In doing so, we observed the same non-monotonic dependence of residual correlation on distance predicted by our theory (Figure 7b; further statistical analysis in Supplementary Note S.3 and Supplementary Figure 5).
In summary, combining theoretical analysis and computer simulations of a multi-layer network reveals a parsimonious model of the sources of shared variability in a visual cortical circuit in vivo. Under this model, positive correlations introduced by shared thalamic inputs to L4C neurons are actively canceled by negative correlations arising from recurrent circuitry so that pairs of L4C neurons at all distances exhibit weak average spike count correlations32,33. Correlations between neurons in L2/3 are introduced by overlapping feedforward inputs from L4C and a low-dimensional source of variability. Correlations arising from overlapping L4C projections are filtered by recurrent circuitry in L2/3 to promote a non-monotonic dependence of correlation on distance. This non- monotonic correlation structure is washed out by low-dimensional latent variability, but can be recovered using GPFA to estimate and remove this variability.
Previous theoretical work on spatially homogeneous balanced networks with dense connectivity shows that they produce very weak spike train correlations 23. We have generalized this theory to account for heterogeneous inputs and distance-dependent connection probability. In this framework we have made two important discoveries.
First, in agreement with the original findings, when lateral synaptic projections are spatially narrower than incoming feedforward projections, correlations are extremely weak on average at all distances. This theoretical finding can explain the weak pairwise correlations observed between neurons in middle layers of macaque primary visual cortex32,33. However, correlations measured in cortical recordings are not always weak25. Our second finding is that networks with broader lateral than feedforward projections produce correlations that do not decay to zero at increasing network size.
In previous studies of balanced networks with spatially homogeneous or clustered connectivity23,42, the asynchrony condition (CSS~O(1/N)) is satisfied and population averaged pairwise correlations vanish in the large network limit. In contrast, spatially extended networks with broad lateral projections violate the asynchrony condition, and consequently the expected pairwise correlations at a specific distance do not vanish. Nonetheless, mean excitatory and inhibitory currents balance and firing rates are moderate even when the asynchrony condition is violated (Supplementary Figure 4a,b and Supplementary Note S.1). This represents a novel solution for balanced networks that, for the first time, formally decouples network-wide asynchrony from excitatory-inhibitory balance.
We focused on the dependence of correlations on distance, but correlations also depend on tuning similarity. Partitioning L2/3 neuron pairs by tuning similarity reveals that correlations are strongest between similarly tuned neurons36(Supplementary Figure 6a,b). Modifying our computational model to capture tuning-dependent correlations produces a non-monotonic dependence of residual correlation on tuning similarity in some parameter regimes, but the relevant parameters have not been measured experimentally (Supplementary Figure 6c–e and Supplementary Note S.4). Nevertheless, the modified theory could explain negative correlations previously observed in computer simulations of networks with tuning-specific connectivity32 and the finding that negative correlations are more frequent between disparately tuned neurons in V143.
As with nearly any computational model, many of the parameters used in our simulations may not reflect their corresponding values in specific cortical areas of specific species. However, our theoretical analysis does not depend on the precise values of these parameters. Our finding that the asynchronous state requires αrec<αffwd is a fundamental property of networks with balanced excitation and inhibition.
We used a simplified model of a visual cortical circuit. In reality, pyramidal neurons in V1 form both local and long-range projections, connection probability in primate V1 depends on both distance and tuning similarity, and these dimensions are coupled44. Moreover, connectivity properties of inhibitory neurons depend on their sub-type 45. We modeled unidirectional connections from L4 to L2/3, but L4 also receives indirect feedback from L2/3 through deeper cortical layers. Spike trains in our model feedforward layer were modeled by homogeneous Poisson processes, in contrast to the oscillatory firing rates evoked by drifting grating stimuli in the data we analyzed. Our model can be extended to account for these additional features without affecting our overall conclusions.
Our findings have important implications for the interpretation of correlations in neural recordings. The average (residual) correlation between cell pairs sampled across a large range of distances could be extremely small, even when nearby pairs are positively correlated with moderate magnitude (Figures 4 and and7).7). Hence, subtracting low-dimensional latent variability and partitioning neuron pairs by distance can reveal correlation structure that would otherwise not be apparent. A previous study31 computed residual correlations as a function of distance in primate V1, but did not report a non-monotonic dependence. While we cannot be certain why their findings differ from ours, the accurate estimation of residual correlations with GPFA depends on the amount of data used to estimate shared variability. Our data are well-suited for this purpose, as they contain over 800 pairs of units per recording on average.
There is a long history of computational models of cortical circuits that consider either networks with spatially dependent coupling1 or balanced excitation and inhibition in spatially homogeneous networks7,23. Only recently has the spatial structure of cortical connectivity been included in networks with balanced excitation and inhibition46,9,37, and guiding theoretical principles are lacking. Our theory has taken this spatial structure into account and produced two core predictions for cortical circuits with long range lateral connections: First, nearby neurons exhibit significant positive correlations; second, the dependence of pairwise correlation on pairwise distance is non-monotonic. These predictions are clearly falsifiable and hence represent strong tests of our theory. The superficial layers of visual cortex have long-range lateral connections34, making them a suitable test bed for our theory of correlations. After accounting for a source of global variability, both of our predictions were verified from population recordings in macaque V1 (Figure 7). Further, a similar noise correlation structure was reported in recordings from mouse V147. The successful validation of our predictions marks our theory as a promising framework for studying the structure of neural variability in cortical circuits. Nevertheless, there are many aspects of cortical dynamics that remain unexplained by balanced networks, such infrequent yet large membrane fluctuations during spontaneous dynamics15,48. Capturing these dynamics in cortical models with balanced architectures remains an open challenge.
We modeled a square of cortex with N neurons, Ne of which are excitatory and Ni inhibitory. The membrane potential of neuron j from the excitatory (a = e) or inhibitory (a = i) population obeyed exponential integrate-and-fire (EIF) dynamics,
each time that exceeds a threshold at Vth, the neuron spikes and the membrane potential is held for a refractory period τref then reset to a fixed value Vre. The leak current is given by IL(V) = −gL(V − EL) and a spike-generating current is defined by f(V) = gLΔTexp[(V −VT)/ΔT]. For excitatory neurons, τm = Cm/gL = 15ms, EL = −60mV, VT = −50mV, Vth = −10mV, ΔT = 2mV, Vre = −65mV and τref = 1.5ms. Inhibitory neurons were the same except τm = 10ms, ΔT = 0.5mV and τref = 0.5ms.
Synaptic input currents were defined by
where Fa(t) is the feedforward input and Ra(t) the recurrent input to neuron j in population a = e, i. The feedforward input was modeled differently for different figures, as described below. The recurrent input was defined by
where is the nth spike time of neuron k in population b = e, i. The 1/√N scaling of synaptic weights is a defining feature of the balanced network formalism and captures the balance between excitatory and inhibitory currents as well as intrinsically generated temporal variability for large N 7. Each term represents the synaptic weight from presynaptic neuron k in population b to postsynaptic neuron j in population a. For all simulations, we modeled synaptic kinetics using ηb(t) = exp(−t/τb)/τb for t > 0 where τe = 6ms and τi = 5ms. All networks are “dense” in the sense that connection probabilities are O(1)23.
For the model in Figure 1, there were N = 20,000 neurons, half of which were excitatory and half inhibitory. For each (presynaptic) neuron, we randomly and uniformly chose 2500 excitatory and 2500 inhibitory postsynaptic neurons in the network. Postsynaptic neurons were sampled with replacement, so that a single presynaptic neuron can make multiple contacts with a postsynaptic neuron. The synaptic weight of each connection depended on the pre- and postsynaptic neuron types (excitatory or inhibitory). Specifically,
where jee = 12.5mV, jie = 20 mV and jii = jei = −50 mV. Note that synaptic weights were scaled by √N= 141 in Eq. (4), so that the actual synaptic weight of each contact was on the order of 0.1mV.
For Figure 1a–d, the feedforward input to each neuron was given by the sum of an input bias and a smoothly varying signal, Here, s(t) is a shared source of smooth, unbiased Gaussian noise defined by its auto-covariance function,
τs = 40ms sets the correlation timescale and σs = 0.1 mV/ms scales the magnitude of the fluctuations. The terms me = 0.015 mV/ms and mi = 0.01 mV/ms introduce a static bias to the input current. The model for Figure 1e–h was identical, except that two independent realizations, s1(t) and s2(t), of the shared input were generated. Half of the neurons received s1(t) and the other half received s2(t). Firing rates for Figure 1a–d were 7.6 Hz on average for excitatory neurons and 3.8 Hz for inhibitory neurons. For Figure 1e–h, average firing rates were 7.4 Hz for excitatory and 3.8 Hz for inhibitory neurons.
To model the spatially extended recurrent network in Figures 3 and and4,4, we arranged Ne = 40,000 excitatory and Ni =10,000 inhibitory EIF model neurons on a uniform grid covering a two-dimensional square domain. The feedforward layer was modeled by a population of NF = 5625 excitatory Poisson-spiking neurons on a uniform grid covering a square that is parallel to the recurrent network. Feedforward input to the recurrent layer was defined by
where is the nth spike time of neuron k in the feedforward layer. Each spike train in the feedforward layer was modeled as independent Poisson processes with rate rF = 5 Hz.
To simplify calculations, we measured distance in units of the side-length of the square domain. In these units, the domain is represented as the unit square, Γ = [0, 1] × [0, 1]. Neurons were connected randomly and the probability that two neurons were connected depended on their distance measured periodically on Γ. The precise algorithm for generating connections is described in Supplementary Note S.5. This algorithm assures that the expected number of synaptic contacts from a presynaptic neuron at coordinates y = (y1, y2) in population b to a postsynaptic neuron at x = (x1, x2) in population a is given by
where g(u; α) is a wrapped Gaussian distribution37. Out-degrees were and . It follows that the network-wide average number of synaptic inputs to excitatory and inhibitory neurons in the recurrent network was K = 3715. Synaptic weights were determined by Eq. (5) where jee = 40 mV, jie = 120 mV, jei = jii = −400 mV, and jeF = jiF = 120 mV. Note again that these terms were divided by √N ≈ 224 as indicated in Eq. (4), so that the actual synaptic weights were between 0.18 mV and 1.8 mV. Excitatory and inhibitory recurrent projection widths were αrec = αe = αi = 0.05 for Figure 3 and αrec = 0.25 for Figure 4. Feedforward connection widths were αffwd = 0.1 in both figures. For the simulations in Figure 3, average firing rates were 3.9 Hz for excitatory and 6.2 Hz for inhibitory neurons. For the simulations in Figure 4, average rates were 4.0 Hz for excitatory and 6.1 Hz for inhibitory neurons.
The first layer (L4C) in the model in Figure 6 was identical to the model in Figure 3 except that αffwd = 0.1, αe = 0.05 and αi = 0.03. The length units used in Figures 6 and and7a7a were determined by interpreting the network domain, Γ, as a 10 mm by 10 mm square. Thus, in physical dimensions, αffwd = 1 mm, αe = 0.5 mm and αi = 0.3 mm. Average firing rates in the L4C layer were 3.7 Hz for excitatory and 6.1 Hz for inhibitory neurons.
Connectivity in the second layer (L2/3) in Figure 6 was identical to the L4C layer except that αe = 0.15 (or 1.5 mm), αffwd = 0.05 (or 0.5 mm). The spike times from Eq. (6) for the L2/3 layer were given by the spike times of neurons in the L4C layer, so that NF = 50,000. The magnitude of feedforward connectivity was also modified by setting KeF = 1406, KiF = 113 and jeF = jiF = 220 mV. A shared gain modulation was implemented by altering feedforward input currents according to
The shared gain modulation, L(t), is a realization of unbiased Gaussian noise defined by its auto-covariance function,
with correlation timescale τL = 40 ms. The dimensionless weight factor, wL(x), depended on the coordinates, x = (x1, x2) Γ, of the neuron and was given by
where c = 0.5 and σL = 0.25 (in physical dimensions, c = 5 mm and σL = 2.5 mm) so that neurons near the center of the network to receive a more intense gain modulation and the strength of the modulation decays slowly with distance from the center. This imparted a long-range distance dependence to the correlations (Figure 6c), similar to that observed in the cortical recordings (Figure 5b). Average firing rates in the L2/3 layer were 5.7 Hz for excitatory and 9.1 Hz for inhibitory neurons.
All simulations and numerical computations were performed on a MacBook Pro running OS X 10.9.5 with a 2.3 GHz Intel Core i7 processor. All simulations were written in a combination of C and Matlab (Matlab R 2014b, Mathworks). The differential equations defining the neuron model were solved using a simple forward Euler method with time step 0.1 ms and all simulations were run for a 22 s duration.
Anesthesia was induced with ketamine (10 mg/kg) and maintained during preparatory surgery with isoflurane (1.0–2.0% in 95% O2). Anesthesia during recordings was maintained with sufentanil citrate (6–18 μg·kg−1·h−1). Vecuronium or pancuronium bromide (0.1–0.15 mg·kg−1·h−1) was used to suppress eye movements. Physiological signs were monitored to ensure adequate anesthesia and animal well-being. Vital signs (EEG, ECG, blood pressure, end-tidal PCO2, temperature, and airway pressure) were monitored continuously. We used supplementary lenses to bring the retinal image into focus. At the end of the recording session, animals were killed and tissue was processed histologically to verify recording locations. All procedures were approved by the Institutional Animal Care and Use Committee of the Albert Einstein College of Medicine.
The pairwise correlations used to make the distributions in Figure 5a were computed using the same methods to those reported previously33 and include 1,613 L2/3 pairs and 469 L4 pairs. The data were recorded from nine anesthetized, adult male macaque monkeys (Macaca fascicularis). Recording were made with a group of five to seven linearly arranged (305 μm spacing) platinum-tungsten electrodes or tetrodes (Thomas Recording), inserted normally to the cortical surface. The electrodes were advanced together through cortex, sampling in 200 μm intervals until all electrodes had exited into white matter. Location of the middle layer sites involved a number of criteria including nominal depth, histology, and CSD analysis, as detailed previously33.
Data for distance-dependent correlations reported in Figures 5b, 5c and and7b7b were from recordings described previously36. The data were recorded from seven hemispheres of four adult male macaque monkeys (Macaca fascicularis). There was one experimental group of normal animals. The array consisted of a 10×10 grid of silicon microelectrodes (1 mm in length) spaced 400 μm apart. We inserted the array partially into cortex, resulting in recordings confined mostly to layers 2–3. In two cases, we recorded simultaneously with a group of seven linearly arranged (2 mm extent) platinum-tungsten microelectrodes or tetrodes (Thomas Recording), positioned so that the nearest electrode was ~5 mm anterior to one edge of the multielectrode array. In this configuration, the distances between these electrodes and the array ranged from ~5 to 10mm. Neuronal receptive fields in all data sets were within 5° of the fovea. Waveform segments that exceeded a user-defined threshold were sorted offline (Plexon Offline Sorter). We quantified sort quality using the signal-to-noise ratio (SNR) of each candidate unit49, keeping units with an SNR of at least 2.3 (for the data in Figure 5a) or at least 2.75 (for Figures 5b and and7b).7b). In both data sets, we eliminated neurons for which the best grating stimulus did not evoke a response of at least 2 spikes/s. Changing the SNR or responsivity threshold did not qualitatively change any of the results described herein.
Visual stimuli were generated with EXPO and displayed on a linearized CRT monitor (mean luminance 40 cd/m2) with a resolution of 1,024 by 768 pixels and a refresh rate of 100 Hz. Stimuli were presented in a circular aperture surrounded by a gray field of average luminance. We mapped the spatial receptive field (RF) of units by presenting small (0.6°) drifting gratings at a range of spatial positions. We then centered our stimuli on the aggregate RF of the recorded units.
Stimuli were viewed binocularly and presented for 1.28 s, separated by 1.5 s intervals of isoluminant gray screen (except in one penetration, for which the interval was 10 s). We presented full-contrast drifting sinusoidal gratings at 8 or 12 orientations spaced equally (22.5 or 30° increments). The spatial frequency (1–1.3 cpd) and temporal frequency (3 or 6.25 Hz) values were chosen to correspond to the typical preference of parafoveal V1 neurons50. The position and size (3.9–10°) of the grating were sufficient to cover the receptive fields of all the neurons. Stimulus orientation was block randomized, each stimulus orientation was repeated 100–200 times per recording session and each of the eight recording sessions yielded 21–68 units (210–2278 pairs) that met our SNR and responsivity thresholds. This yielded a total of 318 units and 6,907 simultaneously recorded pairs.
To compute spike count correlations for Figures 1, ,3,3, ,44 and and66 we first randomly sampled a subset of the excitatory neurons in the recurrent network. The number of neurons sampled is indicated in each figure caption. Neurons with firing rates less than 1 Hz were excluded from the correlation analysis. Pearson correlation coefficients were computed by counting spikes over 250 ms windows then computing the Pearson correlation coefficients between all pairs. The first 2 s of each simulation was excluded from the correlation analysis.
To obtain the residual correlations in Figure 7a, we randomly sampled 500 neurons with firing rates above 2 Hz from the L2/3 layer simulations. After omitting the first 2 s of the simulation, spike counts were computed using 50 ms bins, then partitioned into a sequence of 250 ms windows. To each window, we directly applied the GPFA algorithm used in a previous study31 and downloaded from http://toliaslab.org/publications/ecker-et-al-2014/. This is a modified version of the algorithm introduced previously41. Briefly, GPFA extracts shared fluctuations from a population of spike trains by fitting spike count covariance matrices to a model of the form
Here, c is a column vector of coefficients quantifying the contribution of latent variability to each neuron’s spiking activity. Latent variability is modeled by a one-dimensional Gaussian stochastic process. Hence, ccT is a rank-one matrix of latent covariances and Q quantifies the residual covariability not captured by the latent variable. Residual correlations were computed using the “window” option (so that spike count covariances were computed over 250ms windows) then converting the resulting covariance matrix, Q, to a correlation matrix using the corrcov command in Matlab (Mathworks).
Spike count correlations in Figure 5b were computed across trials, counting spikes over 250ms time bins and subtracting the trial-averaged rate of each unit over each stimulus orientation. Hence, all reported spike count correlations represent correlated trial-to-trial variability (noise correlations) with stimulus correlations removed. Latent covariance in Figure 5c and residual correlation in Figure 7b were computed from the mean-subtracted spike counts using 50ms time bins and 250ms window sizes (treating each 250ms window as a separate trial) then using the same GPFA algorithms and procedures that were used for the simulated data in Figure 7a (see description above). The three p-values reported in Figure 7b were computed using one-sided, unpaired t-tests with t-values 8.1, 2.1 and 2.5 and degrees of freedom 3,937, 2,404 and 2,508 respectively.
The four error bars in Figure 3e were computed from n=880,057, 2,646,645, 4,420,219 and 4,550,579 pairs respectively. The five error bars in Figure 5b were computed from n=1,688, 2,251, 1948, 458 and 562 pairs respectively, and the same for Figures 5c and and7b.7b. The six error bars in Figure 7a were computed from n=3,885, 11,696, 19,493, 27,521, 35,336 and 26,819 pairs respectively. No statistical methods were used to pre-determine sample sizes but our sample sizes are the same as those reported in previous publications33,36. Data collection and analysis were performed blind to the conditions of the experiments. No animals were excluded from the analysis.
Computer code for all simulations and analysis of the resulting data is included in Supplementary Software.
The data that support the findings of this study are available from the corresponding author upon request.
This work was supported by National Science Foundation grants NSF-DMS-1517828 (RR), NSF-DMS-1313225 (BD), NSF-DMS-1517082 (BD), NSF-DMS-1612913 (JER), NSF-DMS-1516288 (JER) and NSF-DMS-1312508 (JER); National Institute of Health grants R01NS070865 (BD, JER), CRCNS-R01DC015139 (BD), R01EY016774 (AK), R01EY022928 (MAS) and P30EY008098 (MAS); two grants from the Simons Foundation collaboration on the global brain (SCGB#325293MC;BD, BD and 364994AK, AK); by the Eye and Ear Foundation of Pittsburgh (MAS); and by Research to Prevent Blindness (AK, MAS). We are grateful to Tai Sing Lee and Anthony Movshon for research support.
Competing financial interests
The authors declare no competing financial interests.