Role of interneuron subtypes in controlling trial-by-trial output variability in the neocortex

Trial-by-trial variability is a ubiquitous property of neuronal activity in vivo which shapes the stimulus response. Computational models have revealed how local network structure and feedforward inputs shape the trial-by-trial variability. However, the role of input statistics and different interneuron subtypes in this process is less understood. To address this, we investigate the dynamics of stimulus response in a cortical microcircuit model with one excitatory and three inhibitory interneuron populations (PV, SST, VIP). Our findings demonstrate that the balance of inputs to different neuron populations and input covariances are the primary determinants of output trial-by-trial variability. The effect of input covariances is contingent on the input balances. In general, the network exhibits smaller output trial-by-trial variability in a PV-dominated regime than in an SST-dominated regime. Importantly, our work reveals mechanisms by which output trial-by-trial variability can be controlled in a context, state, and task-dependent manner.


Introduction
Trial-by-trial variability is a ubiquitous feature of cortical activity in vivo (Arieli et al., 1996;Churchland et al., 2010). Instead of being just noise, trial-by-trial variability varies (typically gets reduced) during the stimulus presentation (Churchland et al., 2010;Oram, 2011), due to attentional shifts (Kanashiro et al., 2017) or external stimulation (De Luna et al., 2017). The presence and change in trial-by-trial variability are not merely a statistical property of 5 neuronal activity as it is necessary for behavior (Waschke et al., 2021) and affects the stimulus-response (Arieli et al., 1996) and behavioral performance (Arazi et al., 2017;Rowland et al., 2021). Given the noisy inputs, stochastic neurons, random connectivity, and unreliable synapses, trial-by-trial variability is not surprising. However, the modulation of trial-by-trial variability is usually attributed to local network structure. Computational studies have suggested that the modulation of trial-by-trial variability, particularly the stimulus-induced decrease, can be attributed to the multiple interneuron types (Kepecs & Fishell, 2014). These neuron types differ not only in their chemical signature but also in the neuron-type specific connectivity (Jiang et al., 2015). The interneuron diversity renders the local 20 network with rich dynamical and computational properties which have not been observed in the E-I network (Lee et al., 2018;Hertäg & Sprekeler, 2019;Hahn et al., 2022). However, it remains unclear when and how different interneurons contribute to trial-by-trial variability.
In the neocortex, different interneurons receive specific inputs and selectively inhibit pyramidal cells (PCs). For instance, PV-expressing interneurons are mainly driven by feedforward inputs, inhibiting perisomatic regions and 25 basal dendrites of PCs, whereas SST-expressing interneurons are targeted by top-down inputs through VIP cells, inhibiting apical dendrites of PCs (Larkum, 2013). Moreover, each interneuron type can also be the target of specific neuromodulators. Therefore, by characterizing the contribution of different types of interneurons, we may uncover new mechanisms by which the brain can control the response variability in a context, state, and task-dependent.
To understand the role of interneurons in variability control, we investigated using a model neocortical layer 2/3 30 consisting of one type of excitatory neuron (Exc) and three types of inhibitory interneuron (PV: parvalbumin, SST: somatostatin, and VIP: vasoactive intestinal polypeptide expressing cells). We refer to it as the EPSV network. First, we characterized the transfer-function of neurons in the EPSV network. Next, we measured trial-by-trial variability at different operative points. In particular, we focused on the input statistics: variance and covariance of input rate to different neuron types. We show that the main determinants of the trial-by-trial variability are the (1) variance 35 ratio of inputs to the different neuron populations and (2) covariance of the inputs. The effect of input covariance is contingent on the ratio of variances. The effect of these variables strongly depends on whether the network is operating in an SST or PV-dominated regime. In general, network connectivity provides a landscape on which input variability could be transformed into output variability by varying the neuron excitability, input connection strengths, and stimulus properties.

Results
To characterize the role of interneurons in the control of trial-by-trial variability, we study how the trial-by-trial input variability is transferred to the output in a model of neocortex (the EPSV model, see Methods). To this end, we first estimated the transfer-function of a typical neuron in the network operating in an asynchronous-irregular and non-oscillatory state.

Neuron transfer-function
Experimentally neuron transfer-function is typically estimated by injecting direct current at different amplitudes.
However, it is more natural to inject spiking inputs. The number of different types of inputs a neuron may receive depends on the number of neuron types in the network. In a simple scenario consider a network with only one type of excitatory neuron with weak connectivity such that the network remains stable and asynchronous. We injected 50 Poisson-type spiking input at different rates into the network and measured the output firing rates. As expected, in this setting, the output firing rate varies monotonically as a function of the input rate ( Figure 1A). If we had injected inhibitory inputs the firing rate would have monotonically decreased from some baseline firing.
However, when a neuron is a part of a network with excitatory and inhibitory neurons (E-I network) it is important to consider both excitatory and inhibitory inputs separately to determine the neuron transfer-function (Kuhn et al., 55 2004). This results in a two-dimensional neuron transfer-function ( Figure 1B) and a neuron can show identical output firing rate (Exc population) for many different combinations of inputs to Exc and Inh populations (shown as iso-firing Three iso-firing rate surfaces with different output rates. (E) Output space of EPSV network with cubic input space as in D. The three-dimensional scatter plot shows the activity of E, PV, and SST neurons. The dot color indicated the firing rate of VIP neurons. (F) 2-D slices of the transfer-function in D to indicate how neuron firing rate changes a function of the input to two populations while the input to the third is fixed. A negative input rate denotes a reduction of input compared to the baseline (see Supplementary Figure S1). rate contours: Figure 1B solid white lines). The neuron transfer-function also reveals how a neuron may respond if the input is varied while maintaining the excitation and inhibition in balance ( Figure 1B dotted lines).
Extending the two-dimensional neuron transfer-function to a neuron in the EPSV circuit ( Figure 1C), we need 60 to consider the combination of four inputs to the four populations in four-dimensional input space, λ in E , λ in P , λ in S , λ in V (see section Stimulus evoked input). Because the VIP population and SST population are strongly mutually coupled to each other, we simplified by keeping λ in V fixed while only changing λ in S as in (Hertäg & Sprekeler, 2019). We estimated this three-dimensional neuron transfer-function by systematically varying the inputs (see section Stimulus evoked input) to the Exc, PV, and SST populations ( Figure 1D left). The negative and positive inputs are relative to the baseline inputs (see section Baseline input). The baseline inputs were chosen to match the network activity to the in vivo firing rates of the four neuron types (Yu et al., 2019). Rate traces and raster plots of the network working in baseline state is shown in Supplementary Figure S1 A.
As a function of the input to Exc and PV neurons, the neuron transfer-function was similar to the two-dimensional transfer-function (compare Figure 1B bottom and F left). However, an increase in SST inputs resulted in a threshold-70 like transition in the output of the Exc population ( Figure 1F middle and right panels). This is because of mutual inhibitory coupling between VIP and SST neurons. Increasing the firing rate of SST neurons reduced the activity of VIP neurons which further amplifies SST neurons' firing rate and effectively silenced the PV neurons. Therefore, at some level of positive input, the network made a transition to the SST-dominated regime characterized by a sharp reduction in the activity of Exc neurons ( Figure 1F middle and right panels).

75
For the three-dimensional neuron transfer-function, we can also define the manifold on which a neuron's output remains constant despite a change in the input (Figure 1 D right). Next, we rendered the output firing rate of all populations together to visualize all possible network states for a range of inputs to these populations. As expected the recurrent connectivity in the network restricted the possible network states (Figure 1 E). More importantly, this figure shows that when E neurons are firing at high rates, either SST or PV neurons are completely silent. This is 80 a consequence of the switching dynamics of the network which has also been reported earlier (Hertäg & Sprekeler, 2019;Hahn et al., 2022). As we show in the next section, the restricted state space, the shape of iso-firing rate surfaces, and the gradient of the neuron transfer-function play a major role in determining how trial-by-trial input variability is transformed into trial-by-trial output variability. (A) A normal distribution of inputs across trials (x-axis) is transformed into a distorted distribution of outputs (y-axis). Different colors denote different distributions. (B) Output trial-by-trial variability (measured as standard deviation) for different mean and variance acrosstrial inputs in a one-dimensional case. (C) Transfer-function of the E-I network with weak recurrency. The four input clouds (denoted by colors) were sampled from distinct trial-by-trial input distributions. (D) Output rate distribution of corresponding input clouds shown in C for the E-I network. (E) Output trial-by-trial variability as a function of mean input to Exc and Inh populations for a fixed covariance matrix. The red(blue) dot denotes the red(blue) distribution in C. (F) Output trial-by-trial variability as a function of covariance σ in EI and variance ratio σ in E/I (slated bars illustrate the orientation of sampling cloud), for a fixed mean input to the two populations. The red(black) dot denotes the red(black) distribution in C.

Trial-by-trial variability: E-I network
No matter how well we control the experimental conditions, the input to a network always has some trial-bytrial variability. How this input variability is transformed into output variance depends on the gain of the neuron transfer-function in the network. If we consider a single neuron population (Figure 2 A), the task-related input follows a one-dimensional normal distribution (for simplicity) characterized by the mean µ in and variance σ in of the input rates across trials. For this case, output variance depends on both µ in , σ in , and the shape of the neuron 90 transfer-function (Figure 2 A,B).
Distinct from the one-dimensional case, there is more freedom in controlling the trial-by-trial output variability in an E-I network. We can visualize the trial-by-trial input variability as a point cloud in a two-dimensional space spanned by inputs to the Exc and Inh populations (Figure 2 B). Each point indicates the input level to Exc and Inh populations in a given trial. It becomes apparent that besides the mean and the variance of the input to Exc 95 and Inh neurons, the shape and orientation of the input point cloud are also important to determine the output variance. When the input point cloud is elongated (elliptical) and aligned to the iso-firing rate lines, output variance (contributed only by Poissonian fluctuation) will be small, irrespective of the individual input variances.
It suggests that, for the two-dimensional case, we need to consider the complete statistics of the inputs to the two populations. For simplicity, we assumed that the inputs across trials follow a bivariate normal distribution. In this setting, we need to define the following to characterize the output variance: We decomposed the distribution of the input point cloud into three factors: the trial-by-trial covariance σ in EI = σ in In addition, we assumed that the total input variance was fixed σ in EE + σ in II = const. First, we fixed the trial-by-trial covariance and variance of inputs to Exc, and Inh populations and systematically varied the Exc, and Inh input means. We used the neuron transfer-function, estimated from network response during and red dots) unless we are operating at very low evoked firing rates. This observation provides a new explanation of the reduction in the trial-by-trial variability during evoked activity (Churchland et al., 2010). However, from the Figure 2 E we can also derive constraints on Exc and Inh inputs which will increase trial-by-trial variability during evoked activity.
Next, we systematically varied covariance between inputs to Exc and Inh neurons and the variance ratio of inputs 115 to Exc and Inh neurons while keeping the mean of inputs constant (Figure 2 F). This analysis is consistent with the idea that output variance is minimal when the input point cloud is elliptical with a positive slope and aligned to the iso-firing rate lines. It shows that when Exc and Inh inputs are anti-correlated ( Figure 2 black distribution), the input point cloud is orthogonal to the iso-firing rate lines resulting in maximal output variance.

Trial-by-trial variability: EPSV network
Next, we used the same approach described in the previous section to isolate the role of different interneurons in shaping trial-by-trial variability. First, we fixed the covariances and variances (assuming a standard normal distribution) and systematically varied the mean of the input to different neurons. Given that network interactions render the neuron transfer-function nonlinear, the trial-by-trial variance was also highly nonlinear as a function of mean input ( Figure 3A left). Unlike the standard E-I network (Figure 2 E), in the EPSV network for most cases, output variability 125 was positively correlated with output rate (Figure 3 A right). However, it was possible to find input configurations when an increase in output firing rate was associated with a decrease in trial-by-trial variability (as is experimentally observed (Churchland et al., 2010)). But such a reduction was only possible for high firing rates when neurons' transfer-function saturated (Figure 3 A black dots).
We also noted that just based on the gradient of the neuron transfer-function, trial-by-trial variability could be 130 roughly divided into three regimes (Figure 3 A right). Above 5 Hz output rate, the slope of output firing rate vs trial-by-trial variability curve was smallest when the network was operating in a PV-dominated regime ( trial-by-trial inputs to both E-I pairs (input cloud with positive slope), whereas the highest variability was in the region of anti-correlated trial-by-trial inputs to both E-I pairs (input cloud with negative slope).
To analyze the effect of covariance, we assumed a fixed ratio of trial-by-trial input variances to the different neuron populations where the input point cloud was in the region of correlated inputs to both E-P, E-S pairs ( in both σ in EP and σ in ES resulted in a decrease in the output trial-by-trial variability (Figure 3 B bottom row). These results show that the covariance between inputs to the three neuron populations affected the trial-by-trial variability 165 in a state-dependent manner. (A left) Output trial-by-trial variability (quantified as variance) of E population as a function of across-trial input means with a fixed covariance matrix (inputs were sampled from an un-correlated tri-variate normal distribution). The color code is in the right panel.
(A right) Output trial-by-trial variability as a function of output rate of E neurons. Each dot corresponds to a specific input (see the left panel). Red, blue, and green dots correspond to the network in a PV-dominated, SST-dominated, and both influencing regime. Black dots refer to a high firing rate region when output variance decreases when increasing the output firing rate. (B) Output trial-by-trial variability as a function of the trial-by-trial covariance of inputs with a fixed variance σ in EE = 1.8, σ in P P = σ in SS = 0.6. From the top row to the bottom row, the network was tuned to operate in PV-dominated (red) or SST-dominated (blue), or PV-SST driven (green) regimes with similar output rates. (see Supplementary Figure S1 A-C for spiking activity rasters.) (C) Output trial-by-trial variability as a function of the trial-by-trial variance ratio with fixed covariance σ in P S = σ in ES = σ in EP = 0.5. Slanted bars indicate the orientation of the input point cloud. Roman numerals in the B and C refer to the four cases simulated in Figure 5 and

Trial-by-trial variability: Network simulation
To confirm the analysis above derived from interpolated neuron transfer-function, we simulated the stimulusresponse of the EPSV network (see Methods section: Trial-by-trial variability of the input). We tuned the EPSV network in a PV-dominated regime and stimulated the three neuron populations with inputs sampled from a multi- First, we simulated the effect of the ratio of trial-by-trial input variance while keeping other input variables fixed ( Figure 4 A top row B). We assumed that the total variance remains constant (i.e. the sum of diagonal in Figure 4 A top row). The results are consistent with the estimate of trial-by-trial variance derived from the neuron transferfunction. For instance, changing the ratio of the trial-by-trial input variance to the E and P neurons (by an increase of trial-by-trial input variance to the excitatory population while decreasing the variance of input to the PV population) resulted in a decrease in the output variability (Figure 4 A I, II). In this example, a change in the input variance ratios altered the alignment of the input point cloud with the iso-firing rate surfaces, therefore, we observed a decrease in 180 the trial-by-trial output variance.
To further illustrate the effect of the orientation of the input point cloud (or the ratio of input variances), we simulated the network response when the input point cloud had a negative slope (Figure 4 B III, IV). In the PVdominated regime, the slope of the E-S input point cloud affected the output trial-by-trial variability by a small amount (Figure 4 C III). However, when the input cloud slope was negative in the E-P dimensions, trial-by-trial  Next, we varied the trial-by-trial covariance between the inputs to the three neuron populations while fixing all other variables ( Figure 5 A top row). The effect of input covariances is best seen when the input cloud is aligned to iso-firing rate surfaces. Therefore, we tuned the input variance ratios accordingly (see diagonal in Figure

Trial-by-trial variability of inhibitory neurons
In our model, the change of trial-by-trial output variability accompanied different correlations between E, P, and S  Figure 4 A columns I and II). In the PV-dominated regime, the trial-by-trial output variability (E population) was mainly controlled by the high-firing level PV neurons. Therefore, an increase or decrease of output 215 variance accompanied a larger or smaller anti-correlation between the outputs of the E-P pair (compare columns in Figure 5). The network connectivity generated negative correlations between the outputs of E-I pairs across trials.
Such a negative correlation could be enhanced or canceled by the distribution of inputs across trials. Consequently, the interneurons (PV and SST neurons) contributed more or less to output variance (E population).
Whether the variability of inhibitory interneurons changes in the same way as the variability of the excitatory 220 population depends on the network connectivity. In our model, the output variability of inhibitory populations depends on how their iso-firing rate manifold is aligned with that of the excitatory population. To illustrate this, we considered a two-population E-I network. This network can operate in two regimes -weak self-coupling (or strong mutual interactions) and strong self-coupling (weak mutual interactions). For the weak self-coupling regime, E-I populations had a strong mutual interaction such that the change of one population influenced the other (Figure 6 B 225 left column). In this regime, E and I iso-firing rate lines were not aligned, therefore a decrease in the variability of the excitatory population will be accompanied by an increase in the variability of the inhibitory population. A similar argument has been made by (Kanashiro et al., 2017). By contrast, in the strong self-coupling regime, the output rate   Figure 5: Simulation of trial-by-trial variability controlled by input covariance in a PV-dominated regime.
(A,B,C) Same arrangement as in Figure 4. From column I to II (III), trial-by-trial input covariance between E and SST (PV) populations was increased. In column IV, covariances between E-P and E-S pairs were both increased.
of a population did not depend on the input from the other population (Figure 6 B right column), and iso-firing rate lines of both E and I populations were aligned. Therefore, a decrease in the variability of the excitatory population 230 will be accompanied by a corresponding decrease in the variability of the inhibitory population.

Discussion
Here we have investigated how the mean, variance, and covariances of trial-by-trial inputs affect the trial-by-trial variability of the neocortical activity. Transfer of input distribution to output is mediated by the neuron trasferfunction. The number of distinct neuron types in a network determines the dimensionality of the neuron transfer-235 function. For the neocortical microcircuit with one excitatory and two or three inhibitory populations, the neuron transfer-function becomes 3-or 4-dimensional. As the dimensionality of the neuron transfer-function increases beyond one, we need to consider not only the variance of the input but also the ratio of input variances and covariance between all pairs of neuron types. Here, we have unraveled how the ratio of variance and covariance of inputs to excitatory, PV, and SST neurons affect the trial-by-trial variance of cortical activity.

240
For the neocortex network model, the non-linearity of steady state neuron transfer-function resulted in non-trivial relation between output rate and output variance ( Figure 3A). Next, the output variance could be significantly altered without any change in the output rate by tuning the input covariance matrix ( (Figure 3 B, C). In general, the shape and orientation of input distribution across trials should align with the iso-firing rate manifold for the corresponding In this configuration iso-firing rate lines between Exc and Inh populations align.
network state to reach low output variability. Finally, the non-linearity of the neuron transfer-function implies that 245 the input covariance matrix may never fully align with the iso-firing rate manifolds hence the trial-by-trial output variability is inevitable. The orientation of the input cloud (ratio of variances) is the main predictor of trial-by-trial variability as it aligned the input point cloud to the iso-firing rate manifold. Once such an alignment is achieved, covariance can be varied to further modulate the trial-by-trial variability.

Control of trial-by-trial variability 250
Trial-by-trial variability is necessary for behavior (Waschke et al., 2021) therefore it is crucial that it can be varied in a context-, behavioral state-, and task-dependent manner. For instance, during the early stages of learning, we would like higher variability to explore the state space but once the task is learned animals should reduce variability in their behavior.
Our model suggests that feedforward connectivity provides a natural way to alter the input statistics. If the 255 two target populations receive independent inputs (Figure 6 A top), there is no correlation between inputs, and the output variability is solely controlled by the trial-by-trial input variance. In such a scenario, feedforward weights can modulate the input variance: the larger the feedforward weights, the larger the trial-by-trial input variance, and vice versa. When the target populations, e.g. E and P populations, share a certain level of input excitation (Figure 6 A middle) due to feedforward divergent connections from the same sources (thalamus), their inputs co-vary and the 260 slope of the input point cloud is positive. In this scenario, the degree of shared input controls the degree of covariance  (Figure 6 A bottom). Thus, the structure and strengths of feedforward connectivity provide a natural control over the distribution of received inputs across trials.

265
Feedforward input synapses can be learned through plasticity mechanisms so that the animal may reach the desired degree of trial-by-trial variability in their activity and thereby in their behavior. Furthermore, modulation of neuronal excitability and synaptic strength by neuromodulators can provide a context-dependent control of the trial-by-trial variability.
In parallel to the feedforward input structure, changes in the local activity dynamics can also modulate trial-by-270 trial variability by changing the gradient of the neuron transfer-function and the curvature of iso-firing rate curves.
At the simplest, this can be achieved by switching the network between SST-dominated and PV-dominated states.
Furthermore, intrinsic excitability and synaptic plasticity, when changing the network between weak (Figure 6 B left) and strong connection (Figure 6 B right) regimes, can modulate trial-by-trial variability by altering the whole landscape of neuron transfer-function.
Here we have further explored the role of feedforward inputs in shaping the trial-by-trial variability in a network with three different interneuron populations. While focusing on the feedforward input rate covariance and variances, we have not ignored the role of recurrent activity dynamics. The neuron transfer-function responsible for the transfer of input variance and covariances is shaped by the recurrent activity. That is, the neuron transfer-function we have 285 used is not the same as we may estimate by current injections in a silent network (e.g. in vitro slices). When the network is operating in different regimes such as an inhibition stabilized network (ISN), the results may change but only to the extent that in an ISN regime neuron transfer may be quite different (Sanzeni et al., 2020). Thus, the approach of using the neuron tranfer-function allows us to combine both input statistics and recurrent activity state.
Thus, in contrast to previous work where the focus was on the input correlation at the spiking activity level, we 290 showed that interneurons contribute to the rate variability in three ways by modulating: (1) the effective transferfunction of the neuron, (2) the operating regime of the network and (3) the trial-by-trial input variance and covariances.

Function of interneurons
Given the diversity of interneurons in the brain, there is an impetus to identify specific functions of each interneuron 295 subtype (Kepecs & Fishell, 2014). Classically, the function of interneurons is thought in terms of control of gain modulation (Silver, 2010;Isaacson & Scanziani, 2011), control of network activity state (Brunel, 2000), decorrelation of network activity (Renart et al., 2010;Tetzlaff et al., 2012), gating of input (Isaacson & Scanziani, 2011;Kremkow et al., 2010) and control of dynamics of oscillations (Lee et al., 2018;Hahn et al., 2022). Moreover, different interneurons such as PV, SST, and VIP have been implicated in specific functions such as layer-specific control of 300 excitation-inhibition balance (Naka et al., 2019) and synchronization of gamma-band oscillations .
Here we have identified a new role of interneurons -in controlling the trial-by-trial variability. In particular, for the first time, we highlight the importance of the structure of task-related feedforward inputs to the interneurons.
Moreover, multiple interneuron types provide more means to control the trial-by-trial variability. Because the cortical networks can switch between SST and PV-dominated regimes depending on the context and attention levels of the 305 animal, interneurons provide mechanisms to modulate trial-by-trial variability in a context-, task-, and behavioral state-dependent manner.

Model limitations
In the model, we have a number of assumptions e.g. all neurons were modeled as point neurons and connectivity was independent of spatial distance among neurons. However, the most crucial assumption we made is that the 310 network is operating in a near asynchronous-irregular state. In such a state, it is straightforward to relate the input variance and covariance to output variance. But when the network is operating in an oscillatory state (Supplementary Figure S1 D), our single neuron transfer-function approach will not be sufficient. Given the oscillations, we will also have to consider the serial correlations in the input spiking activity and input oscillations. Moreover, inputs with serial correlations or oscillations could affect the downstream network by entrainment or resonance (Hahn et al., 2014).

315
While interesting, this is beyond the scope of the current manuscript.
Next, our conclusion that input covariance and variance ratio are crucial for the control of trial-by-trial variability is contingent on the fact that the iso-firing rate surfaces are continuous. It is possible that due to non-linear interactions, iso-firing rate regions in the input space appear as small islands. In such a scenario, the effect of input variance and co-variance also become contingent on the operating point or the mean output firing rate. 320 We have shown that increasing the covariance of inputs to E and PV or SST neurons reduces the trial-by-trial variability. This is true when input variance is small enough such that the input cloud can align with the iso-firing rate manifold. When input variance becomes larger, or the curvature of iso-firing rate surfaces is high, an increase in the covariance can misalign the input point cloud from the iso-firing rate cloud and result in a higher trial-bytrial variability. Moreover, when local activity is high, it might end up altering the input statistics and render our 325 predictions wrong.

Model predictions and model verification
Like any good computational model, we have also made several simplifications in our model (see Methods). However, our model still captures a number of crucial biological details and makes testable predictions. First and foremost, our results suggest that to better understand the modulation of trial-by-trial variability we should measure 330 the variability and co-variability of inputs to different neuron populations. A straightforward prediction of our model is that when inputs to excitatory and inhibitory neurons are correlated, an increase in the input variance may not increase the output trial-by-trial variability. Recently, it has become possible to elicit behavior by optical stimulation of selected neurons which were also activated by an external stimulus (Marshel et al., 2019). Our model predicts that under optogenetic activation/inactivation the trial-by-trial output variability should be different from that observed 335 during sensory stimulation conditions because input structures are qualitatively different in the two experiments.
Our model also predicts that there will be higher trial-by-trial variability when the network is operating in an SSTdriven state. Interestingly, a reduction in SST activity is often seen in experiments, e.g. non-whisking and whisking in the barrel cortex of mice (Yu et al., 2019). The feedforward inhibition circuit motif ensures that excitatory and inhibitory inputs to a neuron are correlated. This kind of correlation in our model would lead to low trial-by-trial 340 variability. So we predict that blocking of feedforward inhibitory motif should not only alter the response mean firing rate but also increase the trial-by-trial variability.
Neuromodulators can alter synaptic strengths and neuronal excitability. In our model, a change in the feedforward synaptic weights implies a change in the input variance and a change in neuronal excitability implies rotation or shift of the neuron iso-firing rate manifolds. Therefore, the effect of any neuromodulator on trial-by-trial variability should be non-monotonic because only for a specific amount of neuromodulator it may be possible to align the input cloud with the iso-firing rate manifold.

Neuron model
The neurons were realized as the adaptive exponential integrate and firing model (Brette & Gerstner, 2005): .

350
Neuron parameters are provided in the Table 1 (Lee et al., 2018).

Synapse model
Neurons were connected using conductance-based synapses. Each incoming spike resulted in a conductance transient which decayed exponentially with a time constant τ syn , e.g.τ e for excitatory synapses and τ i for inhibitory synapses: where t i is the arrival time of i th spike and H is the heaviside step function.

Network model
The network consists of 4800 neurons with 3600 excitatory, 480 PV, 360 SST, and 360 VIP neurons (Rudy 355 et al., 2011). Neuronal connectivity parameters (see Table 2) were taken from (Hertäg & Sprekeler, 2019) with modifications. Given the size of the network, neurons were connected in a distance-independent manner that each pair of neurons has probability p ts to form a connection depending on the types of source and target as in Table 2 left.
For synaptic conductance, we first chose the value g syn for each pair of connections as in Table 2 right. Next, we scaled each synaptic weight by a number randomly drawn from a log-normal distribution (µ = 0, σ = 1) and upper 360 bounded the weights by 0.5 mV (-2.0 mV) for EPSP (IPSP). The distribution of excitatory conductance is shown in Figure 7. In the main text, all the results are shown for a network with log-normal synaptic weight distribution.

Baseline input
To mimic the ongoing activity, each neuron received homogeneous Poisson spike trains. The rates of spike trains to different neuron-type populations were tuned to obtain a baseline firing rate around 2.5 Hz for excitatory and SST 365 neurons, and ≈14 Hz for PV and VIP neurons (spontaneous activity level of free whiskering mice (Yu et al., 2019)).
Note the baseline network state was PV-dominated in our model (Supplementary Figure S1 top row).

Stimulus evoked input
To mimic stimulus-evoked inputs, each population was stimulated by additional inputs with rates λ in E , λ in P , λ in S , λ in V on top of the baseline inputs as shown in Figure 1. Since SST and VIP neurons are mutually coupled, we assumed 370 that modulatory input to VIP neurons is just inverted input to SST neurons as in (Hertäg & Sprekeler, 2019) hence reduced the input dimension to three with λ in V = 0. We measured the steady-state response of four populations, λ out E , λ out P , λ out S , λ out V , to different levels of modulatory inputs covering a cubic input space. The corresponding output space is restricted due to the interaction between excitation and inhibition as shown in Figure 1. Because the firing rate of the excitatory population was taken as the output, the transfer-function was formulated as λ out (Figure 1). Negative rates imply a reduction in the input from the baseline input. To obtain a transfer-function with higher resolution, we interpolated the simulated transfer-function with (tri)linear function (Weiser & Zarantonello, 1988) implemented in SciPy library (Virtanen et al., 2020). The analysis of variability transformation was derived from this interpolated neuron transfer-function.

Trial-by-trial variability of the input 380
Trial-by-trial variability of the modulatory inputs was modeled as a three-dimensional normal distribution (Figure 3) characterized as: For each setting of mean and covariance matrix, the input rates were sampled (10000-point input cloud) from an arbitrary distribution with given − → µ in and σ in σ in σ in and corresponding outputs were extrapolated from the neuron transferfunction.
We investigated the transformation of input distribution to output distribution by systematically changing the mean − → µ in , trial-by-trial balance σ in EE : σ in P P : σ in SS , and across trial covariances σ in EP , σ in ES , σ in P S (assuming σ in EE + σ in P P + σ in SS = const). A negative ratio of trial-by-trial input variance (anti-correlated across trials) was generated by flipping the sign of modulatory inputs: Covariance matrices being not positive semidefinite were discarded. These are marked as white space in the Figure 3 B,C.

385
To confirm the analysis, we simulated the three network states in Figure 3A right with different covariance matrices.
The covariance matrices were chosen to show the trend of variability modulation regarding trial-by-trial covariance and trial-by-trial balance (ratio of input variances): for the former factor, we chose four settings where both pairs, E-P and E-S, had low trial-by-trial input covariance, E-P had large covariance, E-S had large covariance, and both had large covariance (Figure 5 A upper); for the latter factor, we chose four settings where both pairs had low trial-by-trial 390 input variance ratio, both had high ratio, E-S had negative ratio, and E-P had negative ratio (Figure 4 A upper).
For each network state (mean input) and covariance matrix, we simulated 100 trials for 1000 ms each with 250 ms preparation to reach the operating point, 250 ms to reach the stable state after injecting modulatory inputs, and last 500 ms were taken as the steady-state response ( Figure 5C and Figure 4C). Firing rate traces of four populations and corresponding spike raster are illustrated in Supplementary Figure S1. The result of trial-by-trial variability control 395 for the PV-dominated regime is shown in Figure 5 and Figure 4, and the result for the SST-dominated regime is given in Supplementary Figure S4.

Measurement of output trial-by-trial variability
For each simulation, we used the mean firing rate of each population (E, P, and S) in the last 500 ms as the steady-state response rate. Across 100 trials, we calculated the covariance matrix of their steady-state response rate

Simulation and data analysis tools
The simulations were performed using NEST simulator (Jordan et al., 2019). Differential equations were integrated using a fixed timestep of 0.1 ms. The analysis of simulated network activity was done using customized code written 405 in Python. The results were visualized using matplotlib. The code for simulating the network and visual illustrations will be shared online on GitHub. The network can exhibit stochastic oscillation between populations given strong inputs. In our analysis, we avoided such strong inputs and assumed that the network remained in an asynchronous-irregular activity state. trial-by-trial input point clouds sampled from given covariance matrix as columns I to IV in D top (normalized range).
(F) PSTH of excitatory population response. Black lines: individual trial. Red line: average response over 100 trials.

510
The stimulus was given at 250ms, and the output covariance matrix in D bottom is calculated for the last 500ms.
The trend is consistent with Figure