Network States Classification based on Local Field Potential Recordings in the Awake Mouse Neocortex

Abstract Recent studies using intracellular recordings in awake behaving mice revealed that cortical network states, defined based on membrane potential features, modulate sensory responses and perceptual outcomes. Single-cell intracellular recordings are difficult and have low yield compared to extracellular recordings of population signals, such as local field potentials (LFPs). However, it is currently unclear how to identify these behaviorally-relevant network states from the LFP. We used simultaneous LFP and intracellular recordings in the somatosensory cortex of awake mice to design a network state classification from the LFP, the Network State Index (NSI). We used the NSI to analyze the relationship between single-cell (intracellular) and population (LFP) signals over different network states of wakefulness. We found that graded levels of population signal faithfully predicted the levels of single-cell depolarization in nonrhythmic regimes whereas, in δ ([2–4 Hz]) oscillatory regimes, the graded levels of rhythmicity in the LFP mapped into a stereotypical oscillatory pattern of membrane potential. Finally, we showed that the variability of network states, beyond the occurrence of slow oscillatory activity, critically shaped the average correlations between single-cell and population signals. Application of the LFP-based NSI to mouse visual cortex data showed that this index increased with pupil size and during locomotion and had a U-shaped dependence on population firing rates. NSI-based characterization provides a ready-to-use tool to understand from LFP recordings how the modulation of local network dynamics shapes the flexibility of sensory processing during behavior.


Introduction
While there is a long history of using extracellular and scalp signals to characterize brain states (Berger, 1929;Creutzfeldt et al., 1966;Steriade et al., 1990Steriade et al., , 1993, recent studies have shown that the membrane potential of individual neurons, which can be measured via intracellular recording techniques, can serve as a particularly sensitive indicator for neural network activity. Intracellular recordings have been instrumental to characterize network states in unprecedented details, defining them as a set of distinctive dynamical features that include oscillatory patterns and depolarization levels. In particular, in awake rodents these studies have correlated variations of behavioral and physiological markers, such as pupil diameter, whisking activity, and locomotion speed with distinct membrane potential dynamics in sensory cortices (Crochet and Petersen, 2006;Poulet and Petersen, 2008;Okun et al., 2010;Constantinople and Bruno, 2011;Bennett et al., 2013;Polack et al., 2013;Reimer et al., 2014;Schneider et al., 2014;McGinley et al., 2015a;Schiemann et al., 2015;Einstein et al., 2017;Neske et al., 2019;Poulet and Crochet, 2019;Nestvogel and McCormick, 2022). A powerful concept for the modulation of network activity in sensory cortices emerging form these studies is a "Umodel" of network states (McGinley et al., 2015b). This model is based on the observation that, in some sensory detection tasks, performance depends on the arousal level following an inverted U-shape (it is maximal at intermediate arousal levels Neske et al., 2019). This model posits the existence of a continuum of dynamical network states across arousal levels which includes three major and well-documented patterns of membrane potential fluctuations. At low arousal levels (small pupil diameter and absence of motor behavior), membrane potential fluctuations largely exhibit stereotypical d -band oscillations. At moderate arousal (intermediate pupil diameter), single cells are hyperpolarized and display low-amplitude membrane potential fluctuations. At high arousal levels (active motor behavior and/or large pupil dilation), the membrane potential exhibits sustained depolarization with high-frequency fluctuations and high firing activity occurs. Network states identified based on these properties of membrane potential dynamics profoundly modulate perceptual abilities and cortical processing of sensory stimuli Neske et al., 2019).
These results shed light on how the internal state of the animal modulates sensory information processing about external stimuli. However, in many experimental settings using awake animals, extracellular measurement of population-level LFP signals is often preferred over single-cell intracellular recordings, because of low yield and high technical difficulty of intracellular experiments. Although LFPs capture subthreshold and integrative phenomena in a local neuronal population (Buzsáki et al., 2012;Panzeri et al., 2015), it is currently unknown how to identify the variety of network states previously described with membrane potentials from the LFP. Furthermore, while it has been reported that single-cell membrane potentials and population-level LFPs are related and their relation varies considerably across cortical states and behavioral conditions (Poulet and Petersen, 2008;Okun et al., 2010;Neske et al., 2019;Nestvogel and McCormick, 2022), it is not yet fully clear how to predict when they are tightly related and when they are not. Precise classification of network state variability from LFPs and its comparison with network state classification performed on membrane potentials could thus be greatly useful to enhance our understanding of how network states change during behavior and what function they may serve. Moreover, such a classification would enhance our comprehension of the relationship between single-cell and population dynamics.
By combining simultaneous intracellular and extracellular recordings in the somatosensory cortex of awake head-fixed mice with novel analytical methods, here we developed an approach to identify, from the LFP signal alone, low-frequency rhythmic states as well as nonrhythmic network states with different levels of depolarization or hyperpolarization. We first characterized the membrane potential dynamics across different cortical states in awake mice. We then identified the LFP properties that better distinguished network states and we used those LFP properties to derive a method for robust classification of network states. Next, we show that our classification method enables to classify well network states and explains the variability of the relationship between LFPs and membrane potentials observed across recordings of neural activity during wakefulness. Finally, we illustrate the generality of the method by applying the NSI classification on recordings from the "Visual coding -Neuropixels" dataset shared by the Allen Institute (Siegle et al., 2021).

Animals
Experimental procedures involving animals have been approved by the IIT Animal Welfare Body and by the Italian Ministry of Health (authorization # 34/2015-PR and125/2012-B), in accordance with the National legislation (D.Lgs. 26/2014) and the European legislation (European Directive 2010/63/EU). Experiments were performed on young-adult (four to six weeks old, either sex) C57BL/6J mice (Charles River). The animals were housed in a 12/12 h light/dark cycle in singularly ventilated cages, with access to food and water ad libitum.

Experimental design
The experimental procedure for simultaneous extracellular and intracellular recordings in awake head-fixed mice have been previously described (Zucca et al., 2017), and the present dataset was used in a previous study (Zerlaut et al., 2019). Briefly, a custom metal plate was fixed on the skull of n = 4 young (22-24 post-natal days) mice two weeks before the experimental sessions. After a 2-to 3-d recovery period, mice were habituated to sit quietly on a fixed platform (without the ability to freely run) for at least 7-10 d (one session per day and gradually increasing session duration). Mice were not involved in any task, and whisking and pupil size were not measured. The day of the experiment, mice were anesthetized with 2.5% isoflurane and a small craniotomy (0.5 x 0.5 mm) was opened over the somatosensory cortex. A 30-min-long recovery period was provided to the animal before starting recordings. Brain surface was kept moist with a HEPES-buffered artificial CSF (aCSF). Local field potential (LFP) recordings were performed by lowering a glass pipette filled with aCSF into the tissue with the tip placed at ;300 mm from pial surface. Simultaneous current-clamp patch-clamp recordings were carried out on superficial layers (100-350 mm) within the same craniotomy. The distance between the tips of the electrodes was in the 200to 250-mm range. All recorded cells had a regular-spiking response to current pulses (data not shown) and were therefore identified as putative pyramidal neurons (Connors and Gutnick, 1990); 3-6 MV borosilicate glass pipettes (Hilgenberg) were filled with an internal solution containing (in mM): 140 K-gluconate, 1 MgCl 2 , 8 NaCl, 2 Na2ATP, 0.5 Na3GTP, 10 HEPES, and 10 Tris-phosphocreatine to pH 7.2 with KOH. Current-clamp recordings were not corrected for liquid junction potential offset. Electrical signals were grounded at the top of the skull (at the location of the craniotomy) and were acquired using a Multiclamp 700B amplifier, filtered at 10 kHz, digitized at 50 kHz with a Digidata 1440 and stored with pClamp 10 (Axon Instruments). From the previously described recordings, we extracted samples with stable membrane potential periods. The early hyperpolarized period (;5 min) following pipette lowering was discarded from the analysis. Next, periods with action potential peaking below 0 mV or displaying a slow (;1 min) drift in the V m trace were discarded from the analysis. This criterion enabled us to perform the analysis on an absolute scale of membrane potential values. Multiunit activity (MUA) was computed by band-pass filtering (0.3-3 kHz) the extracellular signal and taking the absolute value of the resulting signal (Einevoll et al., 2013).

Wavelet transform
Our signal processing pipeline of the LFP was based on the wavelet transform. We implemented a wavelet transform based on the Morlet wavelet, which has the following equation (Torrence and Compo, 1998): where f is the frequency of the wavelet and d 0 the decay parameter of the envelope. We used a value of d 0 =6 throughout the study. The coefficient C d 0 f ð Þ is the normalization coefficient of the wavelet. Note that, to keep a meaningful link with the physical units of the signal, we did not normalize the wavelet with respect to itself, but with respect to a sinusoid (otherwise the wavelet transform with standard normalization of a sinusoid of frequency f and amplitude 1 has a value greater than 1 at the frequency f). The wavelet normalization coefficient was therefore defined, for a wavelet frequency f and an extent d 0 , as: The transform was implemented by a convolution between the complex conjugate of the Morlet wavelet (Eq. 1) and the signal S t ð Þ, i.e.: where hSi t;T f is the signal average in the window centered at t of extent T f . T f is the frequency-dependent window on which the convolution is performed, it was defined as the extent of the wavelet where its amplitude decays by computed the envelope and the phase at time t of a given frequency f in the signal by taking the norm and argument of the complex number W f; t ð Þ.
Computing the pLFP signal From the LFP time series (Fig. 1), we computed a "processed LFP" (pLFP; Mukovski et al., 2007), which corresponds to the temporally-smoothed high-g envelope variations of the LFP fluctuations (see main text). The pLFP was computed as follows. We consider a frequency band spanning [f 0 =w 0 , f 0 Á w 0 ], where f 0 is a root frequency and w 0 the width parameter of the band. We take a set of N = 5 wavelets uniformly spanning this band (i.e., evenly space from f 0 =w 0 to f 0 Á w 0 ). The pLFP signal was computed as the sum over N of the k wavelet envelopes of fre- This time-varying signal is then smoothed over time with a Gaussian filter to yield the final pLFP signal (see Fig. 2a). The parameters f 0 , w 0 , and the time smoothing width T smooth were set to maximize the correlations between the time course of the membrane potential and the pLFP (see Results) and their value is reported in Table 1.

Computing the Network State Index (NSI)
Then, from the pLFP we computed the NSI, as follows. The pLFP was first downsampled by averaging over bins of 1 ms (keeping the 50 kHz sampling of electrophysiological signals is unnecessary given the much slower time scale of state transitions, see Results). We then computed the distribution of the pLFP over the whole recording and we extract the baseline noise level of the pLFP signal p 0 by taking the value of the first (lowest) percentile of the distribution. We interpret this recording-specific p 0 value as the residual high-g envelope in absence of neural activity (see the systematic depolarization from the activity at pLFP p 0 in Fig. 2a) and we therefore considered it as an estimate of the noise level in the pLFP signal.
We computed the time-varying envelope of the [2,4] Hz band d env of the pLFP signal using the above wavelet transform (Eq. 2). We took a set of 20 wavelets uniformly sampling the [2,4] Hz band and, at every 1-ms time point, we extract the maximum envelope from this band. We constructed a weighted estimate X t ð Þ of the lowfrequency content of the pLFP signal by X t ð Þ ¼ p 0 1a Á d env , where a is the threshold parameter for the rhythmic/nonrhythmic classification (see Results). We also computed a slow average of the pLFP fluctuations with a Gaussian smoothing of time constant T mean = 500 ms, yielding the signal Y t ð Þ. Finally, the NSI was defined from the above computed signals by the following equation: where H is the Heaviside step function. From this time-varying signal, we computed what we termed "validated" network states by running through the time axis in steps of T state /2 = 200 ms and identifying those time periods in which the T state = 400 ms window surrounding each time point does not contain variations of the NSI signal larger than the noise level p 0 . When averaging quantities for a given NSI level 7), we considered only the NSI levels including more than five validated episodes to get a meaningful average (thus discarding a few recordings in the population analysis because of a too low number of rhythmic episodes; see figure legends).
The signal processing steps of this procedure are illustrated on Figure 3a, and all parameters of the analysis are summarized in Table 1 for our dataset.

Analysis of the "Visual coding -Neuropixels" dataset
We analyzed data from the publicly available Allen Institute for Brain Science Brain Observatory data (Siegle et al., 2021). The full data collection methodology as well as the download instructions can be found at the link https://allensdk. readthedocs.io/en/latest/visual_coding_neuropixels.html. We restricted the analysis to recordings in primary visual cortex (V1) of wild-type mice in absence of visual stimulation (i.e., focusing on the 20 min in the center of the grey screen presentation episode in the "functional connectivity" recordings). These constraints resulted in a collection of recordings from n = 11 mice. After selecting the probe located in V1, we defined the LFP signal by choosing a single channel within all available channels in the V1 column (;20 channels per Neuropixels probes). We selected the LFP channel with the strongest mean pLFP envelope in the d range over the whole session. By pooling all the available well-isolated single-unit spikes of the V1 channels, we computed a time-varying population rate by binning spikes in bins of 5 ms and smoothing the resulting trace with a Gaussian smoothing of width 30 ms. We also extracted the running speed and pupil area (p ·width·height after their ellipse fit) from the dataset. All modalities (rate, running speed, pupil) were resampled to the LFP sampling using nearest-neighbour interpolation during the analysis.

Statistical analysis
Experimental data were imported in Python using the neo module (Garcia et al., 2014). All signal processing steps (subsampling, convolution, filtering) were implemented in numpy (Harris et al., 2020). Statistical analysis was performed with the scipy.stats module of SciPy (Oliphant, 2007). We analyzed the linear relationship between continuous samples with a Pearson correlation analysis (function scipy.stats.pearsonr), and we reported the two-tailed p-value of the null correlation hypothesis. For the statistics of samples consisting of an averages over a given recording session (n = 14 recordings sessions), we tested the significance using two-tailed t tests (functions scipy.stats.ttest_rel, scipy.stats.ttest_ind, or scipy. stats.ttest_1samp for paired samples, unpaired samples and single samples respectively). The multiple linear regressions of Figure 6 were performed with the OLS (ordinary least squares) function of the statsmodel module. We analyzed the statistical significance of the single or multicomponent models with an F test, and we report the variance adjusted by the numbers of factors.

Results
Simultaneous intracellular and extracellular dynamics in the somatosensory cortex of awake mice: variability and network states of wakefulness We performed simultaneous closely-located (;200-250 mm) recordings of the LFP and the membrane Values used for the S1 dataset . Note that the p 0 and p thre fluct parameters are data-driven quantities, i.e., varying from recording to recording, deriving from the value of p thre 0 (mean 6 SEM over the n = 14 recordings).
Research Article: Methods/New Tools In the middle panel, the time-varying SD s (LFP) evaluated over a 500-ms sliding window (brown line), the d -envelope d env (purple line), and the g envelope g env (green line) of the fluctuations are shown. c, Sorted histograms of the absolute correlation coefficient between LFP and V m over recordings (n = 14, gray). Blue (recording #2) and orange (recording #11) indicate the example recordings shown in a. d, Frequency spectrum of the V m signal obtained with wavelet-based time-frequency analysis (see Materials and Methods). The power-line frequency was blanked (50 6 2 Hz). We highlighted the d (2-4 Hz) and g (30-80 Hz) bands in purple and green, respectively. e, Histogram of V m d envelope across recordings. Time samples classified as "rhythmic" are shown in purple (see main text). f, Histogram of V m depolarization level for nonrhythmic samples. g, Relationship between mean d envelope of the V m and SD in nonrhythmic episodes per recording for all recordings. h, Same as in d for the extracellular LFP. i, Histogram of the d envelope (d env , left) and the g envelope (g env , right) of the LFP over recordings. In the top inset, we show the histogram of the resulting g -to-d ratio. j, Mean depolarization level (shown as mean 6 SEM over time samples at a given g -to-d level) as a function of the g -to-d ratio for a single recording (recording #11, shown in a, b). We highlight how the g -to-d measure classifies the episodes shown in b (see main text). k, Mean depolarization level as a function of the g -to-d ratio for recording #11. k, Mean depolarization level as a function of the g -to-d ratio over time samples. l, Relationship between MUA (see Materials and Methods) and g -to-d ratio over time samples. For panels d-f, h, i, k, l, we show two example recordings (recording #2 in blue and recording #11 in orange) and the population data as mean 6 SEM. over x-axis levels (n = 14, gray line with shaded area). See Extended Data Figure 1-1 for the relationship between cellular properties and V m -LFP correlations across recordings. potential (V m ) of pyramidal cells (see traces from two example recordings in Fig. 1a) in the superficial layers of the barrel cortex (S1). Recordings were performed in awake head-fixed mice habituated to sit quietly on the recording rig (Materials and Methods; see also Poulet and Petersen, 2008). Recordings had a duration of 5.1 6 2.0 min (n = 14 from 4 mice). Before proceeding to use these data to define indices of network states from the LFP, we document some of its basic properties.
First, we found a notable heterogeneity across recordings (compare recording #2 with #11 in Fig. 1a). In particular, the relationship between extracellular population (LFP) and intracellular (V m ) signals was highly variable (Fig 1c, recordings were sorted by their level of absolute correlation, note the large range of observed correlation values) and this variability could not be explained by cellular properties of the recorded cell (such as the membrane resistance and action potential threshold; see Extended Data Fig. 1-1). Second, intracellular and extracellular dynamics had a rich repertoire of activity patterns (illustrated in Fig. 1b). As previously reported under similar conditions (Poulet and Petersen, 2008;McGinley et al., 2015a;Vinck et al., 2015;Chen et al., 2017;Einstein et al., 2017), both LFP and V m traces displayed epochs of rhythmic activity in the d -band (defined as time samples with high [2,4] Hz envelope, see the V m and LFP spectrums in Fig. 1d,h, respectively, peaks were observed at 3.0 6 0.3 Hz for the V m and 3.3 6 0.5 for the LFP, n = 14 recordings). Those rhythmic epochs presented a high synchronization between LFP and V m (correlation between V m and LFP d envelopes across time samples: 0.55 6 0.11, one-sample t test for positive correlation, p = 2e-10, n = 14). Examples epochs of rhythmic activity are shown in traces number 1 and 2 of recording #11 (Fig. 1b). Next, confirming previous observations Neske et al., 2019;Zerlaut et al. 2019;Nestvogel and McCormick, 2022), we found nonrhythmic epochs (here defined as time samples with a V m d envelope lower than 6 mV; see Fig. 1e) at different depolarization levels (see the population histogram on Fig. 1f ranging from ;À80 mV hyperpolarization levels to ;À45 mV depolarization levels). Example epochs of nonrhythmic activity patterns at different depolarization levels (increasing from epoch 3 to 7) are shown in Figure  1b. Taken together, the epochs 1-7 of Figure 1b recapitulate the different states described by the "U-model" of cortical states (McGinley et al., 2015b), that is rhythmic states of d -band activity and nonrhythmic states at various membrane potential depolarization levels (from hyperpolarized to depolarized). Finally, recordings also differed not only in terms of the average strength of d rhythmicity, but also in terms of the distribution of V m levels over time (Fig. 1f). Recordings with lower average d envelope tended to have lower variability of V m levels over time in nonrhythmic states (Fig. 1g, Pearson correlation between mean V m d envelope per recording and V m SD of nonrhythmic time samples, c = 0.83, p = 2e-4). The diversity of recordings thus filled a continuum between two qualitatively different cases of either recordings displaying mostly nonrhythmic states at an almost constant depolarization level (e.g., recording #2 in Fig. 1a, see V m histogram in Fig. 1f showing a low variability of mean V m over time) and recordings exhibiting overall stronger time-averaged d V m envelope and a much wider variation of the V m depolarization levels (e.g., recording #11 in Fig. 1a, see V m histogram on Fig. 1f). These latter cases exhibited a complex dynamics with the alternation of oscillatory d -band activity together with nonrhythmic activity at very different V m depolarization levels (see example epochs of recording #11 in Fig. 1b). In the next sections, we investigate how to quantitatively classify and differentiate these network states based on the extracellular LFP.
Limitations of the existing spectral analysis of LFP for the characterization of network states with different degrees of membrane potential rhythmicity and depolarization Because of the high impact of different strength of V m rhythmicity and V m depolarization levels on behavior and sensory function, we next considered how to determine a quantitative index of networks states with such V m features from the more easily accessible LFP signal. Ideal properties of this index would include: (1) the ability to predict the rhythmicity and depolarization of V m from the LFP; (2) a U-shaped dependence of the index on depolarization levels of membrane potential and of firing of local neural populations to directly map onto the U-model of network states.
We first evaluated whether existing methods based on simple spectral properties, could be used to characterize in this way, using only the LFP, the diversity of network states observed in the awake neocortex. Previous LFPbased characterization of different network states relied on the ratio between d and g power (Cheng-yu et al., 2009;Saleem et al., 2010). We therefore computed the time-varying d [2,4] Hz and g [30,80] Hz envelope of the LFP, and the g -to-d envelope ratio over time samples (see example single recording and population histograms in Fig. 1i). We investigated the ability of the g -to-d ratio to differentiate between epochs of activity that have dynamical features resembling the network states previously documented with V m and described by the U-model of cortical states (McGinley et al., 2015b). To gain intuition, we first considered the example epochs 1-7, which were sorted according to the g -to-d ratio of their LFP (Fig. 1i, see the corresponding time-varying d and g envelopes for those epochs in Fig. 1b). While this ratio could distinguish the strongly rhythmic epochs (1,2) from the high-g and highly depolarized nonrhythmic epoch 7 (Fig. 1i), it could not distinguish well different V m depolarization levels within the nonrhythmic sets of epochs. Epoch 6 had a mean V m depolarization value .15 mV higher than that of epochs 3, 5, but all these three epochs had similar g -to-d ratios (Fig. 1i). Moreover, a nonrhythmic epoch (4) had similar g -to-d ratio to the two rhythmic epochs (1,2). Overall, when quantifying the dependence of mean depolarization level on the g -to-d ratio across all epochs for either the example recording #11 (Fig. 1j) or across all sessions (Fig. 1k), it was apparent that, using the g -to-d LFP ratio, it would be very difficult to distinguish between rhythmic states and nonrhythmic hyperpolarized states, and to distinguish between hyperpolarized and depolarized states within the nonrhythmic range. Next, because states with nonrhythmic and hyperpolarized V m are accompanied by a reduced level of spiking activity Neske et al., 2019;Zerlaut et al., 2019;Nestvogel and McCormick, 2022), we also analyzed the level of population firing by computing the MUA from the extracellular recordings (see Materials and Methods). Similar to what we observed with the average V m depolarization (Fig. 1k), we found that the g -to-d LFP ratio had a very weak predictive power with regard to the population spiking activity (Fig.  1l). Thus, the g -to-d ratio could not be used to identify states of reduced spiking network activity during nonrhythmic activity. Reduced depolarization levels and spiking activity are important as they identify specific network states which are characterized by different properties of sensory information processing during wakefulness McGinley et al., 2015a;Vinck et al., 2015;Neske et al., 2019).
We concluded that the g -to-d LFP ratio poorly differentiated rhythmic and nonrhythmic states and, more importantly, did not enable to evidence different levels of depolarization within the nonrhythmic states observed in cortical dynamics under awake condition. In the next section, we introduce a processing step of the LFP that allows such characterization.
Using the time-varying high-c envelope of the LFP for a richer network state characterization Because the membrane potential V m is the reference signal for cortical state classification (Poulet and Petersen, 2008;Polack et al., 2013;Reimer et al., 2014;McGinley et al., 2015a;Einstein et al., 2017;Arroyo et al., 2018), and because we wanted to capture from the LFP the finer features of V m dynamics (including the variations in rhythmicity and Figure 2. The time-varying high-g envelope of the LFP displays strong correlations with the membrane potential of pyramidal neurons in awake mice S1. a, Example simultaneous recording of the LFP (top) and V m of a layer 2/3 pyramidal cell (bottom). In the middle, we show the time-varying high-g envelope (brown thin line) and its smoothed fluctuations (brown thick line, the pLFP signal). The p 0 value (brown dotted line) corresponds to the first 100th percentile of the pLFP distribution over the whole recording. b, Cross-correlation between V m and the envelope of the LFP wavelet transform in the frequency band [f/w, f·w] (f is a root frequency and w is a width factor). We show the cross-correlation value after averaging over n = 14 recordings (see the individual values per recording in c). Note the optimal band found for f opt = 72.8 Hz and w opt = 1.83 (brown circle). c, The effect of temporal smoothing on the cross-correlation between the LFP and V m signals. Shown for all recordings [individual recordings are color-coded according to their cc(V m ,pLFP) value, we show the correspondence with the recording index of Fig. 1c in the inset]. At T smoothing = 0 ms, one can see the mean value of b and its variability over recordings (black error bar). d, Cross-correlation between LFP and V m as function of the temporal smoothing parameter plotted after normalizing the raw cross-correlation levels of c by their maximum amplitude and subtraction of their level at T smoothing = 0 ms. With this normalization, a peak is clearly visible at T opt = 42.2 ms. depolarization levels posited by the U model) that cannot be captured by the simple d -to-g ratio, we next investigated whether simple mathematical transformations of the LFP displayed temporal fluctuations more qualitatively similar to those of the membrane potential.
The inverted LFP (-LFP) displays high correlation values (cc ;0.5) with the membrane potential in awake animals (Poulet and Petersen, 2008;Arroyo et al., 2018) and could thus potentially provide a basis signal for the characterization of network states. However, the amplitude of the LFP is strongly dependent on the depth of the recording (Sakata and Harris, 2009;Kajikawa and Schroeder, 2011;Lindén et al., 2011;Smith et al., 2012;Herreras et al., 2015) and is subjected to drifts over short (,10 s, Fig. 2a; epochs i and ii) and long (.1 min) time scales. These factors limit the similarity between V m and -LFP, and they were shown in previous work to prevent robust classification of network state during slow wave (,1 Hz) activity (Mukovski et al., 2007).
Guided by the previous findings in anesthetized animals (Mukovski et al., 2007), we hypothesized that the high-frequency content (f . 40 Hz, including the g band activity) of the LFP would provide a good predictor of the depolarization level V m . We therefore applied a wavelet transform to the extracellular LFP and identified the frequency band maximizing the cross-correlation with the simultaneously recorded membrane potential in our dataset (Fig. 2b). This was performed by independently varying a root frequency f and a width factor w, yielding the frequency band [f/w, f·w]. For each frequency band, we divided the band into 20 evenly spaced wavelet frequencies, we computed the mean over frequencies of the wavelet envelope of the LFP (resulting in the time-varying envelope shown in Fig. 2a), and we analyzed the correlation between this transformed LFP trace and the V m trace after averaging over recordings (see the individual values per recording in Fig. 2c sorted by recording index in the inset). We found that the band maximizing this correlation was achieved for f opt = Figure 3. The Network State Index based on the processed LFP (NSI pLFP ). We define a graded measure of network states based on the mean level (for nonrhythmic activity) or the low frequency envelope (for rhythmic activity) of the time-varying pLFP signal. a, Example epochs of activity at different NSI pLFP levels (see bottom plot), the epochs are identical to those of Figure 1b (recording #11). In the top plot, we superimpose the pLFP fluctuations and the V m fluctuations. In the middle plot, we illustrate the signal processing steps leading to the NSI pLFP measure (see main text and Materials and Methods). Two time-varying quantities derived from the pLFP signal are used to classify network states: a weighted estimate of the low frequency content of the pLFP signal X(t) (purple line, shown for a = 2.87) and the pLFP sliding mean Y(t) (black line). A consistency criterion validates a fraction of those as "validated" network states (brown dots). b, Fraction of rhythmic, nonrhythmic and unclassified states as a function of the parameter a (weighting the propensity to classify as rhythmic states). c, V m envelope of the [2,4] Hz band averaged across all identified rhythmic states (mean 6 SEM over the n = 14 recordings) for different values of the parameter a. We fit the decay with an exponential function (dashed red line) and take its decay parameter as the optimal value a opt for the classification. d, Mean depolarization level (shown as mean 6 SEM over episodes at a given NSI pLFP level) as a function of the NSI pLFP measure for a single recording (recording #11). We highlight how the NSI pLFP measure classifies the episodes shown in a (see main text). e, Mean depolarization level as a function of NSI pLFP for the n = 14 recordings of the dataset. We show the mean relation per recording (color-coded dots) and the mean and variability over recordings after Gaussian smoothing of 5 mV width (black curve and gray area, respectively). f, Relationship between MUA (see Materials and Methods) and NSI pLFP over episodes for all recordings (color-code following e). 72.8 Hz and w opt = 1.83, i.e., the [39.7,133.6] Hz band (see Fig. 2b). We also found that a temporal smoothing of the LFP envelope (with an optimal value T opt = 42.2 ms; see Fig. 2d) enhanced its correlation with the V m signal. We refer in the following to such smoothed high-g envelope as the pLFP (in analogy with the terminology by Mukovski et al., 2007). Following previous literature, we interpreted this quantity as an approximation to the time-varying recruitment of synaptic activity from a local region (diameter: ;100-200 mm) surrounding the extracellular electrode (Ray et al., 2008;Katzner et al., 2009;Lindén et al., 2011;Mazzoni et al., 2011;Ray and Maunsell, 2011;Buzsáki et al., 2012;Gaucher et al., 2012;Einevoll et al., 2013).
After this transformation of the LFP, we observed a qualitative match between the previously reported V m signatures of network states (McGinley et al., 2015b) and specific features of the pLFP signal. We illustrate this finding on the recording shown in Figure 3a. We observed rhythmic activity at different envelope levels (example epochs #1, #2) as well as nonrhythmic fluctuations at various mean levels of pLFP signal (example epochs #3-7). This similarity encouraged us to develop a quantitative NSI based on the pLFP.

Designing a NSI from the pLFP
To better discriminate specific network states of wakefulness from the LFP, we thus developed a quantitative index from the pLFP: the pLFP-based NSI (NSI pLFP ). The rationale and procedure to compute the NSI from the pLFP are described in the following text and is sketched graphically with example data in Figure 3a.
To capture the slow fluctuations of network activity over time, we computed the sliding mean Y(t) of the pLFP over a slow time scale (T mean = 500 ms). Because the pLFP signal had a nonzero value at all points, we quantified the baseline of the raw pLFP signal p 0 (set as the lowest 100th percentile of the pLFP distribution, and a measure of the level of baseline noise in the extracellular signal). We then analyzed pLFP fluctuations relative to this baseline level. We classified the pLFP fluctuations at such slow time scale as either rhythmic or nonrhythmic. Rhythmicity was quantified by the time-varying low-frequency envelope of the pLFP fluctuations d env (t) using a wavelet transform. We observed that establishing the rhythmic condition by only thresholding d env (t) would be misleading. Indeed, the d -band envelope was strongly co-modulated by the mean activity level Y(t) even for nonrhythmic epochs (correlation coefficients between Y(t) and d env (t) in the nonrhythmic epochs defined as NSI pLFP . 0: c = 0.39 6 0.16 across the n = 14 recordings, significance of a positive correlation: p = 3.3e-10, one-sample t test). This indicated that the d envelope could reach high values, and cross an arbitrary threshold in absence of strong rhythmicity. This phenomenon is visible on Figure 3a: the envelope in epoch 7 is equivalent to the envelope in epoch 2 without exhibiting the clear rhythmicity in the pLFP or in the V m signal that characterized epoch 2 (middle and top plots, respectively). We therefore introduced a simple model-based criterion for evaluating rhythmicity based on the following reasoning. In a noiseless, purely rhythmic setting defined by pLFP(t) = p 0 1 d env · (1 1 sin(2·p ·f d ·t)), where d env is the envelope of the oscillation, we have Y(t) = p 0 1d env . If Y(t) has an additional nonrhythmic component, we get Y(t) . p 0 1 d env (t) (i.e., the oscillation alone does not account for the mean level of the signal). We adapted this last relation to build our rhythmicity criterion: a higher signal mean Y(t) than the mean expected from the d component implies nonrhythmicity. Given the nonideal nature of the signal and to compensate for the misestimation of rhythmicity in fluctuating regimes, we rescaled the slow oscillation with a parameter a (see the next section for the determination of a), i.e., we introduced the time-varying quantity X(t) = p 0 1 a·d env (t). We compared the estimate of the rhythmic contribution, X(t), with the mean activity, Y(t), to quantify rhythmicity: if X(t) ! Y(t) activity was set as "rhythmic" (because the slow oscillation pattern is able to account for the mean activity level). Activity was defined as "nonrhythmic" otherwise. Finally, we quantified the amount of the pLFP activity in the two regimes. In the rhythmic regime, the amount of network activity was captured by the amplitude of the oscillation, 2·d env (t). In the nonrhythmic regime, the pLFP deviations from baseline Y(t)-p 0 estimated the level of ongoing activity.
The pLFP-derived NSI was defined as the amount of pLFP activity projected on the negative and positive axis for the rhythmic and nonrhythmic regimes respectively. Such a definition resulted in a continuous index where states with stronger d components had negative values while nonrhythmic states with stronger high-g components had higher positive values (Fig. 3a). As highlighted by the colored area (Fig. 3a, bottom plot, purple and kaki colors for rhythmic and nonrhythmic epochs, respectively), the classification into the two states relied on the sign of the difference between the X(t) and Y(t) signals (Fig. 3a, middle plot) followed by a projection on either the negative part of the axis weighted by the oscillation amplitude 2·d env (t) for rhythmic epochs, or the positive part of the axis weighted by the increase from baseline Y(t)-p 0 nonrhythmic epochs (Fig 3a, bottom plot). The described procedure for the computation of the NSI is formalized in Equation 4 in Materials and Methods.
As the time-varying signal NSI(t) might exhibit fluctuations due to noise in both the Y(t) and X(t) quantities (X and Y are derived from the noisy LFP signal and their difference might amplify noise, see for example the signal jumps in epochs 2, 3 in Fig. 3a), we added a consistency criterion to NSI(t) to obtain a robust state classification of individual epochs. We first defined network state "episodes" with a window of T state = 400 ms and an update every T state /2 = 200 ms. The motivation behind the choice of this time scale T state was that it offered a good compromise between two constraints: T state was long enough to get well defined states (e.g., more than half a cycle for oscillations in the [2,4] Hz range) and it was short enough to catch the fast and frequent switches of network states during wakefulness (McGinley et al., 2015b). The consistency criterion for episode classification required that, within a given time window of duration T state , the fluctuations of the NSI signal remained within a fluctuation threshold, equal to the pLFP noise level p 0 (because this noise level provided an estimate of the amount of signal below which a variation is not a robust signal variation). If this stability condition was met, a network state in this time window was labelled as "validated." The "validated" states are highlighted with brown dots over some sample epochs in Figure 3a. If this stability condition was not met, a network state at a well-defined level could not be assessed and the network state was labelled as "unclassified." The consistency criterion prevented state validations in the presence of strong fluctuations in the time-varying NSI signal (epochs 3 and 4 in Fig. 3a, bottom).

Calibration of the rhythmicity threshold in the NSI definition
The parameter a in the pLFP-based NSI weights the propensity to classify the network activity as rhythmic. Its effect is illustrated in Figure 3b. Increasing a increased the proportion of rhythmic states from ;0% at a , 1 to ;100% at a . 6. We used simultaneous V m recordings to optimize a to ensure that the classification of states as rhythmic using the LFP actually finds states that would be defined as rhythmic based on the membrane potential. In Figure 3c, we show, as a function of a, the average across all episodes classified as rhythmic of the [2,4] Hz d envelope of the membrane potential V m . The d envelope of V m of the states classified as rhythmic decreased exponentially when a increased. We thus set a to a value a opt = 2.87 that was equal to the decay constant of the V m envelope as a function of a. This choice ensures we detect a large enough number of rhythmic states with a genuine amount of V m rhythmicity. When classifying rhythmic states using such a value in the pLFP-based NSI algorithm, we indeed obtained that the states classified as rhythmic had V m d envelope systematically larger than the states classified as nonrhythmic (paired t test, p = 1.2e-7, n = 14 recordings). Importantly, such a a value was found to be very close to the value maximizing the fraction of unclassified states (a = 2.95, Fig. 3b, thin grey dashed line), thus corresponding to the most conservative setting to classify network states.

Electrophysiological signatures of NSI pLFP -defined network states
We then analyzed additional electrophysiological features of network regimes defined by the NSI pLFP measure. We show the relationship between the NSI pLFP level and the mean V m depolarization value for a single recording in Figure 3d and across all recordings in Figure 3e. In Figure  3f, we show the relationship between the NSI pLFP level and the MUA across recordings.
States of robust rhythmicity (high d envelope, e.g., epoch 1 in Fig. 3a) showed strongly negative NSI pLFP values and were associated to intermediate depolarization and MUA levels (see population data in Fig. 3e,f). When the rhythmicity was not present (low values of pLFP d envelope, as, e.g., in epochs 3 and 4 in Fig. 3a), both the depolarization and the MUA levels had values close to their minimum (see population data in Fig. 3e,f). States for which the d component did not significantly contribute to the network activity (quantified by the pLFP) were classified as nonrhythmic (i.e., NSI pLFP . 0) and both their depolarization and MUA levels strongly increased with the mean level of network activity as captured by the NSI pLFP (epochs 4-7 in Fig. 3a and population data in Fig. 3e,f). Importantly, the values of mean depolarization and of mean spiking activity in Figure 3e,f had a much lower SE for any given value of the NSI pLFP than the one that was found when considering the dependence of mean depolarization and of mean spiking activity on the g -to-d ratio (Fig. 1k,l), suggesting that the NSI pLFP is a much tighter predictor of membrane potential and of spiking dynamics than the g -to-d ratio of the LFP.
We conclude that the NSI pLFP , has several strengths, especially when compared to previous indices. The NSI pLFP captured key features of V m -defined network states during wakefulness (McGinley et al., 2015b): it enabled extraction of membrane potential activity regimes ranging from d -band activity, to asynchronous regimes at low activity levels and asynchronous regimes at high activity levels (Fig. 3e,f). The NSI pLFP therefore provided a quantitative measure of network states that allowed extracting the U-shape nature of cortical states from LFP recordings previously documented with intracellular recordings (Fig. 3e,f).

Evaluation of the pLFP-based NSI accuracy in estimating membrane potential-based features
The above considerations suggest that it should be possible to use the NSI pLFP , a measure only based in LFPs, to identify reasonably well states that have either rhythmic or nonrhythmic membrane potential properties, and to identify among the states with nonrhythmic membrane potentials, those that have either depolarization or hyperpolarization of membrane potential. In this section, we quantified the accuracy of such state characterization using the pLFP-based NSI.
To this aim, we first computed the NSI on the V m signal. The lower bound of the signal p 0 was translated into the V m0 value by taking the first percentile of the V m distribution (see Figs. 2a,4a). We computed the sliding mean and the time-varying low-frequency envelope with the same parameters as for the pLFP signal. We derived the X(t) and Y(t) and the "V m -defined NSI" (NSI Vm ) according to Equation 4 (see Materials and Methods). This V m -defined NSI has (by construction) negative values for the states with V m rhythmicity, low positive values for states with hyperpolarized nonrhythmic V m , and high positive values for states with depolarized nonrhythmic V m . Thus, comparing the value of the V m -defined and pLFP-defined NSI during the validated network states enables a simple quantification of how good is the pLFP-defined NSI at identifying states of rhythmic, nonrhythmic, depolarized and hyperpolarized membrane potential. Figure 4a shows the time course over different recording sessions of the NSI computed either on the pLFP or on V m . From these plots, it is apparent that the two NSI indices are remarkably well matched over the time epochs of the recordings, with the occasional presence of episodes in which the two indices were mismatched Figure 4. The pLFP-based NSI across multiple recordings: variability and estimated accuracy. We show seven recordings (from i to vii) covering the whole range of observed values of correlations between the intracellular (V m ) and extracellular (pLFP) signals. a, A 60-s sample of the simultaneous LFP (gray) and V m (black) signals. At the bottom, we show the validated network states with their NSI value both for the "pLFP-defined NSI" (NSI pLFP , brown dots) and "V m -defined NSI" (NSI Vm , black dots). b, Histogram of the NSI pLFP over the whole recording length for each recording. The color code per recording (from red to blue, see also Fig. 2c) represents the value of the correlation coefficient between the V m and pLFP signals. c, Scatter plots of the "V m -defined NSI" (NSI Vm ) and the "pLFP-defined NSI" (NSI pLFP ) values for all validated episodes in a given recording (i.e., extending before and after the recording sample shown in a). We highlight the correct area with a green color and the incorrect area with a red color. The large red circles give examples of the different sort of rejections that may happen during cross-validation (see main text). See Extended Data Figure  4-1 for the accuracy estimate with different dataset segmentation.
(e.g., those marked by red circles in Fig .4a). To assess when the pLFP-based NSI did and did not correctly predict the NSI measured on the membrane potential, we implemented the criterion displayed in Figure 4c. We first determined the scaling factor F between the V mbased NSI and the pLFP-based NSI, by performing the linear fit of the data predicting either both rhythmicity or both nonrhythmicity (i.e., on the lower left or the upper right of the plots in Fig. 4c). This linear relationship yielded a prediction for the NSI Vm value from a NSI pLFP value. Then, at a given time-point t i , the prediction was considered as correct if the difference between the predicted and observed value of NSI Vm lied within a tolerance interval defined by two free parameters, p tol and V tol m . We set the value of p tol as the average noise level of the pLFP signal across recordings, i.e., p noise = 2.85 mV. For the tolerance value V tol m , we took V tol m = 2 mV to be well above the noise level of V m recordings (;0.1 mV) and still have a high resolution in the [0-35] mV range of observed depolarization levels (see Fig. 3e). Correct detection thus occurred when the following conditions were met: F Á NSI V m t i ð Þ1V tol m À Á ,NSI pLFP t i ð Þ1p tol and F Á NSI V m t i ð ÞÀ À V tol m Þ.NSI pLFP t i ð Þ À p tol . The strictness of this criterion is illustrated in Figure 4c. States were taken as incorrect when there was a mismatch between the rhythmic versus nonrhythmic classification (as shown for cases ii, iv-vii in Fig. 4, see the nonmatching events highlighted with a red circle). Importantly, state classification was also taken as incorrect when the graded level of the rhythmicity or the nonrhythmicity was not predicted well enough (see Fig. 4, i, iii). For example, in recording #11 (Fig. 4, i), the NSI Vm value did not display a high enough value to be linearly related to NSI pLFP level according to the relationship F.
Using the tolerance criteria defined in the above paragraph, the accuracy of detection of a V m -based NSI value from the pLFP-based NSI value was 79.7 6 10.2% (mean 6 SEM over the n = 14 recordings and all validated episodes). Decreasing the tolerance parameter values to p tol = 1 mV and V tol m = 1 mV (corresponding to extremely strict matching criteria) led to an accuracy of 57.2 6 15.1%, meaning that, in more than half of the cases, the network state could be identified even with such remarkably high precision. We checked that accuracy was not artificially inflated by overfitting by splitting the dataset into different training and test sets. In this control analysis, we found very similar parameters and we did not find any significant differences in the measured accuracies (see Extended Data Fig. 4-1). Proportions of the misclassifications split into the rhythmic (NSI 0) and nonrhythmic (NSI . 0) cases for both the pLFP and V m signals over all episodes in the n = 14 recordings (related to Fig. 4c). Figure 5. Relationship between pLFP-derived population activity (NSI pLFP ) and single-cell depolarization (V m ) in rhythmic and nonrhythmic regimes. a, Relationship, at the single recording level (shown for recording #1), between the NSI pLFP levels and the properties of the V m fluctuations. For rhythmic states (NSI pLFP 0, purple color), we show the relationship between the NSI level and the V m amplitude of the d -band envelope. For nonrhythmic states (NSI pLFP . 0, kaki color), we show the link between the NSI level and the mean V m depolarization level over a T state = 400 ms window. We highlight with gray dots the values of the single episodes visible on Figure 3a (reproduced on the top inset, V m in black and pLFP in brown). We show the linear regressions for the rhythmic and nonrhythmic data (dashed red line). Note that the plain curve does not reach episode 3 because the minimum number of episodes for averaging is not reached at that level (see Materials and Methods). b, Same than a for recording #11. c, Reproducing the analysis of a, b over all n = 14 recordings (see main text). We show the mean relations over recordings (thin gray lines) and the mean (wide curve) and SD (shaded area) across recordings. Evaluated only for the NSI pLFP levels displayed by multiple recordings (i.e., n = 6 recordings for rhythmic activity and n = 13 recordings for nonrhythmic activity). For each recording, we perform a linear regression with respect to the NSI pLFP levels, we compute the mean across recordings hsi and the probability of a deviation from the 0-slope hypothesis p (paired t test). d, Relationship between the a parameter (that sets the proportion of nonrhythmic episodes; see Fig. 3c) and the mean slope over cells (top, hsi cells ) together with the p-value testing the significance of a non-zero slope (bottom), i.e., reproducing the analysis of c for different a values.
We noted that the misclassifications were not homogeneously distributed across different states (see Table 2). A specific set of network state misclassifications represented 65.3% of the misclassifications over the merged episodes across all recordings (i.e., 13.1% of all episodes). In this set, rhythmic activity predicted from the V m (NSI Vm 0) was misclassified as nonrhythmic activity from the pLFP (NSI pLFP . 0). We analyze in the next section the reasons behind this finding.
Graded aspect of network states in terms of pLFP and V m fluctuations: rhythmic versus nonrhythmic states The above analysis revealed a stronger tendency to classify network states as rhythmic from the V m than from the pLFP fluctuations. This difference in classification suggested that the NSI measures revealed previously unexplored asymmetries between the characteristics of fluctuations of pLFP and V m signal across different network states. We therefore analyzed more closely the correspondence between pLFP and V m fluctuations for different levels of pLFP-based NSI and how this impacted our classification results.
For nonrhythmic activity (NSI pLFP . 0), the level of the NSI pLFP was given by the mean pLFP deflection (Y-p 0 ) in the time window T state = 400 ms. We thus compared the positive NSI values to the mean membrane potential depolarization level in the same window T state . For rhythmic activity (NSI pLFP 0), the level of the NSI was proportional to the d envelope of the pLFP. We thus compared the negative NSI values to the d envelope in the V m signal. We show the relationship between those pLFP-defined and V m -defined levels for two recordings in Figure 5a,b (for recording #1 and recording #11, respectively) and for the population data on Figure 5c.
Overall, we found (Fig. 5) that the membrane depolarization exhibited a strong dependency on the pLFP-based NSI level for nonrhythmic activity (NSI pLFP . 0). We showed single episodes of increasing NSI pLFP levels (numbered 3, 4, 5, 6 in Fig. 5a,b) and their respective episode averages at all NSI pLFP . 0 levels for those two recordings ( Fig. 5a,b, kaki curves). The ;15 mV variability in terms of pLFP based NSI levels corresponded to a ;30 mV variability of depolarization level with a clear monotonic relationship for those two sample recordings. This behavior was confirmed at the population level (Fig.  5c). We standardized the analysis across all recordings by binning the NSI pLFP levels from 0 to 30 mV (a range covering all observed values) in bins of 1 mV. We found that all recordings exhibited depolarizations with a steep dependency on the NSI pLFP level 1.5 6 0.7 mV/mV, significantly deviating from the null hypothesis of a zero slope (p = 1.2e-5, n = 13 recordings, paired t test). On the other hand, we found that the value of pLFP-based NSI for rhythmic activity (NSI pLFP 0) had a much lower impact on the membrane depolarization level (see Fig. 5). The ;10 mV variability in terms of pLFP based NSI levels translated into a weakly modulated V m oscillation with a 5 to 10 mV amplitude (Fig. 5a,b, purple curves). We highlight this weak dependency on the selected samples shown in Figure 5a,b. We extended the analysis to all recordings, after standardizing the data by binning the pLFP-based NSI levels from -30 to 0 mV. At the population level, we confirmed that the mean membrane depolarization had a weak dependency on the pLFP-based NSI level (-0.2 6 0.3 mV mV), which was not significantly deviating from the null hypothesis of a zero slope (p = 0.14, over the n = 6 recordings displaying rhythmic activity within multiple NSI pLFP bins, paired t test). It should be noted that the lack of graded V m variations for rhythmic activity was not related to our "rhythmicity threshold" limiting the set of rhythmic samples to a potentially-biased subset. When varying the rhythmicity-factor a up to a = 5 (where the occurrence of rhythmic NSI pLFP 0 states reaches ;80%; see Fig. 3b), the depolarization level still showed a very weak dependency to the NSI level compared to that for nonrhythmic activity (Fig. 5d, top) with similar statistical significance values (see Fig. 5d, bottom).
We concluded that the graded levels of neural activity captured by the NSI pLFP had a strong correlate in terms of membrane potential depolarization for nonrhythmic activity (NSI pLFP . 0), while the various pLFP-based NSI levels of rhythmic activity (NSI pLFP 0) rather corresponded to a stereotypical V m oscillation with ;5 to 10 mV amplitude.
Those observations explained the results of our crossvalidation analysis ( Table 2). The high precision of the classification in the jointly nonrhythmic case ("NSI pLFP . 0 and NSI Vm . 0" in Table 2, only 5.4% of the misclassifications) resulted from the strong relationship between the pLFP and V m signals during nonrhythmic states (Fig. 5ac, kaki curves). On the other hand, the existence of a few episodes showing a high envelope d oscillation in the V m with a low d envelope in the pLFP (cases such as episode 2 in Fig. 5a,b) created an ambiguous situation for the classifier because rhythmicity was hard to establish from the pLFP signal in those episodes. The cases with mixed rhythmic/nonrhythmic predictions indeed represented 90.3% of the misclassifications ("NSI pLFP . 0 and NSI Vm 0" and "NSI pLFP 0 and NSI Vm . 0" in Table 2). In particular, predicting rhythmicity from the V m and nonrhythmicity from the pLFP was the prevalent misclassification case (65.3% for the case "NSI pLFP 0 and NSI Vm . 0"), consistent with the stronger representation of the d pattern in the V m than in the pLFP (NSI pLFP 0 range in Fig.  5a-c). We confirmed that such misclassification cases originated from episodes of low d envelope in the pLFP signal: the mean d env in misclassified cases was significantly lower than in accurately classified cases (1.5 6 0.3 vs 2.7 6 0.5 mV, p = 3.3e-10 paired t test, in the n = 13 recordings displaying both conditions).
Using the NSI pLFP to quantify network state distributions in the somato-sensory cortex of awake head-fixed mice We analyzed the NSI pLFP state distribution in the whole dataset (n = 14 simultaneous V m and LFP recordings in the somato-sensory cortex of awake head-fixed awake mice). Figure 4 shows eight example recordings representative of the dataset, spanning the full range of correlations between the pLFP and V m fluctuations (see Fig. 2c). We show a 60-s sample of the intracellular and extracellular recordings together with the time-varying NSI pLFP (Fig. 4a) and the histogram of the NSI pLFP levels across the whole recording (Fig. 4b).
The mean absolute NSI pLFP value over rhythmic periods (i.e., NSI pLFP 0) was 5.65 6 1.02 mV (mean 6 SEM over n = 14 recordings) and 6.676 1.34 mV for nonrhythmic periods (i.e., NSI pLFP . 0), yielding a significant increase from rhythmic to nonrhythmic states (p = 3.0e-5, paired t test). The increased amplitude in terms of pLFP signal (i.e., high-frequency content of the LFP) suggested that, as an average, synaptic activity was stronger in the nonrhythmic periods that in the rhythmic ones (see Discussion). This increased range of pLFP level was also true for the maximum level displayed by single recordings with 9.52 6 2.25 versus 12.06 6 3.4 mV (p = 9.8e-3, paired t test) for rhythmic and nonrhythmic activity, respectively. During nonrhythmic periods, we found a mild but significant increase of the variability (SD, from 1.29 6 0.32 to 1.56 6 0.71 mV, p = 8.9e-2, paired t test) and skewness of the pLFP distributions (from 0.35 6 0.33 to 0.7 6 0.41, p = 8.2e-2, paired t test).

Network state variability within individual recordings predicts the average correlation between population and single-cell signals
We next analyzed whether the NSI pLFP state distributions could explain the variability of the correlation cc(V m , pLFP) between the time-varying population signal pLFP(t) and the single-cell signal V m (t) in our recordings (see Fig.  4a,b). We reduced the NSI pLFP distribution per recording to a few components (detailed below) and we used univariate and multivariate linear regressions to analyse how the variability of those components across recordings explained the variability of the observed correlation values.
From top to bottom in Fig. 4b (i.e., from high to low correlation recordings), the distribution of NSI pLFP levels across recordings was quantitatively different. At the top [high correlations, cc(V m , pLFP) . 0.4, panels i-v], distributions were bimodal with two high peaks both in the rhythmic (NSI pLFP 0) and nonrhythmic (NSI pLFP . 0) areas. At the bottom [low correlations, cc(V m , pLFP) , 0.4, panels vi, vii], distributions were dominated by nonrhythmic episodes with rather narrow range of NSI pLFP values within the (NSI pLFP . 0) domain. To quantify these features within the single recording distributions, we decomposed each distribution into the following quantities: (1) the mean NSI pLFP . 0 over the whole recording m NSI ; (2) the variability (SD) of the full NSI distribution s NSI ; (3) the fraction of rhythmic episodes F NSI 0 ; (4) the mean in NSI pLFP values restricted to nonrhythmic episodes m NSI.0 ; (5) the mean in NSI pLFP values restricted to rhythmic episodes m NSI 0 ; (6) the variability in NSI pLFP values restricted to nonrhythmic episodes s NSI.0 ; (7) the variability in NSI pLFP values restricted to rhythmic episodes s NSI 0 .
We show the results of a linear regression analysis for each factor in Figure 6a (ordered by explained variance) and their relationship with the correlation value in Figure  6b. The factor explaining the highest percentage of the variability (55.27% of the full variability) was the SD of the NSI distribution s NSI . Strikingly, the fraction of (synchronous) rhythmic states was only the third most important factor with an explained variance of 33.89%. A more important factor was found to be the variability of network states within nonrhythmic states s NSI.0 with an explained variance of 41.5%. Recording #14 in Figure 4ii provided an example of these observations. It exhibited strong variability in terms of nonrhythmic states and relatively low occurrence of rhythmic activity (e.g., lower than recording #1). However, it still displayed high correlation coefficient [with cc(V m , pLFP) = 0.66]. Other individual factors had weak statistical significance (p . 0.04; see Fig.6a). However, using a multiple linear regression including all factors, we found that different components of the NSI Figure 6. Features in the distribution of the pLFP-based NSI explain the diversity over recordings of the correlation value between the population (pLFP) and single-cell (V m ) signals. a, We perform a linear-regression between the correlation coefficient cc(V m ,pLFP) and the following quantities characterizing the NSI pLFP distribution: (1) the variability s NSI (SD) of the NSI pLFP values across the whole recording, (2) the variability in NSI pLFP values restricted to nonrhythmic episodes s NSI.0 , (3) the fraction of rhythmic episodes F NSI 0 , (4) the mean NSI pLFP over the whole recording m NSI , (5) the variability in NSI pLFP values restricted to rhythmic episodes s NSI 0 , (6) the mean in NSI pLFP values restricted to nonrhythmic episodes m NSI.0 , (7) the mean in NSI values restricted to rhythmic episodes m NSI 0 . We show the explained variance for all quantities on the x-axis and the statistical significance of those linear models (p-values). We also performed multiple linear regressions with all those factors (dark gray bar) and we show the best three component model (light gray bar: s NSI.0, F NSI 0, s NSI 0 ). We report the variance corrected by the number of linear factors (adjusted R 2 ). b, Scatter plot between the value of a given factor and the V m -pLFP correlation coefficient across individual recordings. Shown for the first five factors of the NSI pLFP distribution with the highest percentage of explained variance (see variance explained and p-values in a). The color code of individual recordings matches that of Figure 2c. distributions had complementary contributions in shaping the average correlation per recording. The full linear model indeed yielded an explained variance of 69.17% (after correction by the number of linear factors, f test p = 3.1e-2). Reducing the dimensionality of the linear model (up to three components), we found that combining the significant subcomponents of the NSI variability (s NSI.0 and F NSI 0 ) with the variability within rhythmic states s NSI 0 produced a statistically-significant model (f test, p = 6.6e-3) predicting 59.84% of the variability (corrected by the number of factors).
We concluded that, in the present dataset, the average correlation between a single-cell signal (V m ) and the population signal (pLFP) could be largely explained (up to ;70%) by features of the NSI pLFP distributions. State variability among all NSI pLFP -defined states (s NSI ) was a critical factor in determining such a variability (Fig. 6a, blue). We decomposed this state variability and found that the two most prominent factors are the variability within nonrhythmic states (s NSI.0 ; Fig 6a, orange) followed by the occurrence of synchronous rhythmic activity (F NSI 0 ; Fig. 6a, green).

Computation of NSI pLFP indices in mouse visual cortex
Finally, we demonstrate the general applicability of the method on a publicly available dataset from the Allen Institute for Brain Science (Siegle et al., 2021). We analyzed neural activity recorded with Neuropixels probes (Jun et al., 2017) in the V1 of awake mice on a steering wheel during spontaneous activity (grey screen presentation). We extracted a V1-located LFP signal from the probes (Materials and Methods) and we computed the time-varying population firing rate summing spikes from all well-isolated single units (Materials and Methods). We first re-performed all calibration steps of the method on the new dataset. We found that, with respect to our V m and LFP S1 dataset analyzed in previous figures, the d oscillation in V1 was clearly shifted up in frequency (Fig.  7a, peak at 5.7 6 0.6 Hz, n = 11 sessions). We therefore shifted the d -band to [4,8] Hz, and we decreased the temporal smoothing (T smoothing = 30 ms) to avoid smoothing out the faster temporal dynamics of the d oscillation. We next determined the optimal rhythmicity parameter a by analyzing the amount of d envelope in rhythmic states when varying a (Fig. 7b, analysis as in Fig. 3c). We found an optimal value of a=3.07, close to the a= 2.85 value of the S1 dataset with V m recordings. Optimal parameters for the V1 dataset are summarized in Table 3.
Similarly to our S1 dataset, we observed that network activity in V1 during wakefulness was characterized by a strong variability of network activity patterns, including d oscillation and nonrhythmic episodes at various spiking activity levels (see raw data in Fig. 7f,g and histograms in Fig. 7c). The activity was overall dominated by nonrhythmic episodes with a percentage of 86.7 6 3.3% of all validated episodes (see Fig. 7c), i.e., activity was ;5% less rhythmic than in our dataset. State distributions were less variable across recordings than in our S1 dataset (e.g., rhythmic activity fraction only varies from a minimum of 10.2% to a maximum of 22.5%), likely due to the increased sampling durations (20 min here vs ;5 min in our dataset). We used this dataset to analyze the relationship between NSI pLFP values, behavioral state and population firing rates. If the NSI pLFP index well captures the U-shaped relationship between network states and behavioral states, it should satisfy three properties: (1) high values during running; (2) increasing values with pupil size; (3) a U-shape relationship All other parameters are identical to Table 1. The p 0 parameter is a recordingspecific parameter (reported as mean 6 SEM over the n = 11 sessions) resulting from the setting p thre 0 =1%. As in Table 1, the pLFP threshold for state validation was set to p thre fluct =p 0 for each recording.
continued Figure 7. Application of the NSI analysis to the "Visual coding -Neuropixels" dataset of the Allen Institute (Siegle et al., 2021). a, pLFP envelope in the d range for the Allen dataset in V1 and in our S1 dataset for comparison. b, Mean d envelope in the spiking rate (Materials and Methods) across all pLFP-defined rhythmic episodes as a function of the a parameter. c, Network state variability in the Allen dataset captured by the pLFP-defined NSI. Thin lines represent single session data. d, Histogram of pLFP-defined NSI values depending on running conditions. The "running" condition corresponds to episodes when the absolute speed is .5cm/s and "still" when the absolute speed is ,5 cm/s. e, Relationship between pLFP-defined NSI levels and average pupil area at those levels. f, Example episodes (from left to right) in session 774875821. From top to bottom, The single unit spikes over time, the timevarying rate computed from those spikes (thick transparent line: sliding mean), the raw LFP in the selected channel (see Materials and Methods), the pLFP signal (thick transparent line: sliding mean), and the pLFP-defined NSI with the validated episodes (validated episodes as dots). g, Same as in f for session 768515987. h, Relationship, at the single session level, between the NSI pLFP levels and the rate fluctuations (shown as mean 6 SEM over episodes). For rhythmic states (NSI pLFP 0, purple color), we show the relationship between the NSI level and the rate amplitude of the d -band. For nonrhythmic states (NSI pLFP . 0, kaki color), we show the link between the NSI level and the mean rate level over a T state = 400 ms window. We highlight with black dots the values of the single episodes shown in g. We show the linear regressions for the rhythmic and nonrhythmic data (dashed red line). i, Same as in h for session 774875821. j, Reproducing the analysis of h, i over all n = 11 sessions. Evaluated only for the NSI pLFP levels displayed by multiple recordings (i.e., n = 9 recordings for rhythmic activity and n = 11 recordings for nonrhythmic activity; Materials and Methods). For each recording, we performed a linear regression with respect to the NSI pLFP levels, we computed the mean across recordings hsi and the probability of a deviation from the 0-slope hypothesis p (paired t test). In a,b,c,d,e,j, data are shown as mean 6 SEM (thick line with shaded area) over the sessions of the dataset.
with firing rates. We examined next how well the NSI pLFP index matches these expectations. We first computed the network state distribution in "running" (mean absolute running speed during the episode larger than 5cm/s) and "still" conditions (Fig. 7d). As predicted based on previous studies (Niell and Stryker, 2010;Ayaz et al., 2013;Polack et al., 2013;Vinck et al., 2015), we found a significant shift of the mean network state toward the higher values of NSI pLFP (mean NSI pLFP in running 31.3 6 9.6 vs 18.9 6 4.7 mV in still conditions, p = 2e-4, paired t test, n = 11 sessions). Next, we found a strong and monotonic relationship between NSI and pupil area (Fig. 7e, Pearson correlation, c = 0.85, p = 8e-6) as reported in previous studies McGinley et al., 2015a), suggesting that the NSI pLFP is a good index of the states expected by the U-model (McGinley et al., 2015b). Finally, we investigated the dependence of population firing rates on the NSI pLFP index. We found that the relationship between pLFP and spiking rates significantly deviated from the null hypothesis of 0-slope both in the rhythmic and nonrhythmic regimes (one sample t test with slope values, p , 0.05 in both cases, see values in Fig. 7j), with a positive slope for positive NSI pLFP values and a negative slope for negative NSI pLFP values. Thus, the NSI pLFP had a U shape relationship with the population firing rate.

Discussion
In this study, we developed a method to extract from LFP recordings in the awake mouse cortex network states information that previously could be obtained only with intracellular V m recordings (Poulet and Petersen, 2008;Polack et al., 2013;Reimer et al., 2014;McGinley et al., 2015a;Einstein et al., 2017;Arroyo et al., 2018;Nestvogel and McCormick, 2022). Prolonged membrane potential recordings are difficult to achieve, and most of the times require head-fixation (but see Lee et al., 2006). This strongly limits our ability to describe the complexity of network states and their relationship with behavior. Achieving precise state classification based on LFP recordings, which are technically easier to perform and can be performed in freely moving animals, will greatly increase our ability to understand the cellular and network mechanisms underlying cortical processing during behavior. Previous attempts to classify network states from the LFP have been limited to the d -band activity Chen et al., 2017;Pala and Petersen, 2018). When applied to data gathered in awake animals, the d -to-g state classification used in anesthetized preparations (Cheng-yu et al., 2009;Saleem et al., 2010) was shown to only separate between the two extremes of the spectrum of cortical states: synchronized d -band activity and desynchronized activity at high-g power. In contrast, the new classification method developed in this study, NSI, captured the large spectrum of network states in the awake neocortex. Importantly, it provided quantitative measurement of the "U-model" of cortical states, which was previously developed only based on intracellular membrane potential recordings (McGinley et al., 2015b).
Prior to the NSI classification, we found it essential to apply a preprocessing step to the extracellular LFP. We computed the time-varying envelope after a wavelet transformation in the high-g band (yielding the pLFP signal). In a previous study in cat neocortex (Mukovski et al., 2007), the authors identified active and silent states (up and down, respectively) under anesthesia from the pLFP signal (with a slightly different frequency band of the LFP: 20-100 Hz, instead of the data-driven [39.7,133.6] Hz band used here). We showed that such a signal processing step enabled capturing various states of wakefulness, ranging from d (;3 Hz) oscillatory activity to desynchronized states at various levels of spiking activity.
We then used the NSI to analyze how the distribution of network states within a recording period shaped the average correlations between the single-cell (intracellular) and population (extracellular) signals. Such diverse correlation levels seemed not to be explained by single-cell features (Extended Data Fig. 1-1) and are unlikely to result from a variable distance between electrodes across recordings (Arroyo et al., 2018), given the close and narrow variability setting of our experiments (200 to 250 mm distance between electrodes; see Materials and Methods). In anesthetized preparations, a consolidated view suggests that low-frequency activity is the main source for neural synchrony in cortical networks and thus for cell-to-population correlation (Steriade et al., 1993). Also in our dataset recorded during wakefulness, the fraction of slow (d -band) oscillatory episodes was a factor significantly contributing to the level of correlation between single-cell and population signal (Fig. 6). However, and unexpectedly, we found that this was a weaker factor than the variability in the set of nonrhythmic states NSI values of cortical activity (Fig.  6). This observation was explained by the fact that synchronized d activity represented only a modest fraction of network activity during wakefulness (18.9% in the S1 dataset, 13.3% in the Allen V1 dataset) and by the fact that the diverse nonrhythmic states corresponded to strongly differing levels of both the V m depolarization and the high-g activity McGinley et al., 2015a;Zerlaut et al., 2019;Figs. 2-4). Importantly, the diversity of network state distributions across recordings largely contributed to the high variability in the measured cell-to-population correlations (70% of this variability could be explained by recording-specific features of the NSI pLFP distribution; Fig. 6). This observation further highlights the necessity of network-state monitoring in the interpretation of experimental results in the awake cortex (McGinley et al., 2015b), and it suggests the possible importance of NSI indices to understand and characterize how the degree of coupling between single-cell and population-level activity varies across different behavioral states.
Our results provide possible insights on the circuit dynamics during different cortical states observed during wakefulness. First, nonrhythmic episodes had population activity levels (pLFP) varying over a wide range and seemed to reach up to levels never observed in rhythmic episodes (Fig 3f). Second, in such nonrhythmic states, we observed a tight relationship between the population