Attentional fluctuations induce shared variability in macaque primary visual cortex

Variability in neuronal responses to identical stimuli is frequently correlated across a population. Attention is thought to reduce these correlations by suppressing noisy inputs shared by the population. However, even with precise control of the visual stimulus, the subject’s attentional state varies across trials. While these state fluctuations are bound to induce some degree of correlated variability, it is currently unknown how strong their effect is, as previous studies generally do not dissociate changes in attentional strength from changes in attentional state variability. We designed a novel paradigm that does so and find both a pronounced effect of attentional fluctuations on correlated variability at long timescales and attention-dependent reductions in correlations at short timescales. These effects predominate in layers 2/3, as expected from a feedback signal such as attention. Thus, significant portions of correlated variability can be attributed to fluctuations in internally generated signals, like attention, rather than noise.

N euronal responses to repeated presentations of identical stimuli are highly variable 1 . This trial-to-trial variability can be correlated across populations of neurons [2][3][4] and is often referred to as "noise correlation" 5 . Many studies have investigated the implications of these correlations for population coding 4,[6][7][8][9][10] . However, the origin of these correlations is still not clear. Here we focus on this latter question: what causes noise correlations?
One factor modulating correlations is attention. Studies of population activity in V4 found that attending to a stimulus inside the receptive fields of the recorded neurons reduced correlations in the trial-to-trial variability of the responses of those neurons to identical stimuli, compared to conditions in which attention was directed away from the receptive field 11,12 . These studies concluded that increasing the strength of attention reduces correlated variability by suppressing the shared, noisy input sources thought to give rise to correlated variability in a population 3,4,13 . This perspective on the relationship between correlated variability and attention is illustrated in Fig. 1a.
However, because the subject's state of attention can be controlled only on average but not precisely across trials, the strength and focus of attention may vary from trial to trial even within a given attention condition 14,15 . Here, we refer to such variability as fluctuations in the attentional state. Therefore, shared neuronal variability could also be driven by variability in the state of attention and changes in the level of that variability over time 8 . Indeed, the patterns of shared variability induced by fluctuations in gain-modulating signals such as attention are consistent with experimental data 8,16 if attentional state variability decreases as the strength of attention increases (Fig. 1b).
In other words, correlated variability during attention tasks can be interpreted as evidence for both a suppression of common noise by attention 11,12,17 as well as trial-to-trial fluctuations of attentional state 8,14,15 . Thus, it is unknown to what extent fluctuations in the state of attention indeed contribute to correlated variability in population responses, because the paradigms employed in these studies did not manipulate the level of attentional state variability behaviorally.
Therefore, we developed a novel, cued change-detection task that can dissociate changes in the strength of attention from changes in the variability of the attentional state by manipulating the behavioral relevance of two simultaneously displayed stimuli across task conditions. When only one stimulus is behaviorally relevant, subjects can maximize reward by focusing their attention on a single spatial location over time. However, when two stimuli are relevant, subjects need to attend to both stimuli to some degree. We expect attentional fluctuations to be highest in this latter scenario, if subjects shift the focus of attention between the two stimulus locations, as supported by recent work 18,19 .
Thus, if the dominant factor governing levels of correlated variability is attentional suppression of common noise, we expect correlations to decrease as attentional strength increases, resulting in intermediate levels of correlations when both stimuli need to be attended (Fig. 2a). Alternatively, if fluctuations in attention are the dominant factor modulating correlations, we predict correlations to be highest when both stimuli need to be attended and attentional fluctuations are most pronounced (Fig. 2b) 8 .
We recorded neuronal responses from primary visual cortex of macaque monkeys while they performed this task and find that attention modulates firing rates of V1 neurons. On a timescale of one second, we find that shared variability is highest when both stimuli are behaviorally relevant and lowest in conditions in which only one stimulus is the focus of attention, arguing that, at this timescale, fluctuations in the state of attention, induced by changes in attentional allocation strategies, are an important factor governing shared neuronal variability. On a faster timescale of 200 ms, we find attention-dependent reductions in correlated variability consistent with previous studies. Both effects . Correlated variability is driven by a common noise source (top right), which is suppressed by attention 11,12 . b Hypothesis 2: Attentional gain is increased, but fluctuates from trial to trial 8,14,15 . In this case, we expect attention to switch randomly between the two targets in the "Attend Both" condition, resulting in the highest correlations in this condition predominate in supragranular cortical layers, as expected from a feedback signal such as attention [20][21][22][23] .

Results
Change detection task and manipulation of attention. We trained two rhesus macaque monkeys to perform a cued, orientation-change detection task (Fig. 3a). A trial was initiated when the subject fixated a central fixation spot. Two "noisy" Gabor patches appeared symmetrically in the lower left and lower right visual field 300 ms later. During the zero-coherence period (ZCP), these patches randomly changed their orientation every frame (10 ms per frame; 36 orientations evenly spaced between 0 and 175 degrees). After a random period of time, drawn from an exponential distribution (minimum: 0.01 s, mean: 2.17 s, maximum: 5 s), one of the two stimuli entered the Coherent Period (CP). During the CP one particular orientation, called the "signal" orientation, was shown with a higher probability than the other orientations. By varying this probability, we could control the "coherence" of the stimulus, making the occurrence of the signal orientation more or less obvious over the background orientation noise, to manipulate the difficulty of a trial. The occurrence of this signal orientation was the change the monkey had to detect, which he reported by making a saccade to the changed stimulus within a short reaction time window. On 10% of trials no signal orientation occurred, and the monkey was rewarded for maintaining fixation throughout the trial. We used a cued block design to manipulate the focus of the subject's attentional state (Fig. 3b), where the cue was the color of the fixation spot. Two of these conditions, "Attend In" (AI) and "Attend Out" (AO), were similar to those in typical spatial attention tasks, where the stimulus overlapping the neurons' receptive fields is cued in the AI condition, and the other stimulus is cued in the AO condition. The cues for these conditions (red for AI, blue for AO) were 100% valid, such that the change occurred only at the cued location. In the condition labeled "Attend Both" (AB), indicated by a black fixation spot, either stimulus had an equal probability (50%) of showing the change on a given trial. Our paradigm therefore differs from typical covert attention tasks used to study neuronal variability in two respects. First, during the AI and AO conditions in our task, there are no catch trials with invalid cues 11 or signals in the distractor that need to be ignored 17 . While catch trials are typically used to measure the behavioral shift due to attention, they are likely to induce attentional fluctuations, as they render the cue unreliable and encourage some degree of attentional focus on the non-cued stimulus by rewarding successful performance at that location. As our goal in the AI and AO conditions is to minimize attentional fluctuations, we used 100% reliable cues. In our AB condition, either stimulus was equally likely to change. We used this condition as the baseline to measure the behavioral improvement attributable to attention, analogous to how other paradigms use catch trials.
There were, therefore, three attentional conditions but two attentional strategies that our task engaged. To maximize reward in the AI and AO conditions, attention should be focused on only the cued stimulus. With attention deployed consistently across trials with regard to spatial location, attentional state fluctuations should be minimized. In the AB condition, attention should fluctuate more strongly between the two spatial locations across trials, as ignoring one of the stimuli is no longer a viable strategy for maximizing reward. One way to conceive of this allocation strategy is that the AB condition is comprised of a mixture of the attentional states deployed in the AI and AO conditions. Note, attentional state fluctuations need not be non-existent in the AI and AO conditions but only decreased relative to the AB condition in order to test our hypothesis.
If subjects used the strategies described above, there should be some trials in the AB condition where the subject attended the unchanged stimulus and required a higher coherence level to notice a change in the correct stimulus on that trial. Such occurrences would lead to a rightward shift in the psychometric function and higher detection thresholds in the AB condition. The example session in Fig. 3c Fig. 1).
Overall, our goal was to develop a behavioral paradigm in which attention could fluctuate or shift between two stimulus locations-the AB condition-and remain focused on one location in the other conditions. Recent work suggests that attention is likely to operate in this fashion in the AB condition 18,19 , and our behavioral results, particularly those pertaining to psychophysical threshold, are consistent with this attentional allocation strategy. However, these results are also consistent with a strategy in which attention acts as a zoom lens 24 , widening its focus to encompass both stimuli simultaneously. Note, the fact that detection thresholds are elevated in the AB condition suggests that if attention is allocated to both stimuli simultaneously, the stimuli are not processed to the same degree as they are in the AI or AO conditions. That is, widening the attentional field entails a reduction in attentional strength within the field. As we will see, however, these strategies make different predictions for the patterns of correlated variability we expect to see across our task conditions.
Attentional modulation of neuronal firing rates. While subjects performed the task, we recorded spiking responses from neurons in primary visual cortex using 32-channel silicon probes with a spacing of 60μm between channels (NeuroNexus V1x32-Edge-10 mm-60-177). We recorded 474 single units (15.8 ± 1 units per session) across 30 sessions (N = 7 from Subject B, N = 23 from Subject D) from two male macaque monkeys. The two Gabor stimuli in our task were placed symmetrically in the lower visual field with one stimulus covering the receptive fields of the recorded neuronal population. Given the laminar nature of our recordings, receptive fields overlapped almost completely.
Our highly dynamic stimulus drove neurons strongly, with mean firing rates of 22.4 ± 0.9 spikes/sec across sessions. Consistent with previous studies we found that attention increased firing rates of V1 neurons 25,26 , with on average~31% of single units being significantly modulated by attention in a given session. This modulation was present in both the AI and AB conditions and appeared strongest early in the ZCP (Fig. 4a, b).
Note, our dataset contains fewer trials of long duration, given the exponential distribution of ZCP lengths and a slight tendency of subjects to prematurely abort longer trials (only~40% of valid trials are longer than 1 s, and~15% are longer than 2 s). We thus focused our analyses on the first second after stimulus onset, in which attentional modulation of firing rates was strongest, and on correct trials, where we can have the most confidence that attention was oriented as desired in our task. Additionally, all analyses of firing rates and spike counts were performed during the ZCP, before any changes in stimulus coherence or behavioral responses were made, ensuring that analyses were performed on identical stimuli across conditions.
We first calculated fractional firing rate increases in the AI and AB conditions, relative to the AO condition (Fig. 4c). During this interval, firing rates in the AI and AB conditions were significantly elevated relative to the AO condition (AI: 5.4 ± 1% increase, t(29) = 5.2, p = 0.00001, Bonferroni-corrected t-test, α = 0.0167; AB: 4.1 ± 1%, t(29) = 4.1, p = 0.0003) but not different from each other (t(29) = 1.4, p = 0.17). Amongst the roughly 31% of units showing significant modulation of firing rates by attention, around 32% showed pure gain modulation, around 20% showed pure offset modulation, while the remainder exhibited a mixture of multiplicative and additive modulation. Examples of pure gain-versus pure offset-modulated cells are shown in Fig. 4d. Note, these tuning curves were fit in a manner that assumed preferred orientation and tuning width did not vary as a function of attention condition 25 (see Methods for further details).
Differentiating the effects of attention on shared variability. Our results so far, beyond demonstrating that our task engages attention, are consistent with two different attentional allocation strategies in the AB condition, while we conclude that attention is primarily focused on the single, relevant stimulus in the AI and AO conditions. The first strategy involves widening the focus of attention to encompass both stimuli. In this case, we would expect attentional fluctuations to be negligible. This scenario would support the interpretation that attention suppresses a common noise source 11,12 , and we would expect correlations to be intermediate in the AB condition (Fig. 2a). The second strategy involves shifting the focus of attention randomly between the two stimuli. In this case, we would expect correlations to be highest in the AB condition (Fig. 2b). Note that this scenario does not rule out the possibility that attention suppresses a common noise source, as both mechanisms could be at play. However, given that the same dataset has been interpreted as evidence that attention suppresses noise 11 and that attention fluctuates 14 , it is an important question to quantify to what degree attentional fluctuations induce trial-to-trial variability.
Attentional modulation of shared variability. To measure the degree to which attentional fluctuations induce trial-to-trial variability, we calculated pairwise spike count correlations over repeated presentations of identical ZCP sequences in each attention condition. Our results match the predictions in Fig. 2b and support the hypothesis that fluctuations in the state of attention are the dominant factor inducing shared neuronal response variability in our dataset (Fig. 5a). Spike count correlations were significantly modulated by attention condition (F (2,29) = 15.1, p = 5e-6, rmANOVA), correlations were highest in the AB condition (t(29) = 5.7, p = 4e-6, t-test, see methods), and correlations in the AI and AO conditions were not significantly different from one another (p = 0.8, post-hoc Tukey's test). This relationship held individually for both subjects ( . We believe this result is due to a lack of statistical power, because the expected effect size for Fano factors is smaller than that for the correlation coefficients.
Next, we wanted to investigate the timescale of the correlation effect we found, to better understand its origin. Synaptic processes unfold on the millisecond scale whereas cognitive processes, such as attention, unfold over longer timescales. Behavioral work suggests that voluntarily shifting attention between different stimuli takes on the order of several hundred milliseconds 18,19,27,28 . Thus, if attention is indeed shifting between the two stimulus locations during the AB condition, these psychophysical results provide a lower bound for the timescale over which we expect to see correlations rise in the AB condition.
Using the relationship between spike count correlations and cross-correlograms, described in Bair et al. 3 and modified in Ecker et al. 29 , we calculated spike train cross-correlograms for neuronal pairs in each attention condition and integrated them from 1 to 1000 ms, our maximum counting window. Examining the point at which the resulting correlation levels saturate provides an estimate of the timescale of correlation. The results in Fig. 5c show that correlations in the AB condition began to diverge from the AI and AO conditions after 200 ms, and correlations in the AI and AO condition saturated to similar levels near 400 ms, while AB correlations continued to rise for several hundred milliseconds more. The time course of these results fits well with the estimated time course of changes in attentional state 18,19,27,28 . Interestingly, between 40 and 400 ms, the level of correlations appeared lower in the attended versus unattended conditions (Fig. 5c), consistent with earlier work 11,12,17 , suggesting that attention may indeed suppress common noise at this faster timescale. However, despite being consistent with previous results, this trend was not statistically significant for our overall dataset (F(2,29) = 1.8, p = 0.18 at 200 ms, rmANOVA).
It is worth pointing out here that our analyses in this paper focus on a set of recording sessions in which the two stimuli were horizontally separated from one another by at least 6°(that is, each stimulus was at least 3°from monitor center on the horizontal axis; see Methods for details). We also recorded some sessions in which the stimuli were closer to the vertical meridian. In these sessions, we failed to observe our predicted effect. We reasoned that this lack of effect was likely because the two stimuli were too close to each other, allowing the monkey to attend to both simultaneously. Indeed, the difference between correlations in the AB condition and the average of AI and AO increased as the two stimuli were further separated from one another ( Fig. 5d; To verify that this effect was not a false positive due to post-hoc analysis, we collected an independent 10session dataset at high eccentricities from Subject D, which confirmed the effect (Fig. 5d squares; see Methods for details).
Laminar profile of attention effects. To examine the laminar profile of the attentional modulation of firing rates and shared variability, we calculated the current source density (CSD) 30 across channels for each session from the task-stimulus evoked local field potentials (Fig. 6a). These profiles were quite consistent across sessions, with the most prominent stimulus-evoked sinksource configurations in L5-6 and L1-2/3, largely washing out the earliest sink-source switch typical of the L4-5 boundary (van Kerkoerle et al. 31 report a similar effect). We computed CSDs to aid in the grouping of single units into the supragranular (S), granular (G), or infragranular (I) layers, but we also took advantage of known electrophysiological characteristics of cells in different layers 32 . The most reliable such property was the high spontaneous activity associated with L4C 32 , which was readily discernible from multi-unit activity and was located consistently close to the L4-5 boundary determined from the CSD. Additional factors included the weaker orientation tuning of the deep granular layer and smaller receptive fields (Fig. 6a). The first channel below the L4-5 boundary was our zero-point for relative unit depths. We defined the granular layer as the first 400 μm superficial to the L4-5 boundary, consistent with previous histological 33,34 and recent electrophysiological studies 35,36 . All units above this 400 μm band were labeled supragranular, and all those below it were labeled infragranular. The G-I (L4-5) boundary could be determined most reliably across sessions, but the S-G boundary could not always be determined as precisely. We therefore varied the cut-off boundary between the supragranular and granular groups over a span of nearly 200 μm and recalculated the results presented in Fig. 6. Doing so did not qualitatively affect our results.
Considering the consistency of the finding in previous studies that correlations are reduced in attended conditions, at least at shorter timescales, and the trend we observed at such timescales when not conditioning on laminar position (Fig. 5c), we analyzed correlations at a 200 ms interval by laminar position as well (Fig. 6d). In the supragranular group, correlations were significantly modulated by attention condition (F(2,29) = 3.5, p = 0.036, rmANOVA), and consistent with previous studies, correlations were lower in the AI condition relative to the AO condition (t(29) = 2.9, p = 0.007, t-test). Fixational eye movements cannot account for our results.
Fixational eye movements, also called micro-saccades, have been reported to modulate neuronal activity in the visual system 40,41 , contribute to neuronal response variability 42,43 , and act as an index of the focus of covert spatial attention based on subtle changes in their directionality with attention condition 44 . Given these findings, we considered two means by which microsaccades could account for our results. First, micro-saccade direction may vary as a function of attention condition, differently modulating neuronal firing activity across conditions and potentially generating the pattern of correlated variability we report. However, the direction of micro-saccades did not vary across attention conditions in our task ( Fig. 7a; F Changes in task difficulty cannot account for our results. A further potential confounding variable is task difficulty. Recent work has shown that increasing task difficulty is associated with lower spike count correlations, presumably by modulating the overall level of arousal of the subject 45 . If behavioral conditions in which two stimuli must be monitored for a possible change are more difficult than conditions in which only one stimulus needs monitoring, then correlations should be lowest in the AB condition of our task. In fact, we found correlations to be highest in the AB condition (Fig. 5a), suggesting that increased task difficulty does not account for our results in the AB condition. As noted previously, however, to attempt to balance task difficulty across conditions, we increased coherences by one step in the AB condition. One could argue that this change in coherence may have over-corrected for task difficulty and made the AB condition easier, leading to higher correlations in the AB condition by the converse of the above argument. Several observations argue against this possibility. If the AB condition were easier than the other conditions, we would expect the percentage of changes detected to be higher in the AB condition, which was not the case (Fig. 3e). Additionally, decreased task difficulty in the AB condition cannot account for the positive correlation between stimulus eccentricity and the degree to which correlations are elevated in the AB condition (Fig. 5d), because task difficulty is likely to increase, rather than decrease, with eccentricity.
Finally, exploiting the relationship between task difficulty and arousal level 45 and using pupil size as a measure of the overall arousal level of a subject 46,47 , we assessed whether changes in arousal level across task conditions could account for our results. Because we had not recorded pupil size for the sessions reported above, we collected a new set of behavioral sessions in which we recorded pupil size and for which stimulus parameters were matched to those used in our original dataset. We found no significant difference of pupil sizes between the attention conditions in this new dataset, suggesting that our results cannot be explained by changes in the level of arousal either ( Fig. 7c; F(2, 7) = 2.7, p = 0.11, rmANOVA).
Other potential confounds. Further, our results are not trivially explained by changes in firing rates across conditions, as firing rates in the AI condition were elevated compared to the AO condition (Fig. 4b), but correlation magnitudes were not significantly different in these conditions (Fig. 5a, b). In fact, this dissociation between attentional modulation of firing rates and of spike count correlations is consistent with the predictions of our previously published model of attention 8,48 . Finally, changes in stimulus coherence cannot function as an explanation for elevated correlations in the AB condition, as spike counts were analyzed during the ZCP before any changes in the stimulus coherence occurred.

Discussion
We developed a task to dissociate changes in the strength of attentional modulation from changes in variability in the attentional state by varying the behavioral relevance of two simultaneously presented stimuli and encouraging the use of different attentional allocation strategies across task conditions. We found the effects of attention on correlated variability to differ depending on the timescale analyzed. At a timescale of 1000 ms, levels of shared variability were highest in the condition in which both stimuli were behaviorally relevant, supporting the idea that this condition introduced competition for attentional resources, which increased attentional state variability. In contrast, shared variability was lowest in the conditions in which attention could be focused on only one stimulus, and there was no difference in correlations in the AI and AO conditions at this timescale. These results are consistent with the scenario presented in Fig. 2b, in line with our previous predictions 8 , and support the hypothesis that fluctuations in the state of attention can be a prominent source of shared neuronal response variability. More generally, these results suggest that a significant fraction of shared variability in neuronal populations can be attributed to fluctuations in behaviorally-relevant, internally generated signals, rather than shared sensory noise 8,16,29,[48][49][50][51][52] . Further, at a timescale of 200 ms, we found correlations between neurons in the supragranular cortical layers were lower in the AI relative to the AO condition, consistent with earlier work that considered faster timescales, both in V4 and in V1 11,12,17,53 , and with the scenario depicted in Fig. 2a. Verhoef and Maunsell 54 recently demonstrated how the reduction of correlations under attention could be due to a suppression of (variable) normalizing inputs from the unattended surround 54 , largely consistent with previously hypothesized explanations 11,12 . Taken together, these results suggest that both mechanismssuppression of common noise and attentional fluctuationsimpact levels of correlated variability, but they operate at different timescales.
The importance of timescale could explain why a recent study that employed an attention task with conditions similar to ours, including a neutrally-cued condition akin to our AB condition, found correlations to be intermediate between the attend-in and attend-out conditions at a timescale of 200 ms 55 . Further, both Mayo and Maunsell 55 and Cohen and Maunsell 14 collected data simultaneously from both hemispheres but reported no significant correlation, or anti-correlation as one would expect with a shifting spotlight-like attentional allocation strategy, amongst neurons in opposite hemispheres. Perhaps such a correlation does exist at timescales longer than was analyzed in those studies. Unfortunately, our data cannot resolve this question, as we recorded from only one hemisphere at a time.
Because the impact of variability in the attentional state on correlations manifested on a timescale of individual trials in our task, should we therefore expect that fluctuations in internal signals, in general, only induce correlations on long timescales? Ultimately, this timescale is likely to depend on the mechanism by which such signals impact neuronal populations. Work on orienting of attention and attentional dwell time suggests that voluntarily shifting attention between different stimuli takes on the order of several hundred milliseconds 27,28 . In an experimental paradigm similar to our AB condition, attention was found to alternate between two stimulus locations roughly every 250 ms (4 Hz) 18,19 . This shifting of attention between stimulus locations is the strategy we were hoping to induce in our paradigm and appears to be the likeliest explanation for how attention is allocated across trials in our AB condition, given our behavioral and neurophysiological results. We would, thus, expect that AB correlations should be elevated on a timescale of at least several hundred milliseconds, which is what we found (Fig. 5c).
Note that this line of reasoning stands regardless of whether the shift in attention that occurs involves a narrowly-focused attention field encompassing only one stimulus at a timeresembling the spotlight or narrowly-focused Zoom Lens models 24,56 -or whether some degree of attention is allocated to both stimuli simultaneously, but with one stimulus receiving a greater degree of attention than the other on a given trialresembling the variable precision model of resource allocation 57 . In this latter case, the shift of attention corresponds to alternations in which stimulus receives the greater strength of attentional focus on a given trial. The key, however, is that some change in attentional resources allocated to the receptive field stimulus occurs across trials. Therefore, our results are not consistent with models of attention that suggest that both stimuli are processed simultaneously and that a consistent or uniform degree of attentional processing is distributed across the full field of attention.
Interestingly, we also found a correlation between the horizontal eccentricity of the stimuli and the degree to which correlations in the AB condition were elevated compared to the AO and AI conditions (Fig. 5d). We interpret this finding to suggest that when stimuli are closer to each other, it is easier to attend both simultaneously, resulting in a lower degree of attentional fluctuation in the AB condition. As the stimuli are placed farther apart, attending to both simultaneously becomes increasingly difficult, and subjects are more likely to deploy a switching allocation strategy, leading to more pronounced attentional fluctuations and, thus, higher correlations in the AB condition.
While alternating which stimulus receives the greater strength of attentional processing on a given trial is one means by which attentional state variability increases (across trials), there may be other sources of variability in the attentional state as well. For example, a number of studies have shown that improvements in behavior due to attention, rather than being continuous across time within a trial, appear to exhibit a theta-frequency periodicity, which is related to theta-band cortical oscillations and can occur even with attention focused on only one stimulus [58][59][60] . If attention operates in a periodic manner, as these studies suggest, such oscillations could represent an additional source of variability in the attentional state beyond that induced by alternating attention between stimulus locations. Further studies have suggested that shifts in attention between stimulus locations are also linked to theta-band oscillatory activity 19,59,61 , raising a number of interesting questions. Does attention itself truly operate periodically, or do ongoing cortical oscillations mediate the effects of an otherwise more continuous attention signal, giving the appearance of periodicity? Are shifts in attention only possible at certain phases of these ongoing cortical rhythms? Ultimately, these are important empirical questions that future research should address. To do so will require a combination of behavioral paradigms that allow attention-related performance to be tracked more explicitly over time 18 and multi-electrode array recordings with single-unit-resolution population analyses such as those undertaken in the present study.
Another interesting question is how correlations in an attention task impact behavioral performance. Quantifying precisely how correlations affect the information encoding capacity of a neuronal population in an experimental setting is a challenge because one would have to decode from a large population of simultaneously recorded neurons 9 . Because we do not have such a sufficiently large dataset, we cannot draw any conclusions regarding the impact of correlations on performance. Nonetheless, this is a critical topic for future work to address.
Recent studies have examined the laminar profile of attentional modulation of firing rates 31 or of spike count correlations during passive fixation 35,36 . Only one study has examined the laminar relationship between attentional modulation and shared variability 62 , and ours is the first to do so in V1. Nandy et al. 62 found significant attentional modulation of firing rates in all layers, with the strongest effects in the granular layer. In contrast, van Kerkoerle et al. 31 found the weakest attentional modulation of firing rates in the granular layer of V1. Similar to Nandy et al. 62 , we found significant modulation of firing rates by attention in all layers in the AI condition. However, considering both the AB and AI conditions, our results are in better agreement with those of van Kerkoerle et al. 31 , as we found the strongest attentional modulation of firing rates in the supragranular, followed by the infragranular layers, as expected given the anatomical distribution of feedback cortical connections [20][21][22][23] .
Regarding correlation magnitude across layers, we observed different patterns of results at the two main timescales we analyzed, 200 and 1000 ms. At the 1000 ms interval there was no significant effect of layer on correlation magnitude, whereas at the 200 ms interval, correlations were lowest in the granular layer, consistent with previous laminar studies in V1 35,36 . This 200 ms interval is similar to the window size used in Hansen et al. 35 While Smith et al. 36 found a similar pattern over a 1280 ms interval, they recorded from anesthetized animals where the mechanisms driving correlated fluctuations are likely to be very different from those during wakefulness 29 .
At both timescales, attentional modulation of correlations was confined primarily to the supragranular layers and was not present in the infragranular layers, despite attentional modulation of rates in the AI condition. One reason may be a lack of sufficient statistical power. Most of our isolated single units were from the supragranular layers (just over eight units per session on average), with about half that number isolated in the infragranular layers, and fewer still from the granular layer. The difference could also be attributable to the anatomical and computational characteristics of each layer, which by no means are completely understood 34,63,64 . The infragranular layers additionally receive feedback from and send projections to subcortical regions 65 and such signals may modulate shared variability differently. Ultimately, the finding that attention predominantly modulates correlations in the supragranular layers matches the location where we found the most pronounced attentional modulation of firing rates and accords well with the known anatomy of corticocortical interactions, particularly for feedback signals.
Nandy et al. 62 also found attentional modulation of correlations to be strongest in the same layer in which they found attentional modulation of firing rates to be strongest. Interestingly, this layer was not the supragranular layer but rather the granular layer. As suggested by Nandy et al. 62 , it is possible that the input layer in V4 inherits the correlation pattern from the output (supragranular) layers of V1. Our results at the 200 ms interval in the supragranular layers are consistent with this possibility and match the findings reported by Nandy et al. 62 It is also possible that attention operates somewhat differently in V4 than in V1, with attentional modulation of firing rates typically being stronger overall and occurring earlier in the response period in V4 25,37 .
Overall, correlations in the present study were a bit higher than in our earlier studies with awake fixating animals 49 . The primary difference between these studies is that subjects in the present study perform a demanding task engaging feedback processes such as attention, and our main results demonstrate the effect that fluctuations in such signals have on levels of correlated variability. Although attentional fluctuations are reduced in the focused attention conditions, they are unlikely to be entirely absent, so some elevation in correlation magnitude above zero in these conditions is to be expected. Additionally, correlations are also likely to be somewhat higher given that the highly dynamic stimulus in the current study drives the neurons much more strongly than static or drifting gratings.
Finally, there has been an increasing interest in recent years in leveraging population recording and latent-variable modeling techniques to infer the state of internally-generated, cognitive signals, such as attention, on more behaviorally-relevant timescales, to better understand the nature of these signals and their impact on decision-making and behavior 16,[66][67][68] . To make such inferences, these methods make use of the patterns of covariance in population activity and rely on the assumption that this variability occurs in a low-dimensional space (e.g., the "attention axis" 14 ). A further, but critical, assumption of these techniques is that much of this shared variability is not noise but is attributable to the action of behaviorally-relevant, internally generated signals. However, a clearer demonstration that changes in internal signals indeed contribute significantly to shared neuronal variability was lacking. We presented a paradigm designed specifically to test for such a contribution, and our results provide support for this critical assumption. Additionally, our results demonstrate the subtlety of the effects that internal signals such as attention have on correlated variability, exemplified by the two timescales over which attention modulated correlations.

Methods
Experimental model and subject details. All behavioral and electrophysiological data were obtained from two healthy, male rhesus macaque (Macaca mulatta) monkeys (B and D) aged 12 and 13 years and weighing 11 and 10 kg, respectively, during the time of study. All experimental procedures complied with guidelines of the NIH and were approved by the Baylor College of Medicine Institutional Animal Care and Use Committee (permit number: AN-4367). Animals were housed individually in a large room located adjacent to the training facility, along with around ten other monkeys permitting rich visual, olfactory and auditory interactions, on a 12 h light/dark cycle. Regular veterinary care and monitoring, balanced nutrition and environmental enrichment were provided by the Center for Comparative Medicine of Baylor College of Medicine. Surgical procedures on monkeys were conducted under general anesthesia following standard aseptic techniques. To ameliorate pain after surgery, analgesics were given for 7 days. Animals were not sacrificed after the experiments.
Visual stimuli and behavioral paradigm. Visual stimuli were two Gabor patches (size: diameter of 2-3°depending on eccentricity; spatial frequency: 3-3.5 cycles per degree; contrast: 100% Michelson; eccentricity: 3.7-8.9°) presented on CRT monitors (at a distance of 100 cm; resolution: 1600 × 1200 pixels; refresh rate: 100 Hz) using Psychophysics Toolbox 69 . The monitors were gamma corrected to have a linear luminance response profile. Video cameras (DALSA genie HM640; frame rate 200 Hz) with custom video eye tracking software developed in LabView were used to monitor eye movements.
Monkeys performed a noisy, orientation-change detection task. Trials were initiated by a sound and the appearance of a colored fixation target (~0.15°). Monkeys were required to fixate within a radius of 0.5°-1°, but typically fixated much more accurately, as revealed by offline analysis. After fixating for 300 ms, two Gabor patches were presented symmetrically in the lower left and right visual fields. During what we labeled the Zero-Coherence Period (ZCP), these stimuli changed their orientation pseudo-randomly every 10 ms (uniform distribution over 36 orientations spaced by 5°between 0 and 175°) for a random period of time drawn from an exponential distribution with a minimum of 10 ms, mean of 2170 ms, and maximum of 5000 ms.
After this time one of the two stimuli entered the coherent period (CP), where one particular orientation, called the "signal" orientation, was shown with a higher frequency than the other orientations. The CP lasted 300 ms (30 frames), and from trial to trial the number of frames in the CP showing the signal orientation was selected from a set of five unique "coherences" chosen for that session, which allowed us to vary the difficulty of the trials within a session and compute psychometric functions. After this period, the stimulus returned to the ZCP for a further 200 ms to allow sufficient time for subjects to report whether or not they noticed the presence of the signal orientation by making a saccade to the stimulus showing the change. Subjects were prevented from responding within the first 100 ms of the CP to minimize guessing. Successful identification of the signal orientation was rewarded with a small drop of juice. On 10% of trials in each attention condition no change occurred, and subjects were rewarded for maintaining fixation. Orthogonal signal orientations were used in the left (135°) and right (45°) stimuli.
Note, occurrences of the signal orientation during the CP were not constrained to occur in successive frames. Also note that the left and right stimuli displayed different orientation sequences, so that subjects could not identify a change simply by noticing when the two orientation sequences diverged. Orientation sequences were described as pseudo-random for the following reason. For each trial a random number generator seed was chosen from a set of five such seeds selected for a given recording session. Doing so meant there were five unique stimuli that could be repeated across attention conditions for the purposes of calculating spike count correlations and Fano factors over identical stimuli. Sequences were constrained to show each orientation once before any repetitions were allowed so that the maximum number of signal orientations that could occur by chance in a period of time equal to the CP (300 ms) was two.
Attention was cued in blocks of trials by the color of the fixation spot (Fig. 3b). In the attend out (AO) condition, 100% of the changes occurred in the nonreceptive field stimulus. In the attend in (AI) condition, 100% of changes occurred in the receptive field stimulus. In the attend both (AB) condition, the change was equally likely to occur in either stimulus (50% chance that the change was in the receptive field stimulus). Block transitions occurred after a total of 60 hit and miss trials was achieved (i.e., false alarms did not count). Blocks were randomized in sets of three so that each attention condition was seen before one was allowed to repeat. Coherences were increased by one frame in the AB condition to keep task difficulty approximately constant across conditions. Surgical methods. Our surgical procedures followed a previously established approach 70 . A cranial headpost was first implanted under general anesthesia using aseptic conditions in a dedicated operating room. After premedication with atropine (0.05 mg/kg prior to sedation), animals were sedated with a mixture of ketamine (10 mg/kg) and dexdormitor (0.015 mg/kg). During the surgery anesthesia was maintained using isoflurane (0.5-2%).
After subjects were trained to perform the above described task, they were implanted with a form-fitted titanium recording chamber, designed based on preoperatively obtained anatomical MRI scans, placed at a location over the operculum in V1 determined by stereotactic coordinates 70 . This surgery was performed under identical conditions as described for headpost implantation. The chamber was attached to the skull using orthopedic screws only. We used a small amount of dental cement to seal any openings between the bone and the lower surface of the recording chamber. A custom-made chamber cap was then placed to seal the chamber and prevent infection. A minimum of three weeks was provided for the implant to heal. After healing, small 2-3 mm trephinations could be performed, in aseptic conditions under ketamine (10 mg/kg) sedation with ketoprophen (2 mg/kg) for analgesia and meloxicam (0.2 mg/kg for two days), to enable access for subsequent daily electrophysiological recordings.
Electrophysiology in awake, behaving monkeys. We performed daily electrophysiological recordings beginning 48 h after a craniotomy was performed. Custom-designed 32 channel, linear silicon probes (NeuroNexus V1 x 32-Edge-10mm-60-177) with inter-channel spacing of 60 μm, contact site dimensions of roughly 12 × 15 μm, contact site area of 177 μm 2 and typical impedances around 1 mega-Ohm were mounted in a Narishige microdrive (MO-97) with a nested, stainless steel guide tube composed of one extra-thin walled 23-gauge piece, spanning most of the length of the probe shaft, and a smaller 27-gauge piece (roughly 6 mm long) nested inside such that 4 mm of the smaller tubing protruded beyond the large piece. This design enabled a tight fit around the probe to support it during dural penetrations. We took care during the insertion procedure to ensure that the dura was penetrated only by the probe itself, rather than the guide tube, to minimize damage to the superficial layers of cortex. We alternated lowering the guide tube in steps of 250 μm and extending the probe up to~500 μm beyond the guide tube, retracting and repeating as necessary, until either characteristic changes in the LFP or multi-unit activity, or both, were observed, indicating successful penetration of cortex.
The probe was then lowered in~250 μm steps at <10 μm per second, pausing for several minutes after each step, until activity was seen on all channels. As a result of this procedure there would be variable degrees of tissue compression. Some of this compression was relieved early in the positioning of the probe by retracting the guide tube by~500 μm after the probe was several hundred microns inside the cortex. If compression remained after completely lowering the probe, we could successfully relieve it by slowly retracting the guide tube further. The single most reliable indicator of the position of our probe in cortex before receptive field mapping was a band of high spontaneous activity corresponding to layer 4 C 32 , which could be clearly seen to span roughly 6-7 channels. In general, we found the basic laminar properties described by Snodderly and Gur 32 to be very reliable guidelines. After final positioning of the probe, we allowed between 30 and 60 min for tissue settling and recording stability to become established. The entire insertion procedure typically took around 3-4 h, from penetrating the dura to the start of recording. Receptive field mapping experiments were performed (see Data analysis below for details) to determine where to place one of the two stimuli such that it covered the recorded neurons' receptive fields for that session.
Data acquisition and spike sorting. The methods described below for spike detection and spike sorting were adapted for use with multi-channel silicon probes from our previous methods used for tetrode recordings 29 . Neural signals were digitized at 24 bits using analog acquisition cards with 30 dB of onboard gain (PXI-4498, National Instruments, Austin, TX) and recorded continuously at 32 KHz as broad-band signal (0.5 Hz to 16 kHz). Eye movement traces were sampled at 2 kHz.
Spikes were detected offline when the signal on a given channel crossed a threshold of five times the standard deviation of the corresponding channel. To avoid artificial inflation of the threshold in the presence of a large number of high amplitude spikes, we used a robust estimator of the standard deviation, given by σ = median(|x|)/0.6745 71 . Spikes were aligned to the center of mass of the continuous waveform segment above half the peak amplitude. Code for spike detection is available online at [https://github.com/atlab/spikedetection].
Virtual electrodes consisting of six channels were constructed in a sliding window (stride 2) spanning the length of the probe to aid in the spike sorting process by enabling some degree of triangulation, as with tetrodes. Given a channel spacing of 60 μm, in many cases the waveforms of a single neuron could be detected by several channels. To extract features for spike sorting, we performed principal component analysis on the extracted waveform segments (individually for each channel). This step reduced the data to three dimensions per channel, resulting in an 18-dimensional feature vector. We fit a mixture of t distributions with a Kalman filter on the cluster means to track waveform drift 72 .
The number of clusters was determined based on a penalized average likelihood, where the penalty term was a constant cost per additional cluster. Code for spike sorting is available online at [https://github.com/aecker/moksm]. Following this automatic step, results of the model were examined manually for each virtual electrode and single units were flagged at this time according to degree of cluster isolation, uniqueness of waveforms and size of refractory period. To avoid duplicate single units due to overlapping channel groups used for spike sorting, we included only those single units that had their largest waveform amplitude on one of the two central channels of the group (this was not an issue for the first and last two channels on the probe). . We included recording sessions with at least 10 single units that were visually responsive and significantly orientation tuned in each attention condition. To ensure reliable estimates of neuronal (co-)variability, sessions were also excluded if there were fewer than three (of five possible) valid seed conditions. A seed condition was considered invalid if in any of the three attention conditions there were fewer than three correct trials generated using that seed that had sufficient ZCP length available for spike count analysis. On average for the 1-s analysis window, included sessions had~10 correct trials per seed per attention condition.
After having collected a complete dataset of 13 sessions from Subject B and a dataset of 29 sessions from Subject D, we found that sessions with recording locations close to the vertical meridian did not exhibit our predicted main effect. We reasoned that this lack of effect was likely because the two stimuli were too close to each other, allowing the monkey to attend to both simultaneously. To verify that this result was not a false positive due to post-hoc analysis, we collected an independent 10-session dataset at high eccentricities from Subject D (the termination condition of 10 sessions was set before starting to collect additional data), which confirmed the effect at high eccentricity. The results reported in this paper, except in Fig. 5d, include all sessions with x-axis receptive field eccentricities of at least 3°(representing the median such eccentricities for Subject B), including the separate validation dataset from Subject D.
Data analysis. Data were analyzed in Matlab, using custom Matlab software and the DataJoint processing pipeline 73 .
Trial results were classified as "hits", "misses", "correct rejections" (for successful completion of trials with no change) and "false alarms" (for saccades made to a stimulus before any change occurred). For each session, behavior was analyzed by calculating the fraction of changes detected (hits / [hits + misses]), both conditioned on and marginalized over coherence in each attention condition. Psychometric functions were plotted as the fraction of changes detected versus coherence in each attention condition. Using the psignifit toolbox 74,75 in MATLAB, logistic functions were fit to the attention condition specific curves using the method of maximum likelihood, and 50% performance thresholds were extracted. Reaction times could be calculated using only hit trials and reaction time distributions for each session were quantified by calculating the median deviation for each condition in each session. False alarm rates were calculated using all valid trials ("hits", "misses", "correct rejections", "false alarms").
Prior to starting the main task, we quantitatively mapped receptive fields based on unsorted multi-unit responses using a white noise random dot stimulus. A single square dot of size 0.29°of visual angle was presented on a uniform gray background, changing location and color (black or white) randomly every three frames, or 30 ms, for 1 s. Receptive field profiles were obtained by spike-triggered averaging. Average diameter of multi-unit receptive fields across sessions was 1.14 ± 0.05°.
Our task allowed us to compute orientation tuning curves for each neuron. We binned the spike counts in bins of 10 ms and used linear regression based on a onehot encoding of the 15 stimuli directly preceding the response (i.e., the stimulus is a 36 × 15-dimensional vector, because there were 36 possible stimulus orientations). We defined the optimal latency of each neuron as the time delay that produced the strongest response modulation across orientations (determined by taking the variance of the regression weights across orientations). The optimal latency of most neurons was 50 ms. We then re-estimated the regression using only that single time lag to obtain a tuning curve. Significance of tuning was then tested by projecting the weight vector onto a complex exponential with one cycle, the norm of which was compared to its null distribution calculated by randomly shuffling orientation labels. A p-value was obtained by performing 1000 iterations of the shuffling procedure and using the fraction of runs in which the norm of the shuffled projection was greater than that observed in the real data. Signal correlations were computed for pairs of neurons by calculating the correlation coefficient between the two cells' tuning curves.
For each unit, a von Mises distribution function, parameterized as was fit to the tuning curve obtained across all trials via the method described above. From this fit, the shape and preferred orientation parameters, w 3 and w 4 , were obtained. These parameters were assumed not to change across attention conditions, leaving only the offset, w 1 , and gain, exp(w 2 ), terms to vary across conditions. New von Mises functions were then fit for each attention condition using a linear regression model with a binary indicator variable for attention condition and an interaction term. To illustrate, we write the response y to orientation i as Þ and was obtained from the overall tuning curve as described. Our linear regression model comparing fits in the AO and AI condition, for example, then became: where X i1 ¼ θ i and X i2 2 f0; 1g, with 0 coding the AO condition and 1 coding the AI condition. In this manner we enabled different gain and offset terms to be fit to different attention conditions. We then assessed whether significant attentional modulation was present by performing an F-test comparing the full model above to the reduced model containing only the β 0 and β 1 terms, and when significant, we tested whether the offset and gain parameters differed between conditions with ttests.
Visual responsiveness of neurons was determined by comparing firing rates in the 300 ms fixation interval before stimulus onset to those in the 300 ms immediately following stimulus onset. A t-test was performed to test for a significant change in rate following stimulus onset. Spike density functions (SDFs) were calculated first for a given neuron, across all hit trials grouped by attention condition and stimulus seed, by counting spikes in 50 ms bins relative to stimulus onset and averaging across trials. Averages were then taken across seeds and smoothed with a Gaussian window. To calculate SDFs for a given session, individual neuron SDFs were normalized by the average response in the AO condition, starting from 100 ms after stimulus onset, before averaging across neurons. Fractional firing rate increases were also calculated first at the individual neuronal level, by averaging all available bins from the first second following stimulus onset conditioned on the stimulus seed for each attention condition, and then averaging across seeds. The rates were again normalized by the AO condition rate before averaging across neurons to get a session-level rate modulation for each attention condition. Finally, responses in the AI and AB conditions were converted to fractional changes relative to the AO responses.
Fano factors and spike count correlations were computed on the first 1000 ms of the response. Fano factors were computed as the variance of the spike count divided by its mean. Spike count correlations were computed as the covariance of the two neurons' z-scored responses to identical repetitions of the same stimulus condition (seed). Z-scoring and Fano factor calculations were performed in a block-wise fashion to control for slow fluctuations in firing rate across a recording session. For the analysis of correlation timescale we used the relationship between spike count correlations and cross-correlation functions first described in Bair et al. 3 to compute a cumulative correlation coefficient, r CCG . We compute a spike train cross-correlation function for a pair of neurons j and k, as well as a shift-predictor, which is the cross-correlation function of the spike density functions of neurons j and k. The shift-predictor is subtracted from the cross-correlation function to control for stimulus-induced correlation. This shift-corrected cross-correlation is denoted C jk (τ). The cumulative crosscorrelation is given by Following Ecker et al. 29 , the cumulative correlation coefficient is where T is the last time point in the counting window, in our case 1000 ms. The CSD profile at each time point was calculated as the second spatial derivative of the task-stimulus evoked LFPs across channels, smoothed with a Gaussian kernel to aid visualization 30 . The granular layer was identified according to several criteria used in conjunction. The earliest current sink to source transition (identified by an arrow in Fig. 6a) is one indicator, immediately below which is a complementary source to sink transition in L5. We used additional criteria, described by Snodderly and Gur 32 , to verify this positioning, because there was a prominent current sink to source transition in L6 as well. These criteria included higher spontaneous activity and more poorly defined orientation tuning curves characteristic of the granular layer 32 . Additional reports have described the granular layer to contain smaller receptive fields 76,77 , which we also saw (Fig. 6a). In general across sessions, all of these granular layer features were quite consistent, allowing for confident determination of the L4-5 boundary. The first L5 channel was labeled as the zero-point for depth. Negative depths are more superficial to this point. The granular layer was defined as a roughly 400 μm band just superficial to the zero-point [33][34][35][36] . The supragranular group (L1-3) was defined as everything superficial to the top of the granular layer, and the infragranular group (L5-6) was defined as everything deeper than and including the zero-point.
We identified micro-saccades our subjects made during the ZCP of our task (when spike counts were analyzed) to determine whether our correlation results could be accounted for by an increase in micro-saccade frequency in our AB condition, relative to the AI and AO conditions. Periods of stable gaze were taken to be those intervals during which eye position remained within a 0.1°window, and deviations greater than 0.1°in 10 ms (10°/s velocity) were taken to be microsaccades 78 . The number of micro-saccades during analysis periods was counted for each attention condition in each session and a repeated-measures ANOVA was performed to determine whether micro-saccades differed across conditions. Microsaccades were also grouped according to the direction in which the saccade was made (unit circle divided into 8 equal direction bins) and a two-factor, repeatedmeasures ANOVA was used to assess for effects of direction and condition (the two factors). Pupil size was measured for a set of N = 8 sessions recorded from Subject B using the same camera and software used for eye-tracking described above. Stimulus parameters were matched with those used for the original dataset. Pupil size was determined based on the number of pixels above a threshold brightness value and an effect of attention condition on pupil size was determined using a repeated-measures ANOVA.
Quantification and statistical analysis. Although customary in the field, we did not consider units or pairs as independent samples. Treating units as independent samples ignores the session-to-session variability and leads to underestimated confidence intervals and, consequently, inflated false positive rates. Instead, we first averaged our measurements across observations within a session and then performed all statistical tests across sessions, treating the session averages as independent samples. While this approach sacrifices some statistical power, it leads to conservative estimates of p values.
For statistical analyses involving our attention conditions, repeated-measures ANOVAs were used, with session as the random factor and attention condition as the fixed factor. F-statistic values are reported as F(x,y), where x represents the number of degrees of freedom for the fixed factor of attention condition, and y is the equivalent for the random factor of session. The Tukey-Kramer method was primarily used for post-hoc analyses. To test for significantly elevated AB condition correlations, we performed a one-tailed t-test on a contrast between the AB condition and the average of the AO and AI condition results. This choice is justified by our previously published model 8 , which predicts this effect and its direction and was hypothesized and specified before data collection. Statistics for the t-test are reported as t(x), where x represents the degrees of freedom. Note, in the section discussing laminar results, any reductions in the number of degrees of freedom are due to instances in which insufficient single units were isolated in a particular layer for that session to be included in that particular analysis.
A two-factor, repeated-measures ANOVA was used to test changes in microsaccade direction with attention condition. In this case the F-statistic is reported as F(x,y,z), where x represents the number of degrees of freedom for the factor of attention condition, y represents that for the factor of direction, and z represents that for the random factor of session. For assessments of visual responsiveness and significant increases in fractional firing rates, two-tailed t-tests were used, which, for rate increases, were Bonferroni-corrected for multiple comparisons. Orientation tuning significance was assessed according to the permutation test described above. Statistical comparisons were considered significant at p < 0.05 (p < 0.0167 for Bonferroni-corrected tests for firing rates in association with Fig. 4c, as there were three comparisons; p < 0.025 for those associated with Fig. 6b, given two comparisons). All error bars show the standard error of the mean (SEM; either directly calculated or estimated via ANOVA), except in the Fig. 3c inset, which shows 95% confidence intervals. No blinding was used in the analysis.
Code availability. The code used to process and analyze the data for the current study are available from the corresponding author on reasonable request. Links to some of this code have been provided in the Methods section "Data acquisition and spike sorting." Data availability. The datasets generated during and analyzed during the current study are available from the corresponding author on reasonable request.