Narrow and Broad γ Bands Process Complementary Visual Information in Mouse Primary Visual Cortex

Visual Abstract

The role of g synchronization in the processing of visual stimuli has been investigated for decades . Indeed, visual stimulus features are well known for modulating the power and the central frequency of neocortical oscillations in the g range. Previous studies adopting different animal models (e.g., mice, cats, human and nonhuman primates) have indeed highlighted g power to be critically dependent on orientation Berens et al., 2008;Onorato et al., 2020), size (Gieselmann and Thiele, 2008;Perry et al., 2013), speed Friedman-Hill et al., 2000;Womelsdorf et al., 2006), direction (Liu and Newsome, 2006), and contrast (Logothetis et al., 2001;Henrie and Shapley, 2005;Ray and Maunsell, 2010;Saleem et al., 2017;McAfee et al., 2018;Bartoli et al., 2019) of the visual stimulus.
Such cortical g band oscillations originate from the coordinated interaction of excitation and inhibition (Sohal et al., 2009;Cardin, 2016) and might be detected by local field potential (LFP), although single unit activity is irregular (Buzsáki and Wang, 2012). Many modeling works have indeed described how g synchronization originates in excitatory-inhibitory networks and how it is affected by external stimuli (Wilson and Cowan, 1972;Leung, 1982;Ermentrout and Kopell, 1998;Brunel and Wang, 2003;Börgers and Kopell, 2003;Geisler et al., 2005). Specifically, several studies focused on how such networks could capture contrastinduced modulation of g band in primary visual cortex (V1; Mazzoni et al., 2008Mazzoni et al., , 2010Mazzoni et al., , 2011Battaglia and Hansel, 2011;Jia et al., 2013;Roberts et al., 2013;Barbieri et al., 2014;Lowet et al., 2015;Zachariou et al., 2021). These works captured a variety of properties of g band in the visual cortex of primates, but recent evidence showed that properties and the functional role of g band in the visual cortex of rodents might display some peculiar differences. In particular, a very narrow g band (NB) close to 60 Hz has been observed in the lateral geniculate nucleus (LGN) firing activity of awake mice (Saleem et al., 2017;Storchi et al., 2017) but has not been observed in primates. In a recent work (Saleem et al., 2017) this has been shown to induce a cortical narrow band oscillation at ;60 Hz co-existing with the broad band activity in the 30-to 90-Hz range but with a different functional role. While the cortical broad band power increased with contrast, as in primates, the cortical narrow band power increased with luminance and decreased with contrast (Saleem et al., 2017). This suggests that in mice cortical g band might have two components with specific frequency ranges and different encoding properties, marking a relevant difference from the g band in the visual cortex of primates (see also Discussion, NB and BB in primates).
Neurons highly synchronized at 60 Hz were observed in the LGN (Saleem et al., 2017) and the retina (Storchi et al., 2017), suggesting a subcortical origin for such cortical narrow band. Recent modeling works reproducing in detail the structure of the mouse visual cortex (Arkhipov et al., 2018;Billeh et al., 2020) did not specifically address this issue. Overall, the function of the mouse g narrow band in the V1 and the mechanisms underlying its modulation still need to be properly clarified.
To address this issue, we stimulated head-fixed awake mice (n = 12) with alternating gratings of different contrasts and performed a spectral analysis of the resulting V1 LFPs. We showed that narrow and broad g band (BB) not only have opposite but even complementary sensitivity range. Simulations with recurrent excitatory-inhibitory spiking network accurately reproduced the behavior of both narrow and BB, supporting the origin of the former might be because of a periodic component embedded within the thalamic input.

Experimental design Mice
Experiments were conducted in accordance with the European Community Directive 2010/63/EU and were approved by the Italian Ministry of Health. Animals were housed in a 12/12 h light/dark cycle with food and water available ad libitum. Adult (four to six weeks old) C57BL/ 6J female mice (n = 12) were used in all experiments until contrast level equal 50; only seven of them were instead employed for the maximum contrast level (K = 90) experiments.

Recording implant
Animals (n = 12) were chronically implanted with a custom-made aluminum head post, and a rectangular recording chamber (2 Â 1.5 mm) of dental cement (Ivo-clar Vivadent Inc.) was built over the V1 (i.e., between 0 and 1.5 mm anterior and between 1.5 and 3.5 lateral to the l suture) leaving the skull intact. A ground electrode was placed over the cerebellum. The electrode was connected to a pin socket and secured to the skull by acrylic dental cement. Surgery was conducted under deep avertin anesthesia (7 ml/kg; 20% solution in saline, i.p.; Sigma-Aldrich). Animals were then allowed to recover for 3 d. Following recovery, animals were habituated for 3 d to the head fixation apparatus. A craniotomy overlying V1 was performed 24 h before the first recording session. To preserve the cortical surface, the recording chamber was filled with a layer of agar (Sigma-Aldrich) and the silicone elastomer Kwik (World Precision Instrument) as a protective cap. In order to discard nonvisually evoked neural confounding, animals were restrained from moving while in the head fixation apparatus.

Extracellular recordings in awake mice
Recordings were performed on awake mice (n = 12). Mice were carefully placed in the head fixation apparatus. After removing the protective cap, the recording chamber was filled with sterile saline solution (0.9%) to preserve and moisten the tissue.
A NeuroNexus Technologies 16-channel silicon probe (Fig. 1A) with a single-shank (A1x16-3 mm-50-177) was mounted on a three-axis motorized micromanipulator and slowly lowered into the visual cortex (in the central region of the recording chamber) till the depth of 1000 mm. Before the beginning of the recording, the electrode was allowed to settle for ;5 min. The neurophysiological data were continuously recorded using a 16-channel Omniplex system (Plexon). At the end of the extracellular recording session, the recording chamber was covered with the protective cap as described above.
Each of the twelve animals underwent two recording sessions on two different days (one per day). Each experimental session was made up of at least two visual contrast sweeps. Specifically, each sweep consisted in presenting the mice with five contrast levels [0 10 20 30 50] % (n = 5 animals) or six contrast levels [0 10 20 30 50 90] % (n = 7 animals) in increasing order. Each visual stimulus consisted of 1-Hz alternating gratings presented for 30 s at a single contrast level to the headfixed animals (see below, Visual stimuli). Overall, for each level of contrast below K = 90 we have two recording sessions Â (!2 sweeps) Â 12 animals ! 48 trials of 30 s each. For k = 90, we have two recording sessions Â (!2 sweeps) Â 7 animals ! 28 trials of 30 s each. Few recordings were removed because of artifacts. The actual numbers of trials used for every analysis is reported in the main text.
All the electrophysiological data were processed using custom-written MATLAB codes (The MathWorks) on a Windows 10 Home Dell XPS 15 laptop. The extracellular signals were sampled at 1 kHz and lowpass filtered at 200 Hz (filter type: Bessel, four poles).

Visual stimuli
Visual stimuli were computer-generated using the MATLAB Psychophysics Toolbox with g correction and presented on a display (Sony; 40 Â 9 Â 30 cm; mean luminance 15 cd/m 2 ) placed 25 cm from the head of the mouse (Fig. 1A), covering the center of the visual field. Extracellular signals were recorded in response to abrupt reversals (1 Hz) of vertical square-wave gratings (spatial frequency, 0.06 c/°; contrast levels adopted: [0 10 20 30 50 90] %). Signals were amplified (5000-fold), bandpass filtered (0.5-500 Hz), and fed into a computer for storage and analysis. For each recording, alternating gratings were presented for 30 s at a single contrast level to the head-fixed animals. To ensure consistency across animals (n = 12 for K 50 and n = 7 for K = 90), the orientation of the gratings was always vertical. It is worth mentioning that there is little evidence for columnar organization of orientation-selective neurons in the V1 of the mouse (Ko et al., 2011).

Neurophysiological data analysis Current source density (CSD) analysis and Layer IV identification
For each channel, visual evoked potential waveforms in response to contrast reversals were extracted from the LFPs by signal averaging. For each recording session, CSD was computed by applying a standard algorithm (according to the second spatial derivative estimate of the laminar LFP time series; Haberly and Shepherd, 1973;Freeman and Nicholson, 1975) along with the iCSD toolbox for MATLAB (Pettersen et al., 2006). A value of 0.3 S/m was taken as a measure of cortical conductivity (Gaussian filter: SD = 0.05 mm). We focused our analysis on Layer IV as it is the layer where sensory-induced g oscillations are prominent (Welle and Contreras, 2016). Layer IV was identified in each recording session (two for each of the Square-wave 1-Hz alternating gratings at different contrast levels were used for visual stimulation. A linear 16-channels probe (with 50-mm spacing between electrodes) was inserted into the mouse V1. Right, Mean across animals of current source densities (CSDs) aligned by the earliest current sink. B, left, Mean scalogram for contrast equal to 0 (K=0) within À200 to 900 ms around screen onset (for K = 0). The dashed rectangle depicts the frequency bands ranges: narrow band (NB, middle, black) and broad band (BB, gray). Right, Scalogram magnitude comparison between screen OFF (from À200 to 0 ms) and screen ON condition (from 0 to 900 ms) in high NB (top) and BB (bottom). Each dot represents one stimulus presentation for one mouse (n = 12 animals, 70 experimental points overall). Statistical differences were accounted for by the Wilcoxon's matched pairs signed-rank test. For both g bands, p values were far less than 0.001. C, Examples of filtered local field potential (LFP) recorded in mice V1 while viewing a uniform gray screen (top), or alternating gratings at contrast K = 30 (middle) and K = 90 (bottom). The examples are reported between [À0.5; 1] s around screen onset. Examples were filtered between (1) 10-100 Hz (black traces) just for representative purposes; (2) 45-65 Hz to display the NB; (3) from 20 to 45 Hz and from 65 to 90 Hz to display the BB range. Dashed lines indicate screen onset (black) and the first contrast reversal (gray). Monitors' sketches schematically represent visual contrast. D, LFP modulation of minimal contrast (i.e., K = 0; top) and maximal contrast (i.e., K = 90; bottom) with respect to the power spectral density (PSD) at K = 30 (middle). Modulation is defined as the difference between the power of a frequency at a given contrast level (K = 0 or K = 90 in this case) with the power at reference contrast K = 30, normalized to the latter power. Shaded regions indicate standard error of the mean (SEM). 12 mice) with the channel corresponding to the earliest current sink.

Spectral analysis
From extracellular recordings, we extracted LFPs by low pass filtering at 200 Hz ( Fig. 1C for filtered LFP examples). LFPs were z-scored before spectral analysis. The power spectral density (PSD) of the z-scored LFPs was computed with the Fast Fourier Transform via the Welch method (pwelch function in MATLAB), dividing the time window under investigation into subwindows of 500 ms with 50% overlap. The only exception was when investigating the temporal structure of narrow and broad band PSD (Fig. 2C): z-scored LFPs were segmented in consecutive temporal windows after each contrast reversal. PSDs were, consequently, estimated independently for each of them and averaged across trials. Given the short duration of the time windows (200 and 300 ms) the Welch method was applied with no subwindow.
Since the contrast K = 30 minimized both narrow and broad band g power (see Results), we expressed the spectral response at other contrasts as the spectral modulation relative to the response to this stimulus (Fig. 1D, top and bottom): Where PSD K=30 is the median LFP PSD across recordings and animals when contrast K = 30 is presented (Fig.  1D, middle) and PSD K' is the median LFP PSD across recordings and animals when contrast K = K' is presented. The alternative option of computing the modulation relatively previsual baseline LFP led to similar results (data not shown). . Asterisks indicate significant post hoc pairwise difference (K-W followed by Dunn's test p , 0.05). n.s. indicates nonsignificant statistical difference across contrasts. Individual experimental data points are indicated as gray circles. Lower boxes display the mean difference in PSD, DPSD, between a reference contrast (K = 0 for NB and K = 30 for the BB) and the other contrast levels. Black whiskers indicate the 95% confidence interval obtained through nonparametric bias-corrected bootstrap resampling (Ho et al., 2019). C-E, Same as B for narrow band (C), low broad band (D), and high broad band (E) dividing the response interval in two windows:   We investigated the evolution in time of the LFP by means of wavelet analysis. LFP scalograms were computed using the continuous wavelet transform (cwt function in MATLAB), using the analytic Morse wavelet with symmetry parameter equal to three and the time-bandwidth product equal to 60. Scalograms were separately computed for each experimental recording and then split into 500-ms consecutive time windows: specifically, from À100 to 400 ms around every square-wave contrast reversal. We decided not to include in the analysis of the response to visual contrasts the LFPs evoked by the visual stimulus onset (i.e., the first 500 ms of the recorded LFP after stimulus onset also called screen onset hereafter) since it induced flash-like responses as the prestimulus consisted in a dark screen (Fig. 1B). As for the PSD modulation analysis ( Fig. 2 and Fig. 3), the mean scalogram of the LFPs at K = 30 was the reference time-frequency map on relatively to which we computed the modulation of the other contrast levels (Fig. 4A). The scalogram modulation was computed by averaging the modulated time-frequency maps of every recording at a given contrast level (segmented in 500-ms windows around the grating reversal as described above). We computed then the median of the modulated scalogram of every recording within the narrow or the broad band to compute the time evolutions of these two g bands.
To define the broad and narrow band frequency limits, we computed the modulation of the mean scalogram across trials and animals of K = 90 with respect to K = 0 ( Fig. 2A). We defined broad band (BB) the set of frequencies for which the modulation exceeded the 20%, i.e., the two disjoint ranges  and  Hz. We defined NB, the g band interval between these two ranges, i.e.,  Hz.
Next, we evaluated the onset latency of narrow band activation. We computed, for every recording, the onset time of the intertrial median NB modulation amplitude (Bartoli et al., 2019). First, for each NB signal, we marked the first time point at which the NB amplitude exceeded a threshold (,25th percentile for at least 40 ms in the time window between À100 and 400 ms; ;40% of the trials were discarded as they did not meet this criterion). Next, a 200-ms wide window was extracted around that time point (100 ms before and 100 after). This window was segmented into 50-ms bins with 80% overlap and a linear regression was fit to each bin. The first time point of the bin with the highest slope and smallest residual error was defined as the onset of the narrowband g oscillations.
The same procedure was repeated to estimate broad band onset latency, with a threshold of .75th percentile because of the opposite polarity of the modulation.

Spectral information analysis
To determine how well the PSD modulation of the narrow band encoded the contrast level of visual stimuli, we computed the mutual information (Shannon, 1948) I(K; NB mod ) carried by the narrow band PSD modulation with respect to K = 30 throughout the whole intercontrast-reversal interval, NB mod , about the set of contrasts K: Pðnb mod ; kÞ Pðnb mod Þ ; (2) where P(K) is the probability of contrast K to be presented, P(nb mod ) is the probability distribution of the narrow band modulation over all the contrasts, and P(nb mod , k) is the probability of the NB modulation nb mod to be observed when contrast stimulus K is presented. The probabilities, both marginal and conditional, were computed by discretizing the power modulations into 10 equipopulated bins. Specifically, the discretization boundaries were chosen so that each bin contains the same number of elements. This discretization choice, in place of an "equalwidth" binning procedure, ensures robustness against outlier and maximizes the response entropy (Timme and Lapish, 2018). The same analysis was performed for the broad band. Because of the small number of trials for single animals, probabilities were estimated over all trials and recordings. Similarly, we also computed the mutual information carried by NB and BB when considering the PSD modulation within 200 ms or from 200 to 500 ms following contrast reversal.
Limited dataset bias was accounted for by applying the Panzeri-Treves correction (Panzeri and Treves, 1996).
We further investigated the amount of information carried by the two g bands about low or high contrast levels (low: K , 30 and high: K . 30). Low contrast information was computed as described above, after random permutation of high contrast-response association (average over 500 permutations), and vice-versa. In this case the bias was accounted for by applying bootstrap correction and quadratic extrapolation (Panzeri and Treves, 1996;Magri et al., 2009).
Finally, the synergistic contribution in the mutual information when jointly considering both narrow and broadband g modulation was computed as: where I(K; NB K , BB K ) is the joint information carried by both narrow and broad band modulation, I (K; NB K ) 1 I(K; NB K ) is the sum of the information carried about visual contrast level by narrow and broad band g modulation, respectively, and Red [I (K; NB K ), I (K; BB K )] is the amount of redundancy (i.e., the amount of information overlap existing between the two g bands power about the stimulus levels).
We have employed as a measure of redundancy the one proposed by Williams and Beer (2010), namely the minimum information: redundancy is expressed as the expected value of the minimum information that any response variable carries individually about each one of the stimuli.

Red
IðK; NB K Þ; IðK; Significant information (p = 0.05) was estimated using the 95th percentile of bootstrap information. All information quantities were computed in MATLAB with Information Breakdown Toolbox (Magri et al., 2009).

Network model Spiking network model of mice V1
The simulated network is composed of N = 5000 leaky integrate and fire (LIF) neurons (Tuckwell, 1988): 80% excitatory neurons with AMPA-like synapses, and 20% inhibitory neurons with GABA-like synapses (Braitenberg and Schüz, 1991). The network is sparse and random, the connection probability between any directed pair of cells being 0.2 (Sjöström et al., 2001;Holmgren et al., 2003;Fig. 5A, middle).
The membrane potential V k of each neuron k evolves according to (Brunel and Wang, 2003): where t m is the membrane time constant, g leak is the leak membrane conductance (see Table 1 for values) and I k tot ðtÞ is the total synaptic input current. The latter was given by the sum of all the synaptic inputs entering the kth neuron: Where C jk = 0 if j projects to k, and I k AMPA ðtÞ,I k GABA ðtÞ,I k EXT ðtÞ the different synaptic inputs entering the k-th neuron from recurrent AMPA, GABA, and external AMPA synapses, respectively.
The synaptic inputs currents were modeled as: where g syn are the synaptic conductances (Markram et al., 1997;Gupta et al., 2000;Bartos et al., 2001Bartos et al., , 2002; see Table 1), and V syn are the reversal potential of the synapses (V GABA = À80 mV and V AMPA = 0 mV). The function s syn (t) described the time course of the synaptic currents and depends on both the synapse type and on the kind of neuron receiving the input. Specifically, every time a presynaptic spike occurred at time t*, s syn (t) of the postsynaptic neuron was incremented by an amount described by a delayed difference of exponentials (Brunel and Wang, 2003): where the latency t l , the rise time t r and the decay time t d are listed in Table 1 for each synapse type (Xiang et al., 1998;Zhou and Hablitz, 1998;Angulo et al., 1999;Gupta, 2000;Kraushaar and Jonas, 2000;Bartos et al., 2001). The LFP of the simulated network was estimated as the sum of the absolute value of the GABA and AMPA currents (both external and recurrent) that enter all excitatory neurons (Mazzoni et al., 2015). Simulated LFP was analyzed with the same procedures followed for experimental LFP (see above).

External inputs to the simulated network
For all neurons the external input I k EXT ðtÞ is the sum of two terms: a noisy excitatory external input representing the activity from thalamocortical afferents (Fig. 5A, left) and colored noise mimicking stimulus-unspecific cortical activity ( Fig. 5A, right). This simulated external input was implemented as a series of spike times that activated excitatory synapses with the same kinetics as recurrent AMPA synapses, but different strengths (Table 1). These synapses were activated by independent realizations of random Poisson spike trains, with a time-varying rate identical for all neurons. Based on experimental recordings from mice LGN (Saleem et al., 2017;McAfee et al., 2018), this time-variant rate was given by the superimposition of a constant component and a periodic term at g frequency: v ext ðtÞ ¼ ½SðK; tÞ 1 AðK; tÞ« g ðtÞ 1 # n nðtÞ 1 ; where K is the contrast level, S(K, t) the sustained input rate, A(K, t) the amplitude of g range filtered white constant noise « g (t), and n(t) is the colored noise. The first two terms represent thalamic inputs (Fig. 5A, left). Specifically, « g (t) was obtained by applying a third-order bandpass Butterworth filter of central frequency equal to 57 Hz and bandwidth equal to 10 Hz to white noise. The only exception was when investigating the relationship between cortical and thalamic NB peak frequency (Fig.  5B, middle). In this case, we modulated the central The noise term n(t) is a z-scored colored noise, with the PSD following SðfÞ ¼ 1 f a , with a = 1.5, and an amplitude factor # n ¼ 0:4 sp: ms . Both « g (t) and n(t) were independently generated at every simulation.
to avoid a negative number of spikes which could arise because of the noise terms. In the first part of our work (Figs. 5, 6), we set the external input parameters to be time-invariant, i.e., A(K) and S(K) with a value for each K reported in Table 2 (see parameter selection). In particular, setting A ¼ 0 sp: s and S ¼ 500 sp: s the PSD of the simulated LFPs closely matched the reference experimental spectrum for K = 30 (Fig. 1D, middle; reduced x 2 between experimental and simulated LFPs PSD X 2 r ¼ 0:13). The spectral modulations of the other simulated contrast levels were consequently assessed using this median LFP PSD as a baseline, analogously to the approach adopted for the experimental data.
The values of A(K) and S(K) across K were approximated with two piecewise linear functions (Fig. 6A,B): and SðKÞ ¼ ð384:813:8KÞuðK À 30Þ; where K is the contrast level intended to be simulated and u(K) is the Heaviside step function: u(K) = 0 for u(K) , 0 and u(K) = 1 for u(K) ! 0.
In the second part of the work (Fig. 7), we also took into account the evolution in time of the external input parameters. As for the amplitude of the g range filtered white noise, i.e., A(K, t), we defined it as a function of the time of grating reversal (see above, Visual stimuli): where A 0 (K) is the baseline value, a(K) is the amplitude of the reversal-driven modulation f(t), t* is the time of grating reversal and Dt is 40 ms to mimic the latency of the narrowband observed experimentally. The finite rise time was estimated to ;10 ms ( Fig. 7C,D) and for the sake of simplicity was approximated with a linear growth. The decay time was defined to reproduce the quick decay observed experimentally, with the transient effects ending after ;200 ms ( Fig. 7). The values of A 0 (K) and a(K) are reported in Table 2 (see inputs parameter selection). Similarly, we set the sustained thalamic input S(K, t) to: where S 0 (K) is an additive constant and b (K) is the amplitude of the same reversal-driven modulation f(t) defined in Equation 13. The obtained time course of the thalamocortical afferents spike rate qualitatively matches the one observed experimentally by McAfee et al. (2018). The values of S 0 (K) and b (K) are reported in Table 2 (see below, Input parameters selection).

Input parameters selection
Starting from experimental observations (Figs. 1-4) we defined the synaptic inputs A(K) and S(K) to be complementary.
sulting synthetic LFP spectra closely matched experimental results (see above). For all K . 30 A(K) was set to be zero and for all values K , 30 S(K) was set to be 500 sp: s . To find the optimal parameters AðKÞ for K = [0, 10, 20] we simulated the LFP generated by inputs with values of A(K) ranging from 0 to 100 sp./s in steps of 10 sp./ms. For each input, we estimated the spectral modulation relative to K = 30 and we estimated the agreement with the experimental LFP spectral modulation as follows: where F is the total number of frequencies, mod sim (f, S, A) is the median modulation across simulations of the LFPs when setting the external input parameters to S and A (see Eq. 9), and mod exp (f, K) and s exp mod ðf; KÞ are the median and the standard deviation of the experimental LFPs modulation at a given contrast K across animals and trials. For each K, the value of A(K) minimizing this measure was selected (see Table 2).
We set the optimal values of S(K) for K = [50, 90] (see Table 2) in a similar way, simulating the LFP generated by inputs in which S(K) varied from 500 to 1000 sp./s in steps of 100 sp./s.
Similarly, as we found no significant time structure in the experimental scalogram at K = 30 (Fig. 4A), for the time-variant model we set S 0 (K = 30) = 500 sp./s, A 0 (K = 30) = a(K = 30) = b (K = 30) = 0 sp./s. As we have experimentally observed narrow and broad band to be complementary modulated in two different contrast ranges (Fig. 4), for K , 30 we set b (K) = 0 sp./s and S 0 (K) = 500 sp./s, whereas for K . 30 we set a(K) = A 0 (K) = 0 sp./s. Therefore, to find the optimal values of S 0 (K) and b (K) parameters for K = [50,90], we simulated the LFP generated by inputs with values of S 0 (K) ranging from 500 to 1000 sp./s in steps of 100 sp./ms and of b (K) ranging from 0 to 400 sp./s in eight uniform steps. For each input pair combination, we, therefore, computed the scalogram modulation relative to K = 30, we extracted the BB time evolution and we estimated the agreement with the experimental LFP scalogram modulation as follows: where T is the total number of time points, mod sim BB ðt; S 0 ; b Þ is the median modulation of the simulated time course of the broad band g range when setting the external input parameters to S 0 and b , mod exp BB ðt; KÞ and s exp mod BB ðt; KÞ are the median and the variance modulation of the experimental time course of the broad band g range at a given contrast level K.
Similarly, for K=[0,10,20] we simulated the LFP generated by inputs with values of A 0 (K) ranging from 0 to 100 sp./s in steps of 10 sp./ms and of a(K) ranging from 0 to A 0 (K) in steps of 10 sp./s. For each input pair combination, we, therefore, computed the scalogram modulation relative to K = 30, we extracted the NB time evolution and we estimated the agreement with the experimental LFP scalogram modulation as follows: where T is the total number of time points, mod sim NB ðt; A 0 ; aÞ is the median modulation of the simulated time course of the narrow-band g range when setting the external input parameters to A 0 and a, mod exp NB ðt; KÞ and s exp mod NB ðt; KÞ are the median and the variance modulation of the experimental time course of the narrow-band g range at a given contrast level K.

Phase analysis
To investigate the mechanisms underlying the spectral proprieties of the simulated LFPs, we investigated the relationship between the phase distribution of the LFPs and their corresponding thalamic inputs within the narrow and the BB (Fig. 5B,C). After filtering in the frequency band of interest, LFP and thalamic input phases were estimated via the Hilbert transform. The circular variance of the difference between these phase distributions was computed using the Circstat toolbox in MATLAB (Berens, 2009).

Simulations
Network simulations were performed using a finite difference integration scheme based on the second-order Runge-Kutta algorithm (Hansel et al., 1998;Shelley and Tao, 2001;Press, 2007) with time step Dt = 0.05 ms. To focus on stationary responses, the first 200 ms of every simulation were discarded.
Simulations with time-invariant inputs lasted 10 s and each set of inputs was presented 25 times with different noise terms « g (t) and n(t).
Simulations with time-dependent inputs lasted 2 s: during the first second, the input was set to the baseline values of A ¼ 0 sp: s and S ¼ 500 sp: s ; during the last second, instead, the input was determined by the level of contrast (see Table 2). Each set of input was presented 100 times with different noise terms « g (t) and n(t).
All simulations were conducted with custom made Python scripts within the Brian 2 simulator environment (Goodman, 2008;Stimberg et al., 2019) on a Windows 10 Home Dell XPS 15 laptop.

Statistical analysis
Data processing and statistical analysis were carried on with custom-made MATLAB scripts and available thirdparty data analysis Toolboxes (e.g., for computing mutual information). We employed both custom scripts and built-in data analysis and statistical functions. For each statistical comparison throughout the text, we reported the statistical test and their p values; p values lower than 0.05 were considered significant. All results will be reported as median 6 SEM unless otherwise stated. As for the results presented in Figure 2, we adopted the estimation graphics package presented in (Ho et al., 2019).

Code accessibility
The code described in the paper used for simulating the LFPs is freely available online at https://github.com/ nicolomeneghetti/mice-v1-narrow-broad-gamma-band. The code is available as Extended Data 1 contains the python code for running the simulations described in the manuscript.

Results
Broad and NBs in mice V1 display distinct sensitivity to contrast We analyzed the spectral proprieties of Layer IV V1 LFPs recorded from awake mice presented with alternating gratings with different levels of visual contrast (see Materials and Methods; Fig. 1A). The animals (n = 12) were head-restrained during the recordings and were not allowed to move to avoid movement-driven modulations of the g band (Niell and Stryker, 2010;Lee et al., 2014;Saleem et al., 2017). Coherently with previous results (Saleem et al., 2017), we found that the spectral modulation of the LFP with such stimuli was not uniform over the whole g band The subsequent response (i.e., only to the reversal of contrast gratings) of NB was instead maximal for low contrasts (K = 0; Fig. 1C, top) while for BB was maximal at high contrasts (K = 90; Fig. 1C, bottom). As a consequence, the two bands displayed opposite spectral modulation compared with the intermediate contrast level (K = 30; Fig. 1D). To furtherly highlight this result, in the following we will adopt as reference contrast the value K = 30 (whose PSD is reported in Fig. 1D, middle) on which to compute the spectral modulation relative to the other levels of visual contrasts (see Materials and Methods).
Broad and narrow band are modulated by complementary contrast ranges over different temporal windows To furtherly investigate the differences between the processing of visual information of these two g bands, we quantified the response over time of the narrow band on one hand, and the high ([65-95] Hz) and low ([20-45] Hz) broad band on the other ( Fig. 2A).
First, we tested whether the sensitivity was the same over the whole range of visual contrasts explored when considering these g band responses over the whole intercontrast-reversal interval (i.e., 500 ms). We found that NB at very low contrast (K = 0) had a significantly higher power than at intermediate contrast ( . This indicates that NB is not steadily decreasing for increasing contrast, but it is decreasing when contrast goes from low to intermediate level. Conversely, low BB power was uniform for levels of contrast K 30 (K-W ANOVA1: p = 0.02, post hoc Dunn's test, p . 0.5), while it increased significantly at high contrast (K ! 30; K = 90, 32 recordings across 7 animals: [À28.53 6 0.74] dB/Hz vs K = 30, 79 recordings across 12 animals: [À31.26 6 0.49] dB/Hz, post hoc Dunn's test, p = 0.04; Fig. 2B, middle). High range of BB ([65-95] Hz) was finally not significantly modulated by contrast when we considered the whole intercontrast-reversal interval (K-W ANOVA1: p = 0.46; Fig. 2B, bottom). This suggests that NB and BB (at least its lower frequency component) not only display opposite correlations with contrast, but they are sensitive to complementary ranges of visual contrast levels.
Overall, these results show that: (1) NB is only modulated by the low contrast range and is primarily modulated in the late response (i.e., by the presentation of a static visual contrast); (2) conversely, BB is only modulated by the high visual contrast range in the early phase of the response (i.e., in response to the abrupt reversal of the contrast grating).

Broad and narrow band processing of visual information
The complementary sensitivity to visual contrast of narrow and broad band was furtherly investigated by estimating the information carried about the visual contrast level by the power of these two g bands (for details, see Magri et al., 2009 and Materials and Methods). As high and low broad band have displayed a similar behavior, for this analysis (and in the following ones as well) they will be considered as one BB (as in Saleem et al., 2017).
Considering the whole intercontrast-reversal interval, NB carried significant information (p , 0.05, bootstrap test; Materials and Methods) over the whole range of contrasts and for the low contrast range (K 30), but did not carry significant information about contrasts K . 30 (Fig.  3A). Interestingly, NB did not carry significant information at all when considering only modulations occurring in the first 200 ms after contrast reversal (Fig. 3B) proving that significant NB modulations occur only at a later stage. Indeed, information carried by NB in the late stage (from 200 to 500 ms after contrast reversal) about the visual contrast level was significant for all ranges considered, although the information about low level of contrast was higher (0.11 bits for K , 30 vs 0.08 bits for K . 30; Fig.  3C). This therefore proves that NB can be considered an information channel primarily dedicated to convey information about low levels of static visual contrast.
BB carried instead significant information only about contrasts K . 30 (Fig. 3D) when considering the whole intercontrast-reversal interval, suggesting that BB is instead an information channel dedicated primarily to convey information about high levels of contrast, complementarily to NB. The distribution of information over time of BB modulation was furthermore exactly complementary to NB, with significant information in the early response for all ranges (0.07 bits for K , 30 vs 0.11 bits for K . 30; Fig. 3E) and no significant information carried by late modulations (from 200 to 500 ms after contrast reversal; Fig. 3F).
Coherently with this complementarity, spectral modulations of the two bands displayed a synergistic contribution in the encoding of visual contrasts (0.29 bits over the whole stimulation window and contrast range), such in a way that, when considering the modulations of both bands, the large window carries significant information about all contrast ranges (Fig. 3G). Moreover, both early and late responses carry significant information (Fig. 3H,I). This result depends critically on the complementarity of the sensitivity ranges: if contrast-induced modulations in NB and BB were the opposite, with NB decreasing over the whole range and BB increasing over the whole range, the two bands would carry highly redundant information.
Altogether these results show NB and BB to be complementary information channels, with the former encoding increase in contrast from K = 0 to K = 30 by decreasing power .200 ms after contrast reversal, and the latter encoding increase in contrast from K = 30 to K = 90 by increasing power right after contrast reversal.

Narrow and broad band c dynamics at contrast reversal
As we observed a clear temporal structure of the two g bands, we focused on the temporal evolution of NB and BB modulation (see Materials and Methods). As expected, narrow band modulation at K = 0 showed a sustained activity throughout the whole stimulation with no specific response associated with contrast reversal (Fig.  4B). At null visual contrast level, indeed, the contrast reversal does not produce any change in the visual stimulus. For K = 10, NB power instead decreased after the grating contrast reversal (minimum NB modulation: [À4.5 6 2.3] % reached at [107 6 5] ms; Fig. 4C). This decrease in NB power lasted for ;200 ms (K = 10 offset: [198 6 10] ms).
Of note, NB and BB displayed different temporal dynamics relative to contrast reversal, with the former

Spiking neuronal network captures distinct functional c bands
To model the modulation with contrast of these functional g bands in V1, we took into account the presence of a NB oscillation in rodent thalamic activity (Saleem et al., 2017;McAfee et al., 2018;see Materials and Methods), by incorporating it in the input of a spiking neuronal network for which it was shown the g BB to be modulated by the input intensity (Mazzoni et al., 2008). Briefly, the simulated input was given by the superimposition of two components: a sustained component with an average intensity depending on contrast [S(K)], and an oscillatory component with a fixed frequency of ;60 Hz and a contrast-dependent amplitude [A(K); Fig. 5A, left; Eq. 9].
The oscillatory component of the thalamic input A(K) of the model generated a NB peak in the simulated V1 LFP spectrum (Fig. 5B, left) with a strong correlation coefficient between the cortical NB peak power and the amplitude A of the thalamic NB (Pearson correlation, R = 0.94, p ( 0.001). Note that V1 NB peak frequency is exactly the same as A dominant frequency (Fig. 5B, middle). Moreover, input and output oscillations were tightly phase-locked in the NB but not in the BB range (Fig. 5B,  right). A two-way ANOVA analysis of the circular variance for factors (frequency band, amplitude of thalamic input) found a significant interaction (p = 0.003) between the two factors. The NB band was significantly modulated by the amplitude of the periodic thalamic component A in the LFP NB (K-W test p ( 0.001) but not the LFP BB (K-W test p = 0.8). As the circular variance in the NB diminished as the thalamic periodic input increased in amplitude, the NB circular variance was consequently found to be significantly (p = 0.002) less than the BB circular variance for A . 30, hence indicating a tighter phase locking. This shows that in our model the thalamic periodic input (modulated by the parameter A) generates a phase locked periodic response in V1 with the same peak frequency.
The sustained component of thalamic input S(K) modulated instead to the onset and the strengthening of cortically generated g oscillations, as previously seen in this kind of model (Mazzoni et al., 2008), resulting in an increase of power in frequencies involving both BB and NB (linear correlation coefficient between BB power of the simulated LFPs and the sustained component of the thalamic input equal to 0.98, p = 0.004; Fig. 5C, left). Interestingly, phase-locking between V1 NB ad BB and the associated components of the input decreased when modulating the sustained thalamic input strength (see the  (2) sparse LIF network of excitatory (n = 4000, red) and inhibitory neurons (n = 1000, blue). The size of the arrows represents schematically the strength of single synapses. In addition to recurrent interactions, both populations receive external excitatory inputs; (3) colored noise (see Materials and Methods) modeling ongoing unstructured cortical inputs. B, left, power spectral density (PSD) of simulated LFPs as a function of the amplitude of periodic thalamic input A (legend in the inset indicate A values). Middle, Cortical NB frequency peak as a function of the corresponding frequency peak of the thalamic NB. Dashed gray line indicates the identity line. Right, Circular variance of the phase difference between the simulated LFP and the thalamic input within the NB (black line) and within the BB (gray line) as a function of the amplitude of periodic thalamic input A. Asterisks indicate levels of input for which NB circular variance is significantly less than BB (p = 0.002). C, left, PSD of simulated LFPs for increasing levels of the sustained component of the thalamic input S (legend in the inset indicate S values). Right, Circular variance as a function of the amplitude of the sustained thalamic input S for the phase difference between the simulated LFP and the thalamic input within the NB (black line) and the BB (gray line). Circular variance was significantly modulated by the sustained component of the thalamic input S both within the NB (K-W test p ( 0.001) and the BB (K-W test p ( 0.001). Circular variance was likewise modulated within these two bands [two-way ANOVA: p = 0.99 for the interaction between the value of sustained thalamic input (S) and the frequency band (NB or BB)]. D, Modulation of LFP PSD for K = 0 (left) and K = 10 (right) with respect to K = 30 in experimental (red) and simulated (blue) data. The red shaded area represents SEM of experimental data across animals and recordings. In both simulations, sustained thalamic input was set to S = 500 sp./s, while the periodic input was set to A = 40 sp./s for K = 0 (left) and A = 30 sp./s for K = 10 (right). E, Modulation of LFP PSD for K = 50 (left) and K = 90 (right) with respect to K = 30 in experimental (red) and simulated (blue) data. The red shaded area represents SEM of experimental data across animals and recordings. In both simulations, periodic thalamic input was set to A = 0 sp./s, while the sustained input was S = 600 sp./s for K = 50 (left) and S = 700 sp./s for K = 90 (right). increase of circular variance in Fig. 5C, right), showing an opposite trend compared with the previous case in which the thalamic NB was modulated.
This shows that in our model, the periodic component of the thalamic input generates cortical NB oscillations through entrainment of the cortical activity, while the sustained component of the thalamic input leads to the onset of cortically generated g oscillations encompassing both BB and NB.
By varying A(K) while S(K) was fixed (see Table 2), the model reproduced quantitatively (modulation at K = 0 X 2 r ¼ 0:25; modulation at K = 10 X 2 r ¼ 0:28) the experimentally observed NB modulation for low levels of contrast (Fig. 5D). For the same set of network parameters (see Table 1), the model quantitatively (modulation at K = 50 X 2 r ¼ 0:14; modulation at K = 90 X 2 r ¼ 0:23) reproduced BB peak shape and BB power increase from reference to high contrast (Fig. 5E). This was instead achieved by varying S(K) while A(K) was held fixed (see Table 2). This suggests that the different modulation of these two g bands can be accounted for by the two different neural mechanisms proposed by the model: (1) the sustained thalamic input S(K), determining the intensity of the broad band g through the modulation of endogenous cortical g ; (2) the periodic thalamic input of amplitude A(K), determining the intensity of the cortical narrow band through neural entrainment.

Model of complementary contribution of oscillatory and sustained thalamic input to LFP c activity
We then explored the possibility of capturing in the model the whole range of contrast-dependent g modulation by reflecting in the thalamic input the result of the complementary range of contrast sensitivity we found in the cortex for NB and BB. We defined then a complementary range of contrast sensitivity for the oscillatory input A (K) and sustained input S(K). We set the oscillatory input A (K) to zero in the range K ! 30, while in the low contrasts K , 30 range we selected, for each K, the value of A maximizing the overall similarity between simulated and experimental spectral modulations (see Materials and Methods). We obtained a monotonous negative trend (see Table 2) that could be described by a piecewise decreasing linear function ( Fig. 6A; Eq. 10). When coupled with a fixed value of S(K), this type of tuning of the A(K) parameter led to a decrease of the narrow band with simulated contrast (K-W ANOVA p . 0.9 for K ! 30, K-W ANOVA p = 0.08 for K 20, K-W ANOVA p , 0.01 for K 10; Fig. 6C, top) and negligible fluctuations in the broad band instead (K-W ANOVA p . 0.9; Fig. 6C, bottom).
Conversely, we set the sustained input S(K) to have a fixed baseline value S(K) = 500 sp./s up to K = 30, while in the high contrast range K . 30, for each simulated K value, we selected the value of the sustained input S maximizing the overall similarity between the simulated and experimental spectral modulations. We obtained a monotonous positive trend that could be described by a linear fit (Fig. 6B; Eq. 11). As expected, cortical g band increased with S(K) (linear correlation coefficient equal to 0.98, p , 0.001). This finally led to an increase of both NB (K-W ANOVA p = 0.06 for K 70 and K-W ANOVA p ( 0.001 for K 90; Fig. 6D, top) and BB (K-W ANOVA p = 0.09 for K 60 and K-W ANOVA p ( 0.001 for K 90; Fig. 6D, bottom) for simulated contrast levels K increasing in the high contrast range. This small increase in NB in the high contrast range is in agreement with results in Figures  2B,C, left panel, and 5C, and it is a by-product of the increase of activity in whole g range (see Discussion).
We expect A(K) and S(K) to be present at each time for any level of input in real conditions. To compare experimental and simulation results over the whole range of contrasts we then injected in the network, for each value of K, the thalamic input resulting from the superimposition of both A(K) and S(K) values as defined above. We found that the network reproduced the experimental modulations of both the NB and the BB over the whole set of contrasts (X 2 r ¼ 0:07 and X 2 r ¼ 0:12; Fig. 6E, top and bottom, respectively). Similar results were achieved replacing the values of A(K) and S(K) optimized for each value of K with two piecewise linear functions (X 2 r ¼ 0:14 and X 2 r ¼ 0:23; Fig. 6E, top and bottom, respectively). These results show (1) that the dynamics of g bands in mouse V1 can be reproduced by a standard network architecture by a proper design of thalamic inputs; (2) that very good agreement with experimental data are achieved by modeling thalamic input as the sum of an oscillatory component, sensitive to low levels of contrast, and a sustained component, sensitive, instead, to high levels of contrast.
A broader exploration of input parameters revealed that in general NB is sensitive to the increase in the periodic input A(K) and, to a lesser extent, to the increase in the sustained input S(K) (two-way ANOVA showed p ( 0.001 for parameter A, p ( 0.001 for parameter S and p = 0.3 for their interaction), while BB is only sensitive to the sustained input S(K) (two-way ANOVA showed p = 0.33 for A, p ( 0.001 for S and p = 0.32 for their interaction).

Time-dependent thalamic input model accounts for c band temporal evolution
Modeling results shown in Figures 5, 6 reproduced the average response of the network over time. To reproduce the experimental time-frequency features of V1 LFPs we re-defined the simulated thalamic inputs to be also a function of time (see Eq. 13) by taking into account a reversaldriven modulation at contrast reversal. The function was defined in such a way to be zeroed after 200 ms (as observed experimentally in the LFP scalograms; Fig. 7, last column). We then re-determined A(K,t) and S(K,t) to optimally fit experimental data (see Materials and Methods). This allowed capturing the differences in early and late thalamic stimulation leading to the distinction of early and late cortical responses observed experimentally.
By introducing this temporal dynamic in the thalamic input, we were able to reproduce the experimental LFPs at contrast reversal, both the NB time invariance at K = 0 (Fig. 7A, middle, and X 2 r ¼ 0:03 in Fig. 7A, right), the NB power decrease at K = 10 (Fig. 7B, middle, and X 2 r = 0.02 in Fig. 7B, right. X 2 r = 0.01 for K = 20; data not shown) and the BB power increase for high contrast associated with abrupt contrast reversal (K = 50: Fig. 7C, middle, and X 2 r = 0.07 in Fig. 7C, right, K = 90: Fig. 7D, middle, and X 2 r = 0.21 in Fig. 7D, right).

Discussion
We found that, in V1 of awake mice, narrow and BB displayed an opposite modulation with contrast and were sensitive to complementary contrast ranges with complementary temporal dynamics. Such complementarity of the two g bands is not completely described by a narrow band encoding luminance opposed to a broad band encoding contrast as previously thought (Welle and Contreras, 2017): both bands actually encode contrast, but in a complementary way. To shed light on the network dynamics underlying these experimental results, we developed a simple leaky integrate-and-fire spiking network model of mice V1, which reproduced quantitatively the complementary contrast-driven g modulations of the two bands. The model suggests that cortical g bands complementarity originates in the thalamus, as the sustained and periodic components of the thalamic input are already sensitive to two complementary contrast ranges.

Contrast sensitivity: experimental results
This work focused on the mice V1 processing of visual contrast. The power of the NB was found to decrease with contrast in the low contrast range, reach a plateau for medium contrast, and be weakly modulated by high contrast (Fig. 2B). This result differs from the claim of monotonic decrease of NB observed in (Saleem et al., 2017). This discrepancy of the results might be then because of a lack of a dedicated analysis in their work more than to difference in the experimental set-up: Saleem et al. (2017;see their Fig. 2C) hardly shows a significant decrease in NB power for contrasts above 50 (error bars being SEM).
On the other hand, the BB was found to be sensitive to the variation of high levels of visual contrast. While in (Saleem et al., 2017) the role of BB and NB in providing contrast information seemed redundant, we show that their role is instead complementary. Information theory indeed revealed the two g bands to encode visual contrast in two separated ranges (Fig. 3): the low contrast range (K , 30), in which only the NB is modulated, and the high contrast range (K . 30) in which, instead, only the BB is modulated. These results suggest that in rodents NB  . Temporal evolution of contrast-dependent modulation in simulated data. A, left, Time course of thalamic input model parameters for K = 0. Middle, Scalogram mean modulation at K = 0 with respect to K = 30 for simulated data. Right, Time evolution of narrow band modulation for experimental (red) and simulated (blue) data (shading reflects SEM). B, Same as A for K = 10. C, left, Time course of thalamic input model parameters for K = 50. Middle, Scalogram mean modulation at K = 50 with respect to K = 30 for might have a critical role not only in the encoding of luminance, as previously thought, but also in the encoding of contrast, thanks to a complementary sensitivity to the one of BB.
The NB and the BB are not only modulated by complementary contrast ranges but also primarily modulated in two complementary temporal windows. At low contrasts, NB is indeed desynchronized in the early phase after stimulus presentation, but the prominent effect is a significant reduction of power in the late response compared with baseline (Figs. 3C, 4C). At high contrast, the increase in the BB is instead primarily modulated in the first 200 ms after stimulus presentation (Figs. 3E, 4D).
Circuit dynamics leading to contrast sensitivity: predictions from the model Our model provides a complete candidate description of the way different levels of contrast affect LGN activity and of how this leads to NB and BB modulation in V1.
First, the mechanisms of BB g modulation in rodents V1 is the same that in primates V1, i.e., the increase in the intensity of the thalamic input (i.e., the increase of excitatory external drive to the cortical circuitry) because of higher contrast enhances endogenous g oscillations (Henrie and Shapley, 2005;Mazzoni et al., 2008). This makes BB power proportional to the average value of S(K) (Figs. 5, 6). However, the increase in thalamic input intensity in rodents is not linear but is negligible for low contrasts. Hence V1 BB in rodents is modulated only by high contrast.
Second, there are two different mechanisms originating NB, depending on the level of contrast. At low contrasts, the dominant thalamic input is the periodic one (Fig. 6), which entrains V1 NB oscillation (Fig. 5B). Hence, cortical NB power at low contrast decreases with contrast as it reflects the contrast-driven modulation in power of thalamic NB oscillations (Fig. 6, left column). At high contrast, the dominant thalamic input is instead the sustained one (Fig.  6), which generates cortical g oscillations, encompassing both NB and BB (Fig. 5C). Note that we do not rule out the possibility that LGN NB neurons have the capacity of encoding contrast variations in the high contrast range by furtherly desynchronizing. However, the associated decrease in cortical NB power becomes for high contrasts negligible compared with the increase because of local resonances triggered by the sustained thalamic input. Feed-forward-dependent responses are then overridden by recurrent connectivity-dependent components and cortical NB power at high contrast increases with contrast. Consequently, in the model, we divided the contrast range into two domains: a low one modulating only the periodic thalamic component and a high one modulating only the sustained thalamic component (Fig. 5). In other words, cortical NB power at high contrast increases with contrast as it reflects the increased endogenous g activity because of the increased sustained thalamic drive.
These two ranges correspond also to two different operational modes of V1. In the low contrast range V1 passively mirrors the NB component in the LGN input. In the high contrast range the sustained input drives V1 to become an active intrinsic resonator in the BB (Brunel and Hakim, 1999).
We provide hence a prediction and an explanation for the g band activity in V1 and thalamus given any contrast, providing several testable predictions, such as the complementary piecewise linear modulation of intensity and oscillations in the thalamus as a function of contrast.
Our model also reproduces the temporal evolution of spectral response (Fig. 7), under the assumption that both periodic and sustained thalamic inputs can be decomposed into a static component not varying in time and a component specifically associated to the moment of contrast reversal (see Eqs. 12,14,respectively). It is tempting to speculate that the former is driven by spatial contrast and the latter by temporal contrast. The reversal-locked component might indeed stem from a sudden surge of the thalamus activity right after contrast reversal (maybe because of the peak in temporal contrast; Barbieri et al., 2014), which leads to an increase of the baseline input [i.e., S(K) in Eq. 14] and a temporary disruption of local thalamic NB g oscillation [i.e., A(K) in Eq. 12]. This is compatible with the observed increase of multiunit activity in the early window after contrast reversal as a function of visual contrast level (data not shown).
Our experimental data on the processing of visual contrast are overall coherent with what previously described; however, our analysis adds novel and important insights into the dynamics and the functional role of the NB in encoding this feature.

Luminance sensitivity
One of the main limitations of this study is that we did not explicitly tackle the luminance sensitivity of the two bands as we performed experiments with fixed luminance. However, we observed that stimulation onset (i.e., turning on the screen which is associated to an increase of luminance) was associated with an increase in both NB and BB (Fig. 1B). Previous studies (Saleem et al., 2017;Storchi et al., 2017) found that LGN NB was modulated by light intensity, with a higher luminance being associated with a narrow band with larger power and higher peak frequency (Saleem et al., 2017). Our model predicts that thalamic narrow oscillations entrain the cortical NB with the same peak frequency (Fig. 5B, middle). Hence, we would expect the light-induced increase in LGN g activity observed by Saleem et al. (2017) to drive proportional increases in the cortical NB. As our work highlighted how the narrow band is modulated by both luminance and contrast, future experiments should take into account both these visual features to thoroughly describe their interplay.
continued simulated data. Right, Time evolution of broad band modulation for experimental (red) and simulated (blue) data (shading reflects SEM). D, Same as C for K = 90.

Laminar c in rodents
g Activity is known to spread non homogeneously across layers both when it is spontaneous (Welle and Contreras, 2016;Senzai et al., 2019) and when it is stimuli induced (Xing et al., 2012b;Bastos et al., 2014;van Kerkoerle et al., 2014;Saleem et al., 2017). We decided here to focus on Layer IV in both analysis and modeling (see Materials and Methods), as this layer receives direct thalamic inputs and displays early and strong narrow band g oscillations when presented with visual contrast stimuli (Saleem et al., 2017). However, modeling studies suggest that contrast-dependent interlayer coupling might play a role in the shift between NB and BB at high contrasts. In particular, the temporal decorrelation in population activity at high contrasts leads to g peak broadening (Battaglia and Hansel, 2011). Future works will then investigate interlayer properties of propagation of the two bands, analyzing the activity of different recording sites in laminar electrodes and exploiting the knowledge of interlaminar interactions of advanced models of mice V1 (see below, Models of mouse visual cortex).

Models of mouse visual cortex
Recent models of mouse visual cortex (Billeh et al., 2020; but see also Troyer et al., 1998;Arkhipov et al., 2018), succeeded in reproducing firing rate distribution across cortical layers in response to arbitrary visual inputs with a multilayer network composed by multicompartment generalized LIF neurons. In these models LGN inputs are modeled as a series of spatiotemporal filters of the presented image and lack the g band activity observed in LGN (Saleem et al., 2017;McAfee et al., 2018). g modulation because of visual stimuli is indeed present in the model only with a sensitivity similar to the experimentally observed broad band, coherently with our hypothesis that the anticorrelation between V1 narrow band power and contrast might originate from periodicity in the LGN input. Combining our proposed model of LGN inputs with multilayer models as (Billeh et al., 2020) could unravel the specific interlayer dynamics of the narrow band.
Hippocampal models (Keeley et al., 2017) were able to reproduce the shift between slow and fast g oscillations observed in the area (Colgin et al., 2009). These two rhythms (roughly corresponding to low and high BB in our work) are both generated within the hippocampus because of differences in inhibitory synaptic times (Keeley et al., 2017) rather than by entrainment to external inputs.
A recent interesting modeling study reproduced the effects of LGN input modulation on V1 g responses to visual stimuli (Zachariou et al., 2021) with a single-layer network of conductance-based spiking neurons. However, the model aimed at reproducing the saturation of g power at mid-level of visual contrast that has been observed in human and nonhuman primates  and was not found in our data nor similar murine data (Saleem et al., 2017;McAfee et al., 2018).

NB and BB in primates
Encoding of contrast in the V1 g band is known to be present in primates (Henrie and Shapley, 2005). A comparison of the role of narrow and BB between rodents and primates is however hampered by the fact that nomenclature is different in the two cases. In primates' studies, NB and BB are defined as low versus high g ranges:  versus  Hz (Ray and Maunsell, 2010),  versus  Hz (Bartoli et al., 2019), or even as overlapping  versus  Hz (Hermes et al., 2019).
Nevertheless, many studies have observed that in primates a sudden increase in visual contrast leads to a transient broad band activation followed by a sustained narrow band activity (Swettenham et al., 2009;Ray and Maunsell, 2010;Xing et al., 2012a,b;Murty et al., 2018;Shirhatti and Ray, 2018;Peter et al., 2019). Unlike mice, however, the power and the peak frequency of the sustained g narrow band were found to be positively modulated by visual contrasts and invariant to luminance variation (Peter et al., 2019).
The V1 spiking neurons network we adopted here is not significantly different from previously proposed ones proved to reproduce contrast modulations in primates V1 (Mazzoni et al., 2008(Mazzoni et al., , 2010(Mazzoni et al., , 2011, except for the different thalamic inputs. This implies that we would expect to observe NB as defined in this manuscript also in primates if NB oscillations were present in primates' LGN. However, while the sustained narrow g oscillation has been reported in pre cortical structures such as LGN or even the retina for rodents and cats (Neuenschwander and Singer, 1996;Storchi et al., 2017), simultaneous recording from primate V1 and LGN found no narrow band in V1 and no g oscillation in LGN (Bastos et al., 2014).

Perspectives
Mice are becoming the favorite animal model to study the circuit changes involved in several neurologic disorders. This is because of the availability of sensitive imaging techniques and opto-genetic and chemo-genetic approaches identifying that allow the cellular underpinnings of the disease (Götz et al., 2018;Fagiolini et al., 2020). Mice are also becoming a standard model for the visual cortex (Carandini and Churchland, 2013;Haider et al., 2013;Saleem et al., 2017). In this context, the monitoring of visual responses represents a promising biomarker for preclinical and clinical studies on neurodevelopmental disorders (Lupori et al., 2019). Thus, an effective in silico model reproducing cardinal functions of the mouse visual cortex could lay the ground for a better understanding of the pathogenetic mechanisms underlying functional impairment in neurologic diseases. In future works, we aim at investigating disorders involving the visual cortex with an interplay of experimental and modeling studies starting from the results presented here.