The Dynamical Regime of Sensory Cortex: Stable Dynamics around a Single Stimulus-Tuned Attractor Account for Patterns of Noise Variability

1Neuronal activity throughout cerebral cortex is variable, both temporally during epochs of stationary dynamics and across repeated trials despite constant stimulus or task conditions. Moreover, variability is modulated by a variety of factors, most notably by external sensory stimuli, planning and execution of limb movements, and attention. Modulation of variability occurs at the level of single-neuron activity, e.g., membrane potentials or spike counts, as well as in the patterns of joint activity across populations, as seen in multiunit activity or the local field potential (LFP). Variability modulation shows stereotypical patterns. First, the onset of a stimulus quenches variability overall and, in particular, correlated variability in firing rates that is “shared” across many neurons. Moreover, the degree of variability reduction can depend systematically on the tuning of individual cells. For example, in area MT, variability is quenched more strongly in cells that respond best to the stimulus, and correlations decrease more among neurons with similar stimulus preferences. Although these patterned modulations of variability are increasingly included in quantitative analyses of neural recordings, it is still unclear what they imply about the dynamical regime in which the cortex operates.

2There have been two dynamical mechanisms proposed to explain selected aspects of the modulation of cortical variability by stimuli. In “multi-attractor” models, the network operates in a multi-stable regime in the absence of a stimulus, such that it noisily wanders among multiple possible stable states (“attractors”). This wandering among attractors occurs in a concerted way across the population, resulting in substantial shared variability (Figure 1A, top). Stimuli then suppress this shared variability by pinning fluctuations to the vicinity of one particular attractor (Figure 1A, bottom;,,,,,, ). In chaotic network models, firing rates exhibit strong chaotic fluctuations, and certain types of stimuli can suppress chaos by forcing the dynamical state of the network to follow a specific trajectory, thus quenching across-trial variability (Figure 1B;,,, ). While both the multi-attractor and the chaotic mechanisms can explain the general phenomenon of stimulus-induced reduction of variability, only the former has been proposed to explain the stimulus-tuning of variability reduction. However, even in that case, a considerable fine-tuning of parameters or very strong noise was needed to keep the network in the regime with multiple attractors, such that the system stays near attractors, yet noise can move the system between them.

3Here, we explore a qualitatively different regime of cortical dynamics. We describe activity fluctuations as being driven by noise but shaped by nonlinear, recurrent interactions. In contrast to previous models, our network operates around a single stable point that depends on the stimulus (Figure 1C). Crucially, individual neurons have supralinear (expansive) input/output functions. This causes the gains of neurons, and thus the effective synaptic strengths in the network, to increase with network activation. This is a stochastic generalization of the stabilized supralinear network (SSN) model that has successfully accounted for a range of phenomena related to the stimulus dependence of trial-averaged responses in visual cortex. Introducing stochasticity allows us to model the variability of responses and thus use data on neural variability to identify hallmarks of this regime and distinguish it from previous proposals.

4In our network, stimulus-dependent changes in effective connectivity shape the magnitude and structure of activity fluctuations in the network. Specifically, stimuli change the balance of two opposing effects of recurrent network dynamics on variability: hidden feedforward interactions (“balanced amplification”;, ) and recurrent excitation, which amplify variability and dominate for very weak (spontaneous) inputs; and stabilizing inhibitory feedback, which quenches variability and dominates for stronger inputs.

5By studying this network mechanism in a progression of recurrent architectures with increasingly detailed structure, we find that it naturally and robustly explains the modulation of shared cortical variability by stimuli, including its tuning dependence. We first analyze variability in the simplest instantiation of the model, with two unstructured populations of excitatory (E) and inhibitory (I) cells, and find that an external stimulus can strongly modulate the variability of population activities. In particular, the model predicts stimulus-induced quenching of variability, as well as a reduction of the low-temporal-frequency coherence between local population activity and single-cell responses, as found experimentally. Next, we extend our analysis to a more detailed architecture with structured connectivity to account for the tuning-dependent modulations of Fano factors and noise correlations by stimuli. Critically, these results reveal robust qualitative differences between the predictions of our model and those of previously proposed network mechanisms, based on multi-attractor or chaotic dynamics, for both the spatial patterns and temporal dynamics of variability suppression. We tested these predictions against experimental data and found the SSN model to be the most consistent with previously analyzed data from primary visual cortex (V1) and MT as well as with our novel analyses of published V1 recordings in the awake monkey. Such comparisons of different models are crucial for guiding future experiments that can make targeted measurements to fully resolve the dynamical regime in which the cortex operates—a key first step in identifying the computational strategies underlying perception.

6We used a standard model to study the dynamical evolution of momentary firing rates in a recurrently coupled network of excitatory and inhibitory neurons (Figure 2A;; see also STAR Methods). In this model, neurons integrate their external and recurrent inputs linearly in their membrane potentials,, but their output firing rates, r, are a nonlinear function of the voltage: (Figure 2B). Crucially, we studied variants of this model in which the nonlinearity f is an expansive (supralinear) function (Figure 2B) and in which inhibition was both sufficiently fast and strong and appropriately structured to stabilize the network in the face of recurrent excitation and the supralinear input/output function. This is the stabilized supralinear network (SSN) model. In order to study response variability, we added to this model a stochastic component (slow noise) in the membrane potential dynamics of all cells. Stabilization meant that the network operated around a single steady state, albeit a stimulus-dependent one.

7Real neurons, of course, have an input/output function that ultimately saturates. We focus on an expansive, non-saturating input/output function because V1 cortical neurons show such a relationship between mean voltage and firing rate across their full dynamic range, without saturation even for the strongest visual stimuli. Thus, saturation does not appear to play a role in stabilizing cortical activity, a fact that we capture by using a non-saturating input/output function. Such an expansive input/output function arises in spiking neurons when their firing is driven by voltage fluctuations, with the mean voltage sub- or peri-threshold, a firing regime that produces the highly variable spiking seen in cortical neurons. We assume that the voltage fluctuations giving rise to the expansive input/output function are fast compared to the timescales of variability studied here and do not explicitly model them.

8We focused on analyzing how the intrinsic dynamics of the network shaped fixed input noise to give rise to stimulus-dependent patterns of response variability. We studied a progression of connectivity architectures of increasing complexity, all involving two separate populations of excitatory and inhibitory neurons. We also validated our results in large-scale simulations of spiking neuronal networks.

9We first considered a simple circuit motif: an excitatory (E) unit and an inhibitory (I) unit, recurrently coupled and receiving the same mean external input h as well as their own independent noise (Figure 2A). In this reduced model, the two units represent two randomly connected populations of E and I neurons, a canonical model of cortical networks. Thus, their time-varying activity, and, represents the momentary population-average membrane potential of all the E and I cells, respectively. Despite its simplicity, this architecture accounted well for the overall population response properties in the larger networks, with more detailed connectivity patterns, that we analyzed later.

10Activity in the network exhibited temporal variability due to the stochastic component of the dynamics. We found that this (correlated) variability of and fluctuations, together with their means,, was strongly modulated by the external steady input h (Figures 2C and 2D). When, there was no input to drive the network, and and both hovered around mV, fluctuating virtually independently, with standard deviations essentially matching those that would arise without recurrent connections (gray line in Figure 2D, bottom). For a somewhat larger input, mV, both E and I populations fired at moderate rates (3–4 Hz) (Figure 2D, top), but now also exhibited large and synchronous population fluctuations (Figure 2C, black circle mark). For yet larger inputs ( mV), fluctuations remained highly correlated, but their magnitude was strongly quenched (Figure 2C, green circle mark).

11Figure 2D shows how the temporal (or, equivalently, the across-trial) mean and variability of activities varied over a broad range of input strengths. We observed that population mean increased monotonically with growing external input, first linearly or supralinearly for small inputs, but strongly sublinearly for larger inputs, with growing faster than (Figure 2D, middle;, ). In contrast, variability in both and typically increased for small inputs, peaking around this transition between supralinear and sublinear growth, and then decreased with increasing input (Figure 2D, bottom). Importantly, input modulation of variability required recurrent network interactions. This was revealed by comparing our network to a purely feedforward circuit that exhibited qualitatively different behavior (Figure 2D, gray). In the feedforward circuit, mean remained linear in h, so that mean rates rose quadratically with or h (reflecting the input/output nonlinearity; Figure 2B), and fluctuations in no longer depended on the input strength.

12In order to demonstrate that the overall dynamical regime of the stabilized supralinear network, rather than just a particular instantiation of our model, underlies variability modulation, we used a combination of numerical simulations and analytical results to confirm the robustness of our findings.

13We simulated 1,000 model networks with random parameter values within wide brackets. We found that variability suppression was robust over a broad range of network parameters (connection weights, input strengths and correlations, and the exponent and scale of the firing-rate nonlinearity), as long as they ensured dynamical stability even for strong inputs (Figures S1 and S2). Although the precise amplitude and position of the peak of variance depended on network parameters, the overall non-monotonic shape of variability modulation was largely conserved. In particular, we could show analytically that variability suppression occurs earlier (for smaller input h) in networks with strong connections or, for fixed overall connection strength, in networks that are more dominated by feedback inhibition (Methods S3). More generally, we found that the firing rates at the peak of variability are typically low (2.5 Hz on average over a thousand randomly parameterized stable networks and below 6 Hz for 90% of them; cf. Methods S2). As these rates are comparable to cortical spontaneous firing rates, this predicts that increased sensory drive should generally result in variability quenching in cortical LFPs.

14In order to better understand the robustness of variability suppression in the model, we took advantage of the fact that our network was characterized by a single attractor at each level of the input, h, and analyzed the dynamics of small activity fluctuations,, around this stable state (such that, where is the mean activity in the stable state; STAR Methods). These dynamics are governed by a set of effective connection weights,, that quantify the impact of a small momentary change in the of the presynaptic neuron on the total input to its postsynaptic partner. The dependence of the effective connection weights on the stable state and thus on the external input, h, that determines the stable state is simply given by:where is the strength of the “biophysical” connection from unit j to unit i, and is the slope of the single-neuron firing-rate nonlinearity at the stable state. Importantly, increases with increasing, because f is an expansive, convex nonlinearity (Figure 2B). Thus, in general, effective connectivity increases with increasing h, reflecting the growth of (Figure 2D, middle).

15An increase in effective connectivity can have conflicting effects: it can increase excitatory or driving effects that amplify fluctuations and increase variability, but it can also increase inhibitory feedback, suppressing fluctuations and decreasing variability. Thus, understanding how changes in effective connectivity translate into changes in variability required further analysis (Methods S3). We found that the net behavior of the network indeed included a combination of both effects (Figure S3). As the input grew from zero, variability first rapidly increased, due primarily to the growth of effective feedforward weights (“balanced amplification”; ) but also of recurrent excitatory loops. Then, beginning at firing rates comparable to spontaneous activity as described above, variability steadily decreased with increasing stimulus strength due to increasingly strong inhibitory feedback (Figure 2D, bottom).

16Crucially, we were able to show analytically that variability quenching effects must ultimately dominate, leading to progressively stronger quenching of variability as the input increases. This is due to the faster growth of I activity relative to E activity in the network, which is a robust outcome of dynamic stabilization by feedback inhibition (Figure S1;, ) and which has been observed in rodent S1 and V1. We also found that ignoring the variability-increasing effects, which are characteristic of excitatory-inhibitory dynamics and thus largely absent from models that do not include separate excitatory and inhibitory populations, can fail to capture the full extent of variability modulation and lead to an underestimation of the level of spontaneous variability obtained at zero-to-weak input levels (Figure S4).

17In order to study variability in single neurons and at the level of spike counts, we implemented the two-population architecture of Figure 2A in a network of spiking neurons (Figure 3; STAR Methods). The network consisted of 4,000 E neurons and 1,000 I neurons, randomly connected with low probability and with synaptic weights chosen such that the overall connectivity matched that of the reduced model. Each neuron emitted action potentials stochastically with an instantaneous rate given by Equation 3 (this additional stochasticity accounted for the effects of unmodelled fluctuations in synaptic inputs that occur on timescales faster than the 30 ms effective time resolution of our model; Methods S4). The external input to the network again included a constant term, h, and a noise term that was temporally correlated on a 50 ms timescale with uniform spatial correlations of strength 0.2.

18At the population level, the network behaved as predicted by the reduced model. Neurons fired irregularly (Figure 3A, top), with firing rates that grew superlinearly with small input h but sublinearly with stronger input (Figure S5). Moreover, fluctuations in E and I population activities were strongly synchronized (Figure 3A, upper middle), and LFP variability decreased with increasing h (Figure 3A, lower middle). Importantly, variability quenching also occurred at the level of individual neurons’, accompanied by a reduction of pairwise correlations (Figure 3A, bottom; these required that single neurons shared part of their input noise; Methods S3).

19The model primarily suppressed shared rather than private (to individual neurons) variability (Figure S5), as in experiments. This was because the spatially uniform average connectivity of the network meant that its dynamics were only significantly coupled to patterns of uniform activity across E or across I cells. These patterns were thus the ones affected by stimulus-induced changes in effective connectivity (Figure S3). Correlated noise drove such uniform patterns so that they carried significant variability. Thus, these shared excitatory and inhibitory activity patterns behaved as the activity of the individual units of the previous reduced two-population model, and so variability suppression in the reduced model implied the suppression specifically of shared variability in this more detailed model.

20Our model also accounted for the stimulus-induced modulation of the power spectrum and cross-coherence of LFP and single-cell fluctuations, as observed in V1 of the awake monkey (Figure 3B; ). Strong external input reduced the LFP power at low frequencies, due to enhanced effects of feedback inhibition; increased it at intermediate frequencies, due to the faster timescales associated with relatively enhanced inhibition; and also increased it at high frequencies, due to the larger firing rates, which contributed additional, high-frequency fluctuations in synaptic drive (Figure 3B, top left). This asymmetric modulation of LFP power at low and high frequencies is also seen in experiments (Figure 3B, top right). Moreover, as increasing inputs suppressed variability at the population level, the private noise in the activity of each neuron had a proportionately larger contribution to its overall variability, leading to a drop in pairwise correlations (Figure 3A) and -LFP coherence specifically at low frequencies where the suppression of population variability occurred, as seen in experiments (Figure 3B, bottom).

21Neuronal recordings in visual areas have shown that Fano factors drop at the onset of the stimulus (drifting gratings or plaids) in almost every neuron, which was well accounted for by the randomly connected network we studied above. However, in the experiments, variability did not drop uniformly across cells, but exhibited systematic dependencies on stimulus tuning. This could not be explained by randomly connected architectures, so we extended our model to include tuning dependence in connectivity and input noise correlations.

22We studied an architecture in which the preferred stimulus of E/I neuron pairs varied systematically around a “ring” representing an angular stimulus variable, such as stimulus orientation in V1 or motion direction in MT (Figure 4A; STAR Methods). We describe the case in which the variable is orientation, which ranges from 0 to; identical results describe direction if all angles are doubled. The average input to a cell (either E or I) was composed of a constant baseline, which drove spontaneous activity in the network, and a term that depended on the angular distance between the stimulus orientation and the preferred orientation (PO) of the cell, and that scaled with image contrast, c (Figure 4C). Input noise correlations depended on tuning differences (STAR Methods): cells with more similar tuning received more strongly correlated inputs. The strength of recurrent connections depended on the difference in preferred orientation between pre- and postsynaptic neurons and whether they were excitatory or inhibitory (Figure 4B).

23The bump of stimulus-driven input drove a similar, but narrower, bump of network response (Figures 4D and 4G). Although this architecture appears similar to a form of multi-attractor model that has a continuum of attractors—a bump of activity that (in the absence of stimuli) can be centered at any location (the so-called “ring attractor model”;,, )—our model was actually quite different. While multi-attractor networks show a bump of sustained activity even once the stimulus is removed (leaving only non-specific background excitation), in our network the bump of activity depends on the similar bump of stimulus-driven input. When the stimulus is removed, our network returns within a single membrane time constant to a homogeneous level of baseline activity, driven by the homogeneous baseline input (Figure 4D). As we show below, this dynamical regime is also characterized by fundamentally different patterns of response variability than multi-attractor dynamics.

24We applied this model to study the stimulus dependence of variability quenching in V1 and compared our results to a new analysis we performed of previously published recordings in V1 of the awake monkey. In the absence of visual input (0% contrast), the network exhibited spatially patterned fluctuations in momentary firing rates around a few Hz (Figure 4D) with large across-trial variability in single-cell (Figure 4E). In evoked conditions, the input drove a hill of network activity around the stimulus orientation as in the data (Figures 4D and 4G), resulting in approximately contrast-invariant tuning curves. At stimulus offset, activity rapidly decayed back to spontaneous levels with the cellular time constant (Figure 4D), as observed in cortex when thalamic input is silenced.

25The fluctuating firing rates in spontaneous activity implied super-Poisson variability in spike counts—Fano factors greater than 1 (Figure 4F, top) —given the stochastic spiking mechanism described above (Figure 3). This was consistent with the high level of spontaneous variability in the data (Figure 4F, bottom). Both the model and the data exhibited a pronounced drop in Fano factor following stimulus onset (Figure 4F) and displayed a U-shaped tuning of variability suppression with stimulus orientation (Figure 4G, middle), such that variability suppression was stronger for cells whose preferred orientation was close to the stimulus. The model made similar predictions for variability in membrane potentials: a U-shaped profile of variance suppression in stimulus-evoked conditions relative to spontaneous fluctuations (Figure 4G, bottom).

26Notably, for similar reasons as in the randomly connected network (Figure 3; Figure S5), it was primarily the shared and not the private part of variability that was quenched by stimuli in the model (Figure 4H, left), and this required some degree of spatial correlations in the input noise (Figure S6). This was because the spatially smooth nature of the connectivity meant that only spatially smooth patterns of activity were strongly coupled to the network dynamics. A substantial suppression of shared variability at stimulus onset has been observed across many cortical areas as well as in our analysis of the V1 data (Figure 4H, right;; see also ).

27We again explored a broad range of parameters to show that the tuning of variability suppression was a robust outcome of the model. We found that Fano factor and variance were always most strongly suppressed in the neurons that were most strongly driven by the stimulus (the “dip” of the U shape) consistent with the V1 data (see above). Interestingly, there were some cases when neurons tuned to the opposite stimulus also showed a strong reduction of Fano factor (though not of membrane potential variance; Figure S7)—consistent with recent findings of an M-shaped modulation of Fano factors (and spike count correlations of similarly tuned cells) in area MT of the awake macaque (Figure S7; ). However, while such an M-shaped modulation was previously attributed to marginally stable multi-attractor dynamics, our model still produced this with a single stable attractor: the spike count variability of oppositely tuned cells dropped when input tuning in the model was as narrow as, or narrower than, the tuning of recurrent connections. In this configuration, oppositely tuned cells received so small a net input on average that their membrane potential fluctuations barely crossed the threshold of the firing rate nonlinearity, thus producing very little spiking variability. In turn, this loss of firing rate variance even overcame the effect of dividing by very small firing rates in computing Fano factors for these neurons. Under the same conditions, a similar M shape was apparent for spike count correlations between similarly tuned neurons, as a function of their (common) preferred orientation (Figure S7).

28Next, we analyzed the origin and mechanism of the stimulus-tuning of noise variability in the ring architecture. As mentioned above, for a fixed stimulus, the most prominent feature of population activity was a “bump” of high in the cells with preferred orientations near the stimulus orientation and a lower baseline of activity in the surround (Figure 5A, left and middle). In general, variability in the bump and the baseline captured most of the network’s variance and its suppression with increasing stimulus strength (Figure S8). Here and in the next section we specifically focus on the structure of the quenched noise variability after stimulus onset.

29After stimulus onset, most of the shared variability (87%; Figure S8) arose from variability in the location, μ, and width, σ, of the bump of activity (Figure 5A, middle and right). Notably, fluctuations in bump amplitude and width scaled inversely with one another, as the nonlinear interactions among neurons in our network resulted in strong normalization, preserving overall activity. Each of these small transformations resulted in a characteristic pattern of momentary deviation of network activity from the mean bump (Figure 5B). In turn, these two patterns of momentary fluctuations (Figure 5C) contributed two distinct spatial covariance templates (Figure 5D). For example, sideways motion of the bump increased the firing rates of all the cells with preferred orientations on one side of the stimulus orientation and decreased firing rates for all cells on the other side (Figure 5B, top). This resulted in positive covariances between cells with preferred orientations on the same side of the stimulus orientation and negative covariances for cells on opposite sides (Figure 5D, top: μ-template; ). Conversely, an increase in bump width (and thus a decrease in amplitude) increased the activities of cells on the flanks of the bump, tuned away from the stimulus, while decreasing the activity of cells near the peak, tuned for the stimulus (Figure 5B, bottom). This generated positive covariances within each of these groups and negative covariances between the two groups (Figure 5D, bottom: σ-template).

30Taken together, the ongoing jitter in bump location and width contributed a highly structured pattern of response covariances, which accounted for most of the structure in the full covariance matrix of the network (Figure 5D, compare “bump kinetics” with “full”). In particular, bump kinetics correctly predicted the variances of cells (given by the diagonal of the full covariance matrix indicated by the filled arrow in Figure 5D), showing less variance for cells tuned to the stimulus orientation of than for cells tuned to orthogonal orientations (see Figure 4G, bottom, green), and hence explained the U-shaped modulation of Fano factors (Figure 4G, middle, green). Moreover, the recurrent dynamics generated negative correlations in the fluctuations of cells with orthogonal tuning, despite such pairs receiving positively correlated inputs (Figure 5D, “input” versus “bump kinetics,” secondary diagonal with open arrow).

31For a direct comparison of the dynamical regime of the SSN with previously proposed mechanisms for variability modulation, based on marginally stable or chaotic dynamics, we first studied the predictions of the models for the spatial patterns of spike count noise correlations. Chaotic models have not, and probably can not, predict the tuning of mean responses, let alone that of variability suppression, so we focused on a comparison with a multi-attractor ring model. This model has been suggested to account for stimulus-modulated changes in variability in area MT. We matched it to our model such that it produced similar tuning curves and overall levels of variability (Figure S9).

32While there were several differences apparent in the detailed correlations predicted by the two models (Figures 6A and 6C), many of these could be explained away by trivial factors that neither model captured fully. For example, the average correlation was substantially larger in the SSN than in the attractor network—but this difference could be eliminated by invoking, in the attractor model, an additional (potentially extrinsic) mechanism that adds a single source of shared variability across neurons, resulting in a uniform (possibly stimulus strength-dependent) positive offset to all correlations. As another example, the attractor network always exhibited an M-shaped modulation of correlations, whereas, just as for Fano factors (see above), the SSN mostly showed a U-shaped modulation but could show an M shape for particular parameters (Figure S7).

33Therefore, we focused on distinctions that were robust to model details and followed from a fundamental difference of bump kinetics in the two models: in contrast to the richer patterns of variability generated by the SSN, multi-attractor dynamics showed a more limited repertoire, dominated by sideways motion of the bump with barely any fluctuations in bump width (Figure S9; ). As fluctuations in bump location and width had opposite effects on the correlations between orthogonally tuned cells in the SSN model (Figure 5D insets, gray), their cancellation made these correlations only very weakly modulated by the stimulus (Figure 6A, gray arrow; Figure 6B, gray). In particular, this modulation was much shallower than that for similarly tuned cells (Figure 6A, orange arrow; Figure 6B, orange), for which variability in bump location and width had congruent effects (Figure 5D insets, orange) that added to rather than cancelled each other. In contrast, in the attractor model, there was no such cancellation even for orthogonally tuned cells due to the absence of fluctuations in bump width (Figure S9). This meant that correlations between orthogonally tuned cells were just as deeply modulated as those between similarly tuned cells (Figures 6C and 6D).

34Previous reports on the stimulus-tuning of noise correlations examined only similarly tuned cells and reported mostly M-shaped modulation, which does not distinguish between the models. Therefore, we conducted our own analyses of a previously published dataset of V1 responses in the awake monkey (Figure 6E). The modulation of these correlations by the stimulus could only be accounted for by the SSN. First, we found that correlations between similarly tuned cells were significantly modulated by the stimulus (Figure 6E, orange; repeated-measures ANOVA, ), and this modulation had a U rather than an M shape. More critically, also in agreement with the predictions of the SSN but not of the attractor model, correlations between orthogonally tuned cells were unaffected by the stimulus (Figure 6E, gray; repeated-measures ANOVA, ). While the magnitude of correlations in either model was overall larger than in the data, this simply reflected the relatively small number of neurons in the models (model correlations could be decreased without affecting the shape and extent of their stimulus tuning by substituting each model unit by several neurons with independent spiking noise).

35We hypothesized that the fundamentally different mechanisms responsible for variability modulation in the SSN, the multi-attractor, and the chaotic dynamical regimes (Figure 1) should be revealed in the dynamics of variability suppression at stimulus onset and of variability recovery at stimulus offset. In order to test this, we used the same models for the SSN and multi-attractor models as above, and we implemented the classical chaotic model of (STAR Methods), in which variability suppression had previously been shown to occur. We then measured the across-trial variability (averaged across neurons) following the onset and offset of a step stimulus in each model (Figures 7A–7C, shaded areas), as we parametrically varied the amplitude of the stimulus and therefore the degree of variability suppression (Figures 7A–7C, dark to light colors).

36In the SSN, the timescales on which both suppression and recovery of variability occurred were nearly as fast as the single-neuron time constant (20 ms in these simulations; Figures 7D and 7E, green). In contrast, in chaotic networks, both these timescales were several (4–15) times longer than the single-neuron membrane time constant (Figures 7D and 7E, blue). More importantly, recovery times were much longer than suppression times in the chaotic network and increased with increasing stimulus strength and thus increasing amount of variability suppression during the stimulus period, neither of which was the case in the SSN. In the multi-attractor network, both the dynamics of the network activity and those of variability were much slower than in the SSN (Figures 7D and 7E, red). Moreover, we found that, unlike in the SSN or the chaotic model, variability increased transiently immediately following stimulus onset (before eventually decreasing to its new steady state; Figure 7C). The cause of this behavior was the slow rotation of the activity bump from its random position at the time of stimulus onset to the location where cells’ preferred orientation matched the stimulus orientation (Figure S10). Thus, we expect this behavior to be generic at least to the subclass of multi-attractor models that have a continuous ring of attractors and thus show such rotational response, which likely include those that can address the orientation- or direction-tuning of variability reduction in V1 and MT.

37The timescales of variability suppression and recovery found experimentally in anaesthetized cat V1 and awake monkey MT (Figures 7D and 7E, open square and circle; ) and by our own analysis of awake monkey V1 data (Figures 7D and 7E, dotted square; ) were short and nearly identical. Moreover, recovery times showed little dependence on the amount of variability suppression (comparing across areas), and there was no transient increase in variability at stimulus onset (Figure 4F; ). These results confirm the predictions of the SSN and are at odds with the dynamics of variability modulation as predicted by the multi-attractor and chaotic regimes.

38We studied the modulation of variability in a stochastic, nonlinear model of cortical circuit dynamics. We focused on a simple circuit motif that captured the essence of cortical networks: noisy excitatory and inhibitory populations interacting in a recurrent but stable way despite expansive single-neuron nonlinearities. This stochastic stabilized supralinear network (SSN) reproduced key aspects of variability in the cortex. During spontaneous activity, i.e., for weak external inputs, model neurons showed large and relatively slow synchronous fluctuations in their membrane potentials. These fluctuations were considerably amplified by the network relative to that expected from the input alone and were quickly quenched and decorrelated by stimuli. The model thus explains and unifies a large body of experimental observations made in diverse systems under various conditions. Moreover, the drop in variability was tuned to specific stimulus features in a model of V1/MT, also capturing recent experimental findings as well as our own analyses of a previously published dataset.

39The main insight of our analysis was that in a network of nonlinear neurons with an expansive firing rate nonlinearity, increasing the input increases the effective connection strengths of the network, which in turn modulates the variability of responses. We identified two opposing effects of increasing effective connectivity on variability: the amplification of variability by excitatory-inhibitory interactions (balanced amplification), which dominates at very low (spontaneous) levels of input, and the quenching of variability by increased inhibitory feedback, which dominates for stimulus-driven input. Critically, these network effects preferentially act on smooth patterns of activity that are aligned with the anatomical connectivity of the network, so that it is the shared component of variability that is suppressed and modulated by stimuli. Taken together, we showed that these mechanisms robustly produced experimentally observed spatial and temporal patterns of variability quenching and modulation, whereas the dynamics of the network always remained in the vicinity of a single attractor state, unlike previously proposed mechanisms based on multi-attractor or chaotic dynamical regimes.

40We focused on how the network shapes variability and assumed that the variability originates in correlated noise input to the network; such input correlations could arise due to upstream areas already exhibiting noise correlations (e.g., thalamic input to V1, ) and/or because of feedforward connectivity implying shared inputs (e.g., ). In contrast, other models have focused on how circuits intrinsically generate slow correlated variability. Nevertheless, our model also points to an important mechanism for creating shared variability, namely the strong amplification of the input noise by balanced amplification (see also,, ).

41Although most of our analyses were based on rates, rather than spikes, the effect of fast fluctuations resulting from spiking noise were not ignored, but were incorporated implicitly in the power-law input/output nonlinearity of neurons in the model (Equation 3) and in the stochastic spike-generation mechanism used in our spiking network simulations (Figure 3, STAR Methods). Theoretical work shows that these fast fluctuations are the key factor causing momentary firing rates (on the 30–50 ms timescale of fluctuations considered here) to be a supralinear, power-law function of mean voltages, a critical feature of our model. As experiments, as well as our model, show that only the shared but not the private part of variability is modulated by stimuli, we expect our assumption that the exponent of the threshold power-law nonlinearity can be considered constant (implying that fast private spiking fluctuations are not affected by stimuli) to be valid to a good approximation. We also expect that a more detailed model explicitly including these fast fluctuations would allow a more systematic study of the effects of stimuli on high-frequency (gamma) oscillations, which our current model could only partially account for (Figure 3B).

42While we focused on the sources and modulation of slower, correlated fluctuations, a classical model of cortical variability, the “balanced network”, focused on the origin of fast fluctuations from spiking noise. In that model, very large external and recurrent inputs cancel or “balance” to yield a much smaller net input. This mechanism can self-consistently generate the voltage variability to generate irregular spiking. However, the very strong, very fast inhibitory feedback in the balanced network suppresses correlated rate fluctuations away from the stable state, leaving only fast, private variability due to irregular spiking (though “breaking balance” can restore correlated variability;, ). Because the shared variability is already eliminated, stimuli cannot modulate that variability.

43As opposed to the “tight balance” between excitation and inhibition in the classical balanced network model, the SSN in the stimulus-driven regime is “loosely balanced”: the same mathematical cancellation of external and recurrent input occurs, but in a regime in which inputs are not large and the net input after cancellation is comparable in size to the factors that cancel. This regime is supported by observations that external input is comparable to, rather than very much larger than, the net input received by cortical cells. This loose balance allows correlated variability to persist and be modulated by stimuli. Variability quenching in the stochastic SSN robustly occurred as the input pushed the dynamics to stronger and stronger inhibitory dominance. Consistent with this, with increasing strength of external input, the ratio of inhibition to recurrent excitation received by neurons in the network increases, as observed in layers 2/3 of mouse S1 and V1. In the balanced network, the ratio of inhibitory to excitatory activity would be fixed regardless of the strength of activation. The balanced network also only yields responses that are linear functions of the input (though see ), whereas the loosely balanced regime replicates many nonlinear cortical response properties, including the profound dependence of correlated variability on stimuli. Although our model does not focus on the origins of fast spiking variability, spiking models in the loosely balanced SSN regime can, given noisy inputs (e.g., ), yield the irregular spiking characteristic of cortex (unpublished data).

44We analyzed variability modulation solely as arising from intrinsic network interactions, but other factors may also contribute. External inputs may be modulated; for example, the drop with contrast in Fano factors in the lateral geniculate nucleus (LGN) has been argued to underlie variability decreases in V1 simple cells (; but see ). However, since high-contrast stimuli also cause firing rates to increase in LGN, the total variance of LGN-to-V1 inputs (scaling with the product of the LGN Fano factor and mean rate) is modulated far less by contrast. This provides some justification for our model choice that input variance did not scale with contrast. Changes in input correlations have also been suggested as a potential mechanism underlying variability modulation. However, the proposed mechanism would require a stimulus to specifically increase the correlations of the different inputs onto individual cells (and this increase should be tuned to the stimulus) while leaving the correlation of inputs between cells unchanged. This seems difficult to achieve in cortex, where nearby cells are likely to share a significant amount of input and correlations are generally observed to decrease, rather than increase, with stimulus strength.

45One particular form of external input modulation, that involving changes in brain state, has been proposed to directly contribute to correlated variability in both awake and anesthetized cortex, so that a reduction of state switching would underlie the reduction of shared variability. To the extent that correlated noise in the input to our model is aligned with a uniform activity pattern, this input can also be regarded as having a single scalar “brain state”-like component that is changing in time. However, our analysis suggests that the variability of this component needs not be modulated directly by the stimulus to account for variability quenching in network responses. Instead, our network used its intrinsic mechanisms to quench variability in response to a stimulus. Importantly, these intrinsic mechanisms not only quenched this uniform component of variability (Figure S8), but also produced more complex patterns of variability modulation via “bump” kinetics that a single brain state-dependent mechanism could not account for.

46Cellular factors may also modulate variability. For example, inhibitory reversal potential or spike threshold may set boundaries limiting voltage fluctuations, which would more strongly limit voltage fluctuations in more hyperpolarized or more depolarized states, respectively; conductance increases will reduce voltage fluctuations; and dendritic spikes may contribute more to voltage fluctuations in some states than others. A joint treatment of external input, cellular, and recurrent effects may be needed to explain, for example, why variability appears strongest near the preferred stimulus in anaesthetized cat V1 or why overall variability grows with visual stimulation in some neurons of awake macaque V1.

47Cellular properties may themselves be subject to change over time, thereby causing changes in variability. For example, various mechanisms (e.g., attention, intrinsic and synaptic plasticity, neuromodulators, anesthetics) can change the input/output gain of single neurons and the synaptic efficacies of the network. As all these changes eventually lead to changes in effective connectivity, our work offers a principled approach to study their effects on variability and is thus complementary to previous studies that focused on the consequences of different anatomical connectivity patterns on correlations.

48The nonlinear response properties of our network were crucial for the modulation of variability by stimuli. These nonlinearities had been shown to capture ubiquitous phenomena involving nonlinear response summation to multiple stimuli, including normalization, surround suppression, and their dependencies on stimulus contrast. As such, the SSN reproduces much of the phenomenology of the “normalization model” of cortical responses and provides a circuit substrate for it.

49However, while response normalization has previously been studied for deterministic steady-state responses, our results can be interpreted as showing that it also plays a role in the suppression of ongoing variability by stimuli, as well as shaping the structure of stimulus-evoked noise correlations. Specifically, in the deterministic SSN, steady-state responses to multiple stimuli add sublinearly, and as one stimulus becomes stronger than another, the response to their simultaneous presentation becomes “winner take all,” i.e., dominated by the response to the stronger stimulus alone. This provides an alternative conceptual explanation of why, in the stochastic SSN, a stronger mean input drive relative to the noise input leads to greater suppression of the noise’s contribution to the total network response, thus quenching variability.

50Our results go beyond what could be predicted based on this simple qualitative link between steady-state normalization and variability quenching. First, we found a specific quantitative form of normalization in our network: an approximate conservation of the integrated activity across a bump of activity that forms around cells tuned to the stimulus orientation, despite fluctuations in its width. In turn, this predicted a specific pattern of noise correlations that we found contributed substantially to noise variability in V1 of the awake monkey (Figure 6). Second, we were able to study the dynamics with which variability was suppressed following stimulus onset and recovered following stimulus offset and found a good match to experimental data (Figure 7).

51We found that an increase in inhibitory dominance was necessary for the suppression of variability and correlations in the SSN. In line with that, studied rodent A1 and V1 in various awake and anesthetized brain states and found that desynchronized states with weaker correlations were accompanied by enhanced activity of putative fast-spiking inhibitory neurons. By fitting a recurrent spiking E-I network model to the data, they found that enhanced inhibitory feedback was the key factor capturing the suppression of correlations. However, the enhanced dominance of inhibition with increasing network activation, which suppresses correlations, was artificially incorporated into the model by making the inhibitory conductance an exponential function of the inhibitory spike count. In contrast, our model provides a dynamical mechanism by which inhibition becomes increasingly dominant with increasing network activation.

52proposed a mechanism similar to ours for the top-down suppression of correlated variability by attention, rather than bottom-up suppression by a stimulus. They also proposed that this arises from enhanced inhibitory feedback resulting from increased effective connectivity due to expansive input/output functions. However, their conclusions differed significantly from ours. They found that, for attention to suppress variability, attentional input had to be directed dominantly to inhibitory cells, while for attention to increase the gain of response to stimuli, stimuli had to give input dominantly to excitatory cells. Note that this implies that stimuli would not suppress variability. We have found that neither of these conditions are necessary (Methods S4) and that stimuli robustly suppress variability. In particular, increasing input strength decreased variability across a wide range of relative strengths of input to excitatory versus inhibitory cells (Figure S2). The main reason for these differences in conclusions is the special, non-generic parametrization of the model studied by in which a neuron’s projections to excitatory and to inhibitory neurons were statistically identical, which precluded the SSN regime (Methods S2).

53Two proposals have been made previously to explain quenching of variability by a stimulus: a stimulus may quench multi-attractor dynamics to create single-attractor dynamics, and a stimulus may quench chaotic dynamics to produce non-chaotic dynamics. Our results propose a very different dynamical regime underlying variability quenching, which can be distinguished from the multi-attractor or chaos-suppression models.

54Conceptually, the stochastic SSN differs from previous models of stimulus-driven quenching of shared variability in exhibiting a single stable state in all conditions—spontaneous, weakly driven, strongly driven—whereas the others show this only when strongly driven. Furthermore, quenching of variability and correlations in the SSN is highly robust, arising from two basic properties of cortical circuits: inhibitory stabilization of strong excitatory feedback and supralinear input/output functions in single neurons. In contrast, models of multi-attractor or chaotic dynamics can either account only for the modulation of average pairwise correlations or else require considerable fine tuning of connections to account for more detailed correlation patterns. Moreover, as studied thus far (,,; but see,, ), they typically ignore Dale’s law (the separation of E and I neurons) and its consequences for variability, e.g., balanced amplification. These differences between the SSN and previous models also lead to two main experimentally testable features that we used to distinguish their respective dynamical regimes: the tuning and the timing of variability modulation.

55With respect to the stimulus tuning of spike count Fano factors and noise correlations, we found that multi-attractor networks could only predict an M-shaped modulation while the SSN could produce either M- or U-shaped modulations depending on the tuning width of inputs relative to that of connectivity. Indeed, while most types of stimuli in MT were found to result in an M-shaped modulation, coherent plaids and random moving dots in the macaque as well as moving dot fields and drifting gratings in the marmoset (Selina Solomon, personal communication; ) result in a pronounced U-shaped modulation of Fano factors in MT, and our own analyses of V1 data also revealed a U-shaped modulation. Interestingly, our results also suggested that irrespective of the precise shape of the modulation of spike count statistics, membrane potential variability in the SSN should always exhibit a U-shaped profile (Figure 4), which could be tested in future experiments. Critically, we also identified a rarely analyzed aspect of spatial correlation patterns that could most clearly distinguish between different models: the modulation of correlations between orthogonally tuned cells. The SSN predicted only very weak modulation for such cell pairs, while multi-attractor dynamics resulted in modulations that were as deep as for pairs of similarly tuned cells. We found that data from awake macaque V1 supported the SSN.

56Another distinctive feature of the SSN regime is the speed of its dynamics, and in particular the speed with which variability is modulated as the stimulus is changed. In contrast to multi-attractor and chaotic dynamics, in which variability modulation happens on timescales that are considerably slower than the single neuron time constant, the SSN produces fast variability modulation on a timescale comparable to the neural time constant. The timescales of variability modulation we extracted from data recorded in monkey visual cortical areas were fast, on the order of 20–50 ms, providing further support to the SSN.

57In summary, the SSN robustly captures multiple aspects of stimulus modulation of correlated variability and suggests a dynamical regime that uniquely captures a wide array of behaviors of sensory cortex. In turn, our work suggests a principled approach to use data on cortical variability to identify the dynamical regime in which the cortex operates.

58As Lead Contact, Guillaume Hennequin is responsible for all reagent and resource requests. Please contact Guillaume Hennequin at g.hennequin@eng.cam.ac.uk with requests and inquiries.

59The values of all the parameters mentioned below are listed in the tables below. All differential equations detailed below were integrated using a simple Euler scheme with time step 0.1 ms.

60Our rate-based networks contained excitatory and inhibitory units, yielding a total units. The circuit dynamics were governed by (see also Methods S1):where denotes the of neuron i, is its membrane time constant, is a resting potential, is the (positive or zero) strength of the synaptic connection from neuron j to neuron i, and is the potentially time-varying but deterministic component (the mean) of external input to which a noise term is added (see below, “Input noise”). The momentary firing rate of cell j was given by a threshold-powerlaw function of its membrane potential:

61Experiments support Equation 3 when both membrane potentials and spike counts are averaged in 30 ms time bins. Accordingly, in Equation 2 can be understood as the coarse-grained (low-pass filtered) version of the raw somatic membrane potential; in particular it does not incorporate the action potentials themselves. Thus the effective time resolution of our model was around 30 ms which allowed studying the effects of inputs that did not change significantly on timescales shorter than that. Accordingly, in Equation 2 we assumed that external noise had a time constant ms, in line with membrane potential and spike count autocorrelation timescales found across the cortex.

62Equations 2 and 3 together define the stabilized supralinear network model studied in and, but formulated with voltages rather than rates as the dynamical variables (the two formulations are mathematically equivalent when all neurons have the same time constant, ) and with the crucial addition of noise that enables us to study variability. In all the figures of the main text, the exponent of the power-law nonlinearity was set to (but see Figure S2 for n > 2). Methods S2 explores more general scenarios.

63In the reduced rate model of Figure 2, each unit received the same constant mean input h. In the ring model, the mean input to neuron i was the sum of two components,

64The first term mV is a constant baseline which drives spontaneous activity. The second term models the presence of a stimulus with orientation in the visual field as a circular-Gaussian input bump of “half width” centered around and scaled by a factor c (increasing c represents increasing stimulus contrast), taking values from 0 to 1, times a maximum amplitude. We assumed for simplicity that E and I cells are driven equally strongly by the stimulus, though this could be relaxed.

65The input noise term in Equation 2 was modeled as a multivariate Ornstein-Uhlenbeck process:where is a collection of N independent Wiener processes and is an input covariance matrix (see below). Note that Equation 5 implies.

66In the reduced two-population model (Figure 2), noise terms were chosen to be uncorrelated, i.e., (where if and 0 otherwise), is the E/I type of neuron i, and is the variance of noise fed to population α (see Equation 7 below). In the spiking two-population model (Figure 3), input noise covariance was uniform, such that, with the pairwise correlation coefficient set to (see Figure S5 for the dependence of our results on ρ). In the ring model (Figures 4, 5, 6, and 7), the noise had spatial structure, with correlations among neurons decreasing with the difference in their preferred directions following a circular-Gaussian:where and are the preferred orientations of neurons i and j (exc. or inh.), and is the correlation length (see table “Parameters Used in the SSN Simulations”). The noise amplitude has the natural scalingsuch that, in the absence of recurrent connectivity, the input noise alone would drive fluctuations of standard deviation or, measured in mV, in the E or I cells, respectively. We chose values of that yielded spontaneous Fano factors in the range 1.3-1.5 where appropriate, and chose to make up for the difference in membrane time constants between E and I cells (see table “Parameters Used in the SSN Simulations”).

67The synaptic weight matrix in the reduced model was given bywhere is the magnitude of the connection from the unit of type B (E or I) to that of type A (see table “Parameters Used in the SSN Simulations” for parameter values). In the ring model, connectivity fell off with angular distance on the ring, following a circular-Gaussian profile:where and are the preferred orientations of neurons i and j (exc. or inh.), and sets the length scale over which synaptic weights decay (see table “Parameters Used in the SSN Simulations”). The connectivity matrix was further rescaled in each row and in each quadrant, such that the sum of incoming E and I weights onto each E and I neuron (4 cases) matched the values of,, and in the reduced model. Thus, all connectivity matrices used in the SSN model obeyed Dale’s law.

68To relate the firing rate model to spiking data in Figures 4 and 6, we assumed that action potentials were emitted as inhomogeneous (doubly stochastic) Poisson processes with time-varying rate given by Equation 3. Unlike in the full spiking model (see below), spikes did not “re-enter” the dynamics of Equation 2, according to which neurons influence each other through their firing rates. Spikes were counted in 100 ms time bins and spike count statistics such as Fano factors and pairwise correlations were computed as standard.

69In the spiking model (Figure 3), neuron i emitted spikes stochastically with an instantaneous probability equal to, with time-varying rate given by Equation 3, consistent with how (hypothetical) spikes were modeled in the rate-based case (cf. above). Presynaptic spikes were filtered by synaptic dynamics into exponentially decaying postsynaptic currents (E or I):where the ’s are the firing times of neuron j, ms is the synaptic time constant, and ms is a small axonal transmission delay (which enabled the distribution of the simulations onto multiple compute cores; ). Synaptic currents then contributed to membrane potential dynamics according towhere the synaptic efficacies are described below, and the noise term was modeled exactly as described above.

70For each neuron i, we drew excitatory and inhibitory presynaptic partners, uniformly at random. Connection probabilities were set to and respectively. The corresponding synaptic weights took on values where denote the populations to which neuron i and j belong respectively, and are the connections in the reduced model (see table “Parameters Used in the SSN Simulations”). This choice was such that, for a given set of mean firing rates in the E and I populations, average E and I synaptic inputs to E and I cells matched the corresponding recurrent inputs in the rate-based model. Synapses that were not drawn were obviously set to.

71As a proxy for LFP in Figure 3, we took the momentary population-averaged ( simulated various proxies and, although some proxies were more accurate, they found the average to be reasonably accurate).

72We compared our ring SSN model to a version of the ring attractor model published by. The ring attractor model had a single population with a similar ring topology, and—using the same notation as above—the connectivity took the form (cf. Equation 9)where is the number of neurons, and and are two parameters that control the average connection strength and modulation with tuning dissimilarity, respectively. Note that, in general, this connectivity matrix could violate Dale’s law but with the specific parameters used here it did not (see table “Parameters Used in the Multi-attractor Network Simulations”). Instead, all connections were inhibitory to keep the system in the marginally stable regime (as in ). The dynamics of the network obeyed a similar stochastic differential equation as for the ring SSN (Equation 2), namelywith the momentary firing rate of cell j given by a rectified saturating firing rate nonlinearity (cf. Equation 3):and a noise process identical to the one we used in the SSN (same spatial and temporal correlations, Equations 5 and 6), with a variance adjusted so as to obtain Fano factors of about 1.5 during spontaneous activity (Figure S9B, black). The external input had both a constant baseline, b, and a contrast-dependent modulated component (cf. Equation 4):where A controlled the depth of the modulation, and c represents stimulus strength.

73Note that although the phenomenology and dynamical regime of this model was consistent with that of (Figure S9), the model differed from their original implementation in some of the details: our dynamics were written in voltage form, not in rate form, we had only one unit at each location on the ring (as opposed to small pools of neurons), and our input noise process had spatial correlations to allow for a more direct and consistent comparison with the ring SSN.

74Our analysis of variability in this ring attractor network is presented in Figure S9 in a format identical to that of Figure 5, and shows that shared variability is entirely dominated by the fluctuations in the location of an otherwise very stable bump of activity.

75We also implemented a chaotic rate network of size with the following (deterministic) dynamics (cf. Equations 2 and 13):with an (unrectified) saturating firing rate nonlinearity (cf. Equations 3 and 14)(which could thus go negative as well as positive). Elements of the synaptic weight matrix were sampled i.i.d. from a normal distribution (thus violating Dale’s law, cf. Equations 9 and 12):with, which placed the network in the chaotic regime. The external input was a constant input vector of the form (cf. Equations 4 and 15)where is a phase sampled i.i.d. from a uniform distribution between 0 and, and c represents stimulus strength. See table “Parameters Used in the Chaotic Network Simulations” for all parameter values. As shown in, chaos is suppressed for large enough c.

76We analyzed neural recordings from the V1 of two awake monkeys (Figures 4, 6, and 7). A full description of the experimental protocol and recordings can be found in the original publication. We discarded all cells that were poorly isolated (contamination %), leaving us with 330 cells to analyze. The stimuli consisted of moving gratings of various orientations, all at 100% contrast. We fitted orientation tuning curves (Figure S11; average firing rate in the first 500 ms following stimulus onset, as a function of stimulus orientation) of the form, where θ is the stimulus orientation (thus, we neglected the direction of motion, which could be in either of the two directions orthogonal to the orientation of the grating). The fit was achieved using nonlinear least-squares regression.

77For each neuron, we calculated an orientation tuning index (OTI), defined based on the fitted tuning curve aswhere. As the ring architecture we studied in Figures 4, 5, 6, and 7 only applied to neurons with well-defined tuning curves, we excluded cells that had as well as average evoked rates (measured during the stimulus period) below 1 spike/sec. This left us with 99 well-tuned cells to analyze.

78Our analysis of the stimulus tuning of Fano factors and pairwise spike-count correlations was based on a time window of 100 ms starting at stimulus onset.

79We performed factor analysis of spike counts, either for a single stimulus condition in the model (the model had a natural rotational symmetry), or separately for each stimulus condition (direction) in the V1 dataset, subsequently averaging the reported quantities across conditions. We worked with normalized spike counts, defined as where is the spike count of neuron i in trial k and denotes averaging across trials. Note that the variances of these normalized spike counts are exactly the Fano factors, i.e., the usual measure of spike count variability. This prevented the normalized spike count covariance matrix from being contaminated by a rank-1 pattern of network covariance merely reflecting the tuning of single-neuron firing rates (the “Poisson” part of variability, which tends to scale with the mean count). Factor analysis decomposes as the sum of a rank-k covariance matrix representing k modes of network covariances, and a diagonal matrix. In the rate model, we could near-perfectly estimate the spike count covariance matrix, so we performed factor analysis by direct eigendecomposition of, thus defining whereby the top k eigenvectors of contributed to shared variability in proportion of the corresponding eigenvalues. For factor analysis of the monkey V1 data, we performed direct maximization of the data likelihood, also keeping k factors. In Figure 4, we set, but we observed quenching of shared variability irrespective of k (Figure S12).

80The code used for model simulations and data analysis is available from the Lead Contact, Dr Guillaume Hennequin, upon request.