Modulation of cortical slow oscillations and complexity across anesthesia levels

The ability of diﬀerent groups of cortical neurons to engage in causal interactions that are at once diﬀerentiated and integrated results in complex dynamic patterns. Complexity is low during periods of unconsciousness (deep sleep, anesthesia, unresponsive wakefulness syndrome) in which the brain tends to generate a stereotypical pattern consisting of alternating active and silent periods of neural activity —slow oscillations — and is high during wakefulness. But how is cortical complexity built up? Is it a continuum? An open question is whether cortical complexity can vary within the same brain state. Here we recorded with 32-channel multielectrode arrays from the cortical surface of the mouse and used both spontaneous dynamics (wave propagation entropy and functional complexity) and a perturbational approach (a variation of the perturbation complexity index) to measure complexity at diﬀerent anesthesia levels. Variations in anesthesia level within the bistable regime of slow oscillations (0.1–1.5 Hz) resulted in a modulation of the slow oscillation frequency. Both perturbational and spontaneous complexity increased with decreasing anesthesia levels, in correlation with the decrease in coherence of the underlying network. Changes in complexity level are related to, but not dependent on, changes in excitability. We conclude that cortical complexity can vary within a single brain state dominated by slow oscillations, building up to the higher complexity associated with consciousness.


Introduction
Brain complexity relies on the balance between segregation and integration within cortico-thalamic networks ( Tononi and Edelman, 1998 ).Since changes in brain complexity have been hypothesized to be a key correlate of the brain's capacity to sustain conscious experience ( Tononi and Edelman, 1998 ), a series of recent studies has measured cortical complexity during different brain states (e.g., unresponsive wakefulness syndrome, anesthesia, sleep, wakefulness).Typically, complexity is assessed by two main approaches: either by analyzing spontaneous brain activity ( Andrillon et al., 2016 ;Ferenets et al., 2007, 2006, Schartner et al., 2015, 2017 ) or by quantifying reactivity patterns to cortical perturbations ( Casali et al., 2013 ;Sarasso et al., 2015 ).Measures of cortical complexity provide a good marker of consciousness in unresponsive patients ( Casarotto et al., 2016 ;Sitt et al., 2014 ) and anesthetized humans ( Sarasso et al., 2015 ) therefore understanding their general properties is highly relevant.
Different methods have been developed for estimating cortical complexity.For example, Zhao et al. (2010) suggested a definition of cortical complexity based on a pairwise synchronization measure, defining a low complexity state as an unstructured situation in which the individual oscillators composing the cortical network are weakly coupled.The quantitative value of complexity is computed by evaluating the Shannon entropy of the distribution of the average phase difference values between all pairs of oscillators, providing a measure that reflects the integration-segregation balance in the spontaneous network activity from the viewpoint of oscillatory dynamics.A different approach comes from graph theory ( Bullmore and Sporns, 2009 ) where structural brain networks are described as graphs made of nodes linked by edges representing synaptic and/or axonal projections.In this framework, network complexity is related, for example, with a measure of efficiency -the path length -defined as the minimum number of edges that must be traversed to go from one node of the network to another.
Another approach to quantifying cortical complexity is to use the Perturbational Complexity Index (PCI) ( Casali et al., 2013 ;Comolatti et al., 2019 ) in which neural activity is exogenously perturbed by means of transcranial magnetic stimulation (TMS).PCI measurements in unconscious humans ( Bodart et al., 2017 ;Casali et al., 2013 ;Casarotto et al., 2016 ) show that cortical responses are local and stereotypical and therefore more easily compressible, as quantified by the Lempel-Ziv compression algorithm ( Lempel and Ziv, 1976 ), whereas during conscious states cortical responses are more heterogeneous and widespread across the cortex, and thus less compressible.
Independently of the specific method employed for its assessment, it has been generally reported that complexity is higher during conscious states than during unconsciousness.Hence complexity was found to be lower during non-rapid eye movement (NREM) sleep, when the brain spontaneously generates slow oscillations (SO) ( Andrillon et al., 2016 ;Massimini et al., 2005 ;Schartner et al., 2017 ) and in unresponsive wakefulness syndrome patients ( Casarotto et al., 2016 ).SO are a neural activity pattern characterized by periods of neuronal firing called Up states, interleaved by rather silent periods called Down states ( Sanchez-Vives et al., 2017 ;Steriade et al., 1993 ).The slow oscillation (0.1-1.5 Hz) in itself represents a regular pattern resulting in a simplification of ongoing activity.Most importantly, the fundamental neuronal dynamics underlying the SO, namely the tendency of cortical neurons to fall into a silent period or Down state (also known as an Off-period in other studies; e.g., Rosanova et al. 2018 ) after an initial activation due to adaptation currents ( Compte et al., 2003 ), has been suggested to be a key mechanism that leads to the breakdown of complex causal interactions in cortical circuits.Macroscale measurements in unresponsive wakefulness syndrome patients ( Rosanova et al., 2018 ), mesoscale intracranial recordings in sleeping humans ( Pigorini et al., 2015 ) and microscale measurements in cortical slices ( D'Andola et al., 2017 ) have shown that the occurrence of stereotypical Down states can break the chain of neuronal interactions, thus preventing the emergence of patterns that are both differentiated and widespread.However, whether complexity can vary within a single brain state such as that of slow oscillations is currently unknown.
In the current study we adopt both an observational and a perturbational approach to understand how cortical complexity in vivo varies across different levels of anesthesia, characterized by variable SO patterns.Our aim was to investigate whether the regime of bistability, the dynamic state of the system that oscillates between active Up states and silent Down states (also called slow oscillations), can express different levels of network complexity.Using the experimental model of anesthetized mouse ( Ruiz-Mejias et al., 2011 ), we modulated the emergent SO by modifying the depth of anesthesia.Neuronal activity was measured by means of local field potential (LFP) recordings performed on the surface of the cerebral cortex with 32-channel multielectrode arrays covering almost entirely one brain hemisphere.We found that cortical complexity can vary even within the regime of SO in a manner related to, but not dependent on, the level of excitability.

Animals
All procedures were approved by the Ethics Committee at the Hospital Clinic of Barcelona and were carried out to the standards laid down in Spanish regulatory laws (BOE 34/11370-421, 2013) and in the European Union directive 2010/63/EU.Recordings were performed in eight adult male C57BL/6J mice bred in-house at the University of Barcelona and kept under 12h light/dark cycle with food and water ad libitum.All experiments were performed in the animal's subjective daytime.

Surgical procedures
Anesthesia was induced with an intraperitoneal injection of ketamine (75 mg/kg) and medetomidine (1.3 mg/kg) and maintained by the inhalation of isoflurane in pure oxygen.Atropine (0.3 mg/kg), methylprednisolone (30 mg/kg) and mannitol (0.5 g/kg) were administered subcutaneously to avoid respiratory secretions and edema.Body temperature was constantly monitored and kept at 37°C using a thermal blanket (RWD Life Science, China).Once under the surgical plane of anesthesia, mice were placed on a stereotaxic frame (SR-6M, Narishige, Japan) and a craniotomy and durotomy were performed over either the right or left hemisphere from − 3.0 mm to + 3.0 mm relative to bregma and + 3.0 mm relative to midline.This broad craniotomy and durotomy were chosen to allow access to a large area of the targeted hemisphere, either left or right.
Three levels of anesthesia were reached that were classified according to the provided isoflurane concentrations: Deep = 1.16% ± 0.08 (s.e.m); Medium = 0.34% ± 0.06; Light = 0.1-0%.The volume delivered was 0.8 L/min and the animals were breathing freely via a tracheotomy.Each anesthesia level was maintained for 20-30 minutes and recordings were consistently obtained in a stable slow oscillatory regime (~10 minutes after the change in concentration) and thus, before microarousals appeared due to lighter anesthesia ( Tort-Colet et al., 2019 ).Furthermore, the absence of reflexes was regularly checked.Recordings of spontaneous and stimulated ( "perturbed ") responses were each of 300-500 s duration.With respect to the extreme of the depth of anesthesia, we had already determined in previous experiments the range of anesthetic that could be used to ensure the maintenance of a bistable pattern of activity, while avoiding respiratory depression.

Electrophysiological recordings
Extracellular LFP activity was recorded by means of 32-channel multielectrode arrays ( Pazzini et al., 2018 ) (550 μm spacing between recording points) covering the entire exposed area of the cortex.The signal was amplified by 100 and high-pass filtered above 0.1 Hz (Multichannel Systems, GmbH), digitized at 5 kHz and fed into a computer via a digitizer interface (CED 1401 and Spike2 software, Cambridge Electronic Design, UK).
In each anesthesia level, the first 300-500 seconds of spontaneous activity were recorded and then the same period of stimulation for im-PCI calculation.In this study we will use the term "imPCI " (see below) which stands for intracranial-multiunit PCI.Fifty stimuli were delivered for the calculation of imPCI per anesthesia level.We ensured that all 32 channels were working and recording good signal during the experiments.

Electrical stimulation
For each level of anesthesia, following the spontaneous activity, the responses to stimulation were recorded.Electrical stimulation was delivered to perturb the cortical network through a bipolar electrode (210 μm spacing between tips, FHC Inc., USA) by means of a constant current isolated stimulator (DS3, Digitimer Ltd., UK) controlled by Spike2 software using a CED Power 1401 interface (Cambridge Electronic Design, UK).The tip of the stimulation electrode was aimed at deep cortical layers and located on the anterior pole of the brain at a distance of 500 μm from the anterior edge of the electrode array ( Fig. 1 A, pink circle).Electric pulses had a duration of 1 ms and an intensity of 500 μA and were delivered every 10 s with a jitter of 0.5-1.5 s to avoid entrainment to the stimulation frequency.These stimulation parameters ensured reproducibility of evoked responses across trials and anesthesia levels.No anatomy to identify the exact cortical location of the stimulation electrode tips was done, given that the depth of insertion was constant across experiments, and the intensity of the stimulus aimed at activating the cortical column and to evoke responses in the different areas converted by the array.

Data analysis and statistics
The detection of the Up and Down states was carried out based on three distinctive characteristics of the LFP signal: slow oscillation deflection (SO), gamma rhythm and shooting from the local network (MUA).These three characteristics are reflected in three time sequences of data: the decimated raw signal, the envelope of the variance of the gamma filtered signal ( Mukovski et al., 2007 ) and the estimation of the MUA signal, from the power of the frequencies between 200 and 1500 Hz, in 5 ms windows ( Reig et al., 2010 ;Ruiz-Mejias et al., 2011 ;Sanchez-Vives et al., 2010 ).The spectrum in this frequency band is considered to be a good estimate of the firing of the neuron population, since it is proportional to the density of the Fourier components at high frequencies ( Mattia and Del Giudice, 2002 ).The MUA signal values were logarithmically scaled to compensate for the high fluctuations in the firing of neurons that are very close to the electrode, thus obtaining the log-MUA signal.For each LFP signal a highly processed data time sequence was obtained resulting from the linear combination of the three timesequences described above, using principal component analysis (PCA) to weigh the contribution of each of them.Focusing on the first principal component, a bimodal distribution of the MUA resulted, such that the two peaks of the distribution corresponded to the samples of the activity network belonging to the Up and Down state, respectively.Thus, a threshold value separating the two modes of the distribution was set between the two peaks, such that samples belonged to the Up or to the Down states depending on their position with respect to this threshold.After Up and Down state detection, mean Up and Down state durations were obtained ( Ruiz-Mejias et al., 2011 ).The frequency of the SO was the inverse of the duration of the entire Up-Down cycle.
The reported firing rate (i.e., log(MUA)) is the mean instantaneous firing rate across time (sum of firing rate samples divided by the sampled interval of time).The spatial distribution of the firing rate over the space covered by the electrodes grid was obtained interpolating the values of instantaneous firing rate of each channel with a thin-plate spline method.
We reconstructed the activation waves from the detected Down-to-Up state transitions as in Capone et al. (2019) .Each wave was reconstructed merging together the transitions occurring in multiple electrodes in a reasonable time interval that was iteratively reduced until each wave contained no more than one transition per channel.We rejected the waves occurring in less than eight channels; otherwise, we replaced the missing values with the result of a nearest-neighbor interpolation using the five nearest points.The reconstructed waves were each associated with a vector of relative time lags computed as the difference between the time of occurrence of the wave in each electrode and the average time of occurrence across all the electrodes taking part in the wave propagation.Using the vectors describing the activation waves we constructed, for each anesthesia level, a Time Lag Matrix (TLM) having a number of rows equal to the number of waves and a number of columns equal to the number of recording channels of the experiment, expressing the whole wave collection.The spatiotemporal course of the wavefronts was then obtained by spatially interpolating the time lags without smoothing, using a thin-plate spline method.A PCA was performed and used to group the waves into clusters sharing similar features with a k-means algorithm.The optimal number of clusters to be used to group the waves in each anesthesia state was chosen using a Silhouette method ( Rousseeuw, 1987 ).For each cluster, we obtained the mean wavefront and its spatiotemporal profile by averaging the waves belonging to the same group.Also, the projection of the TLM on the space defined by the first two principal components was used to compute a measure of complexity based on the Shannon entropy of the probability distribution of the waves in the PC1-PC2 plane.
Similarly, the mean wavefront extracted from each cluster was used to compute the velocity and the direction of propagation of the cluster.This was done by calculating the x-and y-derivative of a smoothed version (thin-plate smoothing spline obtained from Matlab's function tpaps) of the vector of relative time lags associated with each cluster.
We computed functional complexity C as the sum of the differences over the m bins between the distribution of a pair-wise correlation matrix p  (r ij ) and the uniform distribution p ′ μ = 1 m as described by Zamora-López et al. ( 2016) , where the diagonal entries of the correlation matrix r ii are discarded from the calculations to focus only on pair-wise interactions and C m = 2 m−1 m represents a normalization factor.Lower levels of complexity are characterized by narrow distributions, which reflect either total independence (r ij ≈ 0) or global synchrony (r ij ≈ 1); while complexity rises in association with broader distributions.
Perturbational cortical complexity was also assessed based on a modified version of the PCI originally described by Casali et al. (2013) and adapted by D'Andola et al. ( 2017) .In the current version of the algorithm, as stated above, we refer to intracranial-multiunit PCI (imPCI), in which the algorithm was modified from that in D' Andola et al. (2017) , with the difference that the MUA signal was low pass filtered at 10 Hz in order to avoid the impact of potential fluctuations in the in-vivo signal.Activity on each LFP channel was aligned to stimulus onset across trials and the estimated firing rate (see above) was calculated during baseline (last 500 ms before stimulus onset) and response (first 2 s after stimulus onset) conditions.The average firing rate across trials was then calculated, normalized to the baseline condition and low-pass filtered below 10 Hz.Finally, statistical differences in firing rate were assessed by a Bootstrap procedure using one-tail above-threshold significance, which retrieved the response binary matrices showing the spatiotemporal distribution of significant activity.Complexity was finally computed over these matrices by applying the Lempel-Ziv algorithm ( Lempel and Ziv, 1976 ).Lempel-Ziv complexity measures the generation rate of new patterns along a digital sequence, being closely related to properties such as entropy and compression ratio ( Amigó et al., 2004 ;Szczepanski et al., 2003 ;Wolff et al., 2019 ).We measure the compressibility of a binary matrix of spatiotemporal significant responses to stimulation, such that the least complex series are the easiest to compress.Statistical significance was assessed at the population level by calculating a Friedman multiple comparison test followed by post-hoc Wilcoxon signed-rank tests using the Benjamini-Hochberg correction for multiple comparisons.A p value of less than 0.05 was used to define significance.All analyses were implemented in Matlab (MathWorks, Natick, USA) using custom-written scripts.
Identification of spatial clusters of electrodes from functional connectivity (pair-wise Pearson's correlation) matrices was performed using the Louvain method to partition complex networks into communities ( Blondel et al., 2008 ).The implementation by Vincent Traag for use with the iGraph package was employed.Since correlation matrices are weighted data rather than binary adjacency matrices, a uniform modularity was employed as a cost function to determine the goodness or quality of the partitions ( Reichardt and Bornholdt, 2006 ).The uniform modularity favors the formation of clusters with internal correlation larger than the mean.Diagonal entries were ignored, the Louvain algorithm was run 1000 times, and the partition giving rise to maximal modularity was chosen.

Data availability
The raw data will be made available in April 2020 in the repositories of the Human Brain Project, which is currently under embargo, due to further analysis and curation.

Modulation of slow oscillations by anesthesia levels
To address whether cortical complexity can be modulated within a brain state dominated by Up and Down states ( Sanchez-Vives et al., 2017 ), we recorded LFP activity by means of 32-channel surface electrode arrays in the cortex of eight isoflurane anesthetized mice during SO activity ( Fig. 1 A).The areas covered by the array are described by De Bonis et al. (2019) .In all cases we varied the level of anesthesia in order to modulate the emergent rhythmic activity without departing from a state characterized by the presence of Down states or Offperiods ( Rosanova et al., 2018 ).We classified the levels of anesthesia as deep, medium, and light ( Fig. 1 B).Our "light " level of anesthesia is the one before the cortex starts expressing other patterns of activity such as micro-arousals ( Tort-Colet et al., 2019 ).Under these conditions, cortical activity fluctuated between periods of high neuronal firing (or Up states) and rather silent periods (or Down states) at SO frequencies ranging between 0.28 ± 0.05 Hz and 0.76 ± 0.10Hz at the deepest and lightest anesthesia levels respectively (Friedman test, p = 0.0022) ( Fig. 1 C).This significant modulation in the frequency of the SO was accompanied by a modulation in the level of excitability of the underlying cortical circuit, which was apparent when analyzing the average high-pass filtered LFP (200-1500 Hz) of the SO signal, which corresponds to an estimate of the multi-unit activity (MUA).Indeed, MUA -or population firing rates -monotonically increased from deep to light anesthesia (Friedman test, p = 0.0015) reflecting an increase in the overall excitability level and an effective change in the state of the underlying network ( Fig. 1 D).This was most apparent when plotting the distributed average MUA as an interpolated matrix representing the area covered by our multi-electrode array ( Fig. 1 E).Indeed, as anesthesia lightened, MUA experienced a monotonic increase in magnitude differently expressed across different cortical areas, showing a maximum along a diagonal axis apparent in the center of the array.Relatively smaller firing rates were instead always observed in the more rostral and caudal regions.

Diversity of wavefront patterns under slow oscillations across anesthesia levels
Slow waves have been described as behaving like global events and the number of propagating wavefronts has been defined as a parameter related to the underlying state of the cortical network ( Muller and Destexhe, 2012 ).Thus, to further characterize the changes in brain state derived from the modulation of anesthesia, we studied the statedependent modulation of SO wavefront patterns.Wavefronts were detected relying on the MUA extracted simultaneously from each of the 32 electrodes of the multielectrode array, following the approach introduced in Capone et al. (2019) .For each anesthesia level, we generated a Time Lag Matrix ( Fig. 2 A).Then the different spatial patterns of generation/propagation were k-means clustered and their frequency and degree of diversity quantified (see Materials and Methods).We found that at deeper anesthesia levels wavefronts could be ascribed to two main antagonistic propagation patterns ( Fig. 2 B, top): one pattern where wavefronts were consistently originated in the most frontal cortical areas covered by our MEA and propagated backwards towards posterior areas of the cortex; and another pattern where the origin site was located in the most posterior areas of the MEA, so slow waves traveled forward.At lighter anesthesia levels origin sites started to diversify ( Fig. 2 B, bottom).Wavefront initiation sites were apparent not only in the most anterior and posterior areas of the cortex, as in deep anesthesia levels, but also in more medial and lateral areas as well, thus inducing the appearance of a plethora of propagation wavefronts traveling in several directions across the cortex.A further quantification of the number of modes of spatiotemporal pattern propagation revealed that they increased from two at deep anesthesia (i.e., anterior to posterior and medial to lateral) to four at light anesthesia (where two medial > lateral modes added up to those above).Furthermore, the average velocity per anesthesia level significantly increased from deep to light anesthesia levels (v = 18.2 ± 1.7 mm/s (deep), v = 23 ± 3.1 mm/s (mid), 30.4 ± 2.9 mm/s (light); p < 0.05 for deep vs. light and mid vs. light, Wilcoxon rank-sum test).
This increment in the diversity of wavefront patterns with the lightening of anesthesia was more apparent when analyzing the distribu-tion of the different propagation modes on the plane of the first two principal components of the time lag matrices ( Figs 2 C and D).The amount of dispersion on this plane -that is, the wavefront entropy -is an indication of the dissimilarity between wavefronts, which predicts the level of intrinsic entropy present in the network.Indeed, wavefront entropy increased with the lightening of anesthesia in accordance with the increment in SO frequency and the higher diversity of waveform patterns ( Fig. 2 E) (Friedman test, p = 0.0076, Wilcoxon signed-rank test with Benjamini-Hochberg corrections; "Deep " with "Medium " ( p = 0.0078) and "Light " ( p = 0.0078); and "Medium with "Light "; p = 0.25).Higher entropy levels have been associated with higher complexity in dynamic systems, thus this intrinsic increment in the diversity of wavefront patterns suggests an effective increase in the overall complexity of the network with the lightening of anesthesia.

Modulation of cortical complexity with anesthesia levels
The change in brain state induced by the level of anesthesia was also reflected in the modulation of the coherent behavior expressed across the cortical network.For instance, a Pearson's correlation analysis of the MUA signals recorded by each of the 32 channels of our multi-electrode array revealed the presence of a highly widespread level of coherence between channels at deep anesthesia levels, reflecting the presence of a highly functionally integrated circuit ( Fig. 3 A).In contrast, when anesthesia lightened, there was an overall breakdown of coherence, as inlets of locally connected clusters of channels started to emerge from the background.This was indicative of a more differentiated circuit ( Fig. 3 A).
The previously described anesthesia-dependent changes in the overall level of coherence and excitability were suggestive of modulations of the dynamical interdependencies between the elements of the cortical network that could potentially modify the complexity of its hosted collective dynamics.To test this possibility, we computed the functional complexity of the cortical network across the three different anesthesia levels following the method described in Zamora-López et al. ( 2016) .Functional complexity measures the time-averaged spatial dynamics of a network based on the pair-wise correlations among its nodes.Thus, high complexity arises when collective dynamics are between total independence and global synchrony.Fig. 3 B shows the grand average of functional complexity calculated in our population of mice at each of the three different anesthesia levels.Overall, functional complexity was significantly modulated by the level of anesthesia (Friedman ( p = 0.0015); Wilcoxon signed-rank test with Benjamini-Hochberg correction, "Deep " with "Medium " ( p = 0.0078) and "Light " ( p = 0.0078); and "Medium with "Light " ( p = 0.0547)), reaching lowest values at deep anesthesia and monotonically increasing at lighter levels.These results revealed that, in our population, collective dynamics were sufficiently modulated within a SO regime to reflect effective changes in the level of complexity of the underlying network.
To further characterize the correlation structure of the cortical network as lightened anesthesia, we performed a clustering analysis of the correlation matrices.The quality of the partitions was evaluated using a uniform modularity.This quality function estimates how well-defined the clusters are, thus representing an indirect measure of segregation.Although the average correlation was stable from deep to light anesthesia, a reorganization of the pair-wise correlations was observed.As mentioned above, at deep levels of anesthesia correlation matrices were characterized by large clusters of electrodes while during light anesthesia the matrices tended to break into few clusters, thus resulting in an increase of partition quality (20-80% increase) as compared with deep anesthesia ( Fig. 3 C).As expected from the outcome of this clustering method, internal correlations within the clusters were stronger than correlations across clusters both in deep and in light anesthesia levels ( Fig. 3 D).However, the separation between within-and crosscorrelation values became more apparent at light anesthesia ( Fig. 3 D), where the distribution of cross-modular correlations broadened towards weaker values as a consequence of an onset of segregation.In summary, we found that the increment of functional complexity in the correlation matrices of the MUA in light anesthesia was associated with an onset of segregation, bringing a more balanced connectivity to the highly integrated networks present in deep anesthesia.
Functional complexity is constrained by the spatial topology of the network, so it focuses more on the spatial aspect of complexity ( Zamora-López et al., 2016 ).In order to also study the temporal domain of complexity we decided to follow a causal approach as in Casali et al. (2013) andD'Andola et al. (2017) .To this end we used our 32-channel surface electrode arrays to record the evoked LFP response induced by perturbing the cortical network with the application of electrical stimulation through a bipolar electrode positioned in the deep layers of the frontal cortex at the same three different anesthesia levels (see Materials and Methods).Under these conditions, electrical stimulation invariably evoked a response characterized by an increase in MUA, followed by a cessation of firing or Down state.Based on the MUA responses we computed a modified version of the PCI originally described by Casali et al. (2013) and adapted by D'Andola et al. ( 2017) (see Materials and Methods), where only the significant component of responses are considered, retrieving binary matrices that are then compressed using the Lempel-Ziv algorithm ( Lempel and Ziv, 1976 ) in order to obtain a single complexity value at each anesthesia level (deep, mid and light).
Fig. 4 shows a representative example of the type of spatiotemporal response obtained after electrical stimulation at the three studied levels of anesthesia.During deep anesthesia ( Fig. 4 A) electrical stim-ulation elicited a rather stereotypical response, which was widespread across channels (high integration), but showed a similar temporal profile across all channels (low differentiation), thus generating an easy-tocompress response (low complexity).In contrast, during light anesthesia ( Fig. 4 C) electrical stimulation elicited a much richer spatiotemporal response profile, which was more diverse across channels in space and time, as each channel became active at different time points and showed a characteristic type of response.This type of response was harder to compress resulting in higher imPCI values as compared with deep anesthesia.Middle anesthesia levels ( Fig. 4 B) showed intermediate behavior; in this case, responses were still quite brief in time and showed a rather similar temporal profile, although some divergence in the type of responses between channels was more apparent than during deep anesthesia.
This pattern in the modulation of complexity according to the level of anesthesia was consistent at the population level.Fig. 4 D shows the grand average of the perturbational complexity obtained after electrical stimulation in our group of mice at each of the three different anesthesia levels.Similar to functional complexity, perturbational complexity showed a clear tendency to increase with the lowering of anesthesia, reaching significantly higher values at lighter levels of anesthesia as compared to medium and deep levels (Friedman p = 0.0046).Similar results were obtained when the perturbational was estimated with the in-vitro sPCI algorithm ( D'Andola et al, 2017 ; Figure S3; see Methods).Thus, these data corroborated our previous results based on functional complexity; that is, collective dynamics are affected by the level of anesthesia and induce an effective change in the state of the underly-ing cortical network, which ultimately induces a net change in the level of complexity.
To assess whether the extent of the perturbation was a limiting factor to revealing changes in complexity concomitant to modulations in the collective dynamics of the network, we next tested whether the variation in the level of complexity followed a similar pattern when considering the intrinsically generated perturbation of the underlying network induced by the onset of spontaneous Up states ( Fig. 4 E).Under these conditions the magnitude of complexity at deep anesthesia levels was kept low at values reminiscent of those obtained under electrical stimulation.However, when anesthesia was lightened, the increment in complexity previously observed during electrical stimulation was absent.For instance, complexity remained low, regardless of the anesthetic regime, at values equivalent to those of deep anesthesia ( Figs 4 D and 4 E).Thus, it seemed that the intrinsically generated perturbation induced by spontaneous Up states was not strong enough to reflect changes in complexity.This was probably due to the ineffectiveness of Up states to cause a departure from the internal dynamics of the cortical network at each anesthesia level and to cause a chain of causal interactions between different cortical areas.

Discussion
Cortical complexity has been found to vary across brain states, including physiological (e.g.sleep vs. wakefulness), pathological (e.g., unresponsive wakefulness syndrome) or pharmacologically induced states (e.g., anesthesia).In this study we have characterized how cortical complexity can be indeed modulated within the same brain state, namely, the slow oscillatory regime.To this end we recorded both spontaneous and evoked cortical LFP activity by means of 32-channel multielectrode arrays on the surface of the cortex of anesthetized mice at different anesthesia levels.We show that cortical complexity, as measured by functional complexity (based on spontaneous activity) and a modified version of the PCI (imPCI; based on "perturbed " or stimulus-triggered activity), varied concomitantly with the level of anesthesia, such that deeper levels of anesthesia correspond to lower levels of cortical network complexity.This modulation occurred within a dynamically stable regime, consisting on Up and Down states organized into slow oscillatory activity, and was accompanied by concurrent changes in the level of excitability and of spatiotemporal integration and segregation.
Cortical complexity has been measured mainly using neuroimaging ( Demertzi et al., 2019 ;Rubinov and Sporns, 2010 ;Sporns, 2011 ) or EEG recordings on the surface of the scalp in humans ( Andrillon et al., 2016 ;Casali et al., 2013 ;Ferenets et al., 2007Ferenets et al., , 2006 ; ;Sarasso et al., 2015 ;Schartner et al., 2015Schartner et al., , 2017 ; ;Wolff et al., 2019 ).However, some recent studies have undertaken a more local and invasive approach by performing intracranial recordings in rats ( Abásolo et al., 2015 ;Comolatti et al., 2019 ) and humans ( Schartner et al., 2017 ), and relying also on voltagesensitive dyes imaging in both monkeys, rats and cats ( Fekete et al., 2018 ).These studies, like the current one, are more informative of the underlying local dynamics of the cortex since they use recordings from probes in close contact with the brain tissue, and are therefore potentially suitable for differentiating between elements of the network.
Complexity has also been computed by different metrics that follow either a spontaneous or perturbational approach, all achieving similar results ( Casali et al., 2013 ;Casarotto et al., 2016 ;Demertzi et al., 2019 ;Sarasso et al., 2015 ;Schartner et al., 2015Schartner et al., , 2017 ) ).That is, complexity has consistently been reported to be higher during conscious as opposed to unconscious states, reaching lowest values when consciousness fades out under both physiological conditions (e.g.NREM sleep-like activity), pathological conditions (e.g.unresponsive wakefulness syndrome), and under different anesthetics.When following a perturbational approach, most studies use transcranial magnetic stimulation (TMS) as an effective non-invasive way to transiently modify the intrinsic dynamics of the human cortical network ( Casali et al., 2013 ;Casarotto et al., 2016 ).More recently, however, it has been shown that direct electrical stimu-lation also represents a valid approach to alter the dynamics of the cortical network in order to measure complexity ( D'Andola et al., 2017 ).Indeed, direct LFP measurements of cortical network dynamics using cortical slices maintained in vitro have effectively demonstrated modulations in cortical complexity under different activity patterns mimicking those present during slow-wave sleep and awake-like conditions in the presence of noradrenergic and cholinergic agonists ( D'Andola et al., 2017 ) or under different levels of excitability controlled by exogenously applied electric fields ( Barbero-Castillo et al., 2019 ).
Here we have followed an invasive approach using the mouse as a model to study complexity at various levels of anesthesia.This allowed us to record activity in close proximity to the cortical tissue by means of multi-electrode arrays, thus having a more direct measurement of the spatiotemporal dynamics of the underlying network.We have also used a spontaneous and a perturbational approach in order to contrast whether both are comparable when measuring brain dynamics at a local scale on the same population of neurons.
The slow oscillatory rhythm (0.1-1.5 Hz) is a low-complexity state that has been proposed to be the default activity pattern of the cortical network ( Sanchez-Vives et al., 2017 ;Sanchez-Vives and Mattia, 2014 ), since it is generated by such networks in different conditions of anatomical or physiological disconnection (slow wave sleep, cortical islands, unresponsive wakefulness syndrome, cortical slices and slabs).There are infraslow oscillations, or rather fluctuations ( He and Raichle, 2009 ) that occur in a lower frequency band ( < 0.1Hz) and that are associated with organization of higher frequencies and with the emergence of consciousness ( Baria et al., 2011 ;Li et al., 2015 ;Northoff, 2017 ) .Here we focus on the band of frequencies covered by the Up/Down oscillatory regime, which ranges between 0.1 and 1.5 Hz.The low complexity values reached during SO are probably the consequence of the synchronized nature of the nonlinear dynamics at the mesoscale of cortical assemblies giving rise to the alternation of highly active (Up) and rather silent (Down) periods of activity.Indeed, Down states (or Off-periods) have been recently reported to hamper the causality of interactions during NREM sleep in humans ( Pigorini et al., 2015 ).Our results are consistent with the results of Hudetz et al. (2015) , who found that the repertoire of large-scale brain states derived from the spatiotemporal dynamics of intrinsic networks is substantially reduced at an anesthetic dose associated with loss of consciousness.But additionally, we show that within a brain state governed by dynamic cortical Up-Down bistability, activity associated with unconsciousness, the complexity may also vary.We show that complexity follows the level of anesthesia under both spontaneous and electrically evoked conditions, being lower at deeper levels of anesthesia and higher at lower levels of anesthesia, but always within a regime characterized by the alternation between Up and Down states.
Our modulation of the level of anesthesia was broad enough to exert changes in the dynamics of the cortical network without exiting the slow oscillatory regime.The three levels of anesthesia were exclusively classified based on the isoflurane concentration provided (see Methods).The frequency of slow oscillations was the reference value to double check that the anesthesia level was changing accordingly.Isoflurane was chosen because it is a common anesthetic widely used in the literature for the induction and maintenance of general anesthesia in small rodents and other animals, since it induces muscle relaxation and decreases pain sensitivity.Further, being a gas anesthetic, it allowed us to precisely adjust the level of anesthesia in a more reliable and rapid manner than injectable anesthetics.The precise mechanisms of action of isoflurane are not yet fully understood; however, it is known that isoflurane acts as an agonist of GABA and glycine receptors, and an antagonist of type 1 glutamate and acetylcholine receptors ( Flood et al., 1997 ;Zhang et al., 2001 ).Enhancing inhibition and blocking excitation results in decreased excitability.Such decreased excitability leads to a decreased SO frequency ( D'Andola et al., 2018 ;Sancristóbal et al., 2016 ) resulting at some point in lower complexity due to lack of activation of the network.On the contrary, we know that the activation of cholinergic (and noradrenergic) receptors induce a transition towards "awake-like " or asynchronous states, thus increasing, complexity in the isolated cortical network ( D'Andola et al., 2017 ).Therefore, antagonizing acetylcholine receptors drives the system into slow oscillations and into a decreased complexity.
It could be argued that the overall change in complexity we observe is simply a consequence of a change in the level of excitability, associated with an unspecific increase in firing rate.The underlying excitability of the network per se has been recently reported to be necessary, but not sufficient, to induce changes in complexity.D 'Andola et al. (2017) showed this in cortical slices by demonstrating that general and unspecific increases in excitability, such as those induced by the administration of kainic acid, are not enough to effectively induce changes in complexity.Furthermore, no relationship was found between the firing rate and the PCI, and situations of hyperexcitability, such as in seizures, result in large responses that lack information content or complexity ( Massimini et al., 2007 ) or in vitro in kainate-induced activity ( D'Andola et al., 2017 ).Instead, the modulation of complexity requires the alteration of excitability by means of specific combinations of neurotransmitters involved in the generation of wakefulness, such as norepinephrine and acetylcholine.This is to some extent replicated in our experiments as we show that the monotonic increase in excitability between deep, medium, and light levels of anesthesia effectively induces significant changes in perturbational complexity only between medium and light levels, where it is more likely needed for ascending pathways involved in awakening to become active ( Franks, 2008 ).
Fine changes in the complexity might also not be fully captured by the perturbational approach as the stimulation-elicited Down states are expected to be rather stereotypical.Indeed, irrespective of the excitability level of the network, the initial phase of the quiescent periods is strongly driven by after-hyperpolarization currents which in turn prevent neurons from firing ( Sanchez-Vives et al., 2010 ).Additionally, the evoked slow waves of activation (Up states) and quiescence (Down states) likely always originate from the stimulation site, which may result in a similar spatiotemporal pattern.However, here we have shown that spontaneous slow-wave activity displays a wide repertoire of propagation patterns whose degree of diversity (i.e., the entropy) increases with the lightening of anesthesia.Hence, such a fine structure of the spontaneous activity could be hidden by strong perturbations which in turn have been proven to be rather effective in uncovering Down states during awake-like brain states ( Barbero-Castillo et al., 2019 ;Pigorini et al., 2015 ).On the other hand, exogenous perturbation may be more efficient at recruiting long-range connections (see below), while spontaneous wave propagation instead relies on local connections.
Inspecting the richness of spontaneous activation wavefronts not only allows us to characterize the complexity changes within SO activity, it also provides additional insights into the way the cortical network reorganizes as the awake state is approached.Indeed, under deep anesthesia we found from our multi-site MUA recordings that two preferred modes of propagation in opposite directions exist, compatible with recent characterizations relying on voltage-sensitive dye (VSD) imaging ( Greenberg et al., 2018 ).This might appear to be in blatant contradiction with the well-established evidence that slow-waves under NREM sleep and anesthesia originate in frontal regions and propagate backward along the rostro-caudal axis ( Massimini et al., 2004 ;Ruiz-Mejias et al., 2011 ;Stroh et al., 2013 ).In fact, our results contribute to the solution of this apparent contradiction by proving that these two physiological conditions can both be expressed within SO activity.Indeed, as anesthesia fades out, the wavefront complexity increases giving rise to propagation modes along an orthogonal direction with respect to the deep-anesthesia waves, likely originating from the medial frontal cortex as in Nir et al. (2011) .As isoflurane potentiates GABAergic synaptic transmission ( Alkire et al., 2008 ), deep anesthesia corresponds to a decreased excitability of the cortical network which in turn makes spontaneous activity more sensitive to the network structure as observed in in-vitro experiments ( Capone et al., 2019 ).However, as anesthesia fades out, a high level of excitability starts to be recovered and propagation modes can override structural constraints, giving rise to more diverse spatiotemporal patterns and faster activation waves.According to this picture we eventually measured an increase of both the Up-state firing rates and the average propagation velocity which changed from 18.2 ± 1.7 mm/s to 30.4 ± 2.9 mm/s as anesthesia faded out from deep to light levels.
Finally, it is also noteworthy that the two types of approaches we have used to measure cortical complexity, namely spontaneous and perturbational, produced similar results.Previous studies examining cortical complexity have used one approach or the other, but the same neural circuits had not been tested with both approaches before.It is interesting to note, however, that when utilizing a perturbational approach, only externally generated electrical stimulation induced an effective change in the level of complexity.In contrast, using the intrinsic generation of Up states as a trigger did not achieve significant results.It should be noted that the method used to detect the Up states as an "endogenous trigger " was the same used to detect them for the propagation analysis.The Up state that is considered the "trigger " is the leading one in a time window, as is the case for propagation detection.With the "endogenous trigger ", complexity remained low at levels comparable to those during deep anesthesia, while electrical stimulation achieved values that were modulated by the level of anesthesia and were concurrent with the level of excitability and correlation, thus following a similar profile to that observed with functional complexity.
This result highlights an interesting difference between observational and perturbational measures of complexity; while the former relies on the pattern that the system is currently displaying, the latter explore the larger range of potential states that a system can visit.This difference may explain why complexity measures relying on strong perturbations, such as those induced by TMS, show higher sensitivity than measures based on the spontaneous EEG in conscious patients.This may be particularly relevant to discriminate between specific states; for example, between minimally conscious and vegetative state patients, or between different levels of anesthesia.Accordingly, during exogenously generated stimulation we observed a broader engagement of the elements of the network that spread across the entire cortex, as reflected by the modulation observed in the LFP in all our recording channels.During spontaneously generated Up states, however, the recruitment of the network is more localized and propagation relies on local rather than long-range connectivity, being unable to trigger causal interactions between areas, since the local mechanisms of control (potassium currents, inhibition, etc.) are effective to impose a restricted activation under physiological conditions.Moreover, at deep anesthesia levels, both electrical and intrinsic stimulation retrieve similar complexity values, which indicates that at this level the network may act as a local sink (a single global attractor) due to the presence of minimum complexity.Still, other measures of cortical complexity based on spontaneous activity remain to be explored further, such as the complexity of wave propagation.In a recent in-vitro cortical study ( Barbero-Castillo et al., 2019 ) it was reported that the network complexity measured as the richness of wave propagation patterns was a more sensitive measure than the slice PCI while excitability was modulated by means of electric fields ( D'Andola et al., 2018 ).
Overall, our results show that changes in collective dynamics are apparent during SO and that these changes are strong enough to allow modulations at the level of complexity.These changes are accompanied by, but not limited to, concurrent modulations in the level of correlation or excitability of the elements of the network, as the outcome for spontaneous complexity shows.The increase in functional complexity goes hand in hand with a transition from a highly integrated state in deep anesthesia to a more balanced state in light anesthesia following an onset of segregation.These results also point to the cortical network of the mouse as a suitable model for studying the neural mechanisms underlying brain complexity by means of genetic tools or transgenic models.

Fig. 1 .
Fig. 1.Modulation of slow oscillations with anesthesia level.Spontaneous local field potential activity was recorded during SO with a superficial 32-channel multielectrode array placed on the cortical surface (scale 500μm).The pink circle indicates the location of the stimulation electrode.(A).The recordings were carried out at three different levels of anesthesia (B) in eight mice.Horizontal black and grey bars in (B) depict Up and Down events respectively exemplified at the deep anesthesia level.The average frequency of the SO cycle (C) and the mean instantaneous firing rate (D) were characterized for each level of anesthesia.A Friedman test revealed a significant effect of the level of anesthesia on SO frequency (p = 0.0022, post-hoc Wilcoxon signed-rank test with Benjamini-Hochberg correction compared "Deep " with "Medium " (p = 0.0234) and "Light " (p = 0.0078); and "Medium with "Light " (p = 0.0156)) and firing rate (p = 0.0015, post-hoc Wilcoxon signedrank test with Benjamini-Hochberg correction compared "Deep " with "Medium " (p = 0.0078) and "Light " (p = 0.0078); and "Medium with "Light " (p = 0.0391)).Box plots represent the first and third quartiles with the median depicted by the horizontal line within the box and extreme values shown by whiskers.Scale bar in panel A represents 500μm and the asterisk shows the location of the bipolar stimulation electrode.(E) Contour plot of the spatial distribution of the firing rate under deep, medium and light anesthesia level.A = anterior, P = posterior, M = medial, L = lateral; FR = firing rate.* p < 0.05, * * p < 0.01, Wilcoxon signed-rank test with Benjamini-Hochberg correction.

Fig. 2 .
Fig. 2. The diversity of slow wave propagation patterns increases with the lightening of anesthesia.(A) Matrix of time lags of Up state onsets simultaneously detected in the 32 electrodes from an example recording.Matrix rows are different Down-to-Up transitions and color code for the temporal delay with respect to the mean wave time computed as the average of the transition time in all the channels.Principal component analysis is used to sort the transitions, bringing together similar modes of slow wave propagation that were grouped into two clusters.(B) Examples of wavefront propagation patterns identified by the k-means clustering algorithm in one experiment under a deep and a light level of anesthesia.(C) An example of the wavefront patterns collected during an experiment under a medium level of anesthesia projected on the plane of the first two principal components (PC1-PC2) extracted from a matrix of time lags containing the waves reconstructed in all the experimental conditions.The grey scale represents the probability distribution of the waves on the PC1-PC2 plane.(D) Three-dimensional representation of the contour plot of the distribution of waves in the PC1-PC2 for three selected examples under different levels of anesthesia.Note the increasing spreading of the distributions at the lightening of anesthesia.(E) Average Shannon entropy values for the three studied anesthesia levels.The entropy of the propagating wavefronts monotonically increased from "Deep " to "Light " anesthesia levels (Friedman test, p = 0.0076, Wilcoxon signed-rank test with Benjamini-Hochberg corrections; "Deep " with "Medium " ( p = 0.0078) and "Light " ( p = 0.0078); and "Medium with "Light " ( p = 0.25).Box plot conventions as outlined in Fig. 1 .* * p < 0.01.

Fig. 3 .
Fig. 3. Functional complexity varies with the level of anesthesia.(A) Top, average correlation matrices of MUA signals recorded by 32-channel surface multielectrode arrays placed on the cortex of anesthetized mice ( N = 8) at three different anesthesia levels.Bottom, histograms showing the distribution of correlation values belonging to the matrices on top.(B) Average functional complexity values at the three studied levels of anesthesia (deep, medium and light) in the population of eight mice.Friedman p = 0.0015; Wilcoxon: Deep-Mid p = 0.0078, Deep-light p = 0.0078, Mid-Light p = 0.054.(C,D) Re-balancing of integration and segregation as anesthesia lightens.Clustering results from correlation matrices out of the MUA signals.(C) Quality function of the clusters (relative difference) increases from deep to light anesthesia for all eight mice.Thin lines represent the values for individual mice, while the tick line shows the average of all mice.(D) Distributions of cross-correlation values within (solid lines) and across (dashed lines) clusters.Results correspond to merged data from all mice.At deep anesthesia there is a narrow separation between the inter-and the cross-modular distributions, evidence for strong overall integration.The distributions separate at light anesthesia with a broadening of the crossmodular correlations towards weaker values reflecting the onset of segregation between clusters.Difference between the medians of inter-modular and cross modular distributions are Δ= 0.094 for deep and Δ= 0.149 for light.Box plot conventions as outlined in Fig. 1 .* p < 0.05, * * p < 0.01, Wilcoxon signed-rank test.

Fig. 4 .
Fig. 4. Modulation of perturbational complexity (imPCI) with the level of anesthesia.The MUA spatiotemporal response to electrical stimulation (ES, black triangle) in one representative mouse is depicted on each matrix for all the 32 channels (y-axis) during the first 2s after stimulus onset (x-axis) during deep (A), medium (B) and light (C) anesthesia levels.The time profile of the MUA responses for all channels is shown to the left of each matrix.Significant responses ( p < 0.05 Wilcoxon signed-rank test) are indicated by the extent of the continuous horizontal bars at the bottom of each channel.The average magnitude of the pre-stimulus activity is represented by the dashed horizontal bar for comparison.The spatial profile of the MUA responses is shown to the bottom of each matrix on a cartoon representation of the recording multielectrode array at three different time points (t1 to t3).(D,E) Overall magnitude of perturbational complexity values under evoked; Friedman p = 0.0046; Wilcoxon Deep-Mid p = 0.19, Deep-Light p = 0.0078, Mid-Deep p = 0.15 (D) and spontaneous; Friedman p = 0.19 " (E) conditions in our population of mice.Box plots represent the first and third quartiles with the median depicted by the horizontal line within the box and extreme values shown by whiskers.* p < 0.05, * * p < 0.01, Wilcoxon signed-rank test.