{'fraction_non_alphanumeric': 0.02959181261670094, 'fraction_numerical': 0.011510553310170413, 'mean_word_length': 4.978521972587255, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Partitioning neuronal variability

Responses of sensory neurons differ across repeated measurements. This variability is usually treated as stochasticity arising within neurons or neural circuits. However, some portion of the variability arises from fluctuations in excitability due to factors that are not purely sensory, such as arousal, attention, and adaptation. To isolate these fluctuations, we developed a model in which spikes are generated by a Poisson process whose rate is the product of a drive that is sensory in origin, and a gain summarizing stimulus-independent modulatory influences on excitability. This model provides an accurate account of response distributions of visual neurons in macaque LGN, V1, V2, and MT, revealing that variability originates in large part from excitability fluctuations which are correlated over time and between neurons, and which increase in strength along the visual pathway. The model provides a parsimonious explanation for observed systematic dependencies of response variability and covariability on firing rate.


Introduction
Neurons transmit information with sequences of action potentials. These responses are variable -repeated measurements under identical experimental conditions give different spike trains -but the origins of this variability are unknown. If spike generation were variable, it might account for response variability, but in vitro measurements indicate that it is highly reliable 1 . Variability in synaptic transmission is another possible source 2 , but its magnitude is also believed to be insufficient to account for the observed variability in spiking responses [3][4] . A more likely explanation is that the variability arises from the accumulation and amplification of small amounts of noise as signals flow through neural circuits 5 . And recent theories propose that the substantial variability in neural responses may arise from the dynamics of recurrent but largely deterministic networks 6-7 . Users may view, print, copy, and download text and data-mine the content in such documents, for the purposes of academic research, subject always to the full Conditions of use:http://www.nature.com/authors/editorial_policies/license.html#terms Regardless of its source, characterizing variability with simple stochastic models has proven useful in understanding the nature of neural coding. The simplest stochastic model is a Poisson process, in which spikes occur independent of one another. A hallmark of the Poisson model is that the variance of the spike count in any time interval is equal to the mean. In visual cortex, the spike count variance typically equals or exceeds the mean, but rarely falls below it 8,9 . This suggests that Poisson-like behavior is a "floor" state of cortical variability, and raises the question of the origin of the excess variance. Arousal, attention, adaptation and other contextual factors are known to modulate sensory responses [10][11][12] . In typical electrophysiological experiments, some of these may be well-controlled, but many are not. The idea that fluctuations in excitability can inflate estimates of neuronal variance has a long history 8,9 , and we wondered whether a more directed analysis of single-neuron responses might reveal the effect of these factors.
We formalize this hypothesis in a doubly stochastic response model in which spikes arise from a Poisson process whose rate is the product of "drive" and "gain" (the "modulated Poisson model", Fig. 1). The drive is a reproducible firing rate response to a sensory stimulus; the gain represents modulatory influences on excitability, and can vary across repeated measurements. Under this model, trial-to-trial variability in spike counts can be partitioned into a sum of Poisson point process variance and variance arising from fluctuations in gain. Likewise, spike count covariation can be partitioned into point process covariance, and covariance arising from correlated gain fluctuations.
We found that this model provides an excellent account of single-neuron response distributions in macaque visual thalamus and cortex. Inferred gain fluctuations are correlated over long timescales (minutes), are larger in anesthetized than in awake animals, and (in anesthetized animals) increase in strength along the visual hierarchy. They are also shared across wide areas of primary visual cortex, while point process variance is more localized. Together, these results suggest that much of the response variability of sensory neurons arises from fluctuations in excitability which are correlated over time and between neurons and which increase in strength along the visual pathway.
where μ is the mean spike rate, and Δt is the duration of the counting window. Assume that the rate arises from the product of two positive-valued signals: (2) where f(S) is some function of the stimulus, and G is a stimulus-independent gain. In this case, the variance and mean of the spike count in a time interval Δt are both equal to f(S)GΔt.
If the gain signal is unobserved, and varies over repeated trials, the observed spike counts over those trials will have a variance larger than predicted by the simple Poisson model. Assume the gain has a mean of one, but fluctuates randomly on a timescale that is slow relative to the interval over which spikes are counted. In this case, the net distribution of spike counts is a mixture of Poisson distributions (Online Methods) whose variance can be decomposed into two components: (3) The first term is simply the variance of a Poisson distribution with mean count f(S)Δt. The second term is the variance of the expected spike count conditioned on the stimulus drive 14 , and is proportional to the square of the first term, with a proportionality factor equal to the variance of the gain signal. Thus, in the presence of gain fluctuations, the spike count variance exceeds the mean by an amount that is proportional to the square of the mean. Note that the spike count variance does not depend on the firing rate or duration of the counting window per se, but on their product, f(S)Δt.
The variance-to-mean relationship expressed in Eq. (3) is inherent to systems with modulatory effects, and depends on gain only through its variance: all gain distributions with the same variance will produce the same variance-to-mean relationship. To fit the model to data, the form of the gain distribution must be specified. If we assume that the gain follows a Gamma distribution, the overall spike count distribution will follow a negative binomial (Online Methods): computed by counting spikes in a 1,000-ms window following response onset. We fit both the Poisson and modulated Poisson models to the responses of each individual neuron.
Consider the actual and model-predicted response distributions for an example V1 cell (Fig.  2a). The mean spike count depends on the direction of motion of the stimulus (center Fig.  2a), and the variance grows with the mean. The variance-to-mean relation is well described by the modulated Poisson model (Fig. 2b). As predicted (eq. 3), this relation does not depend on firing rate or measurement interval (inset Fig. 2b). Comparing histograms of the measured spike counts with the predictions of the models reveals that the modulated Poisson model captures the data much better than the standard Poisson model (Fig. 2a). When the mean spike count is low, both models produce similar response distributions with mode at zero and variance approximately equal to the mean (Fig. 2b). But as the mean spike count grows, model predictions diverge, and the Poisson model fails to account for the shape of the count distributions (Fig. 2a). In contrast, the data are well-described by the modulated Poisson model, suggesting that a substantial part of the response variance arises from fluctuations of modulatory inputs, whose distribution is well-described by a Gamma distribution.
This visual impression is confirmed by statistical analysis of the model fits. To quantify relative goodness of fit of the two models, we performed a cross-validation analysis (CVA), comparing the log-probability of observing the hold-out data under each fitted model (Online Methods). For the example neuron, the value for the modulated Poisson model is much higher than that of the Poisson model. To assess the model fits in absolute terms, we developed a parametric bootstrap test that compares the log-probability of the data with its expectation under the model in question (Online Methods). For the example neuron, the fit of the Poisson model is very poor (P < 0.001, absolute goodness of fit test; Fig. 2c), but the modulated Poisson model cannot be rejected (P = 0.91; Fig. 2c).
In sum, the variable discharge of this V1 cell is well described as originating from three different sources: the stimulus attributes (i.e., direction of motion), a Poisson point process, and Gamma-distributed fluctuations in excitability. To estimate the relative contribution of each source, we used the modulated Poisson model to partition the spike count variance (Online Methods). Surprisingly, Poisson noise accounts for only a small fraction of the total variance (5.5%). The gain fluctuations account for nearly half of the variance (47.5%), a share comparable to the fraction due to variations in the stimulus drive (47%). The latter is dependent on the set of stimuli and the tuning properties of the neuron. To focus our analysis on the variability across repeated measurements, we consider the portion of within-condition variance that is explained by the excitability fluctuations. For the example neuron in Fig. 2, this fraction is 89.6%.
In our model, stronger gain fluctuations lead to a more rapidly accelerating variance-tomean relationship, which deviates more and more from the Poisson expectation as spike count grows (Fig. 2d). Note that this accelerating relationship implies that the ratio of variance to mean (Fano factor) can vary within a single spike train, and over short time scales. As such, a single Fano factor provides an incomplete and potentially biased, measure of neuronal variability.
In the modulated Poisson model, the rate arises from the product of two positive-valued signals, one a function of the stimulus and the other arising elsewhere. An alternative model adds these signals rather than multiplying them. Under this additive model, if the two signals are statistically independent the expected spike count variance is given by the sum of the point process variance and a constant equal to the variance of the drive fluctuations, yielding a variance-to-mean relationship (Fig. 2d) which is unlike that in our data.
Fluctuations in excitability modulate neuronal activity throughout visual thalamus and cortex. Neurons in LGN, V1, V2, and MT all exhibit super-Poisson variability (Fig. 3a). Within each area, the variance-to-mean relationship is consistent with the predictions of the modulated Poisson model (compare Fig. 3a to Fig. 2d), but note that the fitted gain variance increases in strength along the visual processing stream. Specifically, the modulated Poisson model systematically outperforms the standard Poisson model for all areas (Fig. 3b), and the proportion of neurons for which this model is selected increases from LGN to V1 (P < 0.001, t test), from V1 to V2 (P = 0.004) and from V1 to MT (P = 0.004). The strength of fluctuations in excitability, as measured by the coefficient of variation of the gain, grows from LGN to V1 (P < 0.001, Wilcoxon rank sum test), from V1 to V2 (P < 0.001), and from V1 to MT (P < 0.001) (Fig. 3c). As information propagates through the visual hierarchy, fluctuations in excitability not only increase in strength, but also account for a larger share of variance. When stimulated with drifting gratings, within-condition variance primarily reflects Poisson-like noise in LGN (P = 0.02, Wilcoxon signed rank test), but becomes progressively more dominated by excitability fluctuations in cortex (LGN to V1, P < 0.001, Wilcoxon rank sum test; V1 to V2, P < 0.001, V1 to MT, P < 0.001; Fig. 3d).

Response correlations in the modulated Poisson framework
Trial-to-trial response fluctuations are often correlated among simultaneously recorded neurons 15 . Pairwise response correlations can arise when neurons receive shared sensory input 5 , but also when they are subject to correlated modulatory influences that are not sensory in origin 16,17 . Our model provides a vehicle for separating the effects of these two contributions. Specifically, for doubly stochastic processes, the spike count covariance can be decomposed in a manner that is analogous to the variance decomposition introduced in Eq. (3) (Online Methods): (5) where r Pij indicates the point process correlation (assumed to be independent of stimulus S), and r Gij the gain correlation. The first term is the covariance expected for the spike counts of neuron i and j with a constant gain of one. The second component is the covariance of the conditional expectations that arises from correlated gain fluctuations, and can generate spike count correlations even when the two point processes are independent. Note that Eq. (5) reduces to the expression for spike count variance (Eq. (3)) when i=j. The spike count correlation is obtained by dividing this equation by the square root of the product of the spike count variances of neurons i and j, as expressed in Eq. (3), yielding a complicated dependence of spike count correlation on stimulus drive, and on the variance of the gain signals in the two neurons.
Studies of response correlation typically combine normalized responses across conditions to estimate a single spike count correlation for a pair of neurons 18,19 . But equation (5) implies that even if the point process and gain correlations are both stable properties of a given cell pair, measured spike count correlations can vary dramatically with stimulus drive. Specifically, when the geometric mean of the stimulus-driven response of the two neurons is low, the spike count correlation will approximate the point process correlation r Pij In contrast, when this geometric mean is high, the spike count correlation will be dominated by the gain correlation r Gij If r Pij and r Gij differ, the spike count correlation will depend on the stimulus conditions, and any single estimate will provide an incomplete and potentially biased measure of neuronal covariability.
We analyzed the responses of four populations of simultaneously recorded neurons in the superficial layers of macaque primary visual cortex 20 . The recordings were made over 2.5 h, during which gratings drifting in 72 equally spaced directions were presented for 1,280 ms each, interleaved with a 1,280 ms blank screen and repeated 50 times in random order. We analyzed 379 well-isolated units (62, 94, 87, and 136 from the four data sets). We fit the modulated Poisson model to the responses of each neuron, and then fit the model to the spike count correlations of pairwise combinations of neurons in each data set (Online Methods).
Consider the actual and predicted response correlations for three example pairs of V1 cells (Fig. 4). For each neuron, the mean spike count depended on the direction of motion of the stimulus (Fig. 4a-c), and for each pair of neurons, the spike count correlation estimated from the combined normalized responses (r SC ) was small and positive. But estimates of spike count correlation computed separately for different response levels reveal a variety of different behaviors, including a decrease with mean response strength (Fig. 4d), independence of response strength (Fig. 4e), or an increase with response strength (Fig. 4f). The modulated Poisson model can mimic each of these behaviors using particular choices of the underlying point process and excitability correlations, r Pij and r Gij . We estimated both statistics for each cell pair, and Fig. 4g-i illustrate the full spike count correlation surface predicted by the fitted models.

The structure of correlations in primary visual cortex
V1 spike count correlations have been shown to depend on cortical distance and tuning similarity 19 , but as revealed by the examples in Fig. 4, this measure mixes two underlying sources of correlation. We wondered whether these two sources might be differently structured across the neuronal population. In particular, common point process variance might be local in space and time, while joint excitability fluctuations could reflect network changes that affect larger populations of neurons on a slower time scale.
We examined how point process correlation and excitability correlation depend on electrode separation and on the similarity of tuning (Fig. 5, Online Methods). Point process correlations were on average smaller than excitability correlations (P < 0.001, F 1,34 = 24.16, ANCOVA) and fell more rapidly with electrode distance (P < 0.001, F 1,34 = 18.05). Note that neurons recorded on the same electrode (the leftmost datum in Fig. 5a-b) might have inflated the significance of the latter difference since they are vulnerable to spike sorting errors that can alter the measured correlation 15,17,21 . Nevertheless, we find that excluding this datum does not alter the conclusion (P < 0.001, F 1,32 = 41.17). Analysis of tuning similarity effects (Fig. 5c) revealed that the average point process and excitability correlation had a largely similar (linear) dependency on tuning correlation (P = 0.6, F 1,14 = 0.28), differing only in their means (Fig. 5d).
In summary, we find that the relationship between response correlations and mean responses in V1 depends on both inter-neural distance and tuning similarity: Spike count correlations are most likely to decrease with response strength for nearby neurons with similar tuning (Fig. 4d). As inter-neural distance or tuning dissimilarity grows, response correlations will initially tend to be independent of the mean responses (Fig. 4e), and eventually exhibit an increasing relation (Fig. 4f).

The temporal structure of gain fluctuations
We assume that the Poisson point process variance is independent across trials, so all temporal structure in neural responses that is not explained by the stimulus must arise from structure in the gain signal. Consider how normalized spiking activity evolves over 2.5 h for three simultaneously recorded neurons (Fig. 6a). Response strength rises and falls over minutes, suggesting that excitability varies slowly. This is consistent with the slow decay of the autocorrelation of the gain (Fig. 6b), computed from temporal correlations in the trialby-trial responses (Online Methods). At a time lag of 5 sec, 307/379 neurons had positive autocorrelation (t-test, P < 0.05). For longer time lags, this fraction gradually decreased (at a time lag of 2 minutes, r > 0 for 185/379 neurons; at a time lag of 10 minutes, this fraction drops to 81/379). Although the average profiles for the four data sets differ in detail, they share a slow falloff, indicating that changes in excitability persist for many minutes (Fig.  6c).
Trial-to-trial changes in excitability are often correlated across neurons (Fig. 5). It might therefore be expected that the slow component of excitability fluctuations is also shared across neurons. However, across neurons, excitability correlations decrease dramatically at all time lags exceeding 0 sec (Fig. 6d). The slow component of excitability fluctuations thus appears to be local.

Effects of anesthesia
We have shown that neural response variability in visual cortex of anesthetized macaques originates largely from excitability fluctuations that are correlated over time. To what degree does this arise from fluctuations in the state of anesthesia? To address this, we analyzed data recorded in area MT of awake monkeys performing a psychophysical task involving random dot kinematogram stimuli 18,22 . We fit both the Poisson and modulated Poisson model to the responses of 307 individual neurons. As in anesthetized animals, MT neurons in awake monkey typically exhibit super-Poisson variability (Fig. 7a), and the modulated Poisson framework successfully accounts for this variability. The inclusion of gain fluctuations improves goodness of fit for 224/307 neurons (CVA analysis). At a significance level of 5%, the Poisson model is accepted for only 28/307 of neurons, while the modulated Poisson model is accepted for 294/307 of cells (absolute goodness of fit test). Fluctuations in excitability were the main source of within-condition variance (P < 0.001, Wilcoxon signed rank test; Fig. 7b). Finally, the average autocorrelation function reveals that gain fluctuations in awake animals exhibit a slow temporal structure similar in timescale to that seen under anesthesia (Fig. 7c).
Qualitatively, there is thus an excellent agreement between these data sets. However, comparing the coefficient of variation of the excitability signal reveals that its fluctuations are significantly stronger in the anesthetized cortex (compare Fig. 3c with Fig. 7b). It is also notable that the time course of this correlation is substantially shorter in the awake data set, suggesting that the fluctuations in gain seen under the two conditions may have different origins.

Discussion
The distinction between inputs that drive the responses of a neuron and those that modulate its response is well known 23 . But while both factors influence the mean response, response variance is usually treated as if it arises solely from noise in driving inputs. Our analysis suggests instead that variability of modulatory signals may underlie much of the response variability in visual cortex. We have shown that a Poisson spiking model whose rate is multiplied by a fluctuating gain signal can explain why the variance of spike count generally grows faster than the mean (Figs. 2 and 3), and that an additive noise model cannot account for this basic behavior (Fig. 2d). The modulated Poisson model also explains why the covariance of spike counts in pairs of neurons can exhibit a diversity of behaviors depending on their means (Fig. 4). Fitting the model to individual neurons from a variety of visual areas reveals that gain fluctuations account for a substantial fraction of spike count variability, and that in anesthetized animals this fraction increases as one ascends the visual hierarchy (Fig. 3).
Poisson spiking models have been common in neuroscience since the 1950's [24][25][26] , and are implicit in all analyses in which responses are summarized solely with mean spike counts. They provide the simplest statistical description of the data, are readily fit to data, and capture the basic fact that spike count variances grow with the mean. Poisson processes can mimic the spiking behaviors of integrate-and-fire models 4 , as well as the responses of model neurons embedded in cortical networks with balanced configurations of excitatory and inhibitory inputs [5][6][7] . But neuronal responses generally exhibit super-Poisson variability 5,7-9,27-28 , and the model presented here provides a simple but effective means of explaining this behavior. Some authors have reported sub-Poisson variability, especially in brief time bins immediately after a sudden stimulus onset [29][30][31][32] . This is likely a consequence of neuronal refractoriness [33][34] , which could be naturally incorporated into our modeling framework through modulatory spike-feedback terms in a generalized linear model [35][36] .
Our work complements recent work on the temporal dynamics of firing rate variability 14,28 . One study used the mean-matched Fano factor as a proxy for firing rate variability, but did not explore the dependence of this value on the mean 28 . Another estimated firing rate variability by subtracting an estimate of the point process variance from the measured spike count variance 14 . The estimate of point process variance is taken from the minimum of an ensemble of experimental measurements, and is therefore vulnerable to the well-known problems associated with estimating extreme values of sample distributions. In contrast, we found it useful to make explicit distributional assumptions: specifically, that gain is drawn from a Gamma distribution, and that the spikes arise from a Poisson process. The resulting negative binomial spike count distribution (a particular continuous mixture of Poisson distributions) is easily fit to data and validated (Fig. 2d, see also ref. 37), and its parameterization as a modulated Poisson distribution provides a natural interpretation in terms of excitability fluctuations. The spike count distribution bears some resemblance to the discrete mixture-of-Poissons model that has been successfully used to account for super-Poisson variability 38 , but has the added advantage that it has far fewer parameters, and provides an interpretation in terms of modulatory input. It is important to note, on the other hand, that our analyses of joint behaviors are based solely on correlations. A full model for joint neural responses would allow better fitting and validation, but requires an explicit description of both the joint point process distribution, as well as the joint distribution of gain signals.
Our analyses indicate that fluctuations in gain account for a substantial portion of spike count correlation. This has been long-suspected, but it has been difficult to segregate the effects of shared sensory inputs from correlated modulatory influences 16 . Our model provides an explicit method for achieving this (Eq. 5). In V1, this analysis reveals that point process correlation and excitability correlation have different structure. Point process correlations decrease rapidly with cortical distance, while excitability correlations change less. Both point process and excitability correlations increase with tuning similarity (Fig. 5). Together, these patterns suggests that point process correlations are caused by inputs that are shared within local functional circuits, while excitability correlations likely arise from fluctuations in modulatory signals that affect larger populations of neurons (e.g., sensory adaptation, metabolic resource availability, attentional signals, reward signals).
The analysis in Fig. 3 indicates that the portion of spike count variability attributable to modulatory fluctuations increases as one ascends the visual hierarchy. This presumably depends on the context in which responses are measured, including the choice of stimuli, the presence or absence of anesthesia, and the cognitive state of the animal. For example, attentional mechanisms are widely believed to act by increasing the gain of neurons involved in a task. Relative to conditions of uncontrolled attention, we might expect that this would lead to a net increase in the mean, but a decrease in the variance, of the modulatory input. Under the modulated Poisson model, this would result in an increase in spike count, accompanied by a decrease in Fano Factor, as has been reported in area V4 (ref. 39). The corresponding predictions for spike count correlation are more complicated, since the relative contributions of the point process and gain correlations depends on the stimulus drive to each cell, as well as the variance of the gain of each cell (Eq. 5). But if the gain correlation is larger than the point process correlation and the two cells are receiving similar stimulus drive, our model predicts that a reduced modulatory variance would decrease spike count correlation, as has also been observed 39 . The modulated Poisson model has broad implications for the characterization of neuronal response variability. It is common practice to summarize variability with the Fano factor (the ratio of the spike count variance to its mean). As clearly revealed in Figs. 2-3, this measure is highly dependent on the conditions under which it is estimated. For example, the Fano factor for the neuron in Fig. 2b would be near 1 if estimated from trials in which the response was low (less than 1 spike per time bin), but would be significantly higher if estimated from trials with larger responses. Thus, the Fano factor, by itself, does not provide a reliable summary of neural variability. In contrast, spike count variability of the modulated Poisson model may be decomposed into a sum of two distinct values, one corresponding to the Poisson component (which is always equal to the mean) and the second arising from unobserved modulatory influences, which grows as the square of the mean (Eq. 2). The amplitude of this second term, which corresponds to the variance of the modulated gain, provides a compact summary of variability beyond that expected from the basic Poisson model.
An analogous issue arises with the common practice of summarizing the covariation of pairs of neurons with spike count correlation. This measure can depend critically on the conditions under which it is measured (Fig. 4), and thus provides a potentially biased summary of covariability. Indeed, previous studies have noted that neural correlations can increase with firing rate 15 , and that this behavior can be explained by models in which spikes arise when a Gaussian-distributed membrane voltage crosses a threshold 15,40-41 . Our model is quite different but accounts for this phenomenon, as well as the more diverse range of behaviors seen across different cell pairs (Fig. 4), by summarizing covariation with two distinct values. One value represents the point process correlation (which could arise from noise in common sensory inputs 42 ), and the other the correlation of the modulatory influences. At low firing rates, response correlations primarily reflect the former, whereas at high firing rates they reflect the latter. This insight suggests a resolution for a recent controversy in the literature regarding the nature and magnitude of cortical correlations. Despite decades of experimental evidence that cortical cells exhibit modest but significant spike count correlations on the order of 0.1-0.3 (ref. 15), a recent study in macaque V1 found correlations statistically indistinguishable from zero, and concluded that previously reported values arose from experimental confounds 17 . However, average firing rates were unusually low in this study, and measured correlations were significantly positive for the small subset of neuronal pairs with high responses. The interpretation arising from our model and data (Fig. 4) is that the correlation values reported in ref. 17 primarily reflect point process correlations, which are dominant at low firing rates and fall rapidly with cortical distance, whereas most previously reported values reflect correlated modulation, which is only evident at high firing rates, but is generally more substantial and falls more slowly with cortical distance.
Our analysis suggests that fluctuations in gain are correlated over long time scales within but not across neurons (Fig. 6). This implies that the mechanisms underlying slow drifts in response gain differ from the mechanisms that give rise to instantaneous gain correlations. The local nature of these drifts rules out that they result simply from global state changes induced by factors like anesthesia or arousal. That said, comparison of responses in area MT of behaving and anesthetized macaque revealed that cortex is less stable under anesthesia: Fluctuations in gain are significantly stronger (Fig. 3c and 7b) and slower (Fig. 6c and 7c) in the anesthetized cortex.
What are the implications of our findings for understanding the representation of sensory information in the brain? We believe that the gradual increase in the strength of modulatory fluctuations along the visual pathway (Fig. 3) reflects the gradual transformation and combination of visual signals with information from other sources, such as other sensory inputs, top-down signals representing attention, arousal, metabolic state, reward expectations, emotional state, and so forth. Regardless of the source of modulatory variability, its increasing weight as information ascends cortical pathways raises the question of why it does not overwhelm the sensory information encoded in higher-level visual areas 5 , with deleterious effects on the accuracy of behaviors that arise from decoding those neural responses [43][44] . It is perhaps worth noting that information encoded in the relative responses of groups of neurons will not be affected by common modulatory changes, such as those resulting from attentional focus 11 . Resolving this mystery is a critical step in understanding how stable perception of the visual world can arise from apparently volatile neural activity, a step made easier by the principled and proven analysis that we have presented here.

Surgical preparation
We recorded from 18 anesthetized, paralyzed, adult macaque monkeys of either sex (3 Macaca nemestrina, 1 Macaca mulatta, and 14 Macaca cynomolgus). Our standard procedures for the surgical preparation of animals and single-unit recordings have been reported in detail previously 45 . Briefly, experiments typically lasted 5-6 d, during which we maintained anesthesia with infusion of sufentanil citrate (6-30 μg kg -1 h -1 ) and paralysis with infusion of vecuronium bromide (Norcuron; 0.1 mg kg -1 h -1 ) in isotonic dextrose-Normosol solution. We monitored vital signs (heart rate, lung pressure, end-tidal pCO 2 , EEG, body temperature, urine flow, and osmolarity) and maintained them within appropriate physiological ranges. Pupils were dilated with topical atropine. The eyes were protected with gas-permeable contact lenses, and refracted with supplementary lenses chosen through direct ophthalmoscopy. At the conclusion of data collection, the animal was killed with an overdose of sodium pentobarbital. All procedures were conducted in compliance with the National Institute of Health Guide for the Care and Use of Laboratory Animals, and with the approval of the New York University Animal Welfare Committee.

Unit recording
Extracellular recordings were made with quartz-platinum-tungsten microelectrodes (Thomas Recording), advanced mechanically into the brain through a craniotomy and small durotomy. Electrode insertion angle and location varied across experiments, depending on the targeted area. We distinguished V1 from V2 on the basis of depth from the cortical surface and changes in the receptive field location of the recorded units. Area MT was identified from the brisk direction-selective responses of isolated neurons. We made recordings from every single unit with a spike waveform that rose sufficiently above noise to be isolated. Stimuli were presented in random order. Data are reported from every unit for which we completed at least 5 repetitions.

Visual stimulation
We presented visual stimuli on a gamma-corrected CRT monitor (Eizo T966; mean luminance, 33 cd/m 2 ) at a resolution of 1,280 × 960 with a refresh rate of 120 Hz. Stimuli were presented using Expo software (http://corevision.cns.nyu.edu) on an Apple Macintosh computer. For each isolated unit, we first determined its ocular dominance and occluded the non-preferred eye. We presented circularly windowed sinusoidal grating stimuli to map each cell's receptive field, determined its preferred size and speed, and then measured selectivity for orientation or spatial frequency.

Analysis of single electrode recordings
Responses were computed by counting spikes in a 1,000-ms window following response onset. We estimated latency for each cell by choosing the latency that maximized the variance of the tuning curve 46 . The modulated Poisson model describes a doubly-stochastic process (also known as a "Cox process" 47 ), in which spiking responses are generated by a Poisson process whose rate is the product of drive and gain (Eqs. 1 and 2). To make the model identifiable, we assume that the gain, G, is constant within trials, and is distributed across trials according to a Gamma distribution with a mean of one and variance of : with shape parameter , and scale parameter . The mean spike count of the modulated Poisson model is the product of the gain, G, the stimulus driven rate, f(S), and the bin duration, Δt, and is thus also distributed according to a Gamma distribution, with parameters and . The spike count distribution is a Gamma mixture of Poisson distributions, and marginalizing over the gain variable yields a negative binomial distribution (Eq. 4) for the spike count 48 : With parameters r and s described above, this distribution has a mean of rs = f(S)Δt and a variance of as expressed in Eq. (3). The negative binomial can also be derived as a Poly-gamma mixture of Gaussian distributions 49 , a form that allows efficient inference 37 , but lacks the modulatory interpretation provided here.
We fit both the ordinary and modulated Poisson models to the responses of each individual neuron. The maximum likelihood estimator of the parameter μ that characterizes a Poisson distribution is the sample mean, which we estimate separately for each stimulus condition. The modulated Poisson model has two parameters (r and s) that are related to the stimulus drive and the variance of the gain. We used a simplex algorithm (the Matlab function 'fminsearch') to search for the value of {r, s} that minimized the negative log likelihood (NLL) of the full set of observed responses.
To assess the models' goodness-of-fit, we performed a 100-fold cross-validation. We fitted both models to a data-set consisting of all trials except for one randomly chosen trial per stimulus condition and compared the average log-probability of the hold-out data under each model (Fig. 3b). To evaluate the models' absolute goodness-of-fit, we compared the logprobability of the data with that of a large number of simulated data sets drawn from the fitted model. Each simulated data set has the same number of repeats as the observed data. We consider a model fit acceptable if the log-probability of the real data lies within the central 95% of the distribution of log-probability of the simulated data. This test can detect both under-dispersion and over-dispersion.
We used the modulated Poisson model to estimate the fraction of spike count variance that arises from the gain signal. As is standard in ANOVA, one can partition the sum of squares into components arising from variations in the stimulus (S stim ), the gain signal (S gain ), and the point process (S pp ): where N k is the spike count on the kth trial, N̄k is the spike count averaged over those trials in which the presented stimulus was the same as that of the kth trial, and N̄ is the spike count averaged over all trials. The second line follows from the first given the spike count variance of the modulated Poisson model, as expressed in Eq. (3). The fraction of within condition variance that arises from excitability fluctuations is given by .

Analysis of multielectrode array recordings
Full details regarding these data sets may be found in ref. 20. Briefly, an array containing 96 fixed electrodes was used to record from single units in the superficial layers of macaque primary visual cortex. The animal was anesthetized and stimulated with sinusoidal gratings. We included all units that could be tracked throughout the entire experiment and whose mean response exceeded 2 spikes/sec for at least one stimulus condition. We opted to leave out one set of array recordings (data set 3 in ref. 20) because the responses exhibited uncharacteristically strong correlations across time and neurons; we suspect the cortical surface was still recovering from array insertion. In our modeling framework, pairwise response correlations can result from both correlated point processes, as well as correlated gain fluctuations. To separate these two sources, we used the covariance decomposition formula in Eq. (5) and searched for the point process correlation and gain correlation that maximized the likelihood of the observed response relations. Rather than fitting response correlations directly, we performed the optimization on Z-transformed spike count correlations: where r is the spike count correlation, and ln the natural logarithm. The advantages of this Fisher transformation include a quicker convergence to normality and variance stabilization. We used a simplex algorithm to find the [r Pij , r Gij ] values that maximized the likelihood of the observed response correlations under a normal residual distribution. This approach works well for most neuronal pairs (Fig. 4), but the model parameters are not always well constrained. For instance, when the modulation of the direction tuning curve is weak for both neurons, many parameter combinations will yield a similar likelihood. For the population analysis (Fig. 5), we therefore took the uncertainty associated with the parameter estimates into account. We modeled the population distribution of the mean parameter values with a bivariate normal distribution, and searched for the parameter values that maximized the likelihood across all neuronal pairs. Poorly constrained pairs of neurons have a flatter likelihood function and are therefore down-weighted in this analysis. Although we consider this approach preferable, simple averaging of the parameter estimates yielded similar results.
To study the temporal structure of the gain signal, we estimated its autocorrelation function. We first removed the stimulus-induced temporal structure by subtracting the appropriate stimulus-elicited mean response from the observed responses. We then estimated the remaining covariance between pairs of responses separated by different time lags. To obtain the autocorrelation of the excitability signal at a given time lag, this covariance is normalized by the variance of the excitability signal, S gain /(N − 1).

Figure 1.
The modulated Poisson model. Spikes are generated by a Poisson process whose rate is the product of two signals: a stimulus-dependent drive, f(S), that is under experimental control, and a gain signal, G, that summarizes the net effect of stimulus-independent modulatory inputs that are assumed to fluctuate slowly relative to the duration of experimental trials.