Catecholamines alter the intrinsic variability of cortical population activity and perception

The ascending modulatory systems of the brain stem are powerful regulators of global brain state. Disturbances of these systems are implicated in several major neuropsychiatric disorders. Yet, how these systems interact with specific neural computations in the cerebral cortex to shape perception, cognition, and behavior remains poorly understood. Here, we probed into the effect of two such systems, the catecholaminergic (dopaminergic and noradrenergic) and cholinergic systems, on an important aspect of cortical computation: its intrinsic variability. To this end, we combined placebo-controlled pharmacological intervention in humans, recordings of cortical population activity using magnetoencephalography (MEG), and psychophysical measurements of the perception of ambiguous visual input. A low-dose catecholaminergic, but not cholinergic, manipulation altered the rate of spontaneous perceptual fluctuations as well as the temporal structure of “scale-free” population activity of large swaths of the visual and parietal cortices. Computational analyses indicate that both effects were consistent with an increase in excitatory relative to inhibitory activity in the cortical areas underlying visual perceptual inference. We propose that catecholamines regulate the variability of perception and cognition through dynamically changing the cortical excitation–inhibition ratio. The combined readout of fluctuations in perception and cortical activity we established here may prove useful as an efficient and easily accessible marker of altered cortical computation in neuropsychiatric disorders.


Introduction
The modulatory systems of the brain stem send widespread, ascending projections to the specialized circuits of the cerebral cortex that mediate perception, cognition, and goal-directed behavior. These systems regulate ongoing changes in brain state, even during periods of wakefulness [1][2][3][4]. They are recruited during, and in turn, shape cognitive processes such as perceptual inference, learning, and decision-making [5][6][7][8]. Because these systems are implicated in most neuropsychiatric disorders, they are also major targets of the pharmacological therapy of brain disorders [5,9,10]. Taken together, neuromodulatory systems have remarkably specific effects on cognition, despite the widespread nature of their projections to the cortex. An important challenge for neuroscience is to uncover the mechanistic principles by which neuromodulatory systems interact with the cortical computations underlying cognition.
One key parameter of cortical computation that might be under neuromodulatory control is the intrinsic variability-i.e., fluctuations that occur during constant (or absent) sensory input [11,12]. Specifically, it has been proposed that the catecholaminergic neuromodulators noradrenaline and dopamine may shift the cortical computations underlying decision-making from a stable ("exploitative") to a variable ("exploratory") mode [5,13]. A context-dependent adjustment of the variability of cortical computations may also be adaptive for perceptual inference in the face of ambiguous sensory input [14].
Animal work has shown that catecholamines and acetylcholine, another important neuromodulator, alter the intrinsic variability of neural activity [2,[15][16][17][18] through highly selective interactions with specific elements (pyramidal cells and/or inhibitory interneurons) of cortical microcircuits [19,20]. However, it is unknown how these changes at the level of cortical microcircuits relate to the intrinsic variability of perception and cognition.
At the larger scale of cortical mass action that is assessable with noninvasive recordings in humans, activity also fluctuates intrinsically in a spatially and temporally structured manner [21,22]. The temporal structure of these fluctuations is characteristic of so-called "scale-free" behavior: power spectra that scale as a function of frequency according to a power law, P(f) / f β [23,24], indicating long-range temporal autocorrelations [25][26][27][28]. Some studies have linked the spatiotemporal structure of the fluctuations in cortical population activity to specific perceptual and cognitive processes [27,[29][30][31]. However, it is unknown if and how these fluctuations in cortical population activity are dynamically regulated by neuromodulatory systems.
We aimed to close these gaps by systematically quantifying the effects of catecholaminergic and cholinergic neuromodulation on the intrinsic variability in perception and large-scale cortical activity in the healthy human brain. To this end, we combined placebo-controlled, selective pharmacological interventions, psychophysical measurements of fluctuations in perception in the face of a continuously presented and ambiguous visual stimulus, and recordings of fluctuations in cortical population activity using magnetoencephalography (MEG).
Catecholamines, but not acetylcholine, increased both the variability of perception as well as long-range temporal correlations of intrinsic cortical activity in the visual and parietal cortices. Based on previous theoretical and experimental work [32][33][34][35], we interpreted the increase in perceptual variability in terms of an increase in the net ratio between cortical excitation and inhibition in those cortical regions. Simulating a recurrent neural network under synaptic gain modulation enabled us to show that an analogous mechanism may account for the increase of long-range temporal correlations of cortical activity under catecholamines.

Results
We tested for changes in intrinsic fluctuations of perception and cortical population activity under placebo-controlled, within-subjects pharmacological manipulations of catecholamine (using the noradrenaline reuptake inhibitor atomoxetine) and acetylcholine (using the cholinesterase inhibitor donepezil) levels ( Fig 1A, see Materials and methods for details on the pharmacological interventions).
Fluctuations in cortical activity were measured during two steady-state conditions, both of which excluded transients in sensory input or motor output ( Fig 1B): (i) fixation of an otherwise gray screen (a condition termed "Fixation"), as in most studies of human "resting-state" activity [21,22], and (ii) silent counting of the spontaneous alternations in the perceptual interpretation of a continuously presented, ambiguous visual stimulus (dubbed "Task-counting"). In a third condition that was only used for the analysis of perceptual fluctuations, subjects immediately reported the perceptual alternations by button press ("Task-pressing," i.e., associated with movement-related transients in cortical activity). This design capitalized on recent insights into the circuit mechanisms underlying intrinsic perceptual dynamics [32,33,35], which helped constrain the mechanistic interpretation of the results reported below.
The Results section is organized as follows. We first present the effects of the "Atomoxetine" and "Donepezil" conditions (each compared against the "Placebo" condition) on the rate of perceptual fluctuations. These effects were in line with a boost in the relative strength of the excitatory drive of the visual cortex under atomoxetine. We then show how (i) constant sensory and task drive (i.e., Task-counting versus Fixation) and (ii) pharmacological manipulations affect the intrinsic fluctuations in cortical activity. We focus on the (long-range) temporal autocorrelation structure of intrinsic fluctuations in the amplitude of band-limited cortical population activity (see Materials and methods and Fig 2). Control analyses showing In log-log coordinates, <F N > increases approximately linearly as a function of N, with a slope that is the scaling exponent α. (D) Illustration of power spectrum analysis of amplitude envelope. In log-log coordinates, the power spectrum can be approximated by a straight line, with a slope β (power-law exponent) and an AUC (gray) that quantifies the overall variance of the signal. The data can be found at https://figshare.com/articles/Example_Figure_Fig_2_/5787321. AUC, area under the curve; F N , RMSE fluctuation averaged across segments of length N; MEG, magnetoencephalography; RMSE, Root-mean-square error; T, Tesla. the drug effects on other measures of cortical population activity and peripheral physiological signals support the validity and specificity of our conclusions.
We close with simulations of a patch of recurrently connected excitatory and inhibitory integrate-and-fire neurons. The simulations show that the changes in temporal correlations observed in the MEG data can be explained by a modulation of synaptic gain that alters the net ratio between excitatory and inhibitory activity.

Atomoxetine increases the rate of bistable perceptual fluctuations
The ambiguous visual stimulus that was continuously presented during both Task-counting and Task-pressing induced ongoing fluctuations in perception, i.e., spontaneous alternations between two apparent rotation directions of 3D motion (Fig 1B; see S1 Movie), a phenomenon that is referred to as multistable perception. The rate of the perceptual alternations reported by the participants provided a readout of visual cortex circuit state. Current models explain bistable perceptual fluctuations in terms of the interplay between the feed-forward, excitatory drive of stimulus-selective neural populations in the visual cortex; mutual inhibition between them; stimulus-selective adaptation; and neural "noise" [32,33]. Increases in the ratio between feed-forward, excitatory input to and mutual inhibition within the cortical circuit give rise to faster perceptual alternations. This idea is supported by convergent evidence from functional magnetic resonance imaging, magnetic resonance spectroscopy, and pharmacological manipulation of GABAergic transmission [30,35]. We reasoned that neuromodulators such as noradrenaline might dynamically change these parameters [36,37] and thereby alter the rate of perceptual fluctuations.
Atomoxetine increased the rate of perceptual fluctuations compared to both Placebo and Donepezil conditions (Fig 3A; atomoxetine versus placebo: p = 0.007, t = 2.913; atomoxetine versus donepezil: p = 0.001, t = 3.632; donepezil versus placebo: p = 0.966, t = −0.043; all paired t tests, pooled across Task-counting and Task-pressing). The perceptual alternation rates were highly consistent across Task-counting and Task-pressing (S1C Fig), supporting the validity of the counting condition as behavioral readout of bistable perceptual fluctuations. Likewise, the atomoxetine effect on the perceptual fluctuation rate was evident for Task These changes in perceptual fluctuations were not explained by an increase in the rates of eye blinks or fixational eye movements. First, there was no significant increase during atomoxetine compared to placebo in any of five different eye movement parameters measured here (S2 Fig). Second, none of these parameters correlated with the perceptual alternation rate (S2 Fig). Third and most importantly, the effect of atomoxetine on the perceptual dynamics was also significant after removing (via linear regression) the individual eye movement parameters ( Fig 3B).
In sum, atomoxetine had an effect on bistable perceptual fluctuations that was both robust and specific and evident when compared with either placebo or donepezil. This effect was in line with an increase in the strength of the excitatory feed-forward drive of the visual cortex, relative to the strength of mutual inhibition between the neural subpopulations encoding the competing perceptual interpretations of the ambiguous stimulus. Such an effect should have occurred in the motion-sensitive visual cortical areas, in which the visual competition induced by the ambiguous structure-from-motion stimulus is implemented [38,39].

Atomoxetine increases the scaling exponent of fluctuations in cortical population activity
We estimated long-range temporal correlations of band-limited amplitude fluctuations (indicated by the scaling exponent α; see Materials and methods for details) to quantify intrinsic fluctuations in cortical population activity. Our analyses focused on amplitude envelope fluctuations in the 8-12-Hz frequency range ("alpha band") for two reasons. First, as expected from previous work [40], the cortical power spectra exhibited a clearly discernible peak in this frequency range, which robustly modulated with sensory or task drive (suppressed under Task-counting, Fig 1C). Second, previous studies reported robust long-range temporal correlations with peaks in the same frequency range [28].
We first replicated two previously established observations pertaining to the scaling exponent α. First, the average across cortical patches and participants was α = 0.67 (SD = 0.09) during Fixation (Placebo condition only) and α = 0.64 (SD = 0.07) during Task-counting (Placebo condition only), indicative of long-range temporal correlations similar to the ones found in previous work [25,28,41]. Second, the sensory and task drive during Task-counting reliably reduced α compared to Fixation, again as shown in previous work [26,42]. Across all voxels, α was significantly larger during Fixation than during Task-counting (p = 0.0062, t = 2.97, paired t test, Placebo condition only). This difference was significant across pharmacological conditions in large parts of the cortex, including the occipital and parietal regions, that were driven by the motion stimulus (p < 0.05, cluster-based permutation test; Fig 4D).
Having verified the validity of our measurements of α, we then tested for changes in α under the pharmacological conditions (Figs 4B-4E and 5). There was a highly significant increase in α for atomoxetine compared to placebo when collapsing across voxels as well as across Fixation and Task-counting (p = 0.0068, t = 2.93; paired t test, Fig 4B and 4C). This effect was widespread but not homogenous across the cortex, comprising occipital and posterior parietal as well as a number of midline regions, including the thalamus (Fig 4D, p = 0.0022; cluster-based permutation test). Because it is unclear to what extent intrinsic activations from deep sources can be recovered using MEG, we focus our description and conclusions on the effects in cortical regions. Importantly, the atomoxetine effect on α was also present at the level of MEG sensors (S4 Fig) and hence did not depend on the source reconstruction method applied here (see Materials and methods). The effect of atomoxetine on α was subtle, likely due to the low dosage. However, importantly, the effect was highly reproducible across repeated measurements. We assessed reproducibility with two complementary approaches. The first was a region-of-interest (ROI) analysis. We defined an ROI in terms of a significant cluster for Atomoxetine > Placebo (onesided paired t test, p < 0.05, uncorrected) during the first run collected in each session (collapsed across Fixation and Task-counting) and extracted this ROI's mean α from the second run. We then reversed the procedure and so extracted a second, independent ROI-based α and averaged the α-estimates. This approach revealed a strong increase under atomoxetine (p = 0.0023, t = 3.365). The second approach assessed the reproducibility of the spatial pattern of effects across both runs. To this end, we correlated the (nonthresholded) individual maps for the atomoxetine versus placebo difference computed from the first and second run in each session (again pooling across Task-counting and Fixation) and tested the resulting correlation coefficients across participants. The average correlation was significantly different from 0 (mean r = 0.29, p < 0.0001; permutation test against 0).
The atomoxetine-related increases in scaling exponent α were evident during both Fixation and Task (Fig 5C). There was no interaction between the effects of atomoxetine and Task-counting anywhere in the cortex: a direct comparison of the two atomoxetine versus placebo difference maps, from Fixation and from Task-counting, yielded no significant clusters (p > 0.081 for all clusters; cluster-based permutation test). The same cortical regions in which α increased during atomoxetine exhibited decreases during Task-counting: when testing for the task-dependent change in α (Fig 4A), specifically in the regions comprising the conjunction cluster of the atomoxetine effect (Fig 5C), the reduction during Task-counting was also highly significant ( Fig 5D) in all pharmacological conditions.
In contrast to the robust effect of atomoxetine on α, there was no evidence for an effect of donepezil at the dosage used here. The difference between donepezil and placebo (collapsed Taken together, the rich experimental design gave rise to a highly specific and consistent pattern of changes in α under the different experimental conditions, which helped constrain the mechanistic interpretation of the results. The atomoxetine effects were specific and not just due to the application of any drug targeting neurotransmitter systems. It is possible that the absence of detectable donepezil effects on α was due to the low dosage or short administration period used here. However, the control analyses presented in the next section revealed clear effects of donepezil on both cortical activity as well as markers of peripheral nervous system activity.

Control analyses for the drug effects on other features of cortical dynamics or peripheral physiological signals
During Fixation, atomoxetine and donepezil both reduced posterior cortical alpha-band power relative to placebo in both the 8-12-Hz ( Fig 6A; p < 0.05 for all clusters; two-sided cluster-based permutation test) as well as the 2-8-Hz frequency ranges (S7A Fig). This suppression in low-frequency power under cholinergic boost is consistent with previous work in rodents [16,17] and humans [43].
The atomoxetine-induced changes on 8-12-Hz power exhibited a different spatial pattern from the one of corresponding change in the scaling exponent α: within the cluster of the significant main effect of atomoxetine on α (Fig 4D), power did not correlate with the changes in α (group average spatial correlation between pooled difference maps within the cluster; r = 0.073, p = 0.129, BF = 1.065). During Task-counting, neither drug altered MEG power in the low frequencies (8-12 Hz: Fig 6B,  We also controlled for changes in peripheral physiological signals under the drugs as potential confounds of the effect on cortical scaling behavior (Fig 7). As expected, atomoxetine increased average heart rate ( Fig 7A and 7B). Donepezil had no detectable effect on average heart rate during either Fixation (p = 0.8676, t = 0.16; paired t test; BF = 0.8676; Fig 7A) or Task-counting (p = 0.3274, t = 1.0; paired t test; BF = 0.3139; Fig 7B). Both drugs altered heart rate variability, increasing α computed on the time series of inter-heartbeat intervals (see Materials and methods) in both behavioral contexts relative to placebo (Fixation: p = 0.0012, t = 3.62; Task In sum, drug-induced changes in peripheral physiological signals under the drugs, if present, did not account for the atomoxetine-induced changes in the scaling behavior of the fluctuations in cortical activity (Figs 4 and 5). These controls support our interpretation in terms of a specific effect of atomoxetine on cortical variability rather than nonspecific secondary effects due to the systemic drug effects or changes in retinal input due to blinks.

Change in scaling exponent under atomoxetine is consistent with increase in net excitation-inhibition ratio in cortical circuits
Atomoxetine had an effect on perceptual fluctuations that was in line with a relative increase in excitation in the cortical circuits of the occipital and posterior parietal cortices that processed the ambiguous visual motion stimulus. We reasoned that this change in circuit state might have also produced the observed change in the scaling behavior of intrinsic cortical activity fluctuations under atomoxetine. In order to solidify this intuition, we simulated the activity of a neural network model made up of recurrently connected excitatory and inhibitory integrate-and-fire units (Fig 8). In what follows, we use the term "excitation-inhibition [E/I] ratio" to refer to the ratio of excitatory and inhibitory activity across the circuit [44] and "E/I balance" to refer to a specific regime of E/I ratios, in which excitation and inhibition change in a proportional way [45][46][47][48].
We started from a network ( Fig 8A) that was similar to the one developed and analyzed in a previous study [49]. The basic features of the model were as follows. The model was built to generate oscillations of neural mass activity (summed across all units) in the alpha band (8-12 Hz; Fig 8B). The amplitude envelope of these oscillations fluctuated over time, with scale-free long-range temporal correlations. Those scale-free intrinsic fluctuations in cortical activity were sensitive to variations in the percentage of excitatory and inhibitory connections in the circuit (i.e., microstructure). Our previous work [49], which we reproduced here (Fig 8D-8F), showed that such a model accounts for the joint emergence of two scale-free phenomena at different spatial scales (single unit activity versus mass activity) and temporal scales (tens of milliseconds versus hundreds of seconds): (i) neuronal avalanches with an event size distribution following a power law and (ii) long-range temporal correlations of the amplitude envelope fluctuations of the circuits' mass activity. Both phenomena have been established in empirical measurements of cortical population activity [25,50]. Neuronal avalanches are activity deflections (i.e., exceeding a certain threshold) that propagate through the cortical network [50], with an "event size" corresponding to the number of activated units. In line with [51], we quantified the power-law scaling of the size distributions of avalanches in the model with the kappa index (κ): the similarity between the actual event size distribution and a power-law distribution with an exponent of −1.5; a κ of 1 indicates perfect match between the two. We extended this model by means of a multiplicative modulation of synaptic gain [36,52] ( Fig 8B). This allowed us to explore how catecholaminergic effects on neural circuits might change the two phenomena of scale-free neural population activity described above. We first determined the structural connectivity (small squares in Fig 8D-8F) and the timescale The code and the data underlying these plots can be found at https://figshare.com/s/374ab0f973f026535549. EE, excitatory-to-excitatory (recurrent excitation); EI, inhibitory-to-excitatory; IE, excitatory-to-inhibitory; II, inhibitory-to-inhibitory (recurrent inhibition); κ, kappa index; w, synaptic weight. parameters of the model such that the network generated intrinsic alpha-band oscillations ( Fig  8C) with amplitude fluctuations that exhibited neuronal avalanches with scale-free event size distributions ( Fig 8D) as well as long-range temporal correlations (with α~0.85). We then independently modulated specific excitatory or inhibitory connections through the multiplicative scaling of the corresponding synaptic weights in two ways. In the version shown in Fig 8, we modulated only excitatory synapses, but independently on excitatory as well as inhibitory neurons (EE and IE), thus producing asymmetries in the circuits' net E/I ratio, similar to recent modeling work on a cortical circuit for perceptual decision-making [44]. In the second version (S8A Fig), we co-modulated EE and IE and independently modulated inhibitory synapses on excitatory neurons (EI). This was intended to specifically simulate glutamate receptors (AMPA or NMDA) in the former two cases (mediating the effects of excitatory neurons) as opposed to modulations of GABA receptors in the latter case (mediating the effects of inhibitory neurons on others). N EE and N IE were co-modulated by the same factor for simplicity, because we did not assume that excitatory (glutamatergic) synapses would be differentially modulated depending on whether they were situated on excitatory or inhibitory target neurons.
Both types of changes in net E/I ratio robustly altered κ ( Fig 8G and S8B Fig), α ( Fig 8H  and S8C Fig), and the mean firing rate (Fig 8I). The effect of changes in E/I ratio on the scaling exponent α were nonmonotonic and dependent on the starting point: increases in excitation led to increases in α when starting from an inhibition-dominant point but to decreases in α when starting from an excitation-dominant point (Fig 8H, white line). The effects of excitatory and inhibitory gain modulation on the temporal correlation structure of the simulated population activity were qualitatively similar to the effects of changes in the fraction of excitatory and inhibitory synapses simulated (as shown in Fig 8D-8F). The latter simulated individual differences in cortical anatomical microstructure, and the former simulated state-dependent changes in cortical circuit interaction, which occur within an individual brain.
In the model, the scaling exponent α exhibited a nonmonotonic dependence on E/I ratio (see the white diagonal line in Fig 8G-8I and schematic depiction in Fig 9). Consequently, without knowing the baseline state, any change in α was ambiguous with respect to the direction of the change in E/I ratio (i.e., towards excitation or inhibition dominance). Thus, the observed increase in α under atomoxetine during Fixation could have been due to either an increase or a decrease in E/I ratio.
Importantly, insights from animal physiology helped constrain the baseline state during Task-counting: in the awake state, visual drive decreases E/I ratio in primary visual cortex (V1), due to the recruitment of inhibitory mechanisms that outweigh the excitatory sensory drive [53,54]. We assumed that the same held for the Task-counting condition (constant visual stimulation) of our study.
This condition enabled us to infer the change in net E/I ratio under atomoxetine. The rationale is illustrated in Fig 9. The animal physiology results referred to above indicate that the observed decrease in α during Task-counting was due to a shift towards inhibition dominance (yellow point in Fig 9A). Under this assumption, the atomoxetine-induced increase in α was due to an increase in net E/I ratio (Fig 9B). Because the effects of atomoxetine on α were the same during Task-counting and Fixation, it is likely that the same mechanism was at play during Fixation.
In sum, under certain conditions, the simulations provided a mechanistic explanation for the observed MEG effects: effective changes in the cortical E/I ratio due to multiplicative changes of synaptic gain [36] or other mechanisms [19,37]-the same conclusion inferred from the increase in the rate of perceptual alternations above.

Discussion
Neuromodulators regulate ongoing changes in the operating mode of cognitive processes [1,5,6,10,55] as well as of cortical microcircuits [15][16][17]19,20,37]. Here, we unraveled the effect of two major classes of neuromodulators, catecholamines and acetylcholine, on the intrinsic variability of cortical computation, an important parameter shaping the operating mode. We used two separate readouts of this parameter: (i) the rate of fluctuations in perception induced by ambiguous visual input and (ii) the temporal autocorrelation structure of fluctuations of the amplitude envelope of band-limited cortical mass activity. Catecholamines, but not acetylcholine, altered both readouts. Simulations of a recurrent neural network revealed that under wellsupported physiological assumptions, the observed changes in the temporal structure of fluctuations in cortical activity are indicative of an increase in the cortical E/I ratio. Earlier modeling and empirical data [35] show that such an increase in net E/I ratio in the visual cortex is also consistent with the increased rate of perceptual fluctuations under the catecholaminergic boost.

Cortical distribution of atomoxetine effects on cortical activity fluctuations
The atomoxetine effects on the scaling exponent were widespread across the cortex, but not entirely homogenous. They were pronounced across the occipital and parietal cortex but not robust in the frontal cortex (see Fig 5B). This distribution might point to a noradrenergic rather than dopaminergic origin. Atomoxetine increases the levels of both catecholamines, noradrenaline and dopamine, but the cortical projection zones differ substantially between both systems: dopaminergic projections mainly target the prefrontal cortex [56] and only sparsely the occipital cortex [57,58], whereas the noradrenergic projections are more widespread and strong to the occipital and parietal cortices [59]. Alternatively, this distribution may reflect the different receptor composition across cortical regions [59,60]. The relative frequency of the different noradrenaline receptors differs between prefrontal and posterior cortices [59], which might translate a homogenous noradrenaline release into a heterogeneous effect on the activity in these different cortical regions. An important next step will be to investigate the differential role of different noradrenaline receptors and different regional receptor profiles in shaping the cortex-wide effect of noradrenaline on long-range temporal correlations.

Opposite effects of external drive and catecholamines on long-range temporal correlations
Consistent with our current results, previous studies also found a decrease in temporal autocorrelations of cortical activity due to external drive [26,42]. The observation is consistent with the insight from intracellular recordings of cortical neurons in animals that cortical responses to sensory stimulation in the awake state are dominated by inhibition [53,54,61]. One candidate source of this sensory-driven state change is thalamocortical inhibition [62], but intracortical feedback inhibition might also contribute [63]. Modeling work shows that the driven state is associated with shortened temporal autocorrelations as well as a decrease in the entropy of activity states in large-scale cortical networks [64]. Correspondingly, the increase in longrange temporal autocorrelations under catecholaminergic modulation may be associated with an increase in entropy-that is, a tendency to explore a larger set of cortical activity states. It is tempting to link this to the prominent idea that high sustained noradrenaline levels promote an exploratory mode of cortical computation and behavior [5].

Convergent evidence for catecholaminergic increase in cortical E/I ratio
Cortical circuits maintain a tight balance between excitation and inhibition, which is largely preserved across contexts and levels of the cortical hierarchy [46,48]. However, even in the absence of changes in sensory input, neuromodulators such as noradrenaline and acetylcholine can change the cortical E/I ratio [65,66]. The E/I ratio, in turn, shapes the computational properties of cortical circuits [67,68] and thereby the behavior of the organism [36,44,69]. Substantial evidence already points to significant changes in E/I ratio in schizophrenia and autism [70][71][72]. Similar changes might be at play in other brain disorders as well [73].
Our simulations indicated that the temporal autocorrelation structure of neural population activity, as measured with the scaling exponent α, is sensitive to changes in E/I ratio produced through synaptic gain modulation (see the white line in Fig 8H). In both versions of our model, the neuromodulatory effects were not perfectly symmetric (see the deviations of peak scaling exponents from the main diagonal in Fig 8H). While the latter effect was small and may be specific to the details of the model, it remains possible that the subtle changes in scaling exponents we observed were produced through symmetric gain modulations that maintained the net E/I balance (i.e., along the main diagonal). However, two additional lines of evidence converge on our conclusion that catecholamines (in particular, noradrenaline) boosted E/I ratio. First, in the same participants, the catecholaminergic manipulation had a reliable effect on the perceptual switch rate, which is also indicative of cortical E/I ratio [32,33,35]. Second, results from invasive rodent work also point to an increase in cortical E/I ratio under noradrenaline: Noradrenaline was found to decrease spontaneous inhibition in the auditory cortex [66] and mediate a tonic depolarization of visual cortical neurons during locomotion [19].

No evidence for donepezil effects on cortical or perceptual fluctuations
The absence of an effect of donepezil on either perceptual fluctuations or long-range temporal correlations of cortical activity may be due to the small dosage or the single administration of the drug in our study. Even so, our donepezil manipulation was sufficient to robustly change heart rate variability and, more importantly, low-frequency power of cortical activity, an established marker of cholinergic action in the cortex [16,17,43,74]. The lack of the effects of donepezil on perceptual fluctuations and cortical scaling behavior might also be due to the specific properties of cholinergic action on the cortical net E/I ratio. Invasive evidence indicates that acetylcholine can rapidly disinhibit pyramidal cells by activating a chain of two inhibitory interneurons [20], a mechanism that may alter E/I ratio mainly during stimulus-evoked responses [65]. By contrast, noradrenaline also alters the levels of the tonic inhibition of pyramidal cells occurring spontaneously [66]. This might explain the dissociation between the effects of atomoxetine and donepezil under the current steady-state conditions, which excluded (or minimized) stimulus-evoked transients.

Functional consequences of changes in net cortical E/I ratio
We observed a selective increase in the rate of spontaneous perceptual alternations under increased catecholamine levels, adding to evidence that these dynamics are under neuromodulatory control [75]. Such a change could be due to an increase in the intrinsic variability of cortical activity [32]. Future invasive studies should relate catecholaminergic changes in the variability of the spiking activity [76] of neurons contributing directly to the contents of multistable perception.
We suspect that an increase in cortical E/I ratio will have particularly strong effects on behavior when affecting parietal and prefrontal cortical circuits characterized by slow intrinsic timescales [29,31,77] and involved in persistent activity during working memory and the slow accumulation of information over time [69]. It is possible that the catecholaminergic effects on the parietal cortex we observed here reflect an increase in recurrent excitation, which is essential for sustained processes, such as working memory [78], as well as information integration during decision-making [31,79]. Future work should assess this through the use of tasks probing into network reverberation and information accumulation in the association cortex.

A control parameter for critical network dynamics
In our model, long-range temporal correlations in the fluctuations of neural mass activity (i.e., activity summed across the entire local network) [25] and avalanches within the neuronal network [50] jointly emerge at the same E/I ratio. Both phenomena are commonly interpreted as hallmarks of "criticality" [25,28,50,80]-the state of a complex dynamical system poised between order and chaos [81][82][83]. It has been proposed that the cortex operates in a narrow regime around criticality [83,84], potentially optimizing its computational capacities [51,80,[85][86][87]. A number of reports showed that cortical dynamics may continuously vary around the critical state [88][89][90][91], but the source of these fluctuations has so far remained unknown. Here, we have identified catecholaminergic neuromodulation as an endogenous factor controlling these spontaneous variations in critical dynamics.
In complex systems, critical dynamics can emerge in a self-organized fashion [81] or through an external control parameter that fine-tunes the system. The tuning of temperature in the Ising model of spin magnetization [83] is a common example of the latter case. It is tempting to speculate that catecholaminergic tone serves as such a control parameter in the cerebral cortex.

Link between catecholaminergic effects on fluctuations in perception and cortical mass activity
We here used two readouts of catecholaminergic effects, constituting two distinct expressions of the resulting changes in cortical circuit state. The envelope of cortical alpha-band oscillations collapsed across large chunks of cortex is unlikely to encode the contents of perception in the phenomenon studied here. The perceived direction of 3D motion, which fluctuates spontaneously, is encoded in fine-grained patterns of neural population activity within motion-sensitive visual cortical areas [39,92]. The power of alpha-band oscillations is a more global feature of cortical population activity, which is likely insensitive to the fine-grained, within-area patterns of neural population activity. The widespread release of neuromodulators also changes the cortical circuit state, specifically E/I ratio, in a widespread manner. Such changes, in turn, alter the highly specific (fine-grained) interactions between percept-selective populations of visual cortical neurons that give rise to the perceptual dynamics [32,33,35]. Thus, although both readouts likely tap into similar changes in global cortical circuit state, there is no one-to-one mapping between them.

Limitations of the current approach
While our model simulations provided important mechanistic insights, the model has limitations that should be addressed in future work. First, different from the MEG data, the power of alpha-band oscillations behaves similarly to the scaling exponents in the model (S8E Fig). This is because in the model, oscillations emerge from the same recurrent neuronal interactions that also shape the long-range temporal correlations in the amplitude envelopes of these oscillations. By contrast, in the brain, alpha-band power of local cortical mass signals is likely affected by a variety of sources other than local circuits, for instance, alpha-frequency-modulated input from the thalamus [93]. This might lead to dissociations between changes in MEG power and long-range temporal correlations of the power fluctuations, which the model does not capture in its present form. Second, the model lacks long-range excitatory connections, which are prominent in the real cortex and the effects of which on the correlation structure of cortical fluctuations are largely unknown.
Another limitation of our study concerns our mechanistic interpretation in terms of E/I ratio. This interpretation rests on the assumption that sensory drive has the same effect on E/I ratio in humans as previously shown in rodent visual cortex [53,54]. Our study and the previous ones differed in terms of behavioral task, stimulus, and species, all of which might have led to violations of the assumption. However, a number of observations from human electrophysiology point to the notion, consistent with our assumption and the above rodent work, that the awake human cortex operates in an inhibition-dominant regime [94], which is further amplified during active stimulus processing [89,90].

Conclusion
The combined measurement of fluctuations in bistable perception as well as in cortical mass activity under steady-state conditions provides an easily assessable, multilevel readout of pharmacological effects on cortical computation. In our study, this readout supported the idea that catecholamines boost the intrinsic variability of perception and behavior, an effect that might be mediated by an increase in the net E/I ratio in the visual cortical system. This readout may be useful for inferring changes in the cortical E/I ratio in neuropsychiatric disorders or in their pharmacological treatment in future work.

Ethics statement
All participants gave written informed consent before the start of experiment, in accordance with the Declaration of Helsinki. The study was approved by the ethical committee responsible for the University Medical Center Hamburg-Eppendorf (Ethik-Komission der Ä rztekammer Hamburg, ID PV4648).

Pharmacological MEG experiment
Participants. Thirty healthy human participants (16 females, age range 20-36, mean 26.7) participated in the study. Two participants were excluded from analyses, one due to excessive MEG artifacts and the other due to not completing all three recording sessions. Thus, we report results from N = 28 participants (15 females). In one of those participants, one Taskcounting run (during the Atomoxetine condition) was not recorded due to a software problem with the data acquisition computer. In two participants, button presses were not correctly recorded (during the Atomoxetine condition). The behavioral data were dismissed for those participants.
General design. We pharmacologically manipulated the levels of catecholamines (noradrenaline and dopamine) and acetylcholine in a double-blind, randomized, placebocontrolled, crossover experimental design (Fig 1A and 1B). Each participant completed three experimental sessions, consisting of drug (or placebo) intake at two time points, a waiting period of 3 h, and an MEG recording. During each MEG session, participants were seated on a chair inside a magnetically shielded MEG chamber. Each session consisted of 6 runs of different tasks, each of which was 10 min long and followed by breaks of variable duration.
Pharmacological intervention. We used the selective noradrenaline reuptake inhibitor atomoxetine (dose: 40 mg) to boost the levels of catecholamines, specifically noradrenaline and (in the prefrontal cortex) dopamine [10]. We used the cholinesterase inhibitor donepezil (dose: 5 mg) to boost acetylcholine levels. Atomoxetine is a relatively selective inhibitor of the noradrenaline transporter, which is responsible for the natural reuptake of noradrenaline that has been released into the extracellular space. Consequently, atomoxetine acts to increase the extracellular levels of noradrenaline, an effect that has been confirmed experimentally in rat prefrontal cortex [95]. The same study showed that atomoxetine also increases the prefrontal levels of dopamine, which has a molecular structure very similar to the one of noradrenaline and is, in fact, a direct precursor of noradrenaline. Atomoxetine has smaller affinity to the serotonin transporter, and there are discrepant reports about the quantitative relevance of these effects: while one study found no increases in serotonin levels under atomoxetine [95], a recent study reports a significant atomoxetine-related occupancy of the serotonin transporter in nonhuman primates [96] at dosages that would correspond to human dosages of 1.0-1.8 mg/kg. Note that these dosages are substantially higher than the administered dosage in this study (40 mg, independent of body weight). It is therefore unclear to what extent our Atomoxetine condition affected cortical serotonin levels.
Donepezil is a selective inhibitor of the enzyme acetylcholinesterase, which breaks up all the extracellular acetylcholine to terminate its synaptic action. Consequently, donepezil acts to increase the extracellular levels of acetylcholine. Donepezil is also an agonist of the endoplasmatic sigma 1 -receptor, which modulates intracellular calcium signaling.
A mannitol-aerosil mixture was administered as placebo. All substances were encapsulated identically in order to render them visually indistinguishable. Peak plasma concentrations are reached about 3-4 h after administration for donepezil [97] and 1-2 h after administration for atomoxetine [98]. We adopted the following procedure to account for these different pharmacokinetics ( Fig 1A): participants received two pills in each session, one 3 h and another 1.5 h before the start of MEG recording. In the Atomoxetine condition, they first received a placebo pill (t = −3 h) followed by the atomoxetine pill (t = −1.5 h). In the Donepezil condition, they first received the donepezil pill (t = −3 h) followed by placebo (t = −1.5 h). In the Placebo condition, they received a placebo at both time points. The half-life is about 5 h for atomoxetine [98] and about 82 h for donepezil [97]. In order to allow plasma concentration levels to return to baseline, the three recording sessions were scheduled at least 2 weeks apart. This design ensured maximum efficacy of both pharmacological manipulations, while effectively blinding participants as well as experimenters.
Stimuli and behavioral tasks. In each session, participants alternated between three different task conditions (2 runs at 10 min per condition) referred to as Fixation, Task-counting, and Task-pressing in the following (Fig 1B). All conditions entailed overall constant sensory input. Fixation and Task-counting also entailed no overt motor responses and are therefore referred to as "steady-state" conditions in the following. We used these steady-state conditions to quantify intrinsic fluctuations in cortical activity. Task-pressing entailed motor responses and was used for reliable quantification of perceptual dynamics. All instructions and stimuli were projected onto a screen (distance: 60 cm) inside the MEG chamber. The individual conditions are described as follows.
Fixation: Participants were asked to keep their eyes open and fixate a green fixation dot (radius = 0.45˚visual angle) presented in the center of an otherwise gray screen. This is analogous to eyes-open measurements of "resting-state" activity widely used in the literature on intrinsic cortical activity fluctuations.
Task-counting: Participants viewed a seemingly rotating sphere, giving rise to the kinetic depth effect [99,100] (S1 Movie): spontaneous changes in the perceived rotation direction (Fig 1B). The stimulus subtended 21˚of visual angle. It consisted of 1,000 dots (500 black and 500 white dots, radius: 0.18˚of visual angle) arranged on a circular aperture presented on a mean-luminance gray background, with the green fixation dot in the center. In order to minimize tracking eye movements, the sphere rotation was along the horizontal axis, either "forward" (towards the observer) or "backward" (away from the observer), and the dot density decreased along the horizontal axis towards the center of the stimulus. Participants were instructed to count the number of perceived changes in rotation direction and report the total number of perceived transitions at the end of the run. Just like during Fixation, Taskcounting minimized any external (sensory or motor) transients. Subjects silently counted the alternations in perceived rotation direction and verbally reported the total count after the end of the 10-min run.
Task-pressing: This condition was identical to Task-counting, except that participants were instructed to press and hold one of two buttons with their index finger to indicate the perceived rotation direction of the sphere. Thus, each perceptual alternation was accompanied by a motor response leading to a change in the button state. This allowed for a more reliable quantification of participants' perceptual dynamics. In two sessions (Atomoxetine condition), button presses were not registered. Hence, the corresponding analyses were performed on 26 participants.

Data acquisition
MEG was recorded using a whole-head CTF 275 MEG system (CTF Systems, Inc., Canada) at a sampling rate of 1,200 Hz. In addition, eye movements and pupil diameter were recorded with an MEG-compatible EyeLink 1000 Long Range Mount system (SR Research, Osgoode, ON, Canada) at a sampling rate of 1,000 Hz. In addition, electrocardiogram (ECG) as well as vertical, horizontal, and radial electrooculogram (EOG) were acquired using Ag/AgCl electrodes (sampling rate 1,200 Hz).

Data analysis
Eye data. Eye blinks were detected using the manufacturer's standard algorithm with default settings. Saccades and microsaccades were detected using the saccade detection algorithm described in [101], with a minimum saccade duration of 4 samples (equal to 4 ms) and a threshold velocity of 6. For 18 out of 28 participants, only horizontal eye movements were recorded. EOG data. EOG events (blinks and saccades) were extracted using semiautomatic artifact procedures, as implemented in FieldTrip [102]. In short, EOG traces were band-pass filtered using a third-order Butterworth filter (1)(2)(3)(4)(5), and the resulting signal was z-scored. All time points at which the resulting signal exceeded a z-score of 4 were marked as EOG events.
MEG data. Preprocessing: First, all data were cleaned of strong transient muscle artifacts and squid jumps through visual inspection and manual as well as semiautomatic artifact rejection procedures, as implemented in the FieldTrip toolbox for MATLAB [102]. To this end, data segments contaminated by such artifacts (±500 ms) were discarded from the data (across all channels). Subsequently, data were downsampled to 400 Hz, split into low-(2-40 Hz) and high-frequency (>40 Hz) components, using a fourth-order (low-or high-pass) Butterworth filter. Both signal components were separately submitted to independent component analysis [103] using the FastICA algorithm [104]. Artifactual components (eye blinks/movements, muscle artifacts, heartbeat, and other extracranial artifacts) were identified based on three established criteria [105]: power spectrum, fluctuation in signal variance over time (in bins of 1-s length), and topography. Artifact components were reconstructed and subtracted from the raw signal and low and high frequencies were combined into a single dataset. On average, 20 (±14) artifact components were identified for the low frequencies and 13 (±7) artifactual components were identified for the high frequencies.
Spectral analysis: Sensor-level spectral estimates (power spectra and cross-spectral density matrices) were computed by means of the multi-taper method using a sequence of discrete prolate Slepian tapers [106]. For the power spectrum shown in Fig 1C, power spectra were computed using a window length of 5 s and a frequency smoothing of 2 Hz, yielding 19 orthogonal tapers. The focus of this paper was on the fluctuations of the amplitude envelopes rather than on the (oscillatory) fluctuations of the carrier signals per se. The temporal correlation structure of the amplitude envelope fluctuations of cortical activity seems similar across different carrier frequency bands [28]. We focused on amplitude envelope fluctuations in the alpha band because (i) the cortical power spectra exhibited a clearly discernible alpha peak, which robustly modulated with the task, as expected from previous work [40] (Fig 1C), and (ii) the computational model used to study the effect of synaptic gain modulation on cortical activity fluctuations was tuned to produce alpha-band oscillations (see above and [49]).
Source reconstruction: General approach: The cleaned sensor level signals (N sensors) were projected onto a grid consisting of M = 3,000 voxels covering gray matter of the entire brain (mean distance: 6.3 mm) using the exact low-resolution brain electromagnetic tomography (eLORETA; [107]) method. The grid was constructed from the ICBM152 template [108], covering gray matter across the brain. The magnetic leadfield was computed separately for each subject and session using a single shell head model constructed from the individual structural MRI scans and the head position relative to the MEG sensors at the beginning of the run [109]. In case no MRI was available (4 subjects), the leadfield was computed from a standard MNI template brain transformed to an estimate of the individual volume conductor using the measured fiducials (located at the nasion, the left ear, and the right ear).
In order to depict the source-level results, we interpolated the voxel-level results onto the surface of the brain. Activations from structures distant to the surface are not shown and were exponentially attenuated.
Source level estimates of amplitude envelopes and power: For comparing amplitude envelope and power estimates between experimental conditions in source space, we aimed to select a single direction of the spatial filter for each voxel across pharmacological conditions (i.e., MEG sessions), but separately for Fixation and Task-Counting conditions. The rationale was to avoid filter-induced biases in the comparisons between the pharmacological conditions, while allowing that external task drive might systematically change the dipole orientations.
To this end, we first computed the mean source-level cross-spectral density matrix C(r, f) for each frequency band, f, averaged across the three MEG sessions, as follows: where i indicated the MEG session, C i (f) was the (sensor-level) session-and frequency-specific cross-spectral density matrix, and A i was the spatial filter for session i. We then extracted the first eigenvector u 1 (r, f) of the session-average matrix C(r, f) by means of singular value decomposition and computed the unbiased filter selective for the dominant dipole orientation, B i (r, f), as: This procedure ensures that, for each voxel, dipole orientation was chosen such that power was maximized. Please note that this filter was now frequency specific, whereas the previous filters, A i (r), were not. To obtain instantaneous estimates of source-level amplitudes, the sensor-level signal for session i, X i (t), was band-pass filtered (using a finite impulse response filter) and Hilbert transformed, yielding a complex-valued signal H i (f, t) for each frequency band. This signal was projected into source space through multiplication with the unbiased spatial filter, B i (r, f), and the absolute value was taken: where Env i (r, f, t) was the estimated amplitude envelope time course of source location r and frequency f. Next, for each session, unbiased source-level cross-spectral density estimates were obtained from the sensor-level cross-spectral density matrix C i (f) and the frequency-specific, unbiased spatial filter B i (f). The main diagonal of the resulting matrix contains source-level power estimates for all source locations: These computations were repeated separately for the Task-counting and Fixation conditions, session by session. The differences in amplitude envelope fluctuations and power estimates between pharmacological and task conditions reported in this paper were robust, with respect to the specifics of the analysis approach. In particular, we obtained qualitatively similar pharmacological effects in sensor space, as reported in an earlier conference abstract [110].
Detrended fluctuation analysis: The source-level amplitude envelopes Env i (r, f, t) were submitted to detrended fluctuation analysis [111,112] in order to quantify long-range temporal correlations. Detrended fluctuation analysis quantifies the power law scaling of the fluctuation (root-mean-square) of a locally detrended, cumulative signal with time-window length. Different from the analysis of the more widely known autocorrelation function [29,77], the detrended fluctuation analysis provides robust estimates of the autocorrelation structure for stationary and nonstationary time series. The procedure of the detrended fluctuation analysis is illustrated in Fig 2. For simplicity, in the following, we rewrite the amplitude envelope Env i (r, f, t) as x of length T. First, we computed the cumulative sum of the demeaned x ( Fig 2B): where t 0 and t denote single time points up to length T. The cumulative signal X was then cut Fig  2B, top). Within each segment Y i of equal length N, the linear trend Y i_trend (least-squares fit) was subtracted from Y i (Fig 2B, bottom, blue versus red lines), and the root-mean-square fluctuation for a given segment was computed as: where n indicates the individual time points. The fluctuation was computed for all k segments of equal length N and the average fluctuation was obtained through: The procedure was repeated for 15 different logarithmically spaced window lengths N ranging from 3-50 s, which yields a fluctuation function ( Fig 2C). As expected for scale-free time series (103), this fluctuation function follows a power law of the form: The "scaling exponent" α was computed through a linear regression fit in log-log coordinates ( Fig 2C). The longest and shortest window lengths were chosen according to guidelines provided in [112].
A scaling exponent of α~= 0.5 indicates a temporally uncorrelated ("white noise") process. Scaling exponents between 0.5 < α < 1 are indicative of scale-free behavior and long-range temporal correlations [112], whereas exponents of α < 0.5 indicate long-range anticorrelations ("switching behavior") and α > 1 are indicative of an unbounded process [112]. The scaling exponents for alpha-band MEG amplitude envelopes estimated in this study ranged (across experimental conditions, MEG sensors, and participants) between 0.40 and 1.04, with 99.4% of all estimates in the range 0.5-1. This is indicative of scale-free behavior and consistent with previous human MEG work [25][26][27][28]42,113].
Relationship between measures of cortical variability: Scale-free behavior of neural time series has also been quantified via analysis of the power spectrum [23,24]. There is a straightforward relationship between both approaches, which we explain below to help appreciate our results in the context of these previous studies. The power spectrum of the amplitude envelope of cortical activity is typically well approximated by the power law p(f) / f −β , where β is referred to as the power-law exponent (Fig 2D). For power-law decaying autocorrelations, the relationship between the power-law exponent β and the scaling exponent α (estimated through DFA) of a time series is: Analysis of ECG data. ECG data were used to analyze two measures of peripheral autonomic activity: average heart rate and heart rate variability. For both measures, we used an adaptive threshold to detect the R-peak of each QRS-complex in the ECG. Heart rate was then computed by dividing the total number of R-components by time. Heart rate variability was quantified by means of the detrended fluctuations analysis described for MEG above but now applied to the time series of the intervals between successive R-peaks [27,28]. In line with the MEG analyses, we used windows ranging from 3-50 heartbeats (roughly corresponding to 3-50 s).
Statistical tests. Statistical comparisons of all dependent variables between conditions were, unless stated otherwise, performed using paired t tests.
Null effects are difficult to interpret using regular null hypothesis significance testing. The Bayes Factor addresses this problem by quantifying the strength of the support for the null hypothesis over the alternative hypothesis provided by the data, taking effect size into account. Wherever null effects were conceptually important, results obtained from a regular (paired) t test [114] and Pearson correlations [115] were converted into corresponding Bayes Factors.
To map significant changes of scaling exponents α across the brain, we computed a nonparametric permutation test based on spatial clustering [116,117]. This procedure has been shown to reliably control for Type I errors arising from multiple comparisons. First, a paired t test was performed to identify voxels with significant changes (voxel with p < 0.05, uncorrected). Subsequently, significant voxels are combined into clusters based on their spatial adjacency. Here, a voxel was only included into a cluster when it had at least two significant neighbors. Subsequently, the t-values of all voxels comprising a cluster were summed, which yields a cluster statistic (i.e., a cluster t-value) for each identified cluster. Next, a randomization null distribution was computed using a permutation procedure (N = 10.000 permutations). On each permutation, the experimental labels (i.e., the pharmacological conditions) were randomly reassigned within participants and the aforementioned procedure was repeated. For each iteration, the maximum cluster statistic was determined and a distribution of maximum cluster statistics was generated. Eventually, the cluster statistic of all empirical clusters was compared to the values obtained from the permutation procedure. All voxels comprising a cluster with a cluster statistic smaller than 2.5% or larger than 97.5% of the permutation distribution were labeled significant, corresponding to a corrected threshold of α = 0.05 (twosided).

Model simulations
To simulate the effects of synaptic gain modulation on cortical activity fluctuations, we extended a previously described computational model of a local cortical patch [49] by means of multiplicative modulation of synaptic gain. All features of the model were identical to those of the model by [49], unless stated otherwise. The model consisted of 2,500 integrate-and-fire neurons (75% excitatory, 25% inhibitory) with local connectivity within a square (width = 7 units) and a connection probability that decayed exponentially with distance ( Fig 8A). The dynamics of the units were governed by: where subscripts i, j indicated different units, N ij was a multiplicative gain factor, W ij were the connection weights between two units, S j a binary spiking vector representing whether unit j did or did not spike on the previous time step, and I 0 = 0. For all simulations reported in this paper, we optimized the connection weights using Bonesa [118], a parameter-tuning algorithm, such that the network exhibited alpha-band oscillations, long-range temporal correlations, and neuronal avalanches (see below). The optimized values for the connection weights were W EE = 0.0085, W IE = 0.0085, W EI = −0.569 and W II = −2, whereby subscript E indicated excitatory, subscript I indicated inhibitory, and the first and second subscript referred to the receiving and sending unit, respectively. On each time step (dt = 1 ms), I i was updated for each unit i, with the summed input from all other (connected) units j and scaled by a time constant τ i = 9 ms, which was the same for excitatory and inhibitory units. The probability of a unit generating a spike output was given by: with the time constant for excitatory units τ P = 6 ms and for inhibitory τ P = 12 ms. P 0 was the background spiking probability, with P 0 (exc.) = 0.000001 [1/ms] and P 0 (inh.) = 0 [1/ms]. For each time step, it was determined whether a unit did or did not spike. If it did, the probability of that unit spiking was reset to P r (excitatory) = −2 [1/ms] and P r (inhibitory) = −20 [1/ms]. We used this model to analyze the dependency of two quantities on E/I ratio: (i) the powerlaw scaling of the distributions of the sizes of neuronal avalanches [50] estimated in terms of the kappa index, κ, which quantifies the difference between an empirically observed event size distribution and a theoretical reference power-law distribution, with a power-law exponent −1.5 [51], and (ii) the scaling behavior (scaling exponent α) of the amplitude envelope fluctuations of the model's local field potential. To this end, we summed the activity across all (excitatory and inhibitory) neurons to obtain a proxy of the local field potential. We band-pass filtered the local field potential in the alpha band (8)(9)(10)(11)(12) and computed long-range temporal correlations in the alpha-band amplitude envelopes following the procedure described above (see "Detrended fluctuation analysis"), using windows sizes ranging from 5-30 s.
In order to assess the influence of structural excitatory and inhibitory connectivity on network dynamics (Fig 4D-4F), we varied the percentage of units (excitatory and inhibitory) a given excitatory or inhibitory unit connects to within a local area (7 x 7 units; Fig 8A). These percentages were varied independently for excitatory and inhibitory units, with a step size of 2.5%.
The gain factor N ij was the main difference to the model described by [49]. It was introduced to simulate the effects of neuromodulation on synaptic interactions in the cortical network [36]. For this, we kept all the structural parameters fixed (42.5% excitatory connectivity, 75% inhibitory connectivity; small square in Fig 4D-4F) in a range in which the model exhibits both robust long-range temporal correlations as well as neuronal avalanches. Note that any other combination of parameters would yield similar results, as long as the model exhibits these two phenomena. From the chosen starting point, we systematically varied the synaptic gain factors in two different ways. In the first version, we only varied N EE and N IE to dynamically modulate the circuit's net E/I ratio (Fig 8B) in a way consistent with recent modeling of the effects of E/I ratio on a cortical circuit for perceptual decision-making [44]. In the second version, we varied N EE , N IE , and N EI (S8A Fig). Here, N EI was modulated independently from N EE and N IE , which in turn were co-modulated by the same factor.
Per parameter combination, we ran 10 simulations using the Brian2 spiking neural networks simulator [119]. Each simulation was run for 1,000 s, with a random initialization of the network structure and the probabilistic spiking. Neuromodulation was simulated as a gain modulation term multiplied with excitatory (EE and IE) and/or inhibitory (EI only) synaptic weights. (B) κ as a function of excitatory and inhibitory connectivity (with a spacing of 2.5%; means across 10 simulations per cell). The region of κ~1 overlaps with the region of α > 0.5 and splits the phase space into an excitationdominant (κ > 1) and an inhibition-dominant region (κ < 1). (C) Same as (B), but for scaling exponent α. (D) Same as (B) and (C), but for firing rate. In sum, the alternative version of modulation of excitation-inhibition ratio yields comparable results to the version presented in Fig 8. (E) Model power spectra under different levels of synaptic gain modulation (neuromodulation). The code and the data underlying these plots can be found at https://figshare.com/s/ 374ab0f973f026535549. EE, excitatory-to-excitatory (recurrent excitation); EI, inhibitory-toexcitatory; IE, excitatory-to-inhibitory; II, inhibitory-to-inhibitory (recurrent inhibition); κ, kappa index. (TIF) S1 Movie. Movie file of the structure-from-motion ambiguous visual stimulus used in this study. (MP4)