Analysis Pipeline for Extracting Features of Cortical Slow Oscillations

1. Introduction

In a brain manifesting slow-wave activity (SWA), expressed in the cerebral cortex under NREM sleep and deep anesthesia (Steriade et al., 1993), the spiking activity—both single and multi-unit (SUA and MUA), respectively—appears as a regular sequence of Up (high-rate) and Down (almost quiescent) states. Detecting the Up and Down periods and the times of transitions between these states is a long-standing issue for those studies investigating the phenomenon. Depending on the available experimental data, different solutions have been developed in the past. At the microscale, threshold-based detection algorithms, leveraging the alternating polarization of the membrane potential in intracellular recordings, have been proposed (Volgushev et al., 2006; Seamari et al., 2007). At the mesoscale, a similar approach has been proposed that focuses on extracellular field potentials by extracting the power of the signal at specific frequency bands (Mukovski et al., 2006). Besides, multi-unit activity sampled from multi-site recordings at the macroscale led to the development of a relatively more sophisticated algorithm detecting Up and Down state transitions by fitting a Hidden Markov Model (Jercog et al., 2017). Here, we improved an approach that relies on the mesoscale signal provided by the MUA extracted as the high-frequency component of the unfiltered field potential. The analysis pipeline, implemented in MATLAB¹, has already been applied to several experimental data sets, acquired with different setups both in vitro and in vivo from rodents, ferrets, and monkeys (Ruiz-Mejias et al., 2011, 2016; Mattia et al., 2013; Capone et al., 2019). In this work, we describe an extensive revision of the above referenced procedure, including some new features implemented in Python². The novel components focus on labor-saving automatic adjustment of steering parameters, on automated comparison of data coming from different experiments, and on the evaluation of statistical significance, offering a set of software tools that allow going from raw-data to statistical assessment of results, with consequent advantages in terms of reliability, flexibility, and usability. The current stage of development is the starting point for a cooperative effort aiming at the release to the community of an engineered pipeline³.

The enhanced analysis pipeline has been validated on a new set of data collected in vivo using micro-ECoG electrodes (Wang et al., 2009). In particular, the system for brain activity recording employed in this work consists of a 32-electrode array and the test-bench data has been acquired 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 (Pazzini et al., 2018). 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; that is, 11 independent experiments that collectively represent an extensive data sample covering a wide range of biological and unavoidable physiological variability.

FIGURE 1

Figure 1. Representation of the multi-electrode array (MEA) used for the data acquisition (Pazzini et al., 2018), with the numbering introduced in the SWAP. 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 pipeline. On the bottom left, the color legend that identifies the cortical areas on the array surface. On the right, an illustration of the grid of electrodes positioned on the mouse cortex (image credit: Allen Institute); the inspected surface is of the order of 10 mm².

TABLE 1

Table 1. Summary of the data acquisition (DAQ) sessions.

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.

Furthermore, the outcome of the data analysis can be employed to feed the input of a dedicated spiking neural network simulation (as Paolucci et al., 2013; Pastorelli et al., 2018, 2019), 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, as well as constraining simulations of the cognitive effects of sleep (Capone et al., 2019). The analysis pipeline itself, hereafter named Slow Wave 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, 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.

2. The Slow Oscillation Analysis Pipeline

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. Figure 2 details, in terms of phases of action and levels of description, the logic of the SWAP when facing the issue of investigating the features of Up and Down states.

FIGURE 2

Figure 2. The logic of the SWAP when facing the issue of investigating the features of Up and Down states, in terms of phases of action (the sequence of 3 actions listed on the left) and levels of description (in the upper box). The lower box contains information on the statistical treatment presently implemented.

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, 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 (Mattia et al., 2010; Ruiz-Mejias et al., 2011). The logarithm of the MUA is evaluated, and the shape of the log(MUA) 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(MUA) 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.

FIGURE 3

Figure 3. The channel-level description of the single experiment. Purple boxes point out the initial data format and the intermediate outcomes of the workflow (as successive elaborations of the raw data). The black box encompasses the key steps of the process. The procedure is coded in MATLAB^®, the sequence of actions is illustrated in Figure 4. The figure illustrates the main loop (as discussed in section 2.1), at the end of which further checks are carried out on average and median values (full-set level description), in order to identify further anomalies or outliers.

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 (Lidierth, 2009), 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; that is, 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–1,500 Hz⁴, a moving window of 5 ms is adopted⁵. FFT (Fast Fourier Transform) and PSD (Power Spectral Density) are computed using MATLAB functions. Then, for each frequency in the band, the median of the PSD is computed, considering the full set of spectrograms; that is, 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(MUA) identify a spectral content smaller or larger than the median, respectively (Mattia et al., 2010; Ruiz-Mejias et al., 2011).

4. Fit the distribution of log(MUA), isolating with a Gaussian function the peak corresponding to the (dominant) regime of low-rate states (Figure 4A). 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(MUA) 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 4A, inset).

5. Set the UD_THRESHOLD; that is, the level of log(MUA) 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(MUA) 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(MUA) values, and the shape of the distribution (Figure 4A); 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(MUA), 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:

(a) Convert the MUA into a binary sequence (BinaryMUA) that is 1 or 0 depending on the value of log(MUA) above or below the threshold, resulting in the MUA time sequence being tagged as “Up” or “Down” (Figure 4B). Aiming at illustrating the sequence of Up and Down states on different signals, the panel includes log(MUA) and BinaryMUA together with unfiltered and filtered field potential; the latter is obtained applying the MATLAB bandpass function with f_pass = [200,1,500] Hz and default settings. The horizontal axis corresponds to 10 s of recordings; the vertical axis is in arbitrary units (the raw signal is measured in μV, while MUA is dimensionless since obtained from a ratio of homogeneous quantities), signals values are scaled in order to make all the traces fit on the same plot.

(b) Label the transitions as upward or downward, and assign the time of transition (Trigger Time) with a cubic interpolation of the waveform to locate the time at which the level of MUA crosses the threshold that separates the two regimes of low and high firing rate activity.

(c) Study the features of states and transitions, to check the robustness of the algorithm and the stability of the outcome for the different channels. In particular:

i. Histograms of the duration of Down states d_DOWN, Up states d_UP, UD-cycles d_UD. The study of these observables is one of the focuses of the analysis; at the channel-level description of the SWAP, the distributions are superimposed for comparison (Figure 4C).

ii. Raster plots of states and transitions; each transition (Trigger Event) is centered at the trigger time; events can be sorted by their time of occurrence, or by their duration (Figures 4D,E).

iii. Waveform plots, for comparing states and transitions, and plots of the average waveform, obtained considering the full set of transitions. A refined algorithm isolates the transition front, to better investigate the transition dynamics (Figures 4F,G).

FIGURE 4

Figure 4. (A) Distribution of log(MUA) as normalized histogram of values shifted at the position of the first peak, The blue curve is the Gaussian fit of the first peak, the μ position is marked by the blue vertical dashed line. The red dotted line is the “tail,” obtained after subtracting the Gaussian fit from the total distribution (see inset; vertical dashed lines: mode in magenta, median in cyan, mean in black). In this example (Ch. 1, file 01), the tail distribution is far from being Gaussian; the area of the tail is ~27%, therefore expressing a clear bimodality, coherently with the value of the skewness of the distribution, γ = 1.97. (B) 10 s of recordings of log(MUA) (in black), BinaryMUA (in red), unfiltered raw signal (in blue) and filtered raw signal (in green). (C) Histograms of the duration of Up states (red), of Down states (blue) and of the full cycles (UD cycle, in green); dashed lines, mean values. (D,E) Raster plot of upward (Down-to-Up) and downward (Up-to-Down) transitions, sorted by state duration. (F,G) Waveforms (WFs) of upward and downward transitions, respectively. For each plot, the first 5 transitions are shown, centered around the transition time (Trigger Time); the average transition, computed considering the full set of n waveforms, is superimposed (black); the shaded area identifies the profile of the waveform at ±1 standard deviation (SD) off the mean. The inset shows the average transition front; the shaded area identifies the profile of the waveform at ±1 standard error (SEM) of the mean (H,I) As (A,B), for a channel rejected by the pipeline and excluded from further analysis. The channel presents a monomodal log(MUA) with negative asymmetry; as a consequence, no UD_THRESHOLD can be set, and no transitions can be identified, evidenced by the filtered field potential (in green) with no signs of state transition and by the BinaryMUA (in red) that appears constant at low level. Additional details are included in the Supplementary Material.

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(MUA) with a Gaussian function (Figure 4A), 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)⁶. 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 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 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(MUA)—is not a blocking element, since once the discontinuity is removed, the log(MUA) can fulfill the stability requirements. An illustration of a rejected channel is given in Figures 4H,I; further details are included in the Supplementary Material.

2.2. 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(MUA) 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; that is, the σ parameter of the Gaussian (Figure 4A). 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.

FIGURE 5

Figure 5. Stability of the MUA estimate across channels and animals. (A) Stem plot of the SD values obtained as the σ parameter of the Gaussian fit to the Down peak in the log(MUA) distribution (the blue curve in Figure 4A); marker colors identify the different experiments. The distribution of the SD values is presented as histogram (B) and as box plot (C), the presence of outliers is evident; descriptive parameters are reported in the box; the vertical lines in the histograms indicate the positions displaced 2 or 3 standard deviations of the mean. Outlier values are removed, obtaining the symmetric distribution represented in the histogram in (E); descriptive parameters are listed in the box, the sample of channels is reduced by 12 units. The grid in (D) is a synthetic representation of the outlier channels for the test bench data: in red, the experiment outliers; in blue, the ones excluded after the SD Stability Study.

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⁷. 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.

2.3. Description at the Cortical Area-Level

2.3.1. A Thumbnail for the Single Experiment

Once the raw signals have been analyzed for each channel, area-related statistics can be obtained as in Ruiz-Mejias et al. (2011). More specifically, the SWAP provides the following observables:

1. Durations of the Down states, d_Down[s];

2. Durations of the Up states, d_Up[s];

3. Durations of the UD-cycles (i.e., a pair of consecutive Down and Up states), d_UD[s]. The duration of the UD-cycle is an estimate of the SO period;

4. Upward transition slope, $s_{Up}\left[s^{-1}\right]$;

5. Downward transition slope $s_{Down}\left[s^{-1}\right]$;

6. Maximum MUA in the Up state (“peak”), p[a.u.]

7. Frequency, f[Hz], defined as $\frac{1}{d_{UD}}$;

Concerning d_Down, d_Up and d_UD, as already evident from the histograms (Figure 4C) and in agreement with the box plots represented in Figure 6A, 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 6B), 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.

FIGURE 6

Figure 6. Statistics of the measured UD-cycle (d_UD) for an example recording (file 01). (A) Channel-level description; channels are grouped by area, channel numbering is as illustrated in Figure 1; channel 8 is missing because it is tagged as an outlier (Table 1). (B) Area-level description, represented with notches indicating the confidence interval for the median (Q2 = 50%, orange line). The box plot is delimited by quartiles Q1 = 25% and Q3 = 75%; IQR = Q3 − Q1 is the Inter Quartile Range; the lower whisker is at Q1 − 1.5 × IQR; the upper whisker is at Q3 + 1.5 × IQR; values outside the whiskers are marked as outliers; the blue marker indicates the mean of the distribution.

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 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 clear correlation with the animal's phenotype, or with any other parameter measurable during the data acquisition (Brown et al., 2010; Akeju and Brown, 2017). 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.

FIGURE 7

Figure 7. Mean frequency across channels for each experiment; error bars, SEM (standard error of the mean). For each channel, the frequency f[Hz] is computed as $\frac{1}{<d_{UD}>}$, with <> denoting the arithmetic mean of the d_UD[s] values in the channel (< d_UD[s] > is an estimate of the SO period). The text box gives information on the distribution of the observable across the 11 experiments.

2.3.2. 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; that is, 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 8A shows the summary result at the area-level description obtained with normalized data.

FIGURE 8

Figure 8. (A) Box plot representation of the summary results for the observable d_UD at the area-level description, obtained with normalized data. Having 11 subjects, for each cortical area 11 median values are available (the orange line in Figure 6); the arithmetic mean of the 5 medians is the normalization factor for the given experiment. (B) Outcome of the Wilcoxon hypothesis tests on the same data; the test results are represented via a matrix of p-values (with Benjamini-Hochberg correction enabled).

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 8B). 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⁸. To take into account multiple comparisons (the so-called “look-elsewhere effect”⁹) 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 (Seabold and Perktold, 2010).

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_Up, s_Down, p). Slopes are computed considering the average transition (Figures 4F,G), obtained by 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₀ ([t₀ − 0.025, t₀ + 0.010] for downward transitions; [t₀ − 0.010, t₀ + 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₀. 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.

2.4. 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 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_Up, d_UD) or the mean (s_Up, s_Down, p, f) of the observable is taken. This set of values is used to compute an interpolator¹⁰ that is applied to points (xy coordinates) on a mesh-grid, made of 50 steps along x and 90 steps along y, with a Grid Step ~0.05 mm, i.e., 1/10 of the Array Step = 0.550 mm (Figure 9A for the single experiment and for a given observable). The same illustration is applied when representing inter-session normalized data (Figure 9B); 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 9

Figure 9. (A) Contour plot for the single experiment (here, file 01) illustrating the observable d_UD (UDcycleLen) interpolated on the mesh-grid, based on values recorded with the MEA. The array map reflects the DAQ schema of 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 toward occipito-medial regions. These gradients prevail over the borders of cortical areas, suggesting further differentiation within areas and connections among areas.

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¹¹. We use the Wilcoxon rank-sum statistics with BH 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¹².

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 10A).

FIGURE 10

Figure 10. Excitability maps of the probed cortices. (A) Normalized d_UD (as in Figure 9B) 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, i.e., channels are grouped by area (M, motor, S, somatosensory, P, parietal, R, retrosplenial, V, visual). (C,D) Normalized Down state duration d_Down (DownStateLen) and Down-to-Up transition slope s_Up (SlopeUp), averaged across experiments, as in (A).

As a result, the somatosensory cortex emerges as the one having the shortest mean Up/Down cycle, in agreement with what is shown in Figure 8A. In Figure 10B, 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; that is, 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 10C. Indeed, only a few channels display a significant difference in d_Up (not shown) and the p-value matrix for the d_Down mirrors the one shown in Figure 10B. If the local excitability of the cell assemblies underneath the MEA contacts is well represented by the ratio d_Up/d_Down (Ruiz-Mejias et al., 2011; Mattia and Sanchez-Vives, 2012), 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 (Sanchez-Vives et al., 2010), which is significantly steeper in the same cortical locations (Figure 10D).

3. Discussion

Biological data are characterized by richness in detail 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 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, the 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 the SWAP to a large variety of data sets (knock-out mice, subjects in different brain states, data collected with optical techniques), and the stability and reliability of the analysis procedure has so far been confirmed. Moreover, several improvements are ongoing, focusing on the reliability and robustness of the analysis algorithms when increasing the number of channels (electrodes or pixels), as in recently developed multi-electrode arrays and in optical imaging data. The latter in particular is a case where we have already successfully applied our protocol (see Celotto et al., 2018, calcium-imaging data, GCaMP6f model, submitted for publication). In general, the approach we follow can be extended provided that single-channel records allow the reliable disentanglement of Up and Down states of the underlying cortical assemblies. The introduced new features in the analysis pipeline have been coded in Python with the aim of realizing open software tools for the scientific community. The complete transition toward open software is in the list of objectives to fulfill in the near future.

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 (Sanchez-Vives et al., 2010). For instance, a smaller d_Down reveals the case in which excitability translates to 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) acts 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 (Massimini et al., 2004) and rodents (Ruiz-Mejias et al., 2011; Mohajerani et al., 2013; Stroh et al., 2013). In turn, heterogeneous excitability might also underlie the non-global nature of the SWA phenomenon observed when wakefulness is approaching (Nir et al., 2011; Vyazovskiy et al., 2011) or under pathological conditions (Sanchez-Vives et al., 2017).

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 (Sancristóbal et al., 2016; Capone et al., 2019). 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 (Massimini et al., 2004; Ruiz-Mejias et al., 2011; Stroh et al., 2013; Chan et al., 2015), in particular along an oblique axis directed from fronto-lateral toward 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, and 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 using 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.

Data Availability Statement

Data will be available on reasonable request.

Ethics Statement

This study was carried out 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.

Author Contributions

MS-V, MM, MD, PP, and GD conceived and designed the study. MD performed the data acquisition under the supervision of MS-V. MM provided the backbone of the analysis pipeline. GD revised the analysis pipeline, improved it with new functions, applied it to the test bench data, performed the statistical analysis, wrote the first draft of the manuscript. PP 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).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This study was carried out in the framework of the Human Brain Project (HBP¹³), funded under Specific Grant Agreements No. 785907 (HBP SGA2) and No. 720270 (HBP SGA1), in particular within activities of sub-project SP3 (Systems and Cognitive Neuroscience). This manuscript has been released as a Pre-Print at De Bonis et al. (2019).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnsys.2019.00070/full#supplementary-material

Footnotes

1. ^MATLAB^®, The MathWorks, Inc., Natick, Massachusetts, United States.

2. ^Python Software Foundation. Python Language Reference, version 2.7.5. Available online at: http://www.python.org

3. ^A preliminary version of the SWAP, released in the framework of the Human Brain Project, is Available online at: https://github.com/INM-6/wavescalephant

4. ^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.

5. ^For the benchmark data acquired at 5 kHz, the time window contains 25 samples.

6. ^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.

7. ^Experimental outliers, discussed in section 2.1 and listed in Table 1, are excluded from the representation.

8. ^The computation of the Wilcoxon rank-sum statistics is carried out with the statistical function scipy.stats.ranksums of the scipy Python module (Jones et al., 2001)

9. ^https://xkcd.com/882/

10. ^The interpolation is carried out with the scipy Python module (Jones et al., 2001) 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

11. ^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 5D.

12. ^If the number of electrodes in the MEA is n = 32, the number of hypotheses to be checked is $N=\frac{n(n-1)}{2}=496$.

13. ^https://www.humanbrainproject.eu/en/

References

Akeju, O., and Brown, E. N. (2017). Neural oscillations demonstrate that general anesthesia and sedative states are neurophysiologically distinct from sleep. Curr. Opin. Neurobiol. 44, 178–185. doi: 10.1016/j.conb.2017.04.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Brown, E. N., Lydic, R., and Schiff, R. (2010). General anesthesia, sleep, and coma. N. Engl. J. Med. 363, 2638–2650. doi: 10.1056/NEJMra0808281

PubMed Abstract | CrossRef Full Text | Google Scholar

Capone, C., Pastorelli, E., Golosio, B., and Paolucci, P. S. (2019). Sleep-like slow oscillations improve visual classification through synaptic homeostasis and memory association in a thalamo-cortical model. Sci. Rep. 9:8990. doi: 10.1038/s41598-019-45525-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Capone, C., Rebollo, B., Muñoz, A., Illa, X., Del Giudice, P., Sanchez-Vives, M. V., et al. (2019). Slow waves in cortical slices: how spontaneous activity is shaped by laminar structure. Cereb. Cortex 29, 319–335. doi: 10.1093/cercor/bhx326

PubMed Abstract | CrossRef Full Text | Google Scholar

Celotto, M., De Luca, C., Muratore, P., Resta, F., Allegra Mascaro, A. L., Pavone, F. S., et al. (2018). Analysis and model of cortical slow waves acquired with optical techniques. arXiv:1811.11687.

Google Scholar

Chan, A. W., Mohajerani, M. H., LeDue, J. M., Wang, Y. T., and Murphy, T. H. (2015). Mesoscale infraslow spontaneous membrane potential fluctuations recapitulate high-frequency activity cortical motifs. Nat. Comm. 6, 1–12. doi: 10.1038/ncomms8738

PubMed Abstract | CrossRef Full Text | Google Scholar

De Bonis, G., Dasilva, M., Pazienti, A., Sanchez-Vives, M. V., Mattia, M., and Paolucci, P. S. (2019). Slow waves analysis pipeline for extracting the features of the bi-modality from the cerebral cortex of anesthetized mice. arXiv:1902.08599.

Google Scholar

Jercog, D., Roxin, A., Barthó, P., Luczak, A., Compte, A., and de la Rocha, J. (2017). Up-down cortical dynamics reflect state transitions in a bistable network. eLife 6:e22425. doi: 10.7554/eLife.22425

PubMed Abstract | CrossRef Full Text | Google Scholar

Jones, E., Oliphant, T., and Peterson, P. (2001). SciPy: Open Source Scientific Tools for Python. Available online at: http://www.scipy.org/ (accessed November 13, 2019).

Google Scholar

Lidierth, M. (2009). sigTOOL: a MATLAB-based environment for sharing laboratory-developed software to analyze biological signals. J. Neurosci. Methods 178, 188–196. doi: 10.1016/j.jneumeth.2008.11.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Massimini, M., Huber, R., Ferrarelli, F., Hill, S. L., and Tononi, G. (2004). The sleep slow oscillation as a traveling wave. J. Neurosci. 24, 6862–6870. doi: 10.1523/JNEUROSCI.1318-04.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Mattia, M., Ferraina, S., and Del Giudice, P. (2010). Dissociated multi-unit activity and local field potentials: a theory inspired analysis of a motor decision task. Neuroimage 52, 812–823. doi: 10.1016/j.neuroimage.2010.01.063

PubMed Abstract | CrossRef Full Text | Google Scholar

Mattia, M., Pani, P., Mirabella, G., Costa, S., Del Giudice, P., and Ferraina, S. (2013). Heterogeneous attractor cell assemblies for motor planning in premotor cortex. J. Neurosci. 33, 11155–11168. doi: 10.1523/JNEUROSCI.4664-12.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Mattia, M., and Sanchez-Vives, M. V. (2012). Exploring the spectrum of dynamical regimes and timescales in spontaneous cortical activity. Cogn. Neurodyn. 6, 239–250. doi: 10.1007/s11571-011-9179-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Mohajerani, M. H., Chan, A. W., Mohsenvand, M., LeDue, J., Liu, R., McVea, D. A., et al. (2013). Spontaneous cortical activity alternates between motifs defined by regional axonal projections. Nat. Neurosci. 16, 1426–1435. doi: 10.1038/nn.3499

PubMed Abstract | CrossRef Full Text | Google Scholar

Mukovski, M., Chauvette, S., Timofeev, I., and Volgushev, M. (2006). Detection of active and silent states in neocortical neurons from the field potential signal during slow-wave sleep. Cereb. Cortex 17, 400–414. doi: 10.1093/cercor/bhj157

PubMed Abstract | CrossRef Full Text | Google Scholar

Nir, Y., Staba, R. J., Andrillon, T., Vyazovskiy, V. V., Cirelli, C., Fried, I., et al. (2011). Regional slow waves and spindles in human sleep. Neuron 70, 153–169. doi: 10.1016/j.neuron.2011.02.043

PubMed Abstract | CrossRef Full Text | Google Scholar

Paolucci, P. S., Ammendola, R., Biagioni, A., Frezza, O., Lo Cicero, F., Lonardo, A., et al. (2013). Distributed simulation of polychronous and plastic spiking neural networks: strong and weak scaling of a representative mini-application benchmark executed on a small-scale commodity cluster. arXiv:1310.8478.

Google Scholar

Pastorelli, E., Capone, C., Simula, F., Sanchez-Vives, M. V., Del Giudice, P., Mattia, M., et al. (2019). Scaling of a large-scale simulation of synchronous slow-wave and asynchronous awake-like activity of a cortical model with long-range interconnections. Front. Syst. Neurosci. 13:33. doi: 10.3389/fnsys.2019.00033

PubMed Abstract | CrossRef Full Text | Google Scholar

Pastorelli, E., Paolucci, P. S., Simula, F., Biagioni, A., Capuani, F., Cretaro, P., et al. (2018). “Gaussian and exponential lateral connectivity on distributed spiking neural network simulation,” in 2018 26th Euromicro International Conference on Parallel, Distributed and Network-Based Processing (PDP) (Cambridge, UK), 658–665. doi: 10.1109/PDP2018.2018.00110

CrossRef Full Text | Google Scholar

Pazzini, L., Polese, D., Weinert, J. F., Maiolo, L., Maita, F., Marrani, M., et al. (2018). An ultra-compact integrated system for brain activity recording and stimulation validated over cortical slow oscillations in vivo and in vitro. Sci. Rep. 8:16717. doi: 10.1038/s41598-018-34560-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Ruiz-Mejias, M., Ciria-Suarez, L., Mattia, M., and Sanchez-Vives, M. V. (2011). Slow and fast rhythms generated in the cerebral cortex of the anesthetized mouse. J. Neurophysiol. 106, 2910–2921. doi: 10.1152/jn.00440.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Ruiz-Mejias, M., Perez-Mendez, L., Ciria-Suarez, L., Martinez de Lagran, M., Gener, T., Sancristobal, B., et al. (2016). Over-expression of Dyrk1A, a Down syndrome candidate, decreases excitability and impairs gamma oscillations in the pre-frontal Cortex. J. Neurosci. 36, 3648–3659. doi: 10.1523/JNEUROSCI.2517-15.2016

CrossRef Full Text | Google Scholar

Sanchez-Vives, M. V., Massimini, M., and Mattia, M. (2017). Shaping the default activity pattern of the cortical network. Neuron 94, 993–1001. doi: 10.1016/j.neuron.2017.05.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanchez-Vives, M. V., Mattia, M., Compte, A., Perez-Zabalza, M., Winograd, M., Descalzo, V. F., et al. (2010). Inhibitory modulation of cortical up states. J. Neurophysiol. 104, 1314–1324. doi: 10.1152/jn.00178.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Sancristóbal, B., Rebollo, B., Boada, P., Sanchez-Vives, M. V., and Garcia-Ojalvo, J. (2016). Collective stochastic coherence in recurrent neuronal networks. Nat. Phys. 12, 1–8. doi: 10.1038/nphys3739

CrossRef Full Text | Google Scholar

Seabold, S., and Perktold, J. (2010). “Statsmodels: econometric and statistical modeling with python,” in 9th Python in Science Conference (Austin, TX).

Google Scholar

Seamari, Y., Narváez, J. A., Vico, F. J., Lobo, D., and Sanchez-Vives, M. V. (2007). Robust off- and online separation of intracellularly recorded up and down cortical states. PLoS ONE 2:e888. doi: 10.1371/journal.pone.0000888

PubMed Abstract | CrossRef Full Text | Google Scholar

Steriade, M., Nunez, A., and Amzica, F. (1993). A novel slow (<1 Hz) oscillation of neocortical neurons in vivo: depolarizing and hyperpolarizing components. J. Neurosci. 13, 3252–3265.

PubMed Abstract | Google Scholar

Stroh, A., Adelsberger, H., Groh, A., Rühlmann, C., Fischer, S., Schierloh, A., et al. (2013). Making waves: initiation and propagation of corticothalamic Ca2+ waves in vivo. Neuron 77:1136–1150. doi: 10.1016/j.neuron.2013.01.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Volgushev, M., Chauvette, S., Mukovski, M., and Timofeev, I. (2006). Precise long-range synchronization of activity and silence in neocortical neurons during slow-wave sleep. J. Neurosci. 26, 5665–5672. doi: 10.1523/JNEUROSCI.0279-06.2006

CrossRef Full Text | Google Scholar

Vyazovskiy, V. V., Olcese, U., Hanlon, E. C., Nir, Y., Cirelli, C., and Tononi, G. (2011). Local sleep in awake rats. Nature 472, 443–447. doi: 10.1038/nature10009

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W., Degenhart, A. D., Collinger, J. L., Vinjamuri, R., Sudre, G. P., Adelson, P. D., et al. (2009). “Human motor cortical activity recorded with micro-ecog electrodes, during individual finger movements,” in 2009 Annual International Conference of the IEEE Engineering in Medicine and Biology Society (Minneapolis, MN), 586–589. doi: 10.1109/IEMBS.2009.5333704

PubMed Abstract | CrossRef Full Text | Google Scholar