Introduction {#s1}
============

Sensory information is represented in the activity of populations of neurons, but the responses of individual neurons within these populations are not uniquely determined by external stimuli: repeated presentations of the same stimulus elicit different spike trains. Although some of this variability presumably arises from noise in local circuits, a substantial portion appears to be due to fluctuations in neurons' excitability, i.e. their gain ([@bib21]). Since this may result from changes in internal states that modulate sensory responses, such as wakefulness, reward, and expectations, or from changes in top-down signals arriving from other cortical areas, it is likely that a significant component of the gain fluctuations in individual neurons are not "private", but shared across populations ([@bib16]; [@bib60]; [@bib33]). Shared gain fluctuations could have a crucial impact on sensory computations, depending on their structure ([@bib44]). Yet this structure remains largely unknown.

Visual attention provides a well-known example of a process that affects neural gain. When monkeys direct their attention to a visual stimulus, the mean firing rates of specific subsets of visual neurons in striate ([@bib39]; [@bib27]) and extrastriate cortex ([@bib43]; [@bib69]; [@bib70]; [@bib37]; [@bib76]; [@bib8]; [@bib42]) have been found to increase. This motivates the hypothesis that, at the level of the population, attention acts to increase the gain of selected neurons, thereby increasing the signal-to-noise ratio of sensory representations, and hence improving performance on perceptual tasks ([@bib39]; [@bib69]; [@bib8]; [@bib32]; [@bib57]).

This "classical" view of attention has been augmented by recent observations that spatial attention affects more than mean response: in particular, attention reduces normalized measures of spike count variance, as well as stimulus-conditioned spike count correlations (i.e. "noise correlations") across pairs of neurons ([@bib8]; [@bib42]; [@bib26]). These measurements (and similar findings in the context of perceptual learning ([@bib23]; [@bib28]), cognitive challenge ([@bib59]), task engagement ([@bib14]), or wakefulness ([@bib53]) seem to imply that attention does more than just increase neural gain; for instance, it might change the underlying connectivity of the network. But the origin of these effects is difficult to interpret, as changes in pairwise correlations can arise from changes in direct or indirect couplings between neurons, or changes in common modulatory input ([@bib21]; [@bib16]; [@bib4]; [@bib78]). In particular, several groups have suggested that variability in the attentional signal itself might contribute to spike count variance and correlation ([@bib21]; [@bib9]; [@bib25]; [@bib15]).

Here, we develop a functional model for the population activity that accounts for all the above observations. The model includes stochastic modulatory signals that alter the gains of targeted subsets of neurons. We fit the model to spiking data from populations of ∼100 neurons in visual area V4, simultaneously recorded from both hemispheres of a macaque monkey that was performing a change-detection task under directed spatial attention ([@bib8]). The model is fit without specification of the hemispheric origin of the neurons, or the extent to which they were influenced by the various gain factors. The resulting fitted model reveals that the population was predominantly influenced by two independent shared modulatory signals, each operating primarily within one hemisphere, and each targeting the neurons most relevant for the task. The statistics of each of these signals changed significantly depending on the attentional condition, with each modulator exhibiting a decrease in variance when the monkey was cued to attend to the stimuli in the corresponding hemifield. Together with the (classical) increases in mean response, these changes in the statistics of the shared modulatory signals account for the previously-reported decrease in neural variability and noise correlations. Finally, we show that the inferred modulatory signals are correlated with the monkey's behavioral performance on each trial, and are influenced by the reward received on the previous trial. The structure and statistics of shared modulatory fluctuations thus provide a parsimonious account of attentional effects on population coding and behavior. A preliminary account of these results was presented in [@bib54].

Results {#s2}
=======

We analyzed spiking data from populations of single- and multi-units recorded in visual area V4 of two macaques while they performed a visual detection task. The experiment and data are described in detail in [@bib8]. We briefly describe the relevant details here.

Monkeys fixated on a central mark while two oriented grating patches were flashed concurrently on opposite sides of the fixation point. These stimuli (the "standards") were shown for 200 ms, with a delay of 200--400 ms between presentations. After a variable-length sequence of repeated standards, the orientation of one of the two gratings (the "target") was changed. For a small fraction of trials ("catch trials"), no target appeared. The monkey had to detect the target, and to indicate this with a visual saccade to its location. Trials were grouped in blocks of 125, in which the target had an 80% probability of occurring on one side. The identity of this "cued side" was indicated to the monkey through a set of 10 instruction trials preceding each block.

While the monkey performed this task, the spiking responses of populations of neurons in visual area V4 were collected from two microelectrode arrays, one per hemisphere. A total of 50--130 units were concurrently recorded from the two hemispheres per day. Over 36 recording days, 3004 units were collected and deemed suitable for further analysis. We analyzed spiking responses from these populations to a total of ∼10^5^ presentations of the standard (i.e. non-target) stimuli. We excluded responses to the first standard stimulus in each sequence (since those responses exhibit atypical onset transients), and to any stimulus in which the monkey failed to maintain fixation. We examined the total spike count for each 200 ms stimulus period, offsetting the response window by 60 ms to account for stimulus-response latency (as in [@bib8]).

Fitting the shared modulator model {#s2-1}
----------------------------------

We used a computational model to explore the structure of neural population activity in the presence of spatial attention ([Figure 1](#fig1){ref-type="fig"}). We describe the stimulus-driven instantaneous firing rate of each neuron *n* over time as *f~n~*(*s*(*t*)), where the function *f~n~*(·) describes the mapping from stimulus, *s*(*t*), to that neuron's firing rate. In addition, we allow each neuron's gain to be affected by three signals: (1) the current attentional cue, *c*(*t*), (a binary signal); (2) a slowly-varying global drift, *d*(*t*); and (3) a set of shared time-varying modulators, *m*~*k*~(*t*). We assume that there are *K* such modulatory signals, indexed by *k*. The degree to which each of these three signals affects the gain of each neuron is specified by a set of coupling weights, {*u~n~, v~n~, w~n,k~*}, and the net instantaneous firing rate of the neuron is written as:

$$r_{n}\left( t \right) = f_{n}\left( s\left( t \right) \right) \cdot \exp\left\lbrack {u_{n} \cdot c\left( t \right) + v_{n} \cdot d\left( t \right) + \sum\limits_{k = 1}^{K}w_{n,k} \cdot m_{k}\left( t \right)} \right\rbrack$$10.7554/eLife.08998.003Figure 1.Diagram of the modulated population model.Shown are three neurons (yellow boxes), each with a firing rate that is a function of the stimulus multiplied by three time-varying gain signals: the binary attentional cue; a slow global drift; and a set of shared modulators. The influence of each of these signals on each neuron is determined by a coupling weight, indicated by the thickness of the blue and black lines. Only one shared modulator is shown in this schematic, but the model allows for more, each with its own coupling weights. The spike counts of each neuron are conditionally Poisson, given the firing rate.**DOI:** [http://dx.doi.org/10.7554/eLife.08998.003](10.7554/eLife.08998.003)10.7554/eLife.08998.004Figure 1---figure supplement 1.Example slow drifts in spike counts of four simultaneously-recorded units (two from each hemisphere), taken from a recorded population of 77 units.Each point is the average spike count observed over ten consecutive stimulus presentations. The blocked structure of the task (i.e. the alternating cue directions) is indicated with alternating colors. Thick lines indicate the portion of model-estimated firing rate due to the combination of stimulus drive, cue, and the slow drift signal (without the shared modulators).**DOI:** [http://dx.doi.org/10.7554/eLife.08998.004](10.7554/eLife.08998.004)10.7554/eLife.08998.005Figure 1---figure supplement 2.Signatures of modulatory (multiplicative) effects in the neural responses.The shared modulator model is based on the assumption that the stimulus-driven firing rates of sensory neurons are modulated (multiplied) by a set of additional inputs. Evidence for fluctuating modulation of sensory responses was recently provided in ([@bib21]). The statistical argument that these inputs are multiplicative (rather than additive) is based on a consideration of their effects as a function of firing rate: multiplicative noise has its greatest effects on response variance (and covariance) for stimuli evoking high firing rates, while additive noise has its greatest effects on response variance (and covariance) for stimuli evoking low firing rates. This is difficult to assess in the context of the V4 attentional dataset, since the standard stimuli were all identical, and responses to the targets were often interrupted by saccades. Restricting the analysis to target responses that were not interrupted is possible, but such conditioning on a behavioral state (and likely a modulator state) would complicate any interpretation. There was, however, some variability in evoked firing rates arising from adaptation. Many neurons showed a trend of slightly decreasing response to the sequence of standard stimuli, with an average total decrease of ∼5% in firing rate over the ten presentations after the first. This decrease was sufficiently small that including it in the analysis of the main text only marginally improved predictive log-likelihoods (as in [Figure 2a](#fig2){ref-type="fig"}), and did not qualitatively change any of the main results or conclusions. These small changes in firing rate over the stimulus sequence were nevertheless sufficient to examine the hypothesis that neural response variance was due a multiplicative noise source. Following the logic of ([@bib21]), if we assume a multiplicative noise source with variance *σ*^2^, the Fano factor (variance divided by the mean) should increase with firing rate, *μ*:$$FF_{mult}\left( \mu \right)\mspace{600mu} = \mspace{600mu}\frac{\mu\operatorname{}\sigma^{2}\;\mu^{2}}{\mu}\mspace{600mu} = \mspace{600mu} 1 + \sigma^{2}\;\mu$$On the other hand, for an additive noise source of variance *σ*^2^, the Fano factor should *decrease* with firing rate:$$FF_{add}\left( \mu \right)\mspace{600mu} = \mspace{600mu}\frac{\mu\operatorname{}\sigma^{2}}{\mu}\mspace{600mu} = \mspace{600mu} 1 + \frac{\sigma^{2}}{\mu}$$These expected trends are illustrated in panel (**a**). Panel (**b**) shows the mean value of these quantities, estimated for the cue-towards condition across all cells. The data are clearly consistent with a multiplicative noise source. A similar trend is observed in the cue-away condition, albeit with an overall lower mean rate, and higher Fano factor. This analysis assumes that the noise source has constant variance across stimulus conditions. There were some small, non-monotonic changes in the estimated shared modulators' variance over the stimulus sequence. Factoring these in does not change the direction of the predictions or the data shown here. A similar analysis can be performed on pairwise response statistics (as stimuli evoke higher mean rates, a multiplicative model predicts correlations will increase, while an additive model predicts correlations will decrease). But these predictions prove more sensitive to the assumption of stability in *σ*^2^, so we omit them here.**DOI:** [http://dx.doi.org/10.7554/eLife.08998.005](10.7554/eLife.08998.005)

The exponential acts to convert the weighted sum of the three signals, which may take on positive or negative values, into a product of three positive-valued modulatory quantities.

Importantly, only the stimulus *s*(*t*) and cue signal *c*(*t*) are known. The drift, *d*(*t*), and modulators, *m*~*k*~(*t*), as well as all of the coupling weights must be fit to the experimental data. To accomplish this, we assumed that the firing rate of each neuron, *r~n~*(*t*), was constant over the duration of each stimulus presentation (200 ms), and that the observed spike counts arose from a Poisson distribution with that rate. We then fit this probabilistic model to the spike count data by maximizing the posterior probability of the model parameters (see Materials and methods). We describe each of the component signals in turn.

The stimulus and cue signals encapsulate the external factors that are set by the experimenter and available to the monkey. Since we only analyzed responses to the standard stimuli--- which were identical for all trials of a given day---the stimulus-dependent drive in this experiment, f~*n*~(*s*(*t*)), is captured by a single mean firing rate per neuron. For clarity of presentation, we assume that this is the firing rate when the monkey was cued away from this stimulus (i.e. to the opposite side). The response of each neuron is then affected by the cue signal to a different degree, which is captured by an additional free parameter per neuron (the "coupling weight" to the cue signal), *u~n~*. Together, these two factors comprise the "classical" model of attention, wherein the attentional cue results in a change in the mean firing rate of each neuron.

The second modulatory signal, *d*(*t*), is motivated by recent reports of coordinated, global slow fluctuations in neural response. These have been hypothesized to arise from fluctuations in global state variables such as arousal, with different neurons potentially affected by these signals to different degrees ([@bib21]; [@bib16]; [@bib60]; [@bib46]; [@bib40]). The V4 data analyzed here also exhibit slow global drifts in neural gain in each recording, with typical time constants on the order of minutes to tens of minutes ([Figure 1---figure supplement 1](#fig1s1){ref-type="fig"}). We therefore allowed a single global slow signal to affect the population. Its timescale and time course were inferred from the population responses, as well as its coupling weights to each neuron (see Materials and methods). While this substantially improved model predictions (see below), global state fluctuations have been described elsewhere and are not the main focus of this work. Importantly, the inclusion or exclusion of this global drift signal (even allowing more than one) does not qualitatively change the results or interpretation presented here.

Our primary interest was in the *K* fast modulatory signals, *m~k~*(*t*), which we introduced to explain any shared structure in the population activity that remains after the structure of the task and global slow drifts have been accounted for. We assumed that these consisted of a small number of independent sources. The number of sources was itself a free parameter of the model, which we explore below.

To evaluate the model fits, we held aside 20% of the spike counts from the full set of stimulus presentations, randomly chosen across units and time. The accuracy of model predictions for these held-out data improved substantially with the inclusion of each model component ([Figure 2a](#fig2){ref-type="fig"}). Since all scores are *cross-validated*, these improvements indicate that the models are explaining structure in the data rather than overfitting to random response fluctuations. The figure shows that the shared modulators were twice as important as the experimentally-controlled attentional cue signal in making predictions on held-out data. This improvement is remarkable given that the stimulus presented was identical on every trial: these model components are capturing a substantial portion of the trial-to-trial variability in the population responses.10.7554/eLife.08998.006Figure 2.The fitted model explains the observed spiking responses, with estimated modulators that are both anatomically and functionally targeted.(**a**) Performance comparison of various submodels, measured as log-likelihood (LL) of predictions on held-out data. Values are expressed relative to performance of a stimulus-drive-only model (leftmost point), and increase as each model component (cue, slow drift, and different numbers of shared modulators) is incorporated. The grey square shows the predictive LL for a two-modulator model, with each modulator constrained to affect only one hemisphere (i.e. with coupling weights set to zero for neurons in the other hemisphere). This restricted model is used for all results from [Figure 2d](#fig2){ref-type="fig"} onwards, excepting the fine temporal analysis of [Figure 6c](#fig6){ref-type="fig"}. (**b**) Modulators are anatomically selective. Inferred coupling weights for a two-modulator model, fit to a population of units recorded on one day. Each point corresponds to one unit. As the model does not uniquely define the coordinate system (i.e. there is an equivalent model for any rotation of the coordinate system), we align the mean weight for LHS units to lie along the positive x-axis (see Materials and methods). (**c**) Distribution of inferred coupling weights aggregated over all recording days indicates that each shared modulator provides input primarily to cells in one hemisphere. (**d**) Hemispheric modulators are functionally selective. Units which are better able to discriminate standard and target stimuli in the cue-away condition have larger coupling weights (blue line). Discriminability is estimated as the difference in mean spike count between standard and target stimuli, divided by the square root of their average variance (*d′*). Values are averaged over units recorded on all days, subdivided into five groups based on their coupling weights. Shaded area denotes ±1 standard error. Pearson correlation over all units is *r *= 0.42. This relationship is not seen for the weights that couple neurons to the slow global drift signal (gray line, Pearson correlation *r* = 0.00). The relationship between *d′* and cue weight is significant, but weaker than for modulator weight (*r* = 0.24); this is not shown here as the cue weights are differently scaled. (**e**) Same as in (**d**), but with units subdivided into subgroups according to mean firing rate. Each line represents a subpopulation of ∼500 units with similar firing rates (from red to blue: 0--7; 7--12; 12--17; 17--25; 25--35; 35--107 spikes/s). Within each group, the Pearson correlations between *d*′ and coupling weight are between 0.2--0.3, but the correlations between mean rate and coupling weight are weak or negligible.**DOI:** [http://dx.doi.org/10.7554/eLife.08998.006](10.7554/eLife.08998.006)10.7554/eLife.08998.007Figure 2---figure supplement 1.The dataset is sufficient to support the estimation of up to 8 shared modulators.To test whether the number of recovered modulators is limited by insufficiency of the dataset, we simulated data from model neural populations that were under the influence of different numbers of modulators. All simulated datasets were matched in size to the physiological data, and simulated shared gain fluctuations were adjusted in amplitude to produce the same pairwise spike count statistics as the actual data. We then fit the model to these synthetic datasets, and measured how well the model fits recovered the true modulatory structure. More specifically, we extracted default model parameters by fitting each neural population: the stimulus-driven mean firing rates *F*, the cue-dependent gains *C*, and the slow global drifts *D*. Then, for a given number of simulated modulators *K* (from 1 to 12), we sampled a (*T* × *K*) matrix of time-varying modulator values (each value i.i.d. Gaussian), and a (*K* × *N*) matrix of random weights (each weight i.i.d. Gaussian), producing a net modulator matrix *M* as the matrix product of these two. We then sampled spike counts *Y* from the generative model, $\left. Y \right.\sim{Poiss}\left( F1^{T}\; \odot \exp\left( C + D + \lambda M \right) \right)$. For each population and *K*, we chose the scaling parameter *λ* such that the median noise correlation between the simulated neurons matched the median noise correlation between the actual neurons. We next fitted the models to these simulated datasets to see how well they recovered the underlying structure. As for the actual data ([Figure 2](#fig2){ref-type="fig"}), we evaluated the model fits via the predictive log-likelihood on held-out data. Each colored line shows the predictive LLs of the fitted models for a given "ground truth" number of modulators. In comparison, the grey squares show the model performance for the actual data. There are two important patterns here. First, for simulated models containing up to 8 modulators, the predictive LLs are greatest when we fit a model having the same number of modulators as the ground truth number used to simulate the data. This demonstrates that the model is in principle able to recover more modulators than the 4 we fit to the actual data. Second, as the number of simulated modulators increases, the ability of the fitted models to make predictions on held-out data declines. This is because the total energy of the shared gain fluctuations is constrained by the measured noise correlations, and is spread amongst the simulated modulators. In this respect, the model predictions on the actual data are most consistent with simulations of 3 or 4 modulators. Finally, it is worth noting that, in simulation, when fitting more modulators than the ground truth, the predictive performance suffers. This reflects overfitting to the noise in the training set. We do not see as pronounced a decline for the model fits to the actual data: instead, the predictive LLs appear to saturate with the number of modulators. This difference between the actual and synthetic data likely reflects our assumption in the simulations that the modulators were all of equal magnitude. A saturation of predictive LL may arise when there is a small set of dominant modulators, and a number of weaker ones. We do not have the statistical power to explore such a long tail of modulatory influences within this dataset, and focus instead on the strongest components.**DOI:** [http://dx.doi.org/10.7554/eLife.08998.007](10.7554/eLife.08998.007)10.7554/eLife.08998.008Figure 2---figure supplement 2.The structure of the modulators in higher-dimensional modulator models.In the main text, we identify three striking properties of the two dominant shared modulators: (1) they each target one of the two V4 hemispheres; (2) they preferentially target the task-specific neurons within these hemispheres; and (3) their variance decreases under cued attention. Here we show that these features are present within higher-dimensional modulator models. In [Figure 2---figure supplement 3](#fig2s3){ref-type="fig"}, we show that these additional modulatory components do not convey any additional structure in these three domains. It is useful to first view the modulator model as a form of exponential-family Principal Component Analysis (PCA) ([@bib10]; [@bib65];[@bib51]). Standard PCA, like the shared modulator model, uncovers directions in signal space of maximal variation. However, PCA suffers from an identifiability problem: it can uniquely recover the subspace in which a small set of signals lie, but not the coordinate axes. PCA does select a particular orthogonal coordinate system to represent this subspace, but this solution is not unique, is sensitive to noise, and typically reveals little about the underlying generative process. This same identifiability problem is present with the shared modulator model. In the two-modulator case, we are able to resolve the ambiguity in the coordinate system by exploiting anatomical information ([Figure 2b--c](#fig2){ref-type="fig"}; see Materials and methods). However, the problem of identifiability becomes more acute in higher dimensions. Here, we show that the results presented in the main text for the 1-modulator/hemisphere model also hold for the unconstrained 2-modulator model. We also extend the 2-modulator results to the 3- and 4-modulator cases. This is necessary as, unlike standard PCA, the solutions to our equations in lower dimensions do not necessarily lie within subspaces of higher-dimensional solutions. This is because the regularization scheme and algorithm we use create biases that disrupt any strict nesting. We therefore need to explicitly test whether the structures we uncover in the 2-modulator model are also present in the 3- and 4-modulator models. And this needs to be done under the limitations of the identifiability problem, i.e. without choosing a particular coordinate system for the modulation subspace. First column: In [Figure 2b--c](#fig2){ref-type="fig"}, we showed that the vectors of modulator coupling weights for LHS units and RHS units in the 2D modulator model were typically orthogonal. Here we show that this holds in higher dimensions. For each recording day, we measured the angle between the average weight vector for LHS units (${\overline{w}}_{L}$), and the average weight vector for RHS units (${\overline{w}}_{R}$), i.e. the arc cosine of their inner product. The 2-modulator hemisphere-constrained model used in most of the main text has this orthogonality enforced by constraint (top row). For the unconstrained 2-, 3-, and 4-modulator models (remaining rows), the blue histograms show the distribution of these angles across recording days. For comparison, we shuffle the anatomical labels on each unit and repeat the analysis to obtain the red histograms. The clustering of the actual data around *π*/2 indicates near orthogonality of the hemispheric weights. Second column: In [Figure 2d](#fig2){ref-type="fig"}, we showed that neurons which were task-relevant (i.e. had larger *d*′ values) were more strongly coupled to the (1D) hemispheric shared modulators. Here, we show that this holds in higher dimensions. For each recording day, we measured the magnitudes of all units' coupling weight vectors, $\left. ||w \right.||_{2}$. Green histograms show the distribution of magnitudes for the quartile of units with largest *d*′ values; brown histograms show the distribution of magnitudes for the quartile with the smallest *d*′ values. Third column: In [Figure 3b](#fig3){ref-type="fig"}, we showed that the variance of the (1D) hemispheric shared modulators changed according to the attentional cue: specifically, when the cue switched, one hemispheric modulator decreased in variance, while the other increased in variance. To show that this holds in higher dimensions, it is necessary to construct an appropriate metric for this change in second-order statistics that generalizes to higher dimensions, and that also does not depend on a choice of coordinate system. To accomplish this, we measure the effect of the attentional cue as a change in the covariance of the (multivariate) modulator. Considering the change from the cue-right to the cue-left condition, we can measure the effect on the modulator's second-order statistics via the ratio of the two modulator covariances, $C_{cueL}\; C_{cueR}^{- 1}$. The eigenvalues of this matrix then provide a coordinate-system-free measure of how the modulator statistics change. If the largest eigenvalue, *λ~max~*, is significantly greater than 1, then there is a direction in modulation space that became more variable due to the switch in cue. If the smallest eigenvalue, *λ~min~*, is significantly less than 1, then there is a direction in modulation space that became less variable due to the switch in cue. Eigenvalues close to 1 indicate that the variance of modulation in that direction was unchanged by the cue. Thus these two values, *λ~max~* and *λ~min~*, play an analogous role to the ratios of modulator variance examined in [Figure 3b](#fig3){ref-type="fig"}. The scatter plots show the distribution of *λ~max~* and *λ~min~* for the higher-dimensional modulator models. Blue points show these eigenvalues from each recording day; red points show the distributions obtained if we shuffle the cue labels for each trial. Importantly, when *λ~max~* exceeds 1 and *λ~min~* is less than 1 (i.e. when the points lie in the lower right quadrant), then the change in attentional cue is causing an increase in modulator variance in one direction, and a decrease in modulator variance in an orthogonal direction. These effects are clear (and significant, compared with the null distribution in red) in all cases.**DOI:** [http://dx.doi.org/10.7554/eLife.08998.008](10.7554/eLife.08998.008)10.7554/eLife.08998.009Figure 2---figure supplement 3.The modulators\' anatomical, functional, and attentional structure manifests primarily within the dominant two dimensions of modulation.In the main text, we recover a set of anatomical, functional, and attentional properties within a 2-dimensional modulator model. In [Figure 2---figure supplement 2](#fig2s2){ref-type="fig"}, we demonstrate that these properties also manifest within 3- and 4-modulator models. We wondered whether any additional structure can be seen in the extra two dimensions of modulation that we can include beyond the 2-modulator model. To answer this, we partitioned the 4-dimensional modulator space into two 2-dimensional halves. For each recording day, we define a particular 2D subspace of modulation along anatomical grounds. We measured the mean weight vectors for LHS and RHS units, ${\overline{w}}_{L}$ and ${\overline{w}}_{R}$ respectively. As shown in [Figure 2---figure supplement 2](#fig2s2){ref-type="fig"}, these two vectors were always near-orthogonal. We define their 2D span as the "hemispheric subspace of modulation", $H\; = \;{span}\left( {\overline{w}}_{L},{\overline{w}}_{R} \right)$. This is a 2D subspace of the 4D weight vectors, capturing the largest component of hemispheric-specificity in the modulator weights. What remains in the 4D modulation space is the hemispheric subspace's orthogonal complement, $H^{\bot}$. This divides the 4D space into two 2D subspaces, and thus amounts to a partial choice of a coordinate system. We can therefore study the anatomical and functional properties of the coupling weights in $H$ and $H^{\bot}$, and also the attention-dependent statistics within the corresponding 2D spaces of time-varying modulator values. This panel shows that all three properties described in [Figure 2---figure supplement 2](#fig2s2){ref-type="fig"} manifest predominantly in $H$, but not in $H^{\bot}$. In summary, each V4 hemisphere is being driven by a shared modulatory signal, that preferentially affects task-specific neurons, and has statistics that depend on the attentional cue provided to the animal. In addition, there is some evidence that other shared modulatory factors are affecting the population of V4 neurons. However, these latter signals do not share the same properties: their net effects are weaker, they do not appear to have the same anatomical or functional specificity, and they do not appear to be affected by the attentional cue.**DOI:** [http://dx.doi.org/10.7554/eLife.08998.009](10.7554/eLife.08998.009)10.7554/eLife.08998.010Figure 2---figure supplement 4.Units with higher mean firing rates typically had stronger coupling to their respective population modulator (*r*^2^ = 0.21).This observation motivates the control analyses shown in [Figure 2e](#fig2){ref-type="fig"}, [Figure 5---figure supplement 1](#fig5s1){ref-type="fig"} and [Figure 8---figure supplement 1](#fig8s1){ref-type="fig"}.**DOI:** [http://dx.doi.org/10.7554/eLife.08998.010](10.7554/eLife.08998.010)

Dimensionality and connectivity of shared modulation {#s2-2}
----------------------------------------------------

We next asked how many independent shared modulators were needed to explain the recorded responses. In principle, the activity of the recorded populations could be influenced by a very large number of shared modulators, each with its own unique connectivity and temporal patterns. We found that the spike count predictions on held-out data improved substantially with the inclusion of two modulators, and showed a modest additional increase for up to four modulators. Including additional modulators beyond these did not improve predictions.

To verify that this outcome did not reflect insufficiency of the dataset in constraining the model, we generated synthetic datasets by simulating responses of the model with different numbers of modulators, each of equal strength ([Figure 2---figure supplement 1](#fig2s1){ref-type="fig"}). These synthetic datasets had the same number of neurons and trials as the true data, and we scaled the magnitude of the simulated modulators' fluctuations to produce the same average noise correlations as found in the true data. We found that we could accurately recover up to 8 independent modulators from these synthetic datasets -- more than the number recovered from the true data.

We next asked how the estimated modulators, and their associated weights, were structured with respect to the neural populations and the task. These structures are most easily visualized for the two-modulator model, and we therefore restrict our analysis to this case for the remainder of this article. Results for three- and four-modulator models are given in [Figure 2---figure supplement 2](#fig2s2){ref-type="fig"}.

Several striking patterns emerge in the fitted model. First, although the model was not given any information regarding anatomical location or connectivity of the neurons, the estimated coupling weights for each modulator clearly identify the hemisphere in which the corresponding neurons reside ([Figure 2b--c](#fig2){ref-type="fig"}). In a given recording, each modulator had largely positive weights for neurons in one hemisphere (indicating that these neurons' gains were being co-modulated by this signal), and small weights for neurons in the opposite hemisphere. Thus, the weight structure of the estimated model suggests that the two hemispheric V4 subpopulations are modulated by two independent signals. Based on this observation, we examined a restricted model, in which we explicitly enforce the hemispheric assignment of each modulator (from here on, referred to as the LHS and RHS modulators), by setting the weights for neurons in the opposite hemisphere to zero. This enforced assignment reduces the number of weight parameters by a factor of two, and results in a modest improvement in the quality of model predictions on held-out data, reaching the level of the 4-modulator model ([Figure 2a](#fig2){ref-type="fig"}, grey square). For the remainder of this article, we retain this enforced assignment of the modulators to the hemispheres.

In addition to the distinct anatomical connectivity of the two modulators, we found that they also exhibited specific *functional* connectivity. We characterized the task-relevance of each neuron in terms of its ability to discriminate the standard and target stimuli in the cue-away (inattentive) condition. Specifically, for each neuron, we computed the difference in mean spike count for the two stimuli, relative to the standard deviation (known in the perceptual psychology literature as *d*′). Comparison of these values to the modulator coupling weights indicates that the modulators preferentially targeted the most task-informative neurons (*r* = 0.42; [Figure 2d](#fig2){ref-type="fig"}). One might suspect that some portion of this effect is simply due to firing rate---units with higher mean firing rates typically had stronger coupling weights to the modulators, paralleling previous observations ([@bib8]; [@bib42]; [@bib12]; [@bib17]) ([Figure 2---figure supplement 4](#fig2s4){ref-type="fig"}). But the relationship between task informativeness and modulator coupling remains robust even when conditioned on firing rate ([Figure 2e](#fig2){ref-type="fig"}). Finally, this correlation with functional specificity was weaker for the coupling weights to the cue signal (*r* = 0.24), and entirely absent from the coupling weights to the slow global drift signal (*r* \< 0.01; [Figure 2d](#fig2){ref-type="fig"}).

The effect of attention on shared gain {#s2-3}
--------------------------------------

The model thus far reveals the action of two structured modulatory signals, each providing input to task-informative V4 neurons in one of the two hemispheres. Given that these signals were recovered from population activity during a cued attentional task, we next ask whether they exhibit systematic changes across differently cued blocks.

The model recovers estimates of the slow drift signal and each hemisphere's modulator for every stimulus presentation ([Figure 3a](#fig3){ref-type="fig"}). Here, a clear pattern is evident: when the monkey was cued to attend to one visual hemifield, the shared modulator of the corresponding (contralateral) V4 population had a smaller variance ([Figure 3b](#fig3){ref-type="fig"}). The modulators' mean values were unchanged across blocks, as any changes in mean rate are captured by the coupling to the cue signal: 75% of neurons had a positive coupling weight to the cue signal, capturing an increase in their firing rate under cued attention. Thus, in general, attention both *increases* and *stabilizes* the time-varying gain of the corresponding neural population.10.7554/eLife.08998.011Figure 3.Time-varying model signals that determine the gain of units in one V4 hemisphere.(**a**) Example values of the cue signal (imposed by experiment), the slow drift (inferred), and a single hemispheric modulator (inferred) across stimulus presentations for one day and hemisphere. In the model, the gain of each neuron is obtained by exponentiating a weighted sum of these three signals (see [Equation (1)](#equ1){ref-type="disp-formula"}). Histogram in the bottom left shows the distribution of modulator values when the monkey was cued towards the contralateral side (blue), and away from it (red). (**b**) Modulator variance decreases under cued attention. Histogram shows the ratio of modulator variances estimated in the two cue conditions. Averaged across days and hemispheres, cued attention reduces modulator variance by 23%.**DOI:** [http://dx.doi.org/10.7554/eLife.08998.011](10.7554/eLife.08998.011)

These two changes in shared gain, in turn, provide a simple explanation for the observed changes in the statistics of individual and paired neural spike counts. Consider first a simulation of two conditionally-Poisson neurons. In the classical model of attention, these neurons' gains increase when the cue is directed to the appropriate hemifield. This produces two major effects on each neuron's marginal spike count statistics ([Figure 4a](#fig4){ref-type="fig"}): the mean increases, and the variance goes up as well (due to the Poisson mean-variance relationship). The Fano factors (ratio of variance to mean) remain unchanged.10.7554/eLife.08998.012Figure 4.Changes in the statistics of the inferred modulator under cued attention explain the observed changes in spike count statistics.(**a**) The classical model of attention. Simulation of two neurons with positive coupling weights, $u_{1}$ and $u_{2}$, to the cue signal. When the cue is directed to the corresponding spatial location (top), both the mean and variance of the simulated neurons' spike counts increase (bottom). Shaded areas demarcate analytic iso-density contours, i.e. the shape of the joint spike count distributions. (**b**) The effect of the shared modulator. Simulation of two simulated neurons with positive coupling weights, $w_{1}$ and $w_{2}$, to a shared modulator. A decrease in modulator variance leads to a decrease in both the variance and correlation of spike counts (bottom). (**c**) Effects on an example pair of units within the same hemisphere, on one day of recording. The cue increases the gain of both cells (numbers indicate cue coupling weights), and the inferred modulator exhibits a decreased variance in cued trials (again numbers indicate coupling weights; top). The spiking responses of the cells exhibit a combination of the effects simulated in (**a**) in (**b**): increased mean, decreased Fano factor, and decreased correlation (bottom; means from 7.0 to 8.1 and 6.0 to 7.2 spikes/stim respectively; Fano factors from 1.9 to 1.6 and 1.7 to 1.6 respectively; correlation from 0.19 to 0.10). The shaded areas demarcate smoothed iso-density contours estimated from the data.**DOI:** [http://dx.doi.org/10.7554/eLife.08998.012](10.7554/eLife.08998.012)

Next, consider what happens if these simulated neurons are coupled to the shared modulator ([Figure 4b](#fig4){ref-type="fig"}). A decrease in the variance of this shared signal leads to a reduction in the spike count variance of each neuron, without a large change in their mean firing rates ([@bib21]). Consequently, their respective Fano factors decrease. Moreover, since both neurons are coupled to the same modulatory signal, the decrease in modulator variance also causes a decrease in the spike-count correlation of the pair.

The data recorded in the experiment exhibit a combination of these effects. An example is shown in [Figure 4c](#fig4){ref-type="fig"}. Under cued attention, the variance of this V4 hemisphere's shared modulator decreases. When we consider two neurons within this population with strong coupling weights to both the cue signals and the modulators, their marginal and joint spike count statistics are changed by attention as predicted by the simulation: their mean firing rates go up, their Fano factors go down, and their noise correlations decrease.

The behavior shown in [Figure 4](#fig4){ref-type="fig"} crucially depends on the two neurons being coupled to the shared modulator. The model thus makes a very specific prediction: the magnitudes of the Fano factors and noise correlations should increase with the magnitude of the weights with which the neurons are coupled to the hemispheric modulators. We find that this prediction is clearly borne out by the data ([Figure 5a](#fig5){ref-type="fig"}), and is robust when controlled for firing rate ([Figure 5---figure supplement 1a](#fig5s1){ref-type="fig"}). In addition, since the modulator variance decreases under cued attention, the model predicts that those neurons that are more strongly coupled should exhibit a larger attention-induced reduction in these measures. This effect is also apparent in the data ([Figure 5a](#fig5){ref-type="fig"}). Ultimately, the model accounts for the majority (but not all) of the attention-induced changes in Fano factor and noise correlation ([Figure 5b](#fig5){ref-type="fig"}). Thus, our population-level model accounts not only for the single neuron statistics and pairwise correlations, but also for the diversity of effects seen across the population.10.7554/eLife.08998.013Figure 5.Attention-induced changes in neural response statistics are larger for neurons that are more strongly coupled to the shared modulator.(**a**) Observed Fano factor and noise correlations, as a function of model coupling weight. Units from all days are divided into five groups, based on their fitted coupling weight to their respective population modulator (model-based quantities), as in [Figure 2d](#fig2){ref-type="fig"}. Points indicate the average Fano factors and noise correlations (model-free quantities) within each group, when attention was cued towards the associated visual hemifield (blue) and away from it (red). Shaded area denotes ± 1 standard error. Unitwise Spearman correlations: *ρ* = 0.31/0.44/0.40/0.51 (fano cue twds/fano cue away/ncorr cue twds/ncorr cue away). (**b**) Comparison of model-predicted vs. measured decrease in Fano factor and noise correlation. Units are divided into ten groups, based on coupling weights (darker points indicate larger weight). The model accounts for 62% of the cue-induced reduction in Fano factor, and 71% of the reduction in noise correlation. (**c**) Comparison of cue weights and modulator weights. Units that are strongly coupled to the cue signal are typically strongly coupled to the modulator signal, though the relationship is only partial (unitwise Spearman correlation: *ρ* = 0.26). These results are robust when controlled for firing rate ([Figure 5---figure supplement 1](#fig5s1){ref-type="fig"}).**DOI:** [http://dx.doi.org/10.7554/eLife.08998.013](10.7554/eLife.08998.013)10.7554/eLife.08998.014Figure 5---figure supplement 1.Differences in coupling weight explain the observed statistics of single and pairwise firing rates, even when controlled for mean firing rate.The relationship shown in [Figure 5](#fig5){ref-type="fig"} contains a potential confound in that units with higher mean firing rates are typically more strongly coupled to their corresponding shared modulators ([Figure 2---figure supplement 4](#fig2s4){ref-type="fig"}). Here, we perform the same analyses on subpopulations with similar mean firing rate. (**a**) We repeated the analysis of [Figure 5a](#fig5){ref-type="fig"}, but subdivided the total population of units in two ways: first, by mean firing rate into six groups (rows), and then by coupling weight into five subgroups (points on each plot). Each row thus replicates [Figure 5a](#fig5){ref-type="fig"} for a controlled subpopulation of approximately 500 units with similar firing rates. Within each group, the correlation between mean rate and modulator coupling weight was weak or negligible. Nevertheless, the relationships of Fano factor and noise correlation to modulator weight remain. (**b**) We also repeated the analysis of [Figure 5b](#fig5){ref-type="fig"}, subdividing the total population of units by mean firing rate into six groups, as in the rows of (**a**). Again, the relationship between cue and modulator weights remains.**DOI:** [http://dx.doi.org/10.7554/eLife.08998.014](10.7554/eLife.08998.014)

The model does not assume any relationship between the coupling weights for cue and shared modulators. We did observe that neurons which were most strongly coupled to the modulators were also most strongly coupled to the cue signal ([Figure 5c](#fig5){ref-type="fig"}; even when controlled for firing rate [Figure 5---figure supplement 1b](#fig5s1){ref-type="fig"}). But overall, the correlation between the coupling weights to the cue signal and the coupling weights to the hemispheric modulator was modest (Spearman *ρ* = 0.26). We fit a restricted model in which these two signals had identical coupling weights, and found that the predictive (log-likelihood) performance was reduced by 20%. Moreover, while the modulator weights were predictive of the Fano factor and noise correlation effects (as shown in [Figure 5a](#fig5){ref-type="fig"}; correlations *ρ *\> 0.3), the cue weights were much less so (*ρ* \< 0.1). Thus, although the effects of the cue and of the inferred hemispheric modulator overlap, they are not identical, suggesting that they arise (at least in part) from distinct sources.

We made three additional observations regarding the shared modulators. First, the two hemispheres' time-varying modulator values were almost completely uncorrelated ([Figure 6a](#fig6){ref-type="fig"}; *r* = 0.03). The very small positive correlation might result from the influence of a global signal spanning both hemispheres. Second, the modulators exhibited correlations over successive stimulus presentations ([Figure 6b](#fig6){ref-type="fig"}), indicating that gain fluctuations can persist over time scales of seconds. Third, while the data were insufficient to infer the value of the modulators at a sub-trial resolution (i.e. at time scales shorter than 200 ms) for individual trials, we were able to estimate their average sub-trial time course. We found that the shared modulators had only weak effects during response onsets, and exerted their influence on neural gain primarily during the sustained response period ([Figure 6c](#fig6){ref-type="fig"}). This property of the modulators is consistent with previous reports that attentional effects on firing rates (and Fano factors) are greater during sustained periods ([@bib8];[@bib42], and with more general observations that network behavior is less affected by context or state during onset transients ([@bib7]).10.7554/eLife.08998.015Figure 6.Statistical properties of the hemispheric modulators.(**a**) Joint statistics of the two hemispheric modulators. Blue points: simultaneous values of the two modulators aggregated over all days. Thick black ellipse: iso-density contour at one standard deviation of the Gaussian density matching the empirical covariance. Thinner black ellipse: two standard deviations. Dashed lines: principal axes (eigenvectors) of this covariance, with the thicker dashed line indicating the axis with the larger eigenvalue. The vertical elongation of the ellipse shows that the variance of the modulator for the cued side is smaller than the variance of the modulator for the opposite side. The slight clockwise orientation shows that the two modulators have a very small positive correlation ($r = 0.03$, $p$ negligible). (**b**) Autocorrelation of modulators across successive stimulus conditions. Individual lines show the within-block autocorrelation of each estimated modulator; the thick lines shows the average across days and hemispheres. For simplicity of presentation, the targets and the gaps between trials have been ignored. The time constant of this process is on the order of several seconds. (**c**) Average time course of shared modulation within each stimulus presentation. We extended the population response model by allowing the value of the modulator to change over the course of each stimulus presentation. Given limitations of the data at fine temporal resolutions, we assumed that the temporal evolution of the modulator within each stimulus presentation followed some stereotyped pattern (up to a scale factor that could change from one stimulus presentation to the next; see Materials and methods). Fine blue lines: modulators' (normalized) temporal structure extracted for each recording day. Heavy black line: average across days. Grey shaded area: normalized peri-stimulus time histogram (arbitrary units) of spiking responses during presentations of the standard stimuli, averaged across all units, days, and cue conditions, with zero denoting spontaneous rate. Shared modulation predominantly occurs during the sustained period and is nearly absent during the onset transient.**DOI:** [http://dx.doi.org/10.7554/eLife.08998.015](10.7554/eLife.08998.015)

Finally, we note that the relationships between attention and modulator statistics were not only evident in the hemispheric-modulator model, but were also present in models with more modulators ([Figure 2---figure supplement 2](#fig2s2){ref-type="fig"}). Nevertheless, these additional modulators did not reveal any additional anatomical or functional specificity, or attentional dependencies ([Figure 2---figure supplement 3](#fig2s3){ref-type="fig"}).

Relationship of modulation to behavior {#s2-4}
--------------------------------------

We have thus far used the model to expose a set of internal modulatory signals that selectively affect the gains of neurons in the population. We wondered whether these modulatory signals had any effect on behavior. For clarity in describing these results, we use the terms "cued" and "opposite" rather than "left" and "right" to describe the visual hemifields, V4 populations, modulators, and targets. For example, if the monkey was cued to the left during a block, then the left hemifield is cued, the monkey's right (i.e. contralateral) V4 is cued, and the right V4 population's modulator is the cued modulator. If, during a trial in this block, a target is presented on the right, we refer to it as an opposite target.

First, we asked whether the shared modulators had any influence on the monkey's trial-by-trial performance. We found that the values of both the cued and the opposite modulators, during presentations of the standard stimuli, were predictive of whether the monkey would detect the upcoming target stimulus. As may be expected, changes in each hemispheric modulator predicted performance changes for their associated targets: an increase in the cued modulator preceded an increased detection probability for cued targets; while an increase in the opposite modulator preceded higher detection probabilities for opposite targets (top left and bottom right of [Figure 7a](#fig7){ref-type="fig"}; see [Figure 7---figure supplement 1a](#fig7s1){ref-type="fig"} for full psychometrics). Suprisingly, this effect is substantially stronger for the opposite side. This result is similar to a previous study on this dataset aimed at directly decoding the trial-by-trial attentional state of the animal ([@bib9]), though the two results rely on different readouts of the population activity.10.7554/eLife.08998.016Figure 7.Inferred modulatory signals are predictive of behavioral performance, and are influenced by previous reward.(**a**) Average effect of modulator values on subsequent behavioral performance, averaged across all days and difficulty levels. Values show the average change in hit probability for targets on the cued side (left column) and the opposite side (right column) following a unit increase in the cued (top row) and opposite (bottom row) modulators. \*$p < 0.05$, \*\*$p < 0.01$, \*\*\*$p < 0.001$. Full psychometric curves are shown in [Figure 7---figure supplement 1a](#fig7s1){ref-type="fig"}. (**b**) Average effects of previous trial reward on current trial performance. Note that this is a direct comparison of the behavioral data, and does not involve the modulator model. (**c**) Average effects of previous trial reward on the value of the two hemispheric shared modulators.**DOI:** [http://dx.doi.org/10.7554/eLife.08998.016](10.7554/eLife.08998.016)10.7554/eLife.08998.017Figure 7---figure supplement 1.Relationship between modulators and behavior: additional details.(**a**) Psychometric performance, averaged across all days. These plots expand the results of [Figure 7](#fig7){ref-type="fig"} showing the interacting effects of trial difficulty and the two hemispheric modulators on task performance. The hit probability is shown as a function of the orientation change in degrees for trials where the target was on the cued side (black points) or the opposite side (single gray point; opposite side targets were only presented at 12 deg). For each condition, the color gradient shows the effects that values of the cued modulator (left panel) and opposite modulator (right panel) have on the hit probability. We fit a family of psychometric curves to the cued-target conditions, with the two modulator values as regressors; the colored lines in each panel show two of these curves, indicating the biasing effect of $m_{cued} = \pm 3\;\sigma_{cued}$ (left panel), and $m_{opp} = \pm 3\;\sigma_{opp}$ (right panel) on performance. (**b**) Left: average effects of previous trial reward on current trial performance, from [Figure 7b](#fig7){ref-type="fig"}. Note that this is a direct comparison of the behavioral data, and is not dependent on the modulator model. Right: effect of reward on performance predicted by chaining together the effects of reward on modulator ([Figure 7a](#fig7){ref-type="fig"}) and modulator on performance ([Figure 7c](#fig7){ref-type="fig"}). The biasing effect of reward on behavior, as mediated by the V4 modulators, is consistent with the observed data (left), but captures only a relatively small proportion of the total reward bias (∼5--10%). To estimate the total behavioral reward bias, we fitted Bernoulli-GLMs (i.e. GLMs with a Bernoulli observation process) which predict the response (hit/miss), given the previous trial's reward (hit for target on cued side/hit for target on opposite side/other) as regressors. When the current trial's target was cued, we treated the orientation change as an additional regressor, and we included a lapse parameter as behavioral performance typically saturated below 100% correct ([@bib73]). The effect of previous reward in this model manifests as a bias term within the sigmoid (logistic) nonlinearity. To estimate the V4-mediated reward bias, we measured how large these total behavioral reward biases were if they had to pass through the "bottleneck" of the V4 modulators. We thus fitted three GLMs: a Gaussian-GLM which predicts the cued modulator on a trial, given the previous reward (and the previous modulator values); a second Gaussian-GLM which predicts the opposite modulator on a trial in the same way; and a Bernoulli-GLM which predicts the response (hit/miss), with the two modulators on that trial as a regressor. By multiplying these two effects together (the average change in modulators due to each previous reward state in the first and second GLMs, with the modulator coefficients in the third GLM), we obtained the desired quantities.**DOI:** [http://dx.doi.org/10.7554/eLife.08998.017](10.7554/eLife.08998.017)

In addition to this modulator-driven improvement in detection probabilities, we also uncovered a striking deficit in performance on the opposite side. Specifically, an increase in one side's modulator predicted a decrease in detection probability on the other side (top right and bottom left of [Figure 7a](#fig7){ref-type="fig"}). This is not due to anti-correlation of the two modulators (they are nearly uncorrelated; [Figure 6a](#fig6){ref-type="fig"}). Rather, this pattern likely reflects the competitive structure of the task, which requires the animal to make comparisons between the represented stimuli in both hemispheres. The influence of V4 activity on detection probability in this task thus cannot be fully explained by a 1D "axis of attention" ([@bib9]), but depends on the joint (2D) gain modulation across the two hemispheres.

Second, we examined the relationship of the modulator values to the outcome of the previous trial. It is well known that both humans and animals show serial dependence in behavioral tasks: our choices in a task are biased by the percepts, actions and outcomes of previous experience ([@bib61]; [@bib22]; [@bib31]; [@bib6]; [@bib18]). The monkeys in this task also exhibit a serial dependence: after receiving a reward for correctly identifying a cued target, they were more likely to score a hit for a cued target, and a miss for an opposite target, on the next trial. Conversely, after receiving a reward for correctly identifying an opposite target, on the next trial they are substantially more likely to score a hit for an opposite target, and a miss for a cued target ([Figure 7b](#fig7){ref-type="fig"}).

This sequential bias is explicitly captured in the modulatory signals controlling the gains of the V4 populations. A hit (and thus, a reward) on a cued target biased the two modulators on the subsequent trial to shift in favor of the cued population, while a hit on the opposite target biased the two modulators on the subsequent trial to shift in favor of the opposite population ([Figure 7c](#fig7){ref-type="fig"}). Again, the effect is more substantial for the modulator opposite the cue. This effect is consistent (in sign) with the biasing effects of previous reward on the animal's behavior. However, we note that only 5--10% of the total serial dependence can be explained through the V4 modulators ([Figure 7---figure supplement 1b](#fig7s1){ref-type="fig"}). We conclude that reward must also induce additional biases in the activity of downstream neurons involved in the animals' decisions.

In summary, we found that the inferred signals that modulate neural gain in V4 populations are intimately connected with the monkey's previous and subsequent behaviors.

Discussion {#s3}
==========

We have used a statistical model to examine the shared gain fluctuations of a large population of V4 neurons, under the influence of cued spatial attention ([Figure 1](#fig1){ref-type="fig"}). When fit to measured population activity, the model reveals the hidden influence of 2--4 shared modulatory signals ([Figure 2a](#fig2){ref-type="fig"}). The two most significant modulators exhibit pronounced anatomical, functional and attentional structure: they each target neurons in one hemisphere ([Figure 2b--c](#fig2){ref-type="fig"}); they are most strongly coupled to the most task-relevant neurons ([Figure 2d--e](#fig2){ref-type="fig"}); and their temporal statistics change under cued attention ([Figure 3](#fig3){ref-type="fig"}). Specifically, when attention is directed to one hemifield, the modulatory signal associated with the corresponding V4 hemisphere decreases in variance, thus stabilizing the shared gain fluctuations within that subpopulation. This, in turn, provides an interpretable and parsimonious account of the observed effects of attention on V4 responses: a reduction in variability (Fano factor) and pairwise spike count correlation ([Figure 4](#fig4){ref-type="fig"}). The estimated coupling strengths to the modulatory signals predicts the degree to which individual neurons exhibit these effects ([Figure 5](#fig5){ref-type="fig"}). The two hemispheric modulators are nearly uncorrelated with each other, but show temporal structure on both fine (\<100 ms) and coarse (\>1 − 10s) time scales ([Figure 6](#fig6){ref-type="fig"}). Finally, the inferred modulatory signals are correlated with both the rewards received on a previous trial, and successful responses in a current trial ([Figure 7](#fig7){ref-type="fig"}).

Relationship to other results {#s3-1}
-----------------------------

Our findings are consistent with previous reports that, under attention, the activity of single V4 neurons changes in relation to aggregate activity. Specifically, attention reduces the correlation between spiking and concurrent slow fluctuations in local field potential (LFP) ([@bib19]), which themselves are reduced in power ([@bib19]; [@bib64]). These behaviors are predicted by our model as a simple consequence of the attention-driven changes in modulator variance. Since attention reduces the variance of the shared modulator, it would also be expected to reduce the variance of summed activity across the population, such as that measured in LFPs. In addition, since attention reduces the proportion of neural variance due to shared modulation, the correlation between single neuron spike counts and the aggregate LFP activity would also be reduced ([@bib21]).

A primary assumption of our model is that the patterns of shared response variability arise from stochastic signals that modulate the gain of sensory neurons. The multiplicative nature of this interaction is broadly consistent with the patterns of response variability seen throughout visual cortex, which show signatures of multiplicative, rather than additive noise ([@bib21]). Nevertheless, other modeling efforts have found evidence for shared additive noise ([@bib33]), and we cannot fully rule out the possibility that additive noise also contributes to the variability in our data. A strong test of the multiplicative assumption requires analysis of neural responses over a range of stimulus drive levels. In the dataset analyzed here, reliable responses to only one stimulus per neuron are available for analysis (target responses have not been included, since they are corrupted by co-occurring saccades). Nevertheless, the modest reduction in firing rates arising from adaptation to repeated presentations of the standard stimuli provides some opportunity to examine this question. We find that the patterns of response variance along this axis are indeed consistent with the multiplicative hypothesis ([Figure 1---figure supplement 2](#fig1s2){ref-type="fig"}). A more definitive comparison of additive vs. multiplicative interactions could be achieved with more extensive experimental manipulations of stimulus parameters.

One of us recently reported that pairs of neurons within the same hemisphere whose responses provide evidence for opposite behavioral choices can exhibit increased noise correlations under spatial attention ([@bib58]). At first glance, this result seems at odds with the attention-induced reduction of variance in the inferred modulators of our model, which generally leads to a decrease in noise correlations. However, a plausible explanation may arise from considering that populations are not just affected by attention-dependent sources of modulation, but also by additional sources that are attention-independent. In the datasets we study here, we identify two examples: the global slow drifts ([Figure 1---figure supplement 1](#fig1s1){ref-type="fig"}), and the extra, non-hemispheric, fast modulators ([Figure 2---figure supplement 3](#fig2s3){ref-type="fig"}). If the two subpopulations responsible for encoding the two visual stimuli each have their own, separate, attention-dependent modulator, while both populations are subject to an additional set of common, attention-independent modulations, then attention might serve to "unmask" the cross-population correlations arising from these common modulatory signals. Preliminary simulations suggest that such a model could explain these observations.

Finally, the model fits indicate that slow, global drifts and structured, rapid fluctuations in shared gain make substantial contributions to the super-Poisson variability of these neurons. However, the model does not account for all of the observed cue-dependent changes in variability ([Figure 5b](#fig5){ref-type="fig"}), suggesting that the model structure (e.g. Poisson spiking, with rate modulated by an exponentiated sum of shared signals) is too restricted to capture the full extent of the effects. Moreover, the dimensionality and structure of the estimated modulators are limited by the experiment itself ([@bib20]), and thus cannot be taken as a complete account of modulatory activity in V4. The behavioral task is designed to engage two attentional conditions, each associated with stimuli presented over extensive blocks of trials in one hemifield of the visual world. Similarly, the stimuli in the task provide only a minimal characterization of the selectivity of individual neurons (e.g. receptive field locations, tuning properties). We expect that a more complex task, in which attention is precisely focused in terms of visual location, stimulus properties, and/or time, could potentially reveal modulatory signals targeting more specific subpopulations of neurons.

Interpreting the model components {#s3-2}
---------------------------------

Our model is functional, and like other functional models of neural responses, such as those based on Generalized Linear Models ([@bib71]; [@bib52]), its components are phenomenological rather than biophysical. The principal value of these models is in providing a parsimonious quantitative framework that explains the relationship between stimuli, neural responses, and behavior. In this respect, a primary contribution of our work is to show that the patterns of neural variability in a visual cortical population during this attentional task are dominated by a few internal signals. This need not have been the case: the model could have required many more modulators to explain the structure of observed population activity, or the imposed structure of the model might have proven inappropriate or insufficient to provide an account of the data.

Notwithstanding the abstraction of the functional model, the inferred properties of the low-dimensional modulatory signal provide constraints on potential underlying mechanisms. We can envision three broad scenarios for the mechanistic source of shared gain fluctuations: bottom-up (stimulus-driven) input, recurrent activity within the population, and top-down state signals. Each of these mechanisms has previously been proposed as an explanation for excess neural response variance or noise correlations ([@bib21]; [@bib16]; [@bib15]; [@bib17]; [@bib79]; [@bib62]; [@bib63]; [@bib45]; [@bib77]; [@bib3]; [@bib34]). The latter two mechanisms seem compatible with the low-dimensional nature ([Figure 2a](#fig2){ref-type="fig"}), spatial scale ([Figure 2b--c](#fig2){ref-type="fig"}), and temporal scale ([Figure 6](#fig6){ref-type="fig"}) of the modulatory signals we have identified. Top-down and recurrent mechanisms need not operate independently: recent theoretical work proposes a combination of these two mechanisms, in which top-down modulation induces a change in the balance of excitability in a recurrent E-I network ([@bib74]). One potential top-down biological substrate is cholinergic input, which has a multiplicative ([@bib13]) and hemispherically-specific ([@bib41]) effect on neurons in macaque visual cortex, and is known to be enhanced under attention ([@bib27]). However, it has recently been shown that in area V1, cholinergic input mediates attention-dependent increases in mean firing rate, but not changes in response variance or noise correlations, which appear to depend on NMDA receptors ([@bib26]). Consistent with these findings, our model suggests that the attention-mediated increase and stabilization of gain in area V4 might arise from two separate mechanisms: the estimated coupling weights associated with the cue and modulator signals are similar but not identical ([Figure 5c](#fig5){ref-type="fig"}), and when they are forced to be identical, model predictions worsen. A full mechanistic account of the modulators could be developed in future experiments by comparing the estimated modulator signals against other neural or physiological signals (e.g. neurotransmitter concentrations, or pupil dilation ([@bib40]; [@bib56]), and by comparing the weights with detailed anatomical and functional properties of the neurons in the population.

Implications for neural coding {#s3-3}
------------------------------

We have used a model that extracts patterns of shared gain variability within a neural population, encapsulating them as "modulators". In turn, we have examined the modulators' statistical structure, and found that it is consistent with a variety of externally measured or controlled quantities. But we are left with an outstanding question: should we think of these modulators as signal or noise? That is, do they reflect a controlled endogenous process, or are they random fluctuations that arise because of a *lack* of control, and are thus detrimental to the encoding of incoming stimuli?

The prevailing view in the neural coding literature is the latter. There has been extensive debate about how correlated fluctuations in neural activity can confound sensory information, and can be difficult to remove once introduced. From this perspective, any process that reduces such fluctuations would improve neural coding ([@bib8]; [@bib42]; [@bib23]; [@bib28]; [@bib14]; [@bib79]; [@bib1]; [@bib66]; [@bib75]). The benefits of attention would thus be two-fold ([Figure 4](#fig4){ref-type="fig"}): by increasing the mean gain of a relevant neural population, attention would increase the signal-to-noise ratio (since, for Poisson spike counts, the response standard deviation grows only with the square root of the mean); and by simultaneously reducing the variance of the gain, attention would reduce the overall response variance ([@bib21]), and thus the deleterious effects of this nuisance variable.

If this hypothesis is correct, we would expect that gain variability in V4 should follow a general pattern. In the absence of directed attention, all neurons in V4 should be subject to some baseline level of gain variability, due to activity in locally shared circuits. This baseline should presumably be relatively independent of which stimulus features are encoded by which neurons. When attention is directed to a particular neural subpopulation, the shared variability should decrease, improving the coding precision of that subpopulation ([Figure 8a](#fig8){ref-type="fig"}).10.7554/eLife.08998.018Figure 8.Interpreting the role of shared modulation.(**a**) Illustration of how shared gain fluctuations would behave if they were noise, i.e. undesirable random fluctuations. In baseline conditions (red), gain fluctuations would be expected to have similar variance for all neurons in the V4 population. The action of attention would be expected to reduce the variance of gain fluctuations in task-relevant neurons, so as to mitigate their adverse effect on coding precision (see [Figure 4](#fig4){ref-type="fig"}). (**b**) Contrary to this simple "noise" interpretation, the variance of shared gain fluctuations are markedly larger for task-relevant neurons than task-irrelevant neurons in baseline (cued away) conditions. Moreover, although this variance decreases under attentional cueing (cued toward), it remains larger for the task-relevant neurons. Functional relevance for each unit is measured as *d*′ (as in [Figure 2d--e](#fig2){ref-type="fig"}); shared gain variability, $\sigma_{g}^{2}$, is measured as the total variance of model-estimated gain fluctuations (from slow drift and modulators combined). These results are robust when controlled for firing rate ([Figure 8---figure supplement 1](#fig8s1){ref-type="fig"}).**DOI:** [http://dx.doi.org/10.7554/eLife.08998.018](10.7554/eLife.08998.018)10.7554/eLife.08998.019Figure 8---figure supplement 1.Variance of shared gain fluctuations is larger in task-relevant neurons, even when controlled for firing rate.The relationship shown in [Figure 8](#fig8){ref-type="fig"} contains a potential confound in that units with higher mean firing rates are typically more task-relevant, and also exhibit larger gain fluctuations. Here, we perform the same analyses on subpopulations with similar mean firing rate (subpopulations as in [Figure 5---figure supplement 1](#fig5s1){ref-type="fig"}).**DOI:** [http://dx.doi.org/10.7554/eLife.08998.019](10.7554/eLife.08998.019)

This, however, is not the pattern we observe. The gain fluctuations we uncover do not affect every neuron in the population, but specifically target those neurons relevant to the task ([Figure 8b](#fig8){ref-type="fig"}). And while their variance certainly decreases under attention, it remains a puzzle why these gain fluctuations are there at all, given that they are nearly absent from the task-irrelevant neurons (which, presumably, would be relevant for different tasks).

This finding suggests that the fluctuating modulatory signals are not noise, but rather reflect meaningful intrinsic signals that play a role in some ongoing computation in the brain. For instance, one tempting possibility is to identify them as fluctuations in the attentional signal itself (e.g. changes in the spatial locus of attention) ([@bib21]; [@bib9]; [@bib25]; [@bib15]), or, more generally, as the local manifestation of a dynamic resource allocation strategy. This might reflect a shifting belief-state of the animal in the likelihood or utility of each of the two targets, as proposed by a recent model ([@bib24]). Such an account might explain the modulators' targeting of task-relevant neurons, their low-dimensional structure, and their connection to behavior.

Regardless of what we call this endogenous signal, it is worthwhile considering whether the gain fluctuations it produces pose a problem for downstream interpretation of the neural code. We can envisage three scenarios. First, the downstream neurons that decode the V4 activity may be invariant to the shared gain fluctuations \[for example, in a linear decoding framework, the gain fluctuations might lie parallel to the decision boundaries ([@bib44])\]. However, this does not appear to be the case: modulator fluctuations have a direct and immediate effect on the perceptual decisions made by the animal ([Figure 7a](#fig7){ref-type="fig"}).

A second scenario is that downstream neurons might have access to the fluctuating shared gain values, and could thus compensate for their biasing effect on perceptual decisions (e.g. [@bib67]). This does not imply that the gain fluctuations would have no effect: they would still alter the instantaneous signal-to-noise ratio of the neural representation. This would predict that when a modulator increases, detection probability on the corresponding side would improve (as we see in the top left and bottom right of [Figure 7a](#fig7){ref-type="fig"}). However, the observation that the value of a modulator also affects performance on the other side (top right and bottom left of [Figure 7a](#fig7){ref-type="fig"}) is inconsistent with this hypothesis. Any downstream decoder must compare the responses of the two hemispheres in order to perform the task. Thus, the unexpected impact of each hemisphere's modulator on ipsilateral performance demonstrates that the decoder is unable to fully discount the fluctuations in gain.

We are left with a third scenario, in which the downstream neurons cannot perfectly compensate for the shared gain fluctuations in V4. Considered in isolation, shared modulation would thus be detrimental to coding precision. We might speculate that this negative would be offset by some other (as of yet, unknown) benefit of a flexible, varying gain signal.

Overall, we find that the effects of attention on the activity of V4 neurons, while seemingly complex at the level of observed responses, reflect population-level patterns that are simple and low-dimensional. This contributes to an emerging trend in systems neuroscience, wherein the dynamics and statistics of large populations have been found to follow structured, coordinated patterns ([@bib78]; [@bib46]; [@bib38]; [@bib5]; [@bib36]; [@bib68]; [@bib29]; [@bib11]). These findings have been driven both by experimental advances in simultaneous acquisition of responses over populations, as well as the development of statistical modeling tools for isolating low-dimensional latent components of population activity ([@bib30]; [@bib47]; [@bib35]; [@bib72]; [@bib2]). Here, we have shown that these tools may be used to expose and characterize signals that underlie modulatory processes arising under attention.

Materials and methods {#s4}
=====================

Stimulus design and data processing are described in ([@bib8]). We only included data from units with mean firing rates greater than 0.5 spikes/s.

Models {#s4-1}
------

We fitted probabilistic models to population neural responses. For each recording day, we collect the spike counts from the *N* neurons over the *T* stimuli as a (*T* × *N*) matrix, *Y*. We assume that each count was drawn from a Poisson distribution, with the (*T* × *N*) rates obtained from the model response *r~n~*(*t*) as defined in [Equation (1)](#equ1){ref-type="disp-formula"}. It is convenient here to express the rates in matrix form, *R*, such that $\left. Y_{t,n} \right.\sim~{Poiss}\left( R_{t,n} \right)$. We describe the components and fitting procedures below.

Each neuron *n* has a mean firing rate for the stimulus, *f~n~*(*s~t~*). Since we only analyze responses to the standards, this is a single scalar value per unit (though see section "Fine temporal analysis of shared modulation" below). We define this to be the mean rate in the cue-away condition (with the exception of [Figure 2a](#fig2){ref-type="fig"}, where the stimulus-drive-only model uses a mean rate across all conditions). This mean rate is modulated by a set of multiplicative factors. Since these factors have to be positive, we adopt the convention that they are derived from signals that can be positive or negative, which are then transformed elementwise by a nonlinear function *h* that is monotonic and has positive range. Choosing *h* as convex and log-concave simplifies inference further ([@bib48]). The results presented in the main paper use the exponential function exp(·), but we found that using the soft-threshold function log(1 + exp(·)) produces qualitatively consistent results. When using the exponential function, the notation simplifies , and the product of multiplicative factors may be written as a single exponential of a sum, as in Equation (1).

We now describe the gain factors introduced in [Equation (1)](#equ1){ref-type="disp-formula"}. We define the cue signal during the experiment as a length-*T* binary vector, *c*, indicating the cue direction for each trial. Then we can write the elements of the rate matrix as $R_{t,n} = f_{n}\left( s_{t} \right) \cdot \exp\left( c_{t}~u_{n} \right)$, where *u~n~* is a coupling weight of neuron *n* to the cue signal. Here, it is useful to define the cue signal differently for units in the two hemispheres, with *c~t~* = 1 indicating that the cue is directed to a hemisphere's corresponding (contralateral) hemifield, and *c~t~* = 0 otherwise. Thus a positive coupling weight *u~n~* means that a neuron increases its gain under cued attention. In matrix notation, we write *C* = *cu^T^*, and $R = F \odot \exp\left( C \right)$, with $\odot$ the elementwise product, and *F* as the matrix of values $F_{t,n} = f_{n}\left( s_{t} \right)$.

Since the joint activity on any day tended to drift slowly over time, independent of the task blocks, we introduced for each recording a slow drift signal, *d*. These were fitted after discounting the cue-dependent gains for each neuron. We constrained each recording's drift signal *d* to vary slowly over time by placing a Gaussian process prior on it. The time constants of the priors were learned for each recording via evidence optimization ([@bib49]; [@bib50]), and were typically on the order of minutes to tens of minutes ([Figure 1---figure supplement 1](#fig1s1){ref-type="fig"}). The method is described in detail in ([@bib55]). We jointly fitted the slow drift signal *d*, and the coupling weights of each neuron to this signal, *v*, via an expectation-maximization algorithm. Arranging the effects of these slow, global drifts into a (*T *× *N*) matrix *D* = *dv^T^*, this model amounts to $R = F \odot \exp\left( C + D \right)$.

We introduced fast, shared modulators by fitting a (*T *× *N*) gain matrix *M*. The rank of this matrix is a measure of the number of shared modulators: we constrain rank (*M*) = *K* for *K *= 1,2,3,.... In this way, *M* encapsulates the net effect of *K* time-varying modulators *m*^(*k*)^, with coupling weights *w*^(*k*)^, such that $M_{t,n} = \sum_{k = 1}^{K}m_{t}^{(k)}\; w_{n}^{(k)}$. The matrix *M* thus expresses a time-varying latent state of the system ([@bib30]; [@bib47]; [@bib35]; [@bib72]; [@bib2]). We include a zero-mean Gaussian prior on the elements of *M*, with $p\left( M \right) \propto - \frac{\tau}{2}\left. ||M \right.||_{F}^{2}$ to prevent overfitting. We choose the hyperparameter *τ* for each population and rank by cross-validation. We find the MAP estimate of *M* from the model with $R = F \odot \exp\left( C + D + M \right)$, under the low-rank constraint, after fitting *C* and *D* above. Since the low-rank constraint limits the feasible region for the inference problem to a non-convex set, we maximized the log posterior using the Alternating Direction Method of Multipliers (ADMM), in a manner similar to ([@bib51]). We then iterate between solving for *D* and for *C* and *M* (including their respective hyperparameters), though in practice, one or two passes are sufficient for convergence.

We cross-validated model fits as follows. For each stimulus presentation, a random subset of 20% of the spike count observations were set aside as a test set. Thus, during the inference on the training data, the values of the latent variables (in *D* and *M*) for the omitted observations did not contribute to the log likelihood. However, these values were automatically imputed by the priors/constraints on the structure of the latent signals (the drifts *d* to vary slowly in time; the shared modulators *M* to be low rank). After convergence on the training set, the imputed values of these signals were used as predictions on the test set.

Identifiability {#s4-2}
---------------

The procedure above recovers a matrix, *M*, of the shared modulators' effects of neural gain, but it is not identifiable: the results do not uniquely constrain the modulator time series (*m*) and the weights (*w*). We resolve these ambiguities as follows.

The rank-1 case arises when fitting separate modulators to each V4 hemisphere ([Figures 2d--e](#fig2){ref-type="fig"}, [3](#fig3){ref-type="fig"}--[5](#fig5){ref-type="fig"}, [6a--b](#fig6){ref-type="fig"}, [7](#fig7){ref-type="fig"}). Here there is a scale and sign ambiguity: we can write $M = mw^{T} = \left( {\alpha m} \right)\left( {\frac{1}{\alpha}w^{T}} \right)$ for any scalar *α* ≠ 0. We therefore fix Var(*m*) = 1, and resolve the scale factor into the weights. Since almost all weights for a given hemisphere had the same polarity, we resolve the sign ambiguity so that the mean weight is positive.

The rank-2 case arises when fitting two modulators to a whole population ([Figure 2b--c](#fig2){ref-type="fig"} only). In addition to the scale/sign ambiguities, there is also a rotation ambiguity: we can write $M = NW^{T} = \left( {NQ} \right)\left( {Q^{T}W^{T}} \right)$ for any orthogonal matrix *Q*. For each dataset we choose the rotation $Q$ such that the mean weight on the LHS units is $\left( {\overline{w},0} \right)^{T}$ for some value $\overline{w}$. This aligns the brown squares shown in [Figure 2b](#fig2){ref-type="fig"} and brown points in [Figure 2c](#fig2){ref-type="fig"} along the x-axis. We resolve the reflection (sign) ambiguity by requiring that the mean weight across both populations lies in the upper right quadrant. This resolution preserves angles, and thus has no bearing on the observed orthogonality of the weight distributions across paired hemispheres shown in [Figure 2b--c](#fig2){ref-type="fig"} (see [Figure 2---figure supplement 2](#fig2s2){ref-type="fig"}).

Our analyses of higher-dimensional models proceed without resolving these ambiguities ([Figure 2---figure supplements 2--3](#fig2s2){ref-type="fig"}).

Behavioral modeling {#s4-3}
-------------------

To quantify the relationship between modulators and behavior, we aggregate the data across all days, and fit a number of generalized linear models (GLMs) to this ensemble.

We estimated how correct detections for cued targets depend on the modulators ([Figure 7a](#fig7){ref-type="fig"}) by fitting a psychometric curve, with the cued and opposite modulators as regressors. We parameterized the hit probability on trial *t* as:

$$\left. {hit}^{(t)}~ \middle| ~{cued}~{target}~ \sim ~{Bern}~\left( \delta \cdot \sigma~\left( \alpha + \beta^{T}m^{(t)} + \lambda\log\left( \mathrm{\Delta}\theta^{(t)} \right) \right) \right)\mspace{600mu} \right.$$

where the superscript (*t*) indicates trial *t, m* is a vector of the cued and opposite modulators (averaged across the standard stimuli on that trial), *σ* is the logistic inverse-link function, Δ*θ* is the target orientation change, *δ* is a lapse parameter, and the Greek characters *α* (bias), *β* (dependence on the modulators), and *λ* (dependence on task difficulty) are free parameters. Other parameterizations of the dependence on target orientation did not change the main result, nor did omission of the lapse term.

We estimated how correct detections for the opposite targets depend on the modulators by:

$$\left. {hit}^{(t)} \middle| {opposite}\mspace{600mu}~{target} \sim \mspace{600mu}{Bern}\left( \sigma\left( \alpha'~ + ~\beta'^{T}m^{(t)} \right) \right) \right.$$

The values reported in [Figure 7](#fig7){ref-type="fig"} reflect the average change in hit probability from a unit increase in *m~cued~* or *m~opp~* via these two models. We also quantify these effects individually for each cued-target condition in [Figure 7---figure supplement 1](#fig7s1){ref-type="fig"}; this figure also shows the full psychometrics.

We estimated how correct detections for the cued and opposite targets depend on previous rewards ([Figure 7b](#fig7){ref-type="fig"}) by replacing the regressors *m*^(*t*)^ in the above equations with categorical variables for previous reward (hit for target on cued side / hit for target on opposite side / other), as below.

We estimated how previous reward affects the modulator values ([Figure 7c](#fig7){ref-type="fig"}) by fitting a Gaussian-GLM:

$$m^{(t)}\mspace{600mu} \sim N\left( a + \; Br^{(t - 1)} + \; Cm^{(t - 1)},\;\sigma^{2}I \right)$$

where *r*^(*t*−1)^ is the reward from the previous trial, being (1,0)^T^ for a previous hit on a cued target, (0,1)^T^ for a previous hit on an opposite target, and (0,0)^T^ for any other outcome (miss/catch trial/false alarm/invalid trial), and $a$, $B$ and $C$ are free parameters.

In all cases, we assessed the significance of parameter estimates by approximating the posterior (both through a Laplace approximation and MCMC) and estimating its integral above/below zero.

Fine temporal analysis of shared modulation {#s4-4}
-------------------------------------------

With the exception of [Figure 6c](#fig6){ref-type="fig"}, all analyses in the main text assume that the shared modulators are constant in value (and thus have uniform effects) over the course of the response to each stimulus presentation. In [Figure 6c](#fig6){ref-type="fig"}, we sought to quantify the dynamics of the shared modulator at finer time scales.

We extended the population response model presented in [Figure 1](#fig1){ref-type="fig"} and [Equation (1)](#equ1){ref-type="disp-formula"} by allowing the value of the modulator to change over the course of a stimulus presentation. As the data are very limited at fine temporal resolutions, we could not reasonably estimate the modulators' values in small time bins for every stimulus presentation. Instead, we assumed that the temporal evolution of the modulator within each stimulus presentation followed some stereotyped pattern (up to a scale factor that could change from one stimulus presentation to the next).

We extend the response model to capture spike counts within bins of shorter duration (here, 10 ms). We assume that the spike count *Y~t,b,n~* for neuron *n* within bin *b* of stimulus presentation *t* is given by:

$$Y_{t,b,n}\; \sim \;{Poiss}\left( F_{b,n} \cdot exp\left( C_{t,n} + D_{t,n} + M_{t,b,n} \right)\; \right)$$

We thus assume that each neuron has a stimulus-driven mean firing rate that changes from bin to bin (*F~b,n~*), but is identical across repeated stimulus presentations; the cue-dependent gains (*C~t,n~*) and slow global drift (*D~t,n~*) are constant over all bins within each stimulus presentation; and the shared modulators (*M~t,b,n~*) are now free to have structure across stimulus presentations (*t*), neurons (*n*), and bins within each stimulus presentation (*b*). This model is now extremely high-dimensional (the tensor *M* having *TNB* free parameters). To overcome this, we must impose structure on *M*, which we do by assuming that it has low tensor rank.

We model the tensor *M* as an outer product of a rank-2 matrix (with components *M~t,n~*, as in the two-modulator model presented in the remainder of the text), and a vector *ω* indexed over bins, i.e. $M_{t,n,b}\; = \;\left( {\sum_{k = 1}^{2}m_{t}^{(k)}w_{n}^{(k)}} \right)\;\omega_{b}$. The weight vector *ω* thus represents the common temporal evolution of the modulators' effects on neural gain within each stimulus presentation, which is identical (up to scale factors) across neurons and successive stimulus presentations.

We learn the whole low-rank tensor of modulator values, *M*, by maximizing the data likelihood. We perform coordinate descent on its components: we iterate between solving for the matrix with elements *M~t,n~* (via ADMM), then solving for the vector *ω* (via gradient descent). We restrict the fit to the main response period (60 ms to 260 ms after stimulus onset) as the shared fluctuations in spontaneous activity were typically large; the time series *ω* is thus shown in [Figure 6c](#fig6){ref-type="fig"} limited to this period.

Funding Information
===================

This paper was supported by the following grants:

-   http://dx.doi.org/10.13039/100000011Howard Hughes Medical Institute to Neil C Rabinowitz, Robbe L Goris, Eero P Simoncelli.

-   http://dx.doi.org/10.13039/100000002National Institutes of Health 4R00EY020844-03 to Marlene Cohen.

-   http://dx.doi.org/10.13039/100000002National Institutes of Health R01 EY022930 to Marlene Cohen.

-   http://dx.doi.org/10.13039/100006211Klingenstein Third Generation Foundation to Marlene Cohen.

-   http://dx.doi.org/10.13039/100000893Simons Foundation to Marlene Cohen.

-   http://dx.doi.org/10.13039/100000879Alfred P. Sloan Foundation Research Fellowship to Marlene Cohen.

Additional information {#s5}
======================

The author declares that no competing interests exist.

NCR, Conception and design, Analysis and interpretation of data, Drafting or revising the article.

RLG, Conception and design, Drafting or revising the article.

MC, Acquisition of data, Drafting or revising the article.

EPS, Conception and design, Drafting or revising the article.

10.7554/eLife.08998.021

Decision letter

Carandini

Matteo

Reviewing editor

University College London

,

United Kingdom

eLife posts the editorial decision letter and author response on a selection of the published articles (subject to the approval of the authors). An edited version of the letter sent to the authors after peer review is shown, indicating the substantive concerns or comments; minor concerns are not usually shown. Reviewers have the opportunity to discuss the decision before the letter is sent (see [review process](http://elifesciences.org/review-process)). Similarly, the author response typically shows only responses to the major concerns raised by the reviewers.

Thank you for submitting your work entitled \"Attention stabilizes the shared gain of V4 populations\" for peer review at *eLife*. Your submission has been favorably evaluated by Timothy Behrens (Senior editor) and three reviewers, one of whom, Matteo Carandini, is a member of our Board of Reviewing Editors.

The reviewers have discussed the reviews with one another and the Reviewing editor has drafted this decision to help you prepare a revised submission.

Summary:

This paper analyzes the factors that control the activity of populations of neurons, and the impact on those factors of attention and behavior. The paper reanalyzes an impressive data set of population responses (\>100 neurons) measured by Cohen and Maunsell (2009, 2010 and 2011) in area V4 of monkeys that were performing an attention task.

One of the paper\'s achievements is to provide a compact account of the activity of large populations of neurons, and of how this activity is affected by attention. This account is based on a simple model of multiplicative shared variability (similar to those by Goris et al., Lin et al., Ecker et al.). By applying this model, the authors show that \"paying attention\" to one group of neurons corresponds to increasing the average gain of those neurons, while stabilizing (decreasing the variability across time of) that shared gain.

The model assumes a number of gain modulators and the authors determine the fluctuations of these modulators as a function of attentional cueing and a function of the animal\'s behavior (and vice versa, how past behavior affects the modulators). The model identifies up to 4 time-varying shared modulatory signals, which affect the neuronal responses of the recorded neurons. These modulators are also linked to past behavior (reward) and are predictive of upcoming performance.

This simple model explains observations that previously appeared contradictory. One is the well-known fact that attention increases firing rates. Another is the fact that attention reduces the correlation between pairs of neurons (Cohen and Maunsell 2009, Mitchell et al. 2009). The two are apparently contradictory, because increases in firing rate gain, on their own, would have caused the opposite effect on correlations. Yet the model captures both: the first by an increase in mean gain, and the second by a decrease in the variance of the gain (which also explains the well-known reduction in variability caused by attention).

The model allows one to assess trial-by-trial variability in a more useful way than simply measuring single neuron variability or pairwise correlations. Moreover, it is a considerable improvement over prior methods to summarize the effect of attention on a large population. The prior methods involved finding a vector between the mean uncued response and the mean cued response and projecting the data onto that axis (Cohen and Maunsell, 2010). The new approach is more general, and finds the directions in neural state space that best account for trial-to-trial shared neural variability with certain broad features. In this way, variability can be assigned to dimensions that have some interpretation, and the contributions of each dimension can be quantified.

Moreover, the paper presents some key new results. Prior methods reported two dimensions of attentional fluctuations, but couldn\'t quantify how much of the total variability was accounted for by those two dimensions. Here we get a clear and important answer: nearly all of it. Another key new result is that the variability of the attentional fluctuations depends on the attentional cue, being lower in the cued hemisphere than the uncued hemisphere.

The paper also goes some way in clarifying whether the variability in neural activity associated with attention tasks corresponds to \'representation fluctuations\' or to \'attentional fluctuations\'. In most previous studies, the measured correlations were interpreted as fluctuations in the neural representation of the stimulus. Because correlations often impair the ability to accurately decode, these \'representation fluctuations\' were hypothesized to be harmful to the ability to detect changes in the stimulus. Thus, attention would aid detection by increasing signal (by increasing firing rate) and by decreasing noise in the stimulus representation (visible in noise correlations). Cohen and Maunsell 2010 and 2011 described a different form of correlated neural variability: fluctuations in the attentional signal itself. Even when attentional demands are nominally the same across trials, the internal deployment of attention fluctuates, leading to fluctuations in performance: better performance when \'attentional fluctuations\' happen to be higher.

Critically, the previous analyses conjured a complex scenario, in which attentional fluctuations would somehow affect representational fluctuations, so that when the attention signal happens to be higher, the representational fluctuations happen to be lower. The present work, instead, indicates a much simpler solution: the decreases in correlated variability caused by attention can be attributed entirely to decreases in the variability of the attentional signals themselves. There is no need to assume attentional modulation of representational fluctuations. Indeed, it isn\'t even clear if there are large representational fluctuations. Which means we are left with one parsimonious interpretation of the data, which provides a unifying explanation for prior results.

The analysis is well motivated. The paper is very accessible despite its relatively computational approach. It is of interest to a broad audience and is topical.

Essential revisions:

However, the paper also comes with weaknesses, such as the fact that many of its results recapitulate existing results. For example, much of [Figure 2](#fig2){ref-type="fig"} can be inferred from Cohen and Maunsell 2010, including the impact of the cue and two modulatory signals (panel a) and the hemispheric separation (panels b, c). One of the more compelling links to behavior ([Figure 6A](#fig6){ref-type="fig"}) is essentially a reproduction of the effect in Cohen and Maunsell.

Therefore, for this paper to have impact, it is particularly important to ensure that it does \'provide a unified account of attentional effects on population coding and behavior\' (to quote its own words). Yet, the paper seems to fall short of a full commitment to its own findings. Or rather, it commits to the aspects that agree with prior work -- there is modulation in dimensions that seem attentional, as in Cohen and Maunsell 2010 -- while downplaying the aspects that disagree with prior work. Indeed, a conclusion that can be drawn from the paper is that there is no longer any reason to assume that \'representational fluctuations\' are reduced by attention. Rather, this old explanation should be replaced with the new explanation where it is attentional variability itself that is the measured effect. In parts of the Discussion the manuscript commits to that. Yet in other parts, it reverts to explanations of the older variety. For instance, when describing the effect of cholinergic input it reverts to the old model, in which attention causes a reduction in variability. The subsequent paragraph also seems to adhere to the old view. The manuscript seems to avoid undermining prior work, rather than committing to the new unifying interpretation that in other places (including the Abstract and Introduction) it seems willing to endorse.

Ideally, this issue would be tackled not only by changing the Discussion, but by having the model make a new prediction, and look for validation of this prediction in the data. For instance, consider the following idea (provided by one of the reviewers). If the commonly observed decline in \'noise correlations\' is in fact due to a decline in attentional fluctuations (rather than a cue-induced suppression of representational fluctuations) then that makes the following prediction. If the data were divided into \'high attention\' trials and \'low attention\' trials on the basis of the internal signal relevant to that hemisphere, then the gain of the neural response should be higher in the first case, but there should be no difference in the noise correlations. This is a testable prediction.

Another area where the paper should be improved is in explaining which effects are \"explained\" by the model and which ones are implied as a hypothesis. For instance, in the Discussion, the paper reports that neurons that showed stronger reward dependent changes in firing rates (gain) were also more strongly coupled to the modulator. This is similar to statements made earlier in relation to attentional modulation. This argument may be circular: given the fitting and the associated weights, it may be an inevitable consequence of the computational approach, and thus not a real result. The same applies to the argument in the subsection "Interpreting the model components", and as stated to wherever similar arguments are made before. These are potentially major criticisms and they absolutely have to be addressed.

Finally, a major comment involves the analysis in [Figure 6](#fig6){ref-type="fig"}. This analysis employs the difference between *m~c~* and *m~u~*. In fact, this seems like exactly the wrong thing to do in the current study, as it is shown that those two signals are not correlated. Shouldn\'t each hemisphere be analyzed on its own?

\[Editors\' note: further revisions were requested prior to acceptance, as described below.\]

Thank you for resubmitting your work entitled \"Attention stabilizes the shared gain of V4 populations\" for further consideration at *eLife*. Your revised article has been favorably evaluated by Timothy Behrens (Senior editor), by Matteo Carandini (Reviewing editor), and by one reviewer.

The paper has improved considerably and is arguably ready to be published as is. It does a much clearer job at delivering a \'unifying account\' of prior results, in the light of new analyses. Where the prior version hedged its bets, the new version takes a clear stand and paints a compelling picture. However, it would be even stronger if the authors would address two main comments indicated below.

The first major comment is about something that is visible in the new version of [Figure 7A](#fig7){ref-type="fig"}, which shows the impact on behavioral performance of the inferred modulators -- one for the cued side and one for the opposite side. This is a better choice than in the previous submission, where they had expressed it as a function of the difference in these two modulators. However, this new analysis reveals something intriguing. The value of the modulator on the cued side predicts the ability to detect a change on the opposite side, even though it does not impact firing rates on the opposite side (top-right red quadrant of [Figure 7A](#fig7){ref-type="fig"}).

This is remarkable, and seems to completely undercut the theory that the impact of attention on performance is linked to changes in response gain. For example, decreasing the cued-side modulator seems to increase performance on the opposite side, even though this would not increase firing rates in the opposite side, nor decrease their variance. This is because the two modulators are essentially uncorrelated ([Figure 6A](#fig6){ref-type="fig"}). If this reasoning is correct, the cued-side modulator has an impact on performance for opposite-side stimuli without having any detectable impact on the physiology of the opposite V4 neurons. Thus, the impact of the modulator basically can\'t be due to physiological effects at the level of V4.

If this understanding of the data is incorrect, it would be a good idea to dispel it. If instead it is correct, it seems that it should be explicitly discussed. The manuscript stops short of this: it mentions that things cannot be fully explained by a \"1D axis of attention\", and the issue comes up again briefly in Discussion, but it is never really highlighted. Indeed, the finding seems inconsistent with standard models of attention, including the one presented in the paper. And it wouldn\'t be just about 2D vs. 1D: it would be about whether the physiological effects of attention in V4 have anything to do with the behavioral effects.

The second major comment is about the control analysis currently shown in the reply to the reviewers: dividing the data based on the value of the modulator (high \'attention\' versus low \'attention\') and assessing whether variability (Fano factor or noise correlations) is lower for the former (which it would be if attention suppresses other variability such as \'representational fluctuations\' as envisioned in the old view). It may be a good idea to leave this figure out of the paper, as the results are essentially a foregone conclusion if one accepts the results and interpretation of the prior figures (i.e. if one accepts that the modulator captures essentially all the relevant shared variability). However, some readers might still have doubts about whether the modulators truly account for all the shared variability (that conclusion depends on how suitable the model is) and/or may not immediately see the implications of that finding. It may therefore be a good idea to add that analysis to the supplementary figures.

10.7554/eLife.08998.022

Author response

Essential revisions:

However, the paper also comes with weaknesses, such as the fact that many of its results recapitulate existing results. For example, much of [Figure 2](#fig2){ref-type="fig"} can be inferred from Cohen and Maunsell 2010, including the impact of the cue and two modulatory signals (panel a) and the hemispheric separation (panels b, c). One of the more compelling links to behavior ([Figure 6A](#fig6){ref-type="fig"}) is essentially a reproduction of the effect in Cohen and Maunsell.Therefore, for this paper to have impact, it is particularly important to ensure that it does \'provide a unified account of attentional effects on population coding and behavior\' (to quote its own words). Yet, the paper seems to fall short of a full commitment to its own findings. Or rather, it commits to the aspects that agree with prior work -- there is modulation in dimensions that seem attentional, as in Cohen and Maunsell 2010 -- while downplaying the aspects that disagree with prior work. Indeed, a conclusion that can be drawn from the paper is that there is no longer any reason to assume that \'representational fluctuations\' are reduced by attention. Rather, this old explanation should be replaced with the new explanation where it is attentional variability itself that is the measured effect. In parts of the Discussion the manuscript commits to that. Yet in other parts, it reverts to explanations of the older variety. For instance, when describing the effect of cholinergic input it reverts to the old model, in which attention causes a reduction in variability. The subsequent paragraph also seems to adhere to the old view. The manuscript seems to avoid undermining prior work, rather than committing to the new unifying interpretation that in other places (including the Abstract and Introduction) it seems willing to endorse.

We agree that we were a bit "restrained" in our interpretations, and appreciate the encouragement. We spent a long time considering your revision comments (as well as the more general comments above), and our discussions led us to a more precise set of interpretations of the results. Although our conclusions are not exactly in the same language as you have described here, we believe they are in the same spirit, and hope you will find them a significant improvement over the original manuscript.

In particular, we were inspired to undertake an additional analysis of the modulator statistics, shown in the new [Figure 8](#fig8){ref-type="fig"}, which we think provides strong support for our interpretation that the modulators reflect attentional fluctuations (i.e., what we describe as "signal"), as opposed to representational fluctuations (i.e. what we describe as "noise"). We have substantially revised the Discussion to address this (in particular, replacing the old section "Perceptual benefits of attention" with the new section "Modulators and neural coding").

We think of the paragraph on cholinergic input (and the subsequent paragraph) as orthogonal to the interpretation of the modulators. We have tried throughout the text to be agnostic regarding the normative role of fluctuating gain signals, and only address this question in the final section of the Discussion. If we have inadvertently given a different impression elsewhere in the text, please let us know so we can correct it.

As an additional side comment, we appreciate that there are similarities between the structures we uncover (e.g. [Figure 2B-C](#fig2){ref-type="fig"}) and the structures that Cohen and Maunsell, 2010, assume. In this respect, we emphasise that our goal is to discover the patterns of shared modulation from the population data, rather than assume them. The model produces a set of results -- some of which have been previously reported, and some of which are new -- but all of them arise from the same unified population model. As per the reviewers' suggestion below, we have also extended our behavioral results ([Figure 7](#fig7){ref-type="fig"}) to a 2dimensional analysis, which further distinguishes our results from those previously published. This also emphasizes the interpretation we land on, that the modulators are signal rather than noise.

Ideally, this issue would be tackled not only by changing the Discussion, but by having the model make a new prediction, and look for validation of this prediction in the data. For instance, consider the following idea (provided by one of the reviewers). If the commonly observed decline in \'noise correlations\' is in fact due to a decline in attentional fluctuations (rather than a cue-induced suppression of representational fluctuations) then that makes the following prediction. If the data were divided into \'high attention\' trials and \'low attention\' trials on the basis of the internal signal relevant to that hemisphere, then the gain of the neural response should be higher in the first case, but there should be no difference in the noise correlations. This is a testable prediction.

This is an interesting suggestion. At a technical level, the idea is correct, and holds up. We can artificially subdivide the trials (within either cue condition) depending on the value of the modulator. We thus create two subsets of trials with different means. If we make sure that the variance of the modulators is the same within the two groups, the noise correlations are very similar (albeit smaller than in the full, combined group). The results of this analysis are included in [Author response image 1](#fig9){ref-type="fig"}). This validates our argument that the variance of the modulator is a major determinant of the noise correlations.10.7554/eLife.08998.020Author response image 1.Suggested Analysis.Here, we show the results of the reviewers' suggested analysis. For each day and hemisphere, and within a given condition (cue towards, or cue away), we subdivide the data into two based on the values of the modulator (high modulator values, or low modulator values). We choose the threshold for the split so the variance of the modulator is the same for the two conditions. Top row: the mean firing rates are higher when conditioned on high modulator values. Bottom row: the noise correlations in the two subgroups are roughly equal (though still overall lower in the cue towards condition). Note that the subgroups' noise correlations are lower than the full dataset for that cue condition, as we have reduced the (conditional) variance in the shared modulator. While this validates our statement that the noise correlations are the result of the variance of the shared modulator (since artificially reducing the variance of the shared modulator causes a dramatic reduction in the noise correlations), it does not resolve the status of whether the modulator is signal or noise.**DOI:** [http://dx.doi.org/10.7554/eLife.08998.020](10.7554/eLife.08998.020)

At a scientific level, however, we don't believe this provides a strong test of whether the modulators themselves are signal or noise (i.e., whether they capture attentional fluctuations or representational fluctuations), which is the question of interest. Both these hypotheses would provide a complete and selfconsistent description of what we would be doing when we subdivide the trials. However, we were inspired by the spirit of this suggestion to come up with a far stronger test of the two hypotheses, which we present in the new [Figure 8](#fig8){ref-type="fig"}.

Another area where the paper should be improved is in explaining which effects are \"explained\" by the model and which ones are implied as a hypothesis. For instance, in the Discussion, the paper reports that neurons that showed stronger reward dependent changes in firing rates (gain) were also more strongly coupled to the modulator. This is similar to statements made earlier in relation to attentional modulation. This argument may be circular: given the fitting and the associated weights, it may be an inevitable consequence of the computational approach, and thus not a real result. The same applies to the argument in the subsection "Interpreting the model components", and as stated to wherever similar arguments are made before. These are potentially major criticisms and they absolutely have to be addressed.

Thank you for pointing this out. We agree that there were places where we did not properly distinguish assumptions from findings. We've made an attempt to address this throughout the manuscript. With regard to the specific example about the rewardcoupling of individual neurons, we decided that this was a relatively minor point, and unnecessarily confusing. As such, we've removed this point, and feel that the manuscript is better off without it. Please let us know if there are any remaining concerns.

*Finally, a major comment involves the analysis in [Figure 6](#fig6){ref-type="fig"}. This analysis employs the difference between* m~c~ *and* m~u~*. In fact, this seems like exactly the wrong thing to do in the current study, as it is shown that those two signals are not correlated. Shouldn\'t each hemisphere be analyzed on its own?*

Thanks for raising this point. The difference does capture much of the effect, but we went back and reanalyzed the two signals separately, and have now included this in the paper. We also simplified this figure (now [figure 7](#fig7){ref-type="fig"}) as it provided unnecessary detail to the understanding of the major result (this detail is now in [Figure 7--figure supplement 1](#fig7s1){ref-type="fig"}).

\[Editors\' note: further revisions were requested prior to acceptance, as described below.\]

The paper has improved considerably and is arguably ready to be published as is. It does a much clearer job at delivering a \'unifying account\' of prior results, in the light of new analyses. Where the prior version hedged its bets, the new version takes a clear stand and paints a compelling picture. However, it would be even stronger if the authors would address two main comments indicated below.

The first major comment is about something that is visible in the new version of [Figure 7A](#fig7){ref-type="fig"}, which shows the impact on behavioral performance of the inferred modulators -- one for the cued side and one for the opposite side. This is a better choice than in the previous submission, where they had expressed it as a function of the difference in these two modulators. However, this new analysis reveals something intriguing. The value of the modulator on the cued side predicts the ability to detect a change on the opposite side, even though it does not impact firing rates on the opposite side (top-right red quadrant of [Figure 7A](#fig7){ref-type="fig"}).This is remarkable, and seems to completely undercut the theory that the impact of attention on performance is linked to changes in response gain. For example, decreasing the cued-side modulator seems to increase performance on the opposite side, even though this would not increase firing rates in the opposite side, nor decrease their variance. This is because the two modulators are essentially uncorrelated ([Figure 6A](#fig6){ref-type="fig"}). If this reasoning is correct, the cued-side modulator has an impact on performance for opposite-side stimuli without having any detectable impact on the physiology of the opposite V4 neurons. Thus, the impact of the modulator basically can\'t be due to physiological effects at the level of V4.If this understanding of the data is incorrect, it would be a good idea to dispel it. If instead it is correct, it seems that it should be explicitly discussed. The manuscript stops short of this: it mentions that things cannot be fully explained by a \"1D axis of attention\", and the issue comes up again briefly in Discussion, but it is never really highlighted. Indeed, the finding seems inconsistent with standard models of attention, including the one presented in the paper. And it wouldn\'t be just about 2D vs. 1D: it would be about whether the physiological effects of attention in V4 have anything to do with the behavioral effects.

We think this interpretation is incorrect, and we had provided an explanation of this in the last section of the Results in the manuscript. The structure of the task is such that the animal must compare the represented stimuli in both hemispheres: it is a threealternative task (left target, right target, or no target). Thus, even though the modulators are uncorrelated, the value of a modulator on one side can affect performance on the other.

We appreciate that the explanation is a bit subtle, and it was easy to miss in the Results. To make this point more salient in the manuscript, we have added a few sentences to the Discussion (see end of paragraph beginning with "A second scenario...").

The second major comment is about the control analysis currently shown in the reply to the reviewers: dividing the data based on the value of the modulator (high \'attention\' versus low \'attention\') and assessing whether variability (Fano factor or noise correlations) is lower for the former (which it would be if attention suppresses other variability such as \'representational fluctuations\' as envisioned in the old view). It may be a good idea to leave this figure out of the paper, as the results are essentially a foregone conclusion if one accepts the results and interpretation of the prior figures (i.e., if one accepts that the modulator captures essentially all the relevant shared variability). However, some readers might still have doubts about whether the modulators truly account for all the shared variability (that conclusion depends on how suitable the model is) and/or may not immediately see the implications of that finding. It may therefore be a good idea to add that analysis to the supplementary figures.

Thank you for this question. This prompted us to realize that we have somehow failed to show, in a quantitative way, exactly how well the model accounts for the attentioninduced changes in variability (i.e. the Fano factor and the noise correlation). The plot in [Figure 2A](#fig2){ref-type="fig"} does indicate the overall performance of the model in capturing the responses, but this is in units of log likelihood (which are not easily interpretable). [Figure 5](#fig5){ref-type="fig"} demonstrates that we can predict reductions in variability, and that the magnitude of these effects depends on the strength of coupling to the hidden modulator. But what is missing is a quantitative statement of how well the model predicts the reductions in Fano factor and noise correlation.

Given this, we have generated yet another figure ([Figure 5B](#fig5){ref-type="fig"}), and we hope you will find that it satisfies your request. We think this is a more direct demonstration than the figure we sent with our previous responses, and we have included it in the main text.
