Analysis Pipeline for Extracting Features of Cortical Slow Oscillations

Cortical slow oscillations (≲1 Hz) are an emergent property of the cortical network that integrate connectivity and physiological features. This rhythm, highly revealing of the characteristics of the underlying dynamics, is a hallmark of low complexity brain states like sleep, and represents a default activity pattern. Here, we present a methodological approach for quantifying the spatial and temporal properties of this emergent activity. We improved and enriched a robust analysis procedure that has already been successfully applied to both in vitro and in vivo data acquisitions. We tested the new tools of the methodology by analyzing the electrocorticography (ECoG) traces recorded from a custom 32-channel multi-electrode array in wild-type isoflurane-anesthetized mice. The enhanced analysis pipeline, named SWAP (Slow Wave Analysis Pipeline), detects Up and Down states, enables the characterization of the spatial dependency of their statistical properties, and supports the comparison of different subjects. The SWAP is implemented in a data-independent way, allowing its application to other data sets (acquired from different subjects, or with different recording tools), as well as to the outcome of numerical simulations. By using the SWAP, we report statistically significant differences in the observed slow oscillations (SO) across cortical areas and cortical sites. Computing cortical maps by interpolating the features of SO acquired at the electrode positions, we give evidence of gradients at the global scale along an oblique axis directed from fronto-lateral toward occipito-medial regions, further highlighting some heterogeneity within cortical areas. The results obtained using the SWAP will be essential for producing data-driven brain simulations. A spatial characterization of slow oscillations will also trigger a discussion on the role of, and the interplay between, the different regions in the cortex, improving our understanding of the mechanisms of generation and propagation of delta rhythms and, more generally, of cortical properties.


Introduction
In a brain manifesting slow-wave activity (SWA), expressed in the cerebral cortex under NREM sleep and deep anesthesia [1], the spiking activity -both single and multi-unit activity (SUA and MUA respectively) -appears as a regular sequence of Up (high-rate) and Down (almost quiescent) states. In the past 10 years, an accurate procedure for the analysis of electro-neurophysiological data has been developed, aimed at extracting the MUA from raw recordings, identifying alternating Up and Down states associated with SWA and investigating the features of such rhythmic spatio-temporal patterns of activity propagating along the cortex. The analysis pipeline, implemented in MATLAB 1  properly refined over time, and applied to several experimental data sets, acquired with different setups both in vitro and in vivo from rodents, ferrets, and monkeys [2][3][4][5]. The software procedure has been extensively revised and improved, including some new features implemented in Python 2 , to be used with a new set of data collected in vivo using a 32-electrode array from the cerebral cortex of 11 deeply isoflurane-anesthetized wild-type mice. The data used in this study was obtained in accordance with the Spanish regulatory laws (BOE-A-2013-6271) which comply with the European Union guidelines on protection of vertebrates used for experimentation (Directive 2010/63/EU of the European Parliament and the Council of 22 September 2010), and the protocol was approved by the Animal Ethics Committee of the University of Barcelona.
The multi-electrode array (MEA) employed for acquiring the data is described in Figure 1 and covers several cortical areas ranging from sensory (V1, S1), motor (M1) and association (PtA) cortices [6]. The unfiltered field potential (UFP) is sampled from each electrode at a frequency of 5 kHz and each acquisition session lasts at least 300 s (see Table 1), thus ensuring a fine inspection of the signal in both space and time. The 11 recording sessions are each from a different animal, i.e. 11 independent experiments that collectively represent an extensive data sample covering a wide range of biological and unavoidable physiological variability. lateral ←→ medial Figure 1. Representation of the multi-electrode array (MEA) used for the data acquisition [6]. On the top left, a scheme that indicates the scale of the experiment; the reported dimensions have been adopted to define the reference frame in the SWAP. On the bottom left, the color legend that identifies the cortical areas on the array surface. On the right, the grid of electrodes, with the numbering introduced in the SWAP. In the center, an illustration of the MEA positioned on the mouse cortex; the inspected surface is of the order of 10 mm 2 .
The accurate time-and-space sampling together with the richness of the experimental data have driven the development of new analytical tools that aim to characterize the differences between cortical areas when expressing SWA. In addition, given the unavoidable and physiological variability of the data set, particular attention has been also given to the best strategies to adopt in order to perform a thorough comparison of recordings obtained from the 11 different subjects. The guiding principle when combining data was to keep and use the largest amount of signal, avoiding arbitrary removal of noisy channels or the creation of a subset of "golden" cases. Descriptive statements are accompanied by the assessment of statistical significance of the results, taking into account the multiplicity of the hypotheses under testing; the obtained claims can give hints as to the mechanisms underlying SWA in mammals. 2 Python Software Foundation. Python Language Reference, version 2.7. 5 In general, discontinuities are not critical, since the failure corresponds to an interruption of the DAQ for a limited time interval (usually, a couple of discontinuities per channel, lasting from a few hundreds of millisecond to a few seconds); therefore, discontinuities in the DAQ are managed in the SWAP by identifying, for each problematic channel, the time interval at which the data acquisition fails, and removing it from the time sequence of the signal. By contrast, excluded channels are those for which the signal presents several irregularities, usually resulting in a number of identified upward transitions well below the median computed over the full channel set (tagged as "failure in the DAQ"); in addition, SD-outlier channels are also excluded; the tag "negative asymmetry" corresponds to the case of having a strong negative skew.
Furthermore, the outcome of the data analysis can be employed to feed the input of a dedicated spiking neural network simulation (as [7][8][9]), in a data-driven approach of in silico studies of the brain, aimed at computing a less stereotypical and a more accurate reproduction of cerebral rhythms. The analysis pipeline itself, hereafter named Slow Waves Analysis Pipeline, or SWAP, can be adopted to study the output of the simulation, and to define a set of benchmark observables for confronting models and experiments, with the goal of having a reliable and flexible set of tools available for the characterization of the slow-wave signal in a wide set of cases.
The material in this paper is structured as follows. Section 2 offers an overview of the SWAP analysis pipeline, with a collection of results obtained from test-bench data used for illustrating the steps of the procedure; the focus is on the methods implemented for the identification of the Up/Down state alternation in the ECoG traces, on the techniques for "stacking" and comparing different subjects, and on the statistical treatment of data with hypothesis testing. Section 3 is dedicated to the discussion of methods and results, with concluding remarks and suggestions for future research.

The Slow Waves Analysis Pipeline (SWAP)
The study of SWA expressed by the cortex can be tackled starting from a description of the bimodality (i.e. the alternation of Up and Down states), by defining a comprehensive set of observables suited to illustrating the phenomenon of SO. Once this local information is acquired, the second step (not illustrated here) is the characterization of the activity propagation across the brain surface as a wave with delta rhythm. Focusing on the features of Up and Down states, differences between cortical sites can be emphasized.
The analysis procedure consists of two phases. First, each data file is examined separately (Single Experiment), yielding a detailed inspection of the single subject. The results from the different experiments are then combined (inter-session data) to produce Summary Results . Finally, conclusive claims are given with the assessment of statistical significance of the results, taking into account the numerousness of the sample and the multiplicity of the hypotheses under testing.
Three different levels of description can be enabled: (I) the channel level, providing information at each recording site; (II) the area level, stacking up the information of all the electrodes located in the same cortical area; and (III) the full-set level, computing average properties and summary statements for the entire cortex portion under study. The levels of description can be applied at the single experiment, or at the inter-session data; outcomes obtained at the different levels of description are complementary and not strictly separated, and can be superimposed on graphical representations.  2.1. The channel-level description of the single experiment: from the raw data to the estimates of MUA and of transition times The channel-level description of the single experiment is obtained with a set of scripts, coded in MATLAB R , carrying out a loop over the electrodes in the array to extract the MUA from recordings of the raw signal (UFP). For each channel, the Power Spectral Density (PSD) of the signal is computed and the MUA is used as an estimate of the firing rate of the neurons around the electrode tip [2,10]. The logarithm of the MUA is evaluated, and the shape of the log(MU A) distribution can be described as a peak, at low-MUA values corresponding to Down states, and a tail, at high-MUA values corresponding to the Up states. The log(MU A) peak is fitted with a Gaussian function, and the parameters of the fit are used to single out Down-state periods from Up-state periods in the MUA time series. Once the MUA time series is tagged as "Up" and "Down" (binary MUA), a detailed study of the features of such states and of the transitions among them -upward (Down-to-Up) and downward (Up-to-Down) -can be achieved. The workflow is illustrated in Figure 3. In more detail, the initialization phase is carried out with the steering file setParamsAndOptions, which takes into account the specificity of the data acquisition (DAQ) sessions, accommodating the frame of reference (the positions of the electrodes in the array) and loading information about the recordings. Some of the settings concern the capability of the algorithm to identify the state transitions, and are fine-tuned following a heuristic approach. The steering file informs what to check in order to perform the analysis pipeline on the given input data. Once the initialization phase is completed, the main loop starts and the flow in the pipeline is as follows (see Figure 4): 1. Read the analog data (raw signal). For the ECoG data used as a benchmark in this study, input-related actions (open/close the input file and read the input data) are performed by the SONLibrary [11], whose functions provide the access to the stored neurophysiological data.
2. Analyze the raw signal in the frequency domain. A moving window of samples is defined for the calculation of the spectrogram, i.e. the spectral content of the raw signal as a function of time; the extent of the moving widow, together with the sampling frequency, determines the frequency band to be examined. Since the frequency band of interest for the estimate of the MUA is (200 − 1500) Hz 3 , a moving window of 5 ms is adopted 4 . FFT (Fast Fourier Transform) and PSD (Power Spectral Density) are computed using MATLAB R functions. Then, for each frequency in the band, the median of the PSD is computed, considering the full set of spectrograms, i.e. the collection of moving windows (time steps) that constitute the entire acquisition session. The obtained vector of medians, one value for each frequency, is used as a baseline to normalize the PSD.
3. Evaluate the MUA for each time step as the mean amplitude of the normalized PSD. The MUA is intended as an estimate of the collective firing rate r(t) of neurons located at the electrode position. The natural logarithm of the MUA is computed, since logarithmic mapping is adopted to emphasize the bimodality of the distribution; because of the normalization to the median of the PSD, negative values and positive values of log(MU A) identify a spectral content smaller or larger than the median, respectively [2,10].
4. Fit the distribution of log(MU A), isolating with a Gaussian function the peak corresponding to the (dominant) regime of low-rate states (Figure 4.A). A simplified description of the bimodality of the SWA in the cortex would require a bimodal fit, and initially the sum of two Gaussian profiles was adopted to describe the shape of the distribution. Conversely, what was discovered after a thorough inspection of a huge number of channels from several different animals (both physiological and pathological subjects at different anesthesia levels) is that, while the Down-state peak is highly stable despite the large variability in the subjects, the content of the log(MU A) histogram corresponding to the high-rate regime (obtained from the total distribution after subtracting the Down-state peak) expresses over a large span of "shapes", and rather than as a "second peak" (with a definite Gaussian profile) it can be generically appointed as a "tail" ( Figure  4.A, inset).
5. Set the UD_THRESHOLD, i.e. the level of log(MU A) that defines the separation of the two regimes of the bimodality, Up and Down (UD). This is a crucial step in the pipeline, and several options can be adopted to find the optimal criterion, which can depend on the specificity of the data set and on the scope of the analysis. In general, the choice takes into account general settings of the DAQ, the log(MU A) distribution, and the results of the Gaussian fit on the Down-state peak (of the given channel, or considering average properties of the recording session). At the channel level in the main loop, checks are introduced to monitor the convergence of the Gaussian fit, the content of the histogram of log(MU A) values, and the shape of the distribution ( Figure  4.A); alerts are activated to signalize: Weak Bimodality (if the area of the tail is smaller than a given threshold -currently set at 10% of the total); Positive Skewness (if gamma -the coefficient of skewness -of the tail is above 1, γ > 1); Negative Skewness (if γ < −1); Right Peak (if the peak, i.e. the dominant component of log(MU A), is centered at large values, on the right segment of the range, with a tail on the left); Large Threshold (if the selected threshold is larger than the mean of the tail); Few Transitions (if the number of transitions obtained with the selected threshold is below 3, the minimum requisite to isolate at least an Up state and a Down state). Some of these alerts may help in defining the threshold level, others provide an indication of "problematic" channels to be inspected (that may or may not be tagged as outliers and excluded from further processing at the end of the main loop). Finally, some conditions are "blocking", in the sense they require a break in the workflow, as is the case with noisy acquisition channels or spoiled recordings. If no blocking conditions are encountered, meaning that a threshold can be set and Down states and Up states are distinguishable, the following operations are performed: 3 The experience suggests that frequencies in the raw signal outside this range cannot be associated with the electro-physiological signal of the MUA; in particular, the highest frequency components reflect the presence of electrical noise introduced by the acquisition system. 4 For the benchmark data acquired at 5 kHz, the time window contains 25 samples.
(a) Convert the MUA into a binary sequence (BinaryMUA) that is 1 or 0 depending on the value of log(MU A) above or below the threshold, resulting in the MUA time sequence being tagged as "Up" or "Down" (Figure 4 As anticipated above, once the main loop execution has yielded a full description at the channel level, a key parameter to be monitored for the validation of the procedure is the stability of the conditions used for the identification of the two states (low-MUA and high-MUA). More precisely, since a requirement for the separability of Up and Down states is the successful fit of the Down-state peak of the log(MU A) with a Gaussian function (Figure 4.A), a similar value for the standard deviation (SD or σ) of the peak across channels is requested, ensuring comparable SO dynamics of the probed cell assemblies. Indeed, Down states are almost quiescent, and the variability of the MUA is mainly due to the acquisition chain, which has to be the same across channels.
A stable σ allows the application of the same MUA threshold at each recording channel, meaning a unique definition of the Up states, and thus more reliable profiles of traveling wave-fronts and a more sound description of the SWA as a collective phenomenon. Therefore, the choice of a fixed UD_THRESHOLD is a valid option when the goal of the analysis is to study the dynamics of the propagation of the activity as a slow wave across the cortex. On the other end, this means to decide a key parameter a priori, with the burden to be too conservative (increasing the false negative rate) or too tolerant (admitting a larger number of false positives) 5 . Conversely, the option of linking the choice of threshold to "internal agents" (e.g. parameters evaluated during the workflow of the pipeline) can be convenient, for reducing the number of free variables, or for anchoring the false positive rate per channel. A dedicated study of the standard deviation σ has been carried out (Section 2.2) on the test-bench data.
Finally, channels not fulfilling the stability requirements are excluded from the analysis and marked as outliers. Here again, the decision on the stability requirements (which parameters to focus on, which threshold levels, which weight to assign at the different instances, how to define the outliers) is a key-element of the pipeline and may largely depend on specific features of the recording sessions and on the goals of the data analysis. As discussed above, the SWAP pipeline has set in place alerts, warnings and counters based on parameters considered relevant for the test-bench data, but the elements to be monitored can change if conditions vary. Also configurable are the criteria that define outliers; for results presented here (Table 1), excluded channels are those with the Right Peak and with 5 In general, less conservative settings are preferable, since the control of false positives can be addressed by checking the spatio-temporal correlations of the propagating signal across the electrodes grid.  Few Transitions, together with the requirement related to the stability of the Gaussian peak, that leaves out channels with σ above Q3 + 1.5 × IQR (IQR is inter-quartile range Q3 − Q1, with Q1 first quartile and Q3 third quartile). As indicated in the caption of Table 1, the presence of discontinuity in the recording sequence -corresponding to a failure in the data acquisition and a drop in the log(MU A)is not a blocking element, since once the discontinuity is removed, the log(MU A) can fulfil the stability requirements.

Stability of the data sample and channel selection/rejection
A key assumption for comparing acquisitions taken at different sites and with different electrodes (or from different animals) is the comparability of the log(MU A) profile in the low firing rate regime. A quantitative evaluation of such a requirement is obtained by monitoring the width of the peak fitted with a Gaussian function, i.e. the σ parameter of the Gaussian (Figure 4.A). The purpose of this analysis step is to evaluate the selection strategy of the UD_THRESHOLD, which can be fixed or channel-dependent. A dedicated routine has been set in place, operating on a given collection of data (inter-session data from different subjects), aiming at stacking the full set of σ parameters estimated from fitting procedures. In Figure 5 we report the statistical distribution of σ values, and we observe large stability across channels and experiments. In more detail, the plot offers an overview of the range of variability, showing that the behavior of the parameter is pretty stable inside each experiment, and stable in the entire collection as well 6 . Therefore, we adopted a channel-dependent UD_THRESHOLD set at 2σ, corresponding to a fixed false positive rate per channel of about 2.25%. In addition, the results of the Stability Study enable channel selection or rejection based on the entire data sample, giving rise to a further list of outliers, to be added to the one already filled out for each DAQ session at the end of the main loop. Considering both lists, a total of 27 outliers were identified and removed from the analysis when inter-session data were taken into account for producing summary results.

A thumbnail for the single experiment
Once the raw signals have been analyzed for each channel, area-related statistics can be obtained as in [2]. More specifically, SWAP provides the following observables: Concerning d Down , d U p and d UD , as already evident from the histograms (Figure 4.C) and in agreement with the box plots represented in Figure 6.A, the distributions at the channel level are skewed and far from being Gaussian. Therefore, for each channel, the median (and not the mean) has been selected as the representative parameter for the observables. The skewness of the distribution is also noticeable at the area level ( Figure 6.B), where the statistics for each distribution are increased since values of electrodes belonging to the same cortical area are grouped together. Therefore, the median is assumed as the representative parameter also at the area level.
Interestingly, significant differences are apparent when comparing median values of different channels and areas. Differences among areas are consistently observed in all animals, despite the large span of values that each observable expresses considering the total of 11 experiments. Indeed, evaluating each observable for each experiment at the full-set level of description, the large variability of the data sample is clearly seen. This is in agreement with the expected biological variability of the subjects, regardless of the identical surgical treatment they have been undergone, the comparable drug delivery, and the uniform monitoring conditions during the data taking. In other words, each data session is characterized by its own central values for all the observables of interest, with no 6 Experimental outliers, discussed in Section 2.1 and listed in Table 1, are excluded from the representation. clear correlation with the animal's phenotype, or with any other parameter measurable during the data acquisition [12,13]. The spanning of the frequency values across the experiments (Figure 7) is in particular revealing, because frequency is a property immediately linkable with the onset of the bimodality and with the propagation of slow waves along the cerebral cortex.

Assessment of statistical significance at the area-level description for inter-session normalized data
The statistical significance of the effect of differentiation by area, visible for all the observables in the list of interest and for each experiment in the data set, is quantitatively assessed for summary results, i.e. when the outcome of the single experiments are combined to drive comprehensive claims on the phenomena. However, the first obstacle met when trying to compare and aggregate results is the large variability exhibited by the different DAQ sessions, which dominates over any other effect and disguises any similarity or common footprint among subjects. Therefore, to confront the area-level descriptions, the median values of the observables for each area are normalized for each experiment, computing for each observable the arithmetic mean across the cortical areas, and using the obtained mean as a normalization factor for that observable. This procedure highlights any trend expressed at different cortical areas, enabling us to check if different experiments express the same trend. Figure  8.A shows the summary result at the area-level description obtained with normalized data.
The outcome of the hypothesis tests executed to assess the statistical significance of the differences observed between each couple of cortical areas is represented through a matrix of p-values, with a graded scale in a given range of confidence (Figure 8.B). Since no assumptions are made on the model, and because of the already-discussed evidence of non-Gaussianity of the distributions, non-parametric approaches are followed when analyzing the test bench data, so the Wilcoxon rank-sum test is applied 7 . To take into account multiple comparisons (the so-called "look-elsewhere effect" 8 ) and reduce the likelihood of incorrectly rejecting a null hypothesis (type I error) when evaluating a family of simultaneous tests, the robustness of the analysis and the control of the false discovery rate (FDR) is obtained by correcting the p-values with the Benjamini-Hochberg (BH) procedure [15]. Analogously to what was done for state duration, the study of a possible differentiation at the area-level has also been carried out for slopes and maximum MUA (s U p , s Down , p). Slopes are computed considering the average transition (Figure 4.F-G), obtained pooling together the detected Up-to-Down and Down-to-Up transitions. The transition front of the average transition is isolated, considering a 35 ms interval around the transition time t 0 ([t 0 − 0.025, t 0 + 0.010] for downward transitions; [t 0 − 0.010, t 0 + 0.025] for upward transitions); the profile is fitted with a cubic and the slope is obtained as the derivative of the polynomial at t 0 . The average upward transition is also used for the estimation of the maximum MUA of the average waveform in the first 250 ms after the transition. Since slopes and maxima are calculated from the average waveform, for each experiment the observables are represented at the channel level by a single value (and not by a distribution). The area-level description for the single experiment is given by the mean and median of values obtained from all the channels belonging to the specific area. In coherence with the analysis carried out on states duration and in order to apply the same non-parametric Wilcoxon tests, the median is taken as the representative parameter. Summary results are produced after the same normalization adopted for states duration, yielding a similar illustration. A similar approach is followed for frequency f .

Interpolation of the array map
The results obtained with the high spatial resolution probe used for the acquisition of the test bench data can still be improved by interpolating the spaces between the electrodes; the accuracy of 7 The computation of the Wilcoxon rank-sum statistics is carried out with the statistical function scipy.stats.ranksums of the scipy Python module [14] 8 https://xkcd.com/882/. the interpolation is assured by the large surface density of sensors offered by the employed MEA. For each experiment and channel, the median (d Down , d U p , d UD ) or the mean (s U p , s Down , p, f ) of the observable is taken. This set of values is used to compute an interpolator 9 that is applied to points (xy coordinates) on a mesh-grid, made of 50 steps along x and 90 step along y, with a Grid Step ∼ 0.05 mm, i.e. 1/10 of the Array Step = 0.550 mm (Figure 9.A for the single experiment and for a given observable). The same illustration is applied when representing inter-session normalized data (Figure 9.B); here, the size of the marker is inversely proportional to the standard deviation of the distribution of values at the given electrode position, thus measuring the amount of variability registered at that position across the experiments.  Figure 1, with black dots marking the electrode positions and black lines identifying cortical areas; the white circle identifies the outlier channel as resulting from Figure 5 and Table 1. The color bar gives the range of the interpolated variable across the array. (B) Summary results for the observable d UD , represented as a contour plot obtained from normalized data. The interpolation is based on the mean values computed across the experiments, the marker size is inversely proportional to the standard deviation of the distribution of values at each position (the smaller the marker, the less variability measured across the experiments). Since normalized data are used for the plot, the color bar gives indications of trends with respect to the average: for instance, regions colored in blue are characterized by having a duration of the oscillation cycle that is up to 30% less than the average. Comparing (A) (single experiment) and (B) (inter-session data), a similar tendency is visible for the observable: in particular, the contours evidence gradients along the same dominant antero-posterior direction, from fronto-lateral towards occipito-medial regions. These gradients prevail over the borders of cortical areas, suggesting further differentiation within areas and connections among areas. 9 The interpolation is carried out with the scipy Python module [14] using the scipy.interpolate.Rbf class for radial basis function (Rbf) interpolation. More on https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.Rbf.html 2.5. From the array map to the assessment of statistical significance at the channel level for the normalized inter-session data The contour plots traced on the interpolated array maps give qualitative hints on differentiation among inspected cortical sites, enlightening further details within cortical areas. A more quantitative evaluation can be obtained by performing the statistical analysis of the normalized values resulting from inter-session data, considering each pair of electrodes in the MEA and inspecting if there is a statistically significant difference between the distribution of values measured therein. The idea is to compare the normalized values located at two different electrode positions in the map and to test the validity of the null hypothesis, that is the samples are extracted from the same statistical distribution 10 . We use the Wilcoxon rank-sum statistics with B-H correction on p-values; in this case, the FDR correction has a larger impact, since for MEAs the number of tests in the family is of the order of hundreds 11 .
Given the large number of hypotheses, electrodes are sorted to emphasize those expressing a significant difference with the others. More specifically, the top-three electrodes in this sorted list are highlighted as "core nodes", indicating the positions in the cortex displaying the largest significant difference with respect to other positions in the cortex (Figure 10.A).
As a result, the somatosensory cortex emerges as the one having the shortest mean Up/Down cycle, in agreement with what shown in Figure 8.A. In Figure 10.B, the differences in the cycle duration at the channel level show a rather homogeneous rhythm within single areas. It is also apparent that the other primary sensory area, i.e. the visual cortex, is markedly different from frontal and parietal regions, displaying a significantly longer oscillation period. Inspecting the changes in the Up and Down state duration, we found that the modulation of the oscillation cycle is mainly due to a reduction of the d Down , as shown in Figure 10.C. Indeed, only few channels display a significant difference in d U p (not shown) and the p-value matrix for the d Down mirrors the one shown in Figure 10.B. If the local excitability of the cell assemblies underneath the MEA contacts is well represented by the ratio d U p /d Down [2,16], here the leading role of somatosensory cortex is apparent. Such enhanced excitability is further confirmed by inspecting the slope of the Down-to-Up transition in the MUA [17], which is significantly steeper in the same cortical locations (Figure 10.D).

Discussion
Biological data are characterized by richness in details and large variability. The efforts of the data analysis should aim at extracting tendencies and regularities, producing a concise description without hiding or neglecting complexity and details that could convey informative content. This is the guideline followed when developing the Slow Waves Analysis Pipeline (SWAP), starting from a solid backbone that has been deeply revised and enriched with new features. In particular, the opportunity of using the MEA data as a test bench has allowed us to focus on the spatial differentiation of the observables, with the aim of uncovering hints as to the local excitability of the cortical assemblies. The developed methodology is robust and easily re-configurable, flexible and adaptable at different acquisition conditions, and also suitable to be applied to the output of simulations. In this framework, SWAP can be employed in bridge theory, simulations and experiments, providing a set of general tools that allow an effective comparison between heterogeneous data sets. The adoption of a unique analysis procedure is also useful for comparing different simulation engines; the SWAP can be applied to define benchmarks and evaluate the performance of numerical models and implementations. Several studies are ongoing for the application of SWAP to a large variety of datasets (knock-out mice, subjects in 10 The number of samples at each location depends on how many experiments have channels rejected as outliers at that location. For the test bench data set, a maximum of 11 samples can be present at each electrode position, details for each electrode position can be extracted from Figure 5.D. 11 If the number of electrodes in the MEA is n = 32, the number of hypotheses to be checked is N = n(n−1)  Figure 9.B) with the "core-nodes" (i.e the three electrodes with the largest number of significant differences with the other channels). Lines are traced to connect the "core nodes" with their partners in the electrode pair being significantly different (p < 0.05). (B) Matrix of p-values of the statistical test on the d UD differences between channels. The matrix construction is analogous to what is described for Figure 8. The numeric sequence of electrodes in the matrix is as in Figure 6 Concerning the interpretation of the results, the analysis of the large set of data collected from 11 wild-type anesthetized mice with the 32-channel MEA allows us to claim a statistically significant differentiation of cortical areas for several parameters that characterize the onset of SWA along the cerebral cortex. Starting from these observables, the excitability of the cortical tissue expressing SWA can be investigated. Indeed, larger excitability is expected to be associated with faster transition fronts (in particular, upward slopes), shorter duty-cycles (i.e. smaller d UD , dominated by the duration of the Down states) and accordingly larger frequencies [17]. For instance, a smaller d Down reveals the case in which excitability translates in faster Down-to-Up transitions. These features are particularly apparent in the somatosensory area, likely the most excitable cortical region we observed; conversely, the occipital lobe (retrosplenial and visual areas) act as the least excitable. Activation waves traveling across the cortex during SWA are expected to be sensitive to the diverse cortical excitability, as more reactive (i.e. more excitable) areas are expected to display a smaller "inertia" in recruiting neurons involved in the collective phenomena associated with the high-firing Up states. As a consequence, waves might be originating from highly excitable regions such that preferential propagation pathways are expected, as previously highlighted both in humans [18] and rodents [2,19,20]. In turn, heterogeneous excitability might also underlie the non-global nature of the SWA phenomenon observed when wakefulness is approaching [21,22] or under pathological conditions [23].
It must be pointed out that the study here presented is "static", in the sense that the absolute time sequence of the events is not taken into account: states and transitions are time-squashed and stacked regardless of their time of occurrence in the DAQ session, and thus any time-dependent effect, like fatigue and recovery intervals of the neurons, that could affect the excitability and alter the responsiveness of cortical regions, is excluded from this analysis. In other words, the figures presented in this paper can be seen as the average photography of the phenomena investing the monitored cortex during the acquisition period. Related to this, excitability and responsiveness can also be altered by wave propagation dynamics, in particular if different propagation patterns coexist and overlap in the same time interval, for instance when two slow waves originating at distinct sites travel along the cortex at different speeds and each one follows its own path [4,24]. Again, not using the information on the absolute timing of the events at the electrode positions and intending the results presented here to be the average photography of the SWA in the cortex, the included results reflect the dominant propagation pattern, namely the antero-posterior direction [2,18,19,25], in particular along an oblique axis directed from fronto-lateral towards occipito-medial regions, as suggested by values depicted in the contour plots. As in this step of the SWAP we have not focused on time-dependent effects, their impact on the area-differentiation of the bistability modes will be evaluated in further studies.
Finally, although in all the animals involved in this study the level of administered anesthesia was the same, the observed inter-subject variability (Figure 7) could in principle be explained by animals being in different brain states. This suggests the possibility of exploiting the characterization of slow oscillations in the cerebral cortex as a new tool for effective classification of the brain states. Several techniques are currently under study, based for instance on the principal components analysis (PCA) to identify and single out the different sources of variability in the experimental data set. Indeed, a more reliable classification of brain states (i.e. of the DAQ sessions that constitute the statistical sample) would provide a more robust comparison of the experiments, allowing us to overcome the limits derived from the need to use normalized data.
improved it with new functions, applied it to the test bench data, performed the statistical analysis, wrote the first draft of the manuscript. PSP and MM contributed to the data analysis and to the interpretation of the results. AP revised sections of the manuscript. All authors contributed to manuscript revision, and read and approved the submitted version.

Funding
This work has received funding from the European Union Horizon 2020 Research and Innovation Programme under Specific Grant Agreements No. 785907 (HBP SGA2) and No. 720270 (HBP SGA1).