TY - JOUR
AU1 - Capone,, Cristiano
AU2 - Rebollo,, Beatriz
AU3 - Muñoz,, Alberto
AU4 - Illa,, Xavi
AU5 - , Del Giudice, Paolo
AU6 - Sanchez-Vives, Maria, V
AU7 - Mattia,, Maurizio
AB - Abstract Cortical slow oscillations (SO) of neural activity spontaneously emerge and propagate during deep sleep and anesthesia and are also expressed in isolated brain slices and cortical slabs. We lack full understanding of how SO integrate the different structural levels underlying local excitability of cell assemblies and their mutual interaction. Here, we focus on ongoing slow waves (SWs) in cortical slices reconstructed from a 16-electrode array designed to probe the neuronal activity at multiple spatial scales. In spite of the variable propagation patterns observed, we reproducibly found a smooth strip of loci leading the SW fronts, overlapping cortical layers 4 and 5, along which Up states were the longest and displayed the highest firing rate. Propagation modes were uncorrelated in time, signaling a memoryless generation of SWs. All these features could be modeled by a multimodular large-scale network of spiking neurons with a specific balance between local and intermodular connectivity. Modules work as relaxation oscillators with a weakly stable Down state and a peak of local excitability to model layers 4 and 5. These conditions allow for both optimal sensitivity to the network structure and richness of propagation modes, both of which are potential substrates for dynamic flexibility in more general contexts. computer models, cortical layers, cortical networks, in vitro recordings, oscillations, slow waves, wave propagation Introduction Propagating waves are a natural strategy to transfer information across the cortical web of neurons. For example, the processing of sensory stimuli by the brain is the result of the coordinated activity of all its cellular components, which are huge in number and distributed in space. Activation waves traveling across the cortical surface appears to be one of the possible ways in which such orchestration takes place (Ermentrout and Kleinfeld 2001; Sato et al. 2012). Traveling waves emerge in topologically organized media (Cross and Hohenberg 1993), where the nonlinear dynamics of their interacting elements contribute to generating complex and state-dependent responses to external stimulations: a way to implement sensory information processing. A wide body of literature supports this view, as activity waves elicited by sensory stimuli are invariantly found to propagate across the cortex of anesthetized animals (Bringuier et al. 1999; Jancke et al. 2004; Benucci et al. 2007; Xu et al. 2007; Stroh et al. 2013). Of course, under these conditions sensory-evoked waves are expected to be only indirectly related to the way in which information processing is performed by the same cortical network during wakefulness. Nevertheless, qualitatively similar spatiotemporal patterns of cortical activity have been measured in behaving animals (Rubino et al. 2006; Ferezou et al. 2007; Muller et al. 2014). Remarkably, these stimulus-evoked waves are not only a mere echo of the incoming input from the environment. Indeed, stimulus-evoked activations can continuously and cyclically propagate across the visual cortex even when sensory information does not change in time, as in the case of binocular rivalry induced by presenting the 2 eyes with dissimilar patterns (Lee et al. 2007). Moreover, perceptual waves spread well beyond primary sensory cortices giving rise to a self-sustained activity propagation that eventually involves other cortical areas at higher hierarchical levels, and can persist in time for fractions of seconds after stimulus presentation (Ferezou et al. 2007; Xu et al. 2007). Hence, the cortical network, whose structure does not change with varying brain states, naturally expresses activation waves even if their features vary across brain states. It is not by chance then that activation waves are spontaneously expressed by the cortical tissue also when it is kept isolated from the environment, as happens in the intact brain of many species during slow-wave (SW) sleep (Cirelli and Tononi 2008; Siegel 2008) and under deep anesthesia (Alkire et al. 2008). The almost periodic occurrence of such activation waves is a default activity pattern (Sanchez-Vives and Mattia 2014; Sanchez-Vives et al. 2017) that invariably emerges also in brain slices maintained in vitro (Sanchez-Vives and McCormick 2000; Wester and Contreras 2012). Despite the reproducibility of this SW activity, the ongoing spatiotemporal patterns composing this default activity pattern display a high degree of dynamical richness. Such richness is not a mere replica of the huge repertoire of activity modes observed in the awake brain. Indeed, in vivo SWs are endowed with a remarkable balance between stochastic and deterministic components. In this brain state, waves preferentially initiate in the frontal cortex and propagate backward to parietal/occipital cortical areas (Massimini et al. 2004; Mohajerani et al. 2010; Ruiz-Mejias et al. 2011; Stroh et al. 2013; Sheroziya and Timofeev 2014). On the other hand, these propagating wavefronts follow pathways widely distributed across the cortex at variable speeds. This variability is area-specific (Ruiz-Mejias et al. 2011), further highlighting that SW activity is shaped by the underlying structure of the cortex. The relationship between structure and SWs is corroborated also by recent findings showing that wavefront propagation is influenced by past experience (Han et al. 2008) and that activation variability does not wander between random activity patterns but rather it is the stochastic rehearsal of previously encoded sensory responses (Luczak et al. 2009). Here we investigated the balance between global nonlinear dynamics supporting SW activity, the intrinsic richness of ongoing spatiotemporal patterns, and the sensitivity to the structure of the underlying cortical medium. More specifically, we focused on the laminar structure of the cortex, known to shape in vitro SWs (Sanchez-Vives and McCormick 2000; Wester and Contreras 2012). By clustering the phase-locked slow oscillations (SO) between high-firing (Up) and quiescent (Down) states across the cortical slices, we quantitatively characterized the large variability of SO ignition sites and wave propagation modes. Reconstructed traveling waves revealed loci across the slices, reproducibly leading the propagation of Up state onsets. These loci were distributed along a smooth strip that, upon checking with anatomical reconstruction, largely overlapped with layers 4 and 5. By matching semiquantitatively large-scale simulations of a multimodular slice model to the measured SWs, we found that cell assemblies in this smooth excitable strip must work as relaxation oscillators, with finite-size fluctuations destabilizing an otherwise stable Down state. We also found that simulations reproduce the observed overlap between the excitable strip and the region with maximum firing rate and longest Up state durations, only if the intramodule and intermodule coupling strengths are optimally balanced, thereby further strengthening the relationship between spontaneous activity and network structure. In summary, novelties of the present work include both an experimental improvement, in that the electrophysiological signals are collected from 16 electrodes unevenly positioned on the slice to provide a multiscale access to the network dynamics, and a detailed characterization of the different layers’ excitability, consistent with theoretical interpretations supported by realistic simulations. Despite the relatively small number of electrodes, our method of analysis allowed an effective fine-grained description of the SW dynamics, which we a posteriori validated experimentally. Materials and Methods Slice Preparation and Histology In vitro experiments were done on 0.4 mm thick visual cortical slices (areas 17, 18, and 19) from 5 to 10 months old male ferrets (Mustela putorius furo) that were anesthetized with sodium pentobarbital (40 mg/kg) and decapitated. Their brains were removed and placed in cold (4–10 °C) cutting solution. Ferrets were cared for in accordance with the European Union guidelines on protection of vertebrates used for experimentation (2010/63/eu directive of the European Parliament and the Council of 22 September 2010) as well as approved by the local ethical committee. While cutting the slices on a vibratome, the tissue was placed in a solution in which the NaCl was substituted with sucrose. Following the transfer of the tissue to the interface-style recording chamber (Fine Sciences Tools, Foster City, CA), the slices were incubated in “traditional” slice solution containing 126 mM NaCl, 2.5 mM KCl, 2 mM MgSO4, 1 mM NaHPO4, 2 mM CaCl2, 26 mM NaHCO3, and 10 mM dextrose, aerated with 95% O2, 5% CO2 to a final pH of 7.4. After approximately 1 h, the slice solution was modified to contain the same composition except for 1 mM MgCl2, 1 mM CaCl2, and 4 mM KCl (Sanchez-Vives and McCormick 2000). Bath temperature was maintained at 34–36 °C. Extracellular multiple unit recordings were obtained with flexible arrays of 16 electrodes arranged in columns as in Figure 1A (Illa et al. 2015). The multielectrode array (MEA) covered most of the area occupied by a cortical slice. It consisted of 6 groups of electrodes positioned to record electrophysiological activity from supragranular and infragranular layers (692 μm apart) and from what should correspond to 3 different cortical columns (1500 μm apart). The unfiltered field potential (raw signal) was acquired at 5 kHz with a Multichannel System amplifier and digitized with a 1401 CED acquisition board and Spike2 software (Cambridge Electronic Design). Figure 1. View largeDownload slide Multielectrode recordings of slow Up/Down oscillations traveling across ferret cortical slices. (A) Matrix of 16 electrodes (filled circles) superimposed to the histological slice of an example experiment. Gray circles, unused channels. (B) Long-tailed distribution of the MUA measured from a representative electrode (histogram), fitted to the sum of 2 Gaussian distributions (red and blue) representing the activity during the Down and Up state, respectively. Vertical dashed lines, mean of the Gaussian distributions. (C) Histograms of the detected Up and Down state durations (red and blue, respectively). Vertical dashed lines, average durations. (D) Rastergram of MUA (bottom) centered around the detection time of Down-to-Up transitions (0 ms). Up states are sorted by duration. Top, average MUA. (E) Representative time course of the MUA simultaneously estimated from used channels, that is, those displaying Down-to-Up transitions. Horizontal dotted line, MUA threshold to detect Up state onset. Diamonds, crossing time of MUA threshold. Colors are as in panel A, and they characterize different electrodes whose numbers are shown on top of the horizontal color bar (top). Vertical dashed line, average time of Up state onsets composing the same wavefront. The time lags of Down-to-Up transitions from this wavefront centroid are represented as a horizontal color bar (top). Figure 1. View largeDownload slide Multielectrode recordings of slow Up/Down oscillations traveling across ferret cortical slices. (A) Matrix of 16 electrodes (filled circles) superimposed to the histological slice of an example experiment. Gray circles, unused channels. (B) Long-tailed distribution of the MUA measured from a representative electrode (histogram), fitted to the sum of 2 Gaussian distributions (red and blue) representing the activity during the Down and Up state, respectively. Vertical dashed lines, mean of the Gaussian distributions. (C) Histograms of the detected Up and Down state durations (red and blue, respectively). Vertical dashed lines, average durations. (D) Rastergram of MUA (bottom) centered around the detection time of Down-to-Up transitions (0 ms). Up states are sorted by duration. Top, average MUA. (E) Representative time course of the MUA simultaneously estimated from used channels, that is, those displaying Down-to-Up transitions. Horizontal dotted line, MUA threshold to detect Up state onset. Diamonds, crossing time of MUA threshold. Colors are as in panel A, and they characterize different electrodes whose numbers are shown on top of the horizontal color bar (top). Vertical dashed line, average time of Up state onsets composing the same wavefront. The time lags of Down-to-Up transitions from this wavefront centroid are represented as a horizontal color bar (top). At the end of each experiment, slices were fixed in paraformaldehyde (4%), cryoprotected in 30% sucrose solution in phosphate buffer (0.1 M, pH 7.4), frozen in dry ice and cut in a Thermo Scientific MICROM HM 450 freezing sliding microtome. The sections (50 μm thick) were mounted on glass slides and Nissl-stained, dehydrated, coverslipped, and photographed. Layer limits were drawn according to the size and density of observed cells, and the MEA position was superimposed on the basis of the identification of marks performed at the end of the recording session. To check the reliability of the interpolation procedure introduced to investigate the relationship between spontaneous electrophysiological activity and network structure in the slice regions in between the electrodes, we performed the same experiments on 5 additional cortical slices including the simultaneous recordings with an additional extra electrode placed sequentially in 5 different spots on each slice (see Supplementary Material for details). MUA Estimate and State Detection Up and Down states and multiunit activity (MUA) were estimated from the recorded raw signals (Reig et al. 2010; Sanchez-Vives et al. 2010). Briefly, the power spectra were computed from 5-ms sliding windows of the raw signal. MUA was estimated as the relative change of the power in the (0.2, 1.5) kHz frequency band. Average MUA following Up-to-Down transitions was used as reference value and set to 1. Such spectral estimate of the MUA was not affected by the electrode filtering properties, and it provided a good estimate of the relative firing rate of the pool of neurons nearby the electrode tip. We then smoothed the log(MUA) by performing a moving average with a sliding window of 80 ms. From the long-tailed histogram of log(MUA) (Fig. 1B), an optimal threshold separating Up and Down activity states was set at 4 times the standard deviation (SD) of the Down peak. To avoid biases in the identification of wavefronts, the same threshold was used for each recording in the same slice, as all having Down peaks with similar SD. This was done to avoid introducing biases across electrodes in the characterization of activation and silencing wavefronts. The threshold was also the same for detecting both Down-to-Up and Up-to-Down transitions. SW Detection and Reconstruction Activation waves were detected when a Down-to-Up transition from multiple electrodes occurred within a time window ΔT initially set to 1 s. Then ΔT was iteratively reduced by a factor of 0.75, until each detected wave contained no more than one state transition per electrode. In order to recompose waves involving a subset of columns of electrodes wrongly detected as separated waves, a minimum interwave interval (IWI) was chosen and waves occurring within the IWI were considered to compose a single wave. The initial IWI was set to half the maximum time lag between initially detected waves. Then the IWI was iteratively reduced by 25% if the collected waves included the same electrode more than once or if the number of time lags composing the wave size was larger than the number of electrodes. After each IWI reduction, wave recollection was recomputed. This iterative process ended once the IWI was neither reduced nor increased, or when the number of waves including the same electrode more than once was less than 10% of the total number of detected waves. The iteration was stopped if the number of recomposed waves in the first step of iteration was higher than a reasonable value of expected waves, that is the median number of transition for each channel. Otherwise we may have had an artefactual amplification of small waves. Waves including the same electrode more than once and waves with state transitions in less than 4 electrodes were rejected. Such constrain was dictated by the necessity to fully sample the area covered by the MEA when an activation wave crossed the whole slice. Indeed, in this case the state transitions had to be detected in each of the 3 electrode columns (see below). Having at least one detected Up state onset per column allowed having reliable predictions about ongoing network activity also relatively far from the position of the electrodes. Each wave was associated with an array Tti of relative time lags (where, t = 1, … W and i = 1, … N are the wave and the electrode indices, respectively), resulting from the difference between the detected state transition time in each electrode, and their average across electrodes. Such arrays of relative time lags in turn composed the rows of the time lag matrix (TLM) shown in Figure 2A. TLM columns were grouped into 3 sets of electrodes, each related to a different electrode column of the MEA. Partial and full waves were those involving 2 and 3 columns, respectively. Unless otherwise specified, only full waves were considered for the following analysis. In order to infer the wavefront time lags T(x,y) in the slice regions not covered by the electrodes, we interpolated without smoothing such times with a thin-plate spline. For each wave only electrodes with detected MUA state transitions were taken into account. From T(x,y) the speed V(x,y) was computed as follows: V(x,y)=1(∂T(x,y)∂x)2+(∂T(x,y)∂y)2 Figure 2. View largeDownload slide Variability of propagation modes of Up state onsets. (A) Time lag matrix of detected Up state onsets from single traveling waves of a representative experiment. Rows are time lags as the top horizontal bar in Figure 1E. Gray pixels, failed detections of Down-to-Up transition. Dashed vertical lines separate channels arranged in different vertical columns of the MEA (see panel C). Wavefronts are sorted by number of columns involved. Full and partial waves are those with 3 and 2 columns involved, respectively. Right, number of channels detecting the Up onset in each wave. (B) Time lag matrix of full waves with failed detections of Up onsets replaced by estimates from similar wavefronts (see text). Waves are pooled in 10 groups by k-means clustering, separated by horizontal dashed lines. Gray-scale bar distinguishes different clusters. (C) Time course of average wavefronts from 2 example clusters. Color code as in panels B and C. Consecutive wavefronts are separated by 50 ms. Time lags across the depicted area are interpolated from the average Up onset times across a wave cluster, relying on thin-plate splines (see text). (D) Box plots of the fraction of full and partial waves (all = partial + full) detected across experiments (n = 12 slices), together with the column of electrodes from where waves originated (side: from first and third column, i.e., from left to right and from right to left; center: from second column, i.e., propagating from the center in 2 opposite directions). (E) Box plots of the average absolute velocity of wave clusters across experiments (n = 12) along the horizontal (X) and vertical (Y) directions of the electrode array, parallel and perpendicular to the cortical surface, respectively. In each box plot, whiskers (dashed lines) are extreme values of the distribution, edges are the first and third quartiles and the central mark is the median. Where shown, circles are average values of the distribution. (F) Example of interpolated activation wavefronts obtained including or not data from an extra electrode (solid and dashed lines, respectively). Red circle represent the location of the extra electrode. (G) Scatterplot for a representative slice/electrode position of Down-to-Up transition times measured from the extra electrode (δtm) versus the interpolated ones estimated without taking into account this electrode (δti). (H) Histogram of the average absolute error estimated for each slice/electrode position (n = 24). Figure 2. View largeDownload slide Variability of propagation modes of Up state onsets. (A) Time lag matrix of detected Up state onsets from single traveling waves of a representative experiment. Rows are time lags as the top horizontal bar in Figure 1E. Gray pixels, failed detections of Down-to-Up transition. Dashed vertical lines separate channels arranged in different vertical columns of the MEA (see panel C). Wavefronts are sorted by number of columns involved. Full and partial waves are those with 3 and 2 columns involved, respectively. Right, number of channels detecting the Up onset in each wave. (B) Time lag matrix of full waves with failed detections of Up onsets replaced by estimates from similar wavefronts (see text). Waves are pooled in 10 groups by k-means clustering, separated by horizontal dashed lines. Gray-scale bar distinguishes different clusters. (C) Time course of average wavefronts from 2 example clusters. Color code as in panels B and C. Consecutive wavefronts are separated by 50 ms. Time lags across the depicted area are interpolated from the average Up onset times across a wave cluster, relying on thin-plate splines (see text). (D) Box plots of the fraction of full and partial waves (all = partial + full) detected across experiments (n = 12 slices), together with the column of electrodes from where waves originated (side: from first and third column, i.e., from left to right and from right to left; center: from second column, i.e., propagating from the center in 2 opposite directions). (E) Box plots of the average absolute velocity of wave clusters across experiments (n = 12) along the horizontal (X) and vertical (Y) directions of the electrode array, parallel and perpendicular to the cortical surface, respectively. In each box plot, whiskers (dashed lines) are extreme values of the distribution, edges are the first and third quartiles and the central mark is the median. Where shown, circles are average values of the distribution. (F) Example of interpolated activation wavefronts obtained including or not data from an extra electrode (solid and dashed lines, respectively). Red circle represent the location of the extra electrode. (G) Scatterplot for a representative slice/electrode position of Down-to-Up transition times measured from the extra electrode (δtm) versus the interpolated ones estimated without taking into account this electrode (δti). (H) Histogram of the average absolute error estimated for each slice/electrode position (n = 24). SW Feature Characterization We collected wave vectors describing full waves and we clustered them using a k-means algorithm to separate different modes of propagations characterized by different origin, speed and direction. In order to perform k-means clustering, we needed complete transition time vectors. We filled each empty entry in the vectors, corresponding to a nonrecorded transition during a wave, with the average transition time recorded by the same electrode in the 5 most similar waves. The median error in the time lag of the “filled” elements estimated in this way was 30.0 ± 8.0 ms, that is, one order of magnitude smaller than the characteristic time scale of wave propagation, which typically took more than 600 ms to cross the whole slice (see Supplementary Material). For each recording we grouped waves into 10 clusters. We rejected clusters, and relative waves, with a mean speed lower than 2 mm/s. This threshold value has been chosen starting from the observation that the distribution of the speed over all experiments (not shown) is well fitted by the sum of 2 distributions. These distributions peak just above and below such threshold, respectively. Of these 2 components, the one corresponding to the highest speeds is related to propagating modes described in the Results. The other component with speed always smaller than 0.6 mm/s seemed not to be associated with traveling activation wavefronts and for this reason we did not analyze them. On average, 3 clusters for each experiment were rejected and the accepted ones were composed of 15 waves. For each wave from an accepted cluster, we computed the horizontal (Vx) and vertical (Vy) components of the wavefront velocity from T(x,y) as follows: {Vx=1∂T(x,y)∂xVy=1∂T(x,y)∂y. In order to characterize the spatiotemporal complexity of the identified waves we considered only the full waves described by the time series {Tti} and evaluated the autocorrelation ACn of the wavefront shapes at wave lag n as follows: ACn=1W−n∑i=1N∑t=1W−n(Tt+ni−Ti¯)(Tti−Ti¯), with N being the number of electrodes, W the number of detected waves and Ti¯=1/W∑t=1WTti ⁠. Note that lag n did not explicitly relate to the time but rather to the ordinal number of the detected wave. In order to avoid spurious effects leading to possible lack of correlation due to the presence of partial waves between full waves, we took into account only the full waves with an IWI lower than 1.5 of the average Up/Down cycle measured for each specific experiment/slice. In each experiment only the ACn measured from at least N = 20 couples of waves were taken into account to compute averages (Fig. 3D, E) in order to have a reliable enough estimate of the autocorrelation at lag n. As a result, the experiments/slices contributing to the median and SD of ACn were 10, 9, and 4 out of 12 for lag 1, 2, and 3, respectively. Figure 3. View largeDownload slide Past modes of propagation do not influence next Up wavefronts. (A) First 2 principal components (PC1 and PC2) of the time lag matrix in Figure 2B. Each circle is an Up onset wavefront (an array of time lags), whose gray shading indicates the wave cluster it belongs to (gray levels as in Fig. 2B). (B) Box plot of variance explained by PC1 and by the first 2 PCs across all experiments. (C) Time course of the phase ϕ in the (PC1, PC2) plane of full (filled circles) and partial (dots) waves. Gray filling and experiment as in A. (D) Autocorrelation ACn of the time lag arrays (wavefront shapes) of full waves at lag n (experiments with at least 20 couples of waves are n = 10, 9, and 4 out of 12 for lag 1, 2, and 3, respectively; see Materials and Methods). Red line, average ACn for a specific outlier experiment. Black line and gray shading, median and SEM of ACn across experiments not including the (red) outlier, respectively. AC1 of the outlier is significantly different from the median of the other n = 9 slides (P < 0.001, Wilcoxon signed-rank test). (E) Scatterplot of AC1 wavefront shapes versus the average duration of the Up states from a reference channel (the one with the highest MUA). Experiments included as in D for wave lag n = 1. Figure 3. View largeDownload slide Past modes of propagation do not influence next Up wavefronts. (A) First 2 principal components (PC1 and PC2) of the time lag matrix in Figure 2B. Each circle is an Up onset wavefront (an array of time lags), whose gray shading indicates the wave cluster it belongs to (gray levels as in Fig. 2B). (B) Box plot of variance explained by PC1 and by the first 2 PCs across all experiments. (C) Time course of the phase ϕ in the (PC1, PC2) plane of full (filled circles) and partial (dots) waves. Gray filling and experiment as in A. (D) Autocorrelation ACn of the time lag arrays (wavefront shapes) of full waves at lag n (experiments with at least 20 couples of waves are n = 10, 9, and 4 out of 12 for lag 1, 2, and 3, respectively; see Materials and Methods). Red line, average ACn for a specific outlier experiment. Black line and gray shading, median and SEM of ACn across experiments not including the (red) outlier, respectively. AC1 of the outlier is significantly different from the median of the other n = 9 slides (P < 0.001, Wilcoxon signed-rank test). (E) Scatterplot of AC1 wavefront shapes versus the average duration of the Up states from a reference channel (the one with the highest MUA). Experiments included as in D for wave lag n = 1. Quantitative Characterization of SWs For each wave an early propagation strip (EPS) was identified as the set of points in which the wavefront tip at different times first occurred in the horizontal direction. In other words, EPS(x) was the vertical position y where the time-lag surface T(x,y) displayed a minimum at fixed x: EPS(x)=argminT(x,y) ⁠. EPSs were identified for all waves from accepted clusters. All the strips were convolved with a Gaussian kernel along the vertical direction and averaged to obtain a smoothed probability density function p(x,y). A cut-off to this density was further adopted by setting p(x,y) = 0 for p(x,y) < 0.75 pmax(x), where pmax(x) was the maximum of p(x,y) at a given x. After renormalizing this cut density, the average EPS(x)=∫yp(x,y)dy and its SD ∫y2p(x,y)dy−EPS(x)2 were finally computed. For each electrode both the average Up state duration, and the maximum MUA during the Up state were computed for each wave in an accepted cluster. As for the time-lag surface T(x,y), these values were interpolated by a thin-plate spline allowing to obtain TUP(x,y) and MUA(x,y), respectively. From these surfaces, the curves where maximum Up duration and MUA can be found at different horizontal positions x (⁠ argmaxyTUp(x,y) and argmaxyMUA(x,y) ⁠, respectively) were computed (Fig. 4C–E). The distances between these curves and EPS(x) at the horizontal locations of the electrodes in the array (n = 8) were computed to estimate histograms in Figure 4F, G. Figure 4. View largeDownload slide Global Up wavefront properties match local dynamical features of SO. (A) Two examples of Up wavefronts propagating in opposite directions as in Figure 2C, from another experiment. Gray curves with arrows, early propagation strip (EPS). (B) Density of EPS from all detected waves. Brown thick line, weighted average of EPS density along the X direction. Gray lines, EPSs of each cluster of traveling waves computed as in Figure 2. (C) Average Up state duration (TUp) interpolated across the slice area covered by the matrix of electrodes (circles). Green thick line, maximum Up duration found at different X. (D) Average maximum MUA during Up state interpolated as in panel C. Cyan thick line, maximum MUA found at different X. (E) Superimposed weighted average of EPS, maximum Up state duration and MUA from panels B–D (brown, green and blue curves, respectively). (F, G) Histograms of the vertical distance across all experiments (n = 12) between EPS and line of maximal Up duration (F) or MUA (G). (H) EPS estimated with and without the extra electrode (dashed and solid lines, respectively). Two examples of different positions of the extra electrode are shown (top and bottom panels, see Supplementary Material for details). (I) Distributions of the average vertical error in interpolating the EPS across slices and extra electrode locations (n = 24). (J) Correlation between speed difference and dTUp/dX (the change per unit length of Up state duration along the EPS) in the same experiment in panels A–E. Speed difference is between median velocities of waves propagating in opposite directions (see Materials and Methods for details). (K) Histogram of the Pearson correlation, ρ, as in panel H for all experiments (n = 12). To the right of the red dashed line, correlations ρ are significant (P < 0.05). Figure 4. View largeDownload slide Global Up wavefront properties match local dynamical features of SO. (A) Two examples of Up wavefronts propagating in opposite directions as in Figure 2C, from another experiment. Gray curves with arrows, early propagation strip (EPS). (B) Density of EPS from all detected waves. Brown thick line, weighted average of EPS density along the X direction. Gray lines, EPSs of each cluster of traveling waves computed as in Figure 2. (C) Average Up state duration (TUp) interpolated across the slice area covered by the matrix of electrodes (circles). Green thick line, maximum Up duration found at different X. (D) Average maximum MUA during Up state interpolated as in panel C. Cyan thick line, maximum MUA found at different X. (E) Superimposed weighted average of EPS, maximum Up state duration and MUA from panels B–D (brown, green and blue curves, respectively). (F, G) Histograms of the vertical distance across all experiments (n = 12) between EPS and line of maximal Up duration (F) or MUA (G). (H) EPS estimated with and without the extra electrode (dashed and solid lines, respectively). Two examples of different positions of the extra electrode are shown (top and bottom panels, see Supplementary Material for details). (I) Distributions of the average vertical error in interpolating the EPS across slices and extra electrode locations (n = 24). (J) Correlation between speed difference and dTUp/dX (the change per unit length of Up state duration along the EPS) in the same experiment in panels A–E. Speed difference is between median velocities of waves propagating in opposite directions (see Materials and Methods for details). (K) Histogram of the Pearson correlation, ρ, as in panel H for all experiments (n = 12). To the right of the red dashed line, correlations ρ are significant (P < 0.05). We repeated the same analysis in an additional set of similar experiments (5 slices) using an extra electrode placed sequentially in 5 different positions far from other electrodes of the MEA (see Fig. 2F and Supplementary Material), collecting a total amount of 25 recordings. One of these recordings was discarded due to data acquisition instability, resulting in n = 24 samples in the data set. For each SW, the propagating wavefront was computed relying on the thin-plate spline-based interpolation procedure described above. Such wavefronts were estimated both with and without the extra electrode (see Fig. 2F and Supplementary Materials) to check the reliability of the method by inspecting the resulting differences, that is, the measure of their misalignment in time. The same approach was used to estimate the measure error due to relying or not on the interpolation procedure both in the identification of the EPS, and in the strips in the slices where the Up state duration and the MUA were maximal (see Supplementary Material for further details). To study the wave speed dependence of propagation modes (Fig. 4H, I), we computed separately the average speeds VL(S) and VR(S) along the average EPS curvilinear abscissa S in the leftward and rightward direction, respectively. We finally computed an asymmetry index as the speed difference ΔV(S) = VR(S) − VL(S). Similarly, we also computed the gradient of Up state durations dTUp(S)/dS ⁠. All these measures were sub-sampled in order to avoid spurious correlation effects. Decimated points were selected at steps ΔS = 0,21 mm, which corresponds to a length always greater than the decay constant of the autocorrelation of TUp(S). To test the overlap between EPS and some specific laminar structures of the cortical slices we considered the region where the EPS density p(x,y) > 0.05. Thus we measured the portion of its area included in each layer (Fig. 7). Where not specified averages across the n = 12 experiments are reported as mean ± SD. All off-line analyses were performed using MATLAB (The MathWorks, Natick, MA). Theory, Models, and Simulations We used both network models of integrate-and-fire (IF) neurons and a minimal rate model to set the simulation parameters. The in silico slice is a multimodular interconnected network of spiking neurons. Similarly to Mattia et al. (2013), cortical modules were composed of 350 leaky integrate-and-fire (LIF) excitatory (E, 80%), and inhibitory (I, 20%) neurons with spike frequency adaptation (SFA). Membrane potential V(t) of LIF neurons evolved according to dV(t)dt=−V(t)τα+Isyn(t)−IAHP(t) ⁠, where Isyn(t) was the synaptic incoming current and τα was the membrane decay constant (τE = 20 ms and τI = 10 ms). Spikes were emitted when V(t) crossed a 20 mV threshold, after which a 15 mV reset potential was set for an absolute refractory period of 10 ms (6 ms) for excitatory (inhibitory) neurons. IAHP(t) was the activity-dependent after-hyperpolarizing potassium current acting as a fatigue mechanism modeling SFA for excitatory neurons: dIAHP(t)dt=−IAHP(t)τAHP+gAHP∑kδ(t−tk) ⁠, with τAHP = 1 s and gAHP = 0.06 mV/s. The δ(t − tk) were the spikes emitted by the neuron. Synaptic transmission was instantaneous, and Isyn(t)=∑jkJjδ(t−tkj+δj)+∑kJext,kδ(t−text,k) ⁠. The k-th spike emitted at t = tjk by the local presynaptic neuron j affected the postsynaptic membrane potential with synaptic efficacy Jj after a transmission delay δj. Synaptic efficacies were randomly chosen from a Gaussian distribution with mean Jαβ and SD ΔJαβ depending on the type of presynaptic (β = {E, I}) and postsynaptic (α = {E, I}) neurons. We initially set JEE = 0.73 mV, JIE = 0.95 mV, JEI = 2.55 mV, and JII = 2.55 mV, whereas ΔJαβ = 0.25 Jαβ for any α and β. Intramodular connectivity was cloc=80% (probability to have synaptic contact between 2 neurons of the same module) for any α and β, unless otherwise specified. Excitatory and inhibitory synaptic transmission delays were drawn from an exponential distribution with average 13.3 and 3.3 ms, respectively. Each neuron also received spikes from neurons outside the cortical module it belonged to, collectively modeled as a Poisson process with average spike frequency νext chosen to have under mean-field approximation average firing rates νE = 0.7 Hz and νI = 1.75 Hz. Synaptic efficacies Jext,k were randomly chosen from a Gaussian distribution with the same moments as the local excitatory connections. The neurons in the excitatory population were grouped into 2 subpopulations: 25% in a foreground population (F) displaying Up/Down SO, and the remaining 75% in a background population (B) always firing at a relatively low-firing rate (Fig. 5B). Oscillatory dynamics in the foreground population was obtained by potentiating synaptic self-excitation JFF=w+JEE ⁠, with a potentiation level w+=1.8 ⁠. To keep νF = 0.7 Hz and νB = 0.89 Hz fixed, JBF=JFB were decreased accordingly. Figure 5. View largeDownload slide Slow-wave activity in a spiking neuronal network model of a cortical slice. (A) Simulated networks are organized as a lattice of interacting modules (small black circles). Top, each module is composed of 2 pools of excitatory (E) and inhibitory (I) leaky integrate-and-fire (LIF) neurons. Excitatory neurons receive an additional state-dependent self-inhibition (magenta link) modeling spike frequency adaptation (SFA). (B) Example of slow oscillation in a simulated module. Top raster plot displays a subset of emitted spikes. Recurrent connectivity between excitatory neurons is shaped in order to have only a subset of them firing during Up states (depicted in panel A, top as a dashed circle in the excitatory pool). Center, firing rates ν(t) of different pools, and fatigue level c(t) proportional to self-inhibition inducing SFA. Bottom, raster of spikes emitted by neurons from different excitatory modules along the X axis of the modeled slice. (C) Nullclines (curves) and fixed-points (circles) of single module dynamics in the (ν,c) plane, under mean-field approximation. Green curve, example trajectory of the module in the (ν,c) plane, with initial condition on the right of the ν-nullcline. (D) Mean-field nullclines for the same module in (C) receiving 2 different inputs obtained by increasing the external rate of spikes by 7% (dark gray) and 20% (light gray). Superimposed are example trajectories for these 2 conditions, colored to encode time (see top color bar). (E) MUA in logarithmic scale estimated from the firing rate ν(t) of the example cortical module in (B). Inset, histogram of Up and Down state durations detected from the modeled MUA with the same analytical tools used for the in vitro recordings, as in Figure 1C. Figure 5. View largeDownload slide Slow-wave activity in a spiking neuronal network model of a cortical slice. (A) Simulated networks are organized as a lattice of interacting modules (small black circles). Top, each module is composed of 2 pools of excitatory (E) and inhibitory (I) leaky integrate-and-fire (LIF) neurons. Excitatory neurons receive an additional state-dependent self-inhibition (magenta link) modeling spike frequency adaptation (SFA). (B) Example of slow oscillation in a simulated module. Top raster plot displays a subset of emitted spikes. Recurrent connectivity between excitatory neurons is shaped in order to have only a subset of them firing during Up states (depicted in panel A, top as a dashed circle in the excitatory pool). Center, firing rates ν(t) of different pools, and fatigue level c(t) proportional to self-inhibition inducing SFA. Bottom, raster of spikes emitted by neurons from different excitatory modules along the X axis of the modeled slice. (C) Nullclines (curves) and fixed-points (circles) of single module dynamics in the (ν,c) plane, under mean-field approximation. Green curve, example trajectory of the module in the (ν,c) plane, with initial condition on the right of the ν-nullcline. (D) Mean-field nullclines for the same module in (C) receiving 2 different inputs obtained by increasing the external rate of spikes by 7% (dark gray) and 20% (light gray). Superimposed are example trajectories for these 2 conditions, colored to encode time (see top color bar). (E) MUA in logarithmic scale estimated from the firing rate ν(t) of the example cortical module in (B). Inset, histogram of Up and Down state durations detected from the modeled MUA with the same analytical tools used for the in vitro recordings, as in Figure 1C. In the in silico slice, modules were spatially arranged on a 2-dimensional grid composed of H × B sites (generally H = 9, B = 20). The intermodular connectivity cglob was set only between pairs of subpopulations with presynaptic excitatory neurons, modeling corticocortical long-range synaptic connections. When not otherwise specified, cglob=0.2 ⁠. Intermodule spike delays were sampled from an exponential distribution with a 21 ms average. Connectivity decayed following a Gaussian function of the distance with SD λ = 0.7 IMD (intermodules distance). The resulting networks had a total size of Nneu = 63 000 neurons and Nsyn = 20 millions of synapses. The modulation of in silico slice excitability (Fig. 7D) was performed by changing the spike rate νExt of external spikes (νExt was multiplied by a factor ranging from 0.89 to 1.00). The layered structure of the in silico slices (Fig. 6) was incorporated by creating an excitable strip of modules with higher connectivity, such that modules far from the strip had the lower values of connectivity: the intermodular and intramodular connectivity (⁠ cglob and cloc ⁠, respectively) were reduced by a factor 0.8. This factor was linearly increased as the distance from the maximum connectivity strip (MCS) where cglob and cloc reach their maximum value, independently chosen for intermodular and intramodular connectivity. Numerical simulations of in silico slices were performed relying on an event-based approach (Mattia and Del Giudice 2000), spanning network times of 100 s. Figure 6. View largeDownload slide Spiking cortical model with a nonhomogeneous structure. (A) Simulated networks have the same multimodular arrangement as in Figure 5A, with the only difference that modules have a different degree of intramodular (local) and intermodular (global) connectivity (see sketch above). Darkest circles (thickest synaptic link) are those with the highest connectivity. (B) Nullclines (curves) and fixed-points (circles) of single module dynamics in the (ν,c) plane, varying the connectivity level under mean-field approximation. Gray levels code for connectivity, same as in (A). (C) Energy landscapes (top) around the fixed-points varying connectivity as in (B). Bottom, a (rotated) zoom of the (ν,c) plane in (B). (D) Snapshots at different times of the propagation of an Up wavefront across the cortical slice model. White dashed line is the maximum connectivity strip (MCS) corresponding to the region of the most excitable modules as in (A). Color coded is the MUA from simulations. (E) Two examples of Up wavefronts propagating in opposite directions in the same simulation, as in Figure 4A. (F) Comparison between estimated EPS (brown), and strips of maximum Up state duration (green) and MUA (cyan), together with the exact position of the MCS set in the modeled slice (black dashed line). Figure 6. View largeDownload slide Spiking cortical model with a nonhomogeneous structure. (A) Simulated networks have the same multimodular arrangement as in Figure 5A, with the only difference that modules have a different degree of intramodular (local) and intermodular (global) connectivity (see sketch above). Darkest circles (thickest synaptic link) are those with the highest connectivity. (B) Nullclines (curves) and fixed-points (circles) of single module dynamics in the (ν,c) plane, varying the connectivity level under mean-field approximation. Gray levels code for connectivity, same as in (A). (C) Energy landscapes (top) around the fixed-points varying connectivity as in (B). Bottom, a (rotated) zoom of the (ν,c) plane in (B). (D) Snapshots at different times of the propagation of an Up wavefront across the cortical slice model. White dashed line is the maximum connectivity strip (MCS) corresponding to the region of the most excitable modules as in (A). Color coded is the MUA from simulations. (E) Two examples of Up wavefronts propagating in opposite directions in the same simulation, as in Figure 4A. (F) Comparison between estimated EPS (brown), and strips of maximum Up state duration (green) and MUA (cyan), together with the exact position of the MCS set in the modeled slice (black dashed line). In order to use the same analytical approach developed for in vitro recordings in our in silico slices, we estimated the instantaneous firing rates ν(t) of the simulated population of excitatory neurons in each module, and used them as MUA after adding a relatively weak white noise to emulate unspecific background fluctuations observed in the experimental recordings. The optimal choice for maximum cglob and cloc (Fig. 7) was found by performing different simulations with different values of these parameters and looking for the minimum value of the identification error between the EPS and the MCS. The identification error was estimated as the average vertical distance between the 2 strips. Simulations were evaluated on a grid of 9 × 9 points in the (max. cloc ⁠, max. cglob ⁠) plane, ranging from 0.8 to 1.1. For each point 4 simulations were performed. The resulting surfaces were interpolated with a thin-plate spline with smoothing parameter p = 0.95. Figure 7. View largeDownload slide In silico optimal balance between global and local connectivity. (A) Speed of traveling wavefronts in cortical slice models varying local (⁠ cloc ⁠) and global (⁠ cglob ⁠) connectivity in the MCS. Connectivity of modules outside the MCS are adjusted accordingly (see text for details). Dashed line, set of network parameters yielding a constant speed of waves of 9 mm/s compatible with experiments. (B) Identification error of MCS from EPS for the same parameter set as in (A). (C) Comparison between EPS and MCS for 3 different example slice models corresponding to the filled circles in (B). (D) Average Up state duration versus autocorrelation AC1 between time lag arrays of consecutive Up wavefronts, computed as in Figure 3. Each circle is a simulated network with a different external excitatory input Iext to each excitatory neuron, which modulates the module excitability. Filling color codes for the changes ΔIext of the reference external input Iext in the network of panels (A, B) with cloc=cglob=1 ⁠. (E) Autocorrelation ACn of the time lag arrays at different wave lags, as in Figure 3D. Black line and gray shading depicts the average and SD over all simulations, respectively. Red line, ACn of an outlier slice model with the highest Iext ⁠. Figure 7. View largeDownload slide In silico optimal balance between global and local connectivity. (A) Speed of traveling wavefronts in cortical slice models varying local (⁠ cloc ⁠) and global (⁠ cglob ⁠) connectivity in the MCS. Connectivity of modules outside the MCS are adjusted accordingly (see text for details). Dashed line, set of network parameters yielding a constant speed of waves of 9 mm/s compatible with experiments. (B) Identification error of MCS from EPS for the same parameter set as in (A). (C) Comparison between EPS and MCS for 3 different example slice models corresponding to the filled circles in (B). (D) Average Up state duration versus autocorrelation AC1 between time lag arrays of consecutive Up wavefronts, computed as in Figure 3. Each circle is a simulated network with a different external excitatory input Iext to each excitatory neuron, which modulates the module excitability. Filling color codes for the changes ΔIext of the reference external input Iext in the network of panels (A, B) with cloc=cglob=1 ⁠. (E) Autocorrelation ACn of the time lag arrays at different wave lags, as in Figure 3D. Black line and gray shading depicts the average and SD over all simulations, respectively. Red line, ACn of an outlier slice model with the highest Iext ⁠. The dynamical features of each module were set relying on mean-field theory. The mean-field (minimal rate) model had dynamics determined by the gain function Φ (Wilson and Cowan 1972) as follows: {dνidt=ϕi(ν→,c→)−νiτdcidt=−ciτc+giνi where ν→={νi} and c→={ci} ⁠, with i = {F,B,I}, pointing out the excitatory foreground, excitatory background and inhibitory populations composing the single model, respectively. Since the single module included multiple interacting neural populations, the mean-field dynamics would be described by a multidimensional gain function. Our gain function Φ was computed as an “effective” gain function ϕ˜, along the lines previously proposed (Mascaro and Amit 1999). This approximation allows reducing the multidimensional mean-field problem (Amit and Brunel 1997) to a 2-dimensional one corresponding to the dynamics of the firing rate (ν) and of the fatigue level (c) of the foreground population of interest (Mattia and Sanchez-Vives 2012): {dνdt=ϕ˜(ν,c)−ντdcdt=−cτc+gν For this reduced dynamics of a single module, nullclines (⁠ ν̇=0 and ċ=0 ⁠) can be derived whose intersections represent fixed points. Modules were set to have a weakly stable Down state (meaning that the low-firing rate intersection was on the lower branch of the ν nullcline, rather close to the point where this branch ended, see Fig. 5C). Energy landscapes near the fixed point were defined by the integral of the function (ϕ˜(ν,c)−ν) ⁠, assuming for c its value in the same fixed point. Results We probed spontaneous Up/Down SO traveling across 12 cortical slices from 5 ferrets by measuring the extracellular raw signals simultaneously recorded from a matrix of 16 electrodes with a multiscale arrangement (Fig. 1A, see Materials and Methods). MUA of each electrode was estimated as the relative power change of high-frequency components of the raw signals in sliding windows of 5 ms (Reig et al. 2010; Sanchez-Vives et al. 2010). MUA estimated in this way represents the spike rate emitted by the population of neurons nearby an electrode (Mattia et al. 2013). The large majority of the recording channels displayed a long-tailed distribution of the logarithmically scaled MUA (Fig. 1B, 80% of the 156 used electrodes) with high- and low-firing rate peaks corresponding to Up and Down states, respectively. For each slice, an optimal MUA threshold was selected relying on the Down peak SD, to detect transition times between Up and Down states (see Materials and Methods). From these we measured the SO statistical features including state duration distributions of representative channels (Fig. 1C) and the average oscillation frequency of 0.31 ± 0.12 Hz (n = 12 slices), similarly to what has been previously found (Sanchez-Vives and McCormick 2000). MUA time course around state transitions (Fig. 1D) was qualitatively similar to that of conventional extracellular recordings (Reig et al. 2010; Sanchez-Vives et al. 2010), further proving the reliability of the new surface electrodes employed here in measuring MUA from cortical slices. Diverse Modes of Wave Propagation Across Cortical Slices Spontaneous Up state onsets at the single channel level did not occur as isolated events. As shown in Figure 1E, Up states traveled across the cortical slice (Sanchez-Vives and McCormick 2000; Sanchez-Vives et al. 2010; Wester and Contreras 2012), giving rise to waves identifiable by the set of time lags between each MUA onset occurring across channels (Fig. 1E, top, see Materials and Methods). Wavefronts were estimated by interpolating between such time lags as explained below. The average time of such local activations was used as the reference time for each wave. In each experiment we found a wide variability of the propagation patterns (Fig. 2A), the large majority of which involved at least 2 of the 3 vertical columns of channels composing the MEA (96.7 ± 3.4% of detected waves, average across all n = 12 slices). The remaining 3.3% of events involving only one channel column were not considered as activation waves and were thus discarded from the analysis. In order to uncover whether such diversity hid some spatiotemporal organization of spontaneous SO, we analyzed with a k-means clustering the subset of full waves that propagated across all the 3 channel columns of the MEA—3 mm apart—(which constituted an average of 80 ± 19% of the detected waves per slice, n = 12; Fig. 2B, see Materials and Methods for details). Using a fixed number of 10 wave clusters, for each we estimated the time course of average activation wavefronts across the rectangular field covered by the MEA by performing a spatial interpolation of the averaged time lags per channel (Fig. 2C, see Materials and Methods for details). From these wave profiles it was apparent that not only different directions of the horizontal propagation spontaneously occurred, but also that across experiments there was almost no preference in the ignition site (Fig. 2D): on average 59 ± 22% of the waves appeared first from the lateral columns of electrodes, while 34 ± 25% originated close to the central column of the MEA (n = 12). Another piece of evidence that added complexity to this phenomenon was the fact that more than one ignition site could contribute to the slice activation, as it is apparent from the first wave cluster in Fig. 2B (bottom) where an early Up state onset (blue matrix elements) occurred in the 2 most lateral electrode columns. Besides, whatever the mode of propagation, average wavefronts always displayed a shape that was squeezed along the vertical axis, perpendicular to the pia of the cortical slice. This highlighted a slower velocity along the horizontal (lateral) direction with respect to the vertical (deep) one (Fig. 2E): this difference was confirmed at the population level, since the average velocities were 5.8 ± 1.8 and 8.9 ± 3.4 mm/s, respectively (n = 12), compatible with previous measures in similar experiments (Sanchez-Vives and McCormick 2000; Sanchez-Vives et al. 2010). As mentioned above, these results rely on an interpolation procedure to obtain estimates from the whole area covered by the array. We tested the precision of such interpolation by performing an additional set of experiments including an extra electrode placed in different locations in between the MEA electrodes to check for the matching between the actual wavefronts and the interpolated ones (Fig. 2F, see Materials and Methods, and Supplementary Material). Up state onset times measured including the extra electrode were directly compared with the interpolated ones and yielded a high degree of correlation (Fig. 2G). At the population level, the average of the absolute error between the real and the extrapolated time lags resulted, across recordings, in a rather narrow distribution (Fig. 2H), with a mean error of 51.0 ± 6.5 ms (n = 24 from different slices and extra electrode positions) which corresponds to a shift of approximately 10% of the time taken by a wave to travel across the whole slice. This is a result per se, as it confirms the smoothness of SW propagation across cortical slices although the dynamics underlying their generation is expected to be intimately nonlinear. Memoryless Generation of Spontaneous Waves Traveling wavefronts of Up state onsets belong to the realm of phenomena occurring in noisy excitable media, known to spontaneously express spatiotemporal patterns (Cross and Hohenberg 1993; Sagués et al. 2007). Experimental evidence supports this view, as local excitation in space elicits propagating self-sustained wavefronts of in vitro activity (Wu et al. 1999; Wester and Contreras 2012). This phenomenon in slices displays a low degree of complexity, as spontaneous propagation modes are unidirectional (Golomb and Amitai 1997; Wu et al. 1999; Sanchez-Vives and McCormick 2000; Wester and Contreras 2012), that is, wavefronts mainly propagate in a direction parallel to the cortical surface. Here, we characterized the dimensionality of the Up wavefronts by performing a principal component analysis (PCA) of the time lag matrices (Fig. 3A). Detected spatiotemporal patterns had low dimensionality (maximal dimensionality was 15 as the electrodes in the MEA were 16 and we set to 0 the average time lag per wave), as the first 2 principal components (PC) explained a large part of the time lag variance (Fig. 3B, average variance explained 82% ± 16%, n = 12). This is consistent with previous reports on propagation patterns measured with optical imaging (Wu et al. 1999; Wester and Contreras 2012). Moreover, in the (PC1, PC2) plane where each full wave corresponded to a circle (Fig. 3A), 2 main clusters were often recognizable, representing the preferred modes of propagations, hence further decreasing the degree of complexity of spontaneous activity. This indeed corresponded to a rather stereotyped activation of the slices when they belonged to the same mode of propagation (wave cluster; see Supplementary Material). In the specific example shown in Figure 3A, the 2 preferred modes corresponded to those shown in Figure 2C for the same experiment, although other 2 smaller clusters were apparent in the first and third quadrant representing waves with different number and position of ignition sites. To uncover possible complexity in the time domain, we inspected the degree of similarity between wavefronts occurring at different distances in time (see Materials and Methods). Visual inspection of the time course of ϕ (Fig. 3C), that is, the angle (phase) associated in the plane (PC1, PC2) to each circle/wave in Figure 3A, highlighted a rather irregular nature in the selection of the next wave onset between the 2 preferred propagation modes shown in Figures 3A and 2C. This qualitative picture was confirmed at the population level; the average autocorrelation ACn of the time lag arrays characterizing the waves at lag n (see Materials and Methods) was close to 0 (Fig. 3D). Also at the single slice level, wave selection was memoryless (ACn = 0), with only one exception (Fig. 3D, red plot). In this experiment we found a significantly positive autocorrelation (AC1 = 0.04; Wilcoxon signed-rank test, P < 0.01). Intriguingly, there we measured the longest average Up state duration (540 ± 90 ms, more than 200 ms longer than the average duration of 310 ± 110 ms across experiments), supporting the hypothesis that this cortical slice was in a relatively more excitable state with respect to the others. The correlation between AC1 and Up state duration (Fig. 3E; Pearson correlation, ρ = 0.80; P < 0.01) confirmed that the memory degree in the next wave selection was tightly related to the excitability of the network, which in turn is expected to be associated with longer high-firing states. From this perspective, spontaneous waves in our basal excitability levels and frequencies between 0.11 and 0.40 Hz did not take into account past events as expected in a subexcitable state (Sagués et al. 2007). Under this condition, endogenous fluctuations of the activity in the Down state were expected to elicit in random locations the onset of different Up wavefronts explaining the memoryless renewal origin of spontaneous waves in our in vitro model of SO. Structure has a Role in Shaping Up State and Onset Wavefronts If SO in cortical slices arose as in a subexcitable medium, we might expect that activity reverberation did not completely overwhelm the underlying intracortical connectivity structure. Under this hypothesis, spontaneous activity should be (at least partially) shaped by the structure and excitability distribution of the cortical network, such that wave propagation and SO in single electrodes should retain and share some of this information. To investigate such possible interplay between structure and spontaneous activity, we looked for features of the excitation waves that were preserved across different modes of propagation. More specifically, following each traveling Up onset wavefront we identified an EPS (see Materials and Methods), corresponding to the most advanced tip of Up onset times along the X (horizontal) direction of the MEA (Fig. 4A, gray lines with arrows). This would likely represent the most excitable part of the slice, where neuronal pools reacted first to the synaptic input provided by nearby active cortical regions. Intriguingly, EPS extracted from excitation waves were remarkably similar in a given experiment, irrespective of their mode of propagation (Fig. 4B), further supporting the hypothesis that the underlying intracortical connectivity establishes preferred routes across the slice, thus modeling the shape of wavefronts independently of ignition site and direction of propagation. A natural hypothesis is that EPS corresponds to sites of maximal excitability, which in turn led us to predict that also local SO properties would be different along the EPS. To test this hypothesis, we computed the average Up state durations TUp across the area covered by the MEA, by interpolating the TUp measured from the different MEA channels (Fig. 4C), similarly to what was done for working out the traveling wavefronts (see Materials and Methods). Another mean feature of SO related to the local excitability of neuronal assemblies was the maximum MUA during Up states (Fig. 4D), which was computed as above by resorting to interpolation. The most excitable pools were expected to have the longest Up duration and maximum MUA, corresponding to a maximally stable Up state (Mattia and Sanchez-Vives 2012). We found these pools at sites composing continuous strips (green and cyan curves in Fig. 4C, D, respectively), which in turn displayed a remarkable overlap with the EPS estimated from wave propagation (Fig. 4E). These overlaps of excitable strips estimated from different features of SW activity in the slice were apparent across all experiments, as the vertical distance between EPS (YEPS) and the position of maximum MUA (Ymax. MUA in Fig. 4F) and maximum TUp (Ymax. Tup in Fig. 4G) were not significantly different from 0 (average distances were 47 ± 140 and 10 ± 90 μm, respectively; Wilcoxon sign-rank test, P > 0.05; see Materials and Methods for details). As also the EPS and the strips of maximum MUA and TUp relied on the interpolation method, we investigated the potentially introduced error by measuring the differences when the extra electrode was or was not taken into account (see Materials and Methods, and Supplementary Material). When computed as an average across all waves, EPSs estimated with and without the extra electrode (solid and dashed line in Fig. 4H, respectively) widely overlapped. The average vertical distance had a narrow distribution across slices and extra electrode position (Fig. 4I), with a relatively small mean EPS location error of 39.5 ± 7.3 μm, corresponding to approximately 5% of the MEA vertical extension. Similar results were obtained for the strips of maximum MUA and TUp (see Supplementary Material). Finally, it is important to remark that all population analyses were computed at the horizontal center of the electrode columns where the misalignment between interpolation and real measures was the smallest. What has been shown here highlights a nonmonotonic excitability across the cortical depth (vertical Y direction of MEAs), displaying a maximum at intermediate layers where local neuronal populations emitted spikes at higher firing rates during Up states, which is compatible with the leading role of Layer 5 in the generation of SO found in previous in vitro studies (Sanchez-Vives and McCormick 2000; Wester and Contreras 2012). On the other hand, a heterogeneous excitability was also apparent along the direction parallel to the cortical surface (horizontal X direction of MEAs), as shown in Figure 4C, D, as both maximum MUA and TUp increased from left to right. A gradient of maximum MUA and Up duration (TUp) along the excitable strip (solid lines in Fig. 4C, D) would intuitively produce an acceleration of the waves traveling in the same direction, while a deceleration would be expected for those waves propagating towards less excitable regions of the slice. From a theoretical perspective, such inhomogeneity of the medium has been recently highlighted to underlie a hysteresis phenomenon (Capone and Mattia 2017). Indeed, a different speed profile should result depending on the propagation mode, as acceleration in one way is expected to be higher than the absolute value of the deceleration in the opposite direction. In fact, accelerating waves would find a more excitable tissue while, on the other way, wavefronts would be pushed ahead by less excitable cortical assemblies. According to this, a correlation would be expected between the gradient of the excitability (the change of TUp per unit horizontal length dTUp/dX, see Materials and Methods) and the speed difference measured in the same position between the leftward and the rightward propagating wave (Fig. 4J). This correlation, related to the interplay between the spatiotemporal patterns and the horizontal circuits of the cortex, was found in the majority of the experiments (Fig. 4K, 7 out of 12 experiments had significant correlations, P < 0.05 with a Pearson correlation, ρ > 0.5), further supporting the intimate relationship between spontaneous activity and the underlying cortical organization. Weakly Stable Down States Behind SO The analysis described so far shows that cortical slices spontaneously express waves generated by nonlinear self-excitation; wave propagation exhibits a rich repertoire of spatiotemporal patterns, which we have shown to be related to the structural organization of the slice. Taken together, these observations suggest that the cortical network behaves as a nonlinear excitable system that is balanced such that it avoids collapsing into stereotyped excited collective states (e.g., epileptic) while keeping sensitivity to the spatial structure of the connectivity. We therefore undertook to understand, in a theoretical model, the determinants of such interplay between excitability and sensitivity to structural features. For this purpose, we set up a large-scale model of a cortical slice as a network of modules arranged uniformly on a lattice, each of them exciting nearby modules with an intensity that decreased with distance (Fig. 5A, see Materials and Methods). Modules were composed of 350 LIF excitatory (80%) and inhibitory (20%) neurons each. SO occurred as the interplay between recurrent synaptic excitation and history-dependent self-inhibition (Fig. 5B), analogously to previous works (Bazhenov et al. 2002; Compte et al. 2003; Holcman and Tsodyks 2006; Mattia and Sanchez-Vives 2012). During Down states, activity fluctuations can drive the module through recurrent excitation towards a high-firing Up state, self-sustained in time by local synaptic reverberation involving a subset of both excitatory and inhibitory neurons (Fig. 5B, top). Each excitatory neuron possesses an activity-dependent inhibitory current modeling SFA due to the influx of a hyperpolarizing potassium current proportional to a fatigue level c(t), which increases each time a spike is emitted (Koch 1999). During periods of high-firing rate ν(t), fatigue c(t) accumulates (Fig. 5B, middle, purple line), eventually crossing a threshold level that destabilizes the Up state and determines the Up-to-Down transition and the beginning of the recovery phase during the Down state. Excitatory intermodular connectivity promotes a fast chain reaction following a local Down-to-Up transition, recruiting nearby inactive modules after recovery from previous active states, eventually generating global SO as in Figure 5B, bottom. A synthetic MUA from the firing rate of each module is computed by adding a Gaussian white noise unavoidably affecting electrophysiological recordings (Fig. 5E, see Materials and Methods). This allows analyzing and comparing experiments and simulations with the same tools as those used to work out state duration distribution (Fig. 5E, inset). In order to quantitatively reproduce the statistics of experimental SO, 2 minimal requirements had to be met. Firstly, Down states should be weakly stable, such that Down-to-Up transitions can occur with a non-negligible probability due only to fluctuations of local activity ν(t); alternatively, Down-to-Up transitions can be primed by the additional input due to the activation of nearby modules. Secondly, modules must be excitable, such that small suprathreshold inputs can elicit the onset of an Up state. Cortical modules with these features can be designed by relying on mean-field theory, in which the module nonlinear dynamics can be fully depicted in the phase plane (ν,c) (see Materials and Methods). Figure 5C shows the nullclines (loci of zero time derivative) for ν (black) and c (red), for a parameter choice meeting the above requirements: the (stable) Down state is the only intersection between the nullclines (i.e., the only point attractor of the dynamics), and the green curve describes an example trajectory followed by the system for an initial condition to the right of the ν nullcline. As the Down state is very close to the knee of the ν nullcline, its stability is weak and noise can promote an escape from this state. If endogenous fluctuations are taken into account, while the system approaches the Down fixed point, it can fluctuate for a long time around it (Fig. 5D, colored curve), until a large enough fluctuation makes it jump towards the upper branch of the ν nullcline (the Up state). If the synaptic input from other modules is increased without changing the other parameters, the ν nullcline shifts to the left (Fig. 5D, gray curve), and an almost deterministic Down-to-Up transition is elicited. The 2 illustrated scenarios address the requirements above, and would represent (1) the modules which due to fluctuations prime the activation wave onset (black nullcline) and (2) those which receiving input from nearby active modules contribute to the wave propagation (gray nullcline). Heterogeneous Excitability Shapes Spatiotemporal Activity Patterns As the next step, the cortical slice model was employed to investigate the mentioned interplay between local excitability and structural features suited to reproduce the spatiotemporal properties of propagating waves found in our experiments (Figs 2, 3 and 4). For this purpose, we differentially modulated the excitability level of the cortical modules in the slice model, with the most excitable modules in a strip (dark inverted U strip in Fig. 6A, see Materials and Methods). Following (Mattia and Sanchez-Vives 2012), the excitability of the cortical modules was regulated by changing the local connectivity cloc (probability to establish a synaptic contact between 2 excitatory neurons in the same module) and/or the external input from other modules, here dependent on the global connectivity cglob (probability of a connection with an excitatory neuron from nearby modules). Increasing cloc raised the firing rate of the Up states and strengthened its stability (upper branch of the ν nullcline became higher and wider, gray to black curves in Fig. 6B). At the same time, Down state stability was weakened as the corresponding attracting energy valley became shallower (in Fig. 6C, top, gray to black curves). In this cortical slice model, synchronized SO across the lattice of modules collectively emerged as an activation wave (Fig. 6D). The spatiotemporal pattern started with the Up state onset in a neuronal assembly close to the most excitable strip (Fig. 6D, white dashed line), where modules with maximum cloc were located (MCS). From there, the activation wave propagated horizontally following the maximally excitable modules. Maximum MUA (dark red) occurred along the MCS as the tip of the Up onset wavefront, consistent with what is shown in Figure 4. More specifically, we characterized activation waves by resorting to the cluster analysis used for the experiments, recognizing also in simulation different modes of propagation (Fig. 6E). Similar to what we found in vitro, the in silico cortical slice displayed a remarkable overlap between the MCS (dashed line), the EPS estimated from wave propagation and the strip where maximum MUA and TUp were measured from SO of the single modules (Fig. 6F). This overlap did not depend on the particular shape of the MCS, as we tested by performing simulations with different strip shapes (not shown). An Optimal Balance Between Short- and Long-Range Connectivity The capability to produce traveling waves of Up state onsets by assigning a pivotal role to a subset of maximally excitable neuronal assemblies, which in turn shape these spatiotemporal patterns, may appear as a rather straightforward expectation. With this in mind, we explored the sensitivity of the wave features expressed by the in silico slice to several key parameters of the network. First, we investigated the differential role of local and global connectivity. As expected, by increasing the number of excitatory synaptic connections (both local and global), an increase of the wave average speed was measured (Fig. 7A, see section Materials and Methods for details), as both the excitability of single modules and the coupling strength between modules raised accordingly. However, looking at the difference between the EPS and the MCS (“identification error”, see Materials and Methods), we found that it attained its minimum value for a specific connectivity configuration (Fig. 7B, green circle). For low maximum Cloc and Cglob (Fig. 7B, red circle), EPS displayed a wide distribution (Fig. 7C-top), on an average covering range of 100 mm from exact MCS (here an oblique straight line, dashed in Fig. 7C). This was the result of the reduced heterogeneity of the slice excitability, yielding more modules with similar Down-state stability, and hence capable to win the competition in leading the onset of the activation wave. The identification error was also high for relatively large maximum (Fig. 7B, blue circle). This weak relationship between structure and spontaneous activity occurred for a different reason: highly excitable modules had weakly stable Down states (Fig. 6B, shallower energy valleys), such that SO occurred also without the incoming input provided by the activation of a nearby cortical module. As result, the in silico slice under this condition was composed of oscillators which were more independent, and producing less stereotyped wavefronts. The hot spot with maximum Cloc ≅ Cglob ≅ 1 (Fig. 7B, green circle, same network parameters used in Fig. 6) was the only one to display a remarkable overlap between EPS distribution and MCS (Fig. 7C-middle). The existence of an optimal network minimizing the distance between EPS and MCS was apparent also when model slice was constrained to produce spontaneous waves at the same constant speed (see black dashed line in Fig. 7A, B corresponding to a speed of 9 mm/s). The uniqueness of the dynamical regime needed to reproduce the collected in vitro experimental evidence described in Figure 5 emerged also by testing the model network under the modulation of another critical parameter known to modulate the nonlinear dynamics of networks expressing SO (Mattia and Sanchez-Vives 2012): the intensity of an external excitatory input Iext modeling the incoming background activity produced outside the modeled neuronal field (see Materials and Methods). By increasing Iext ⁠, SO in single modules displayed Up states with longer durations (Fig. 7D), which in turn positively correlated (Pearson correlation, ρ = 0.56; P < 0.001) with the autocorrelation AC1 between the time lag arrays representing wave propagation, that is, the same as the one measured in the experiments (Fig. 3): the larger the excitability, the longer the Up state and the higher the correlation between the mode of propagation of successive waves (wave lag 1). Similarly to what we found in the experiments, on average no memory was retained about the way waves propagated previously in time (Fig. 7E, black line at ACn≅0 ⁠). A positive autocorrelation of the time lag arrays occurred only for simulations with a large enough Iext (red circles in Fig. 7D) representing the most excitable model networks. In this case, the average correlation (Fig. 7E, red line) was significantly different from 0 (Wilcoxon rank sum test, P < 0.01), in agreement with what we found in vitro (Fig. 3D). Layers 4 and 5 Lead Wave Propagation and SO Our matching observations from experiments and in silico slices strongly suggest that the cortical slice structure determines the features of in vitro spontaneous SO and their propagation across the tissue. A natural question follows up: does the most excitable region of the cortical slices inferred from their spontaneous activity correspond to a specific anatomical structure? Previous in vitro studies have highlighted a leading role of Layer 5 (L5) in the generation of SO (Sanchez-Vives and McCormick 2000). Neurons in this layer were the first to activate with respect to the others located at different depths in the same cortical column. However, the columnar interactions between supragranular and infragranular layers (L2/3 and L5) have been found to underlie the activity propagation of in vitro SO (Wester and Contreras 2012). Here, we addressed the question by matching the laminar organization characterized histologically and the excitable regions inferred from the spontaneous activity of the slices (Fig. 8A). More specifically, in the cortical slices in which histology was available (n = 9 out of 12 experiments), we computed the area covered by the distribution of EPS estimated from recorded waves in different layers (see Materials and Methods for details). As a result, we found a major overlap between L4 and L5 and EPS distribution (Fig. 8B), and hence with the slice region where maximum MUA and Up duration were also found. Fractions of EPS in L2/3 and L6 were significantly lower than those in L4 and L5 (Wilcoxon sign-rank test, P < 0.001). Although EPS area overlapping with L4 was on average lower than the one overlapping with L5, the difference was not significant, such that the functional identification of the most excitable region of the cortical slice highlighted a shared role between these 2 layers. Figure 8. View largeDownload slide L4 and L5 overlap with the excitable cortical region predicted from spontaneous activity propagation of in vitro SO. (A) EPS density and area covered by the MEA superimposed to the image of the corresponding cortical slice stained after the electrophysiological recordings. Dashed lines, boundaries between layers. EPS density below 0.04 is not shown. (B) Box plot of the area covered by EPS density higher than 0.04 in each layer of the cortical slices for which histological characterization was performed (n = 9 out of 12). Figure 8. View largeDownload slide L4 and L5 overlap with the excitable cortical region predicted from spontaneous activity propagation of in vitro SO. (A) EPS density and area covered by the MEA superimposed to the image of the corresponding cortical slice stained after the electrophysiological recordings. Dashed lines, boundaries between layers. EPS density below 0.04 is not shown. (B) Box plot of the area covered by EPS density higher than 0.04 in each layer of the cortical slices for which histological characterization was performed (n = 9 out of 12). Discussion Our study on the propagation of spontaneous, alternating Up and Down states across cortical slices provides a novel perspective on the dynamical origin of this multiscale phenomenon. We relied on an experimental setup in which a multielectrode array arranged in small electrode clusters acted as a multiresolution probe to investigate simultaneously local and global network dynamics. By this means, we reproducibly found a rich repertoire of propagation modes all of which shared 3 distinctive features. Firstly, activation waves systematically propagated parallel to the cortical surface, in agreement with previous works (Wu et al. 1999; Sanchez-Vives and McCormick 2000; Wester and Contreras 2012), although a wide distribution of speeds and ignition sites was found within each single experiment. Almost no local activity bumps were observed. Wavefronts of Up state onset were mildly curved, corresponding to a measured vertical speed along the columnar orientation, that was twice the horizontal one. Secondly, irrespective of their variability, wavefronts revealed loci across the slices that reproducibly led SO propagation and that composed a longitudinal smooth strip. Anatomical reconstruction showed that this strip largely overlapped with L4/L5, compatible with the evidence that the most excitable cell assemblies reside in L5 both in vitro (Sanchez-Vives and McCormick 2000; Wester and Contreras 2012; Krause et al. 2014) and in vivo (Sakata and Harris 2009; Chauvette et al. 2010; Amigó et al. 2015). Finally, wave onset was an almost memoryless process, that is, no correlations were found between propagation modes of successive waves. Such composite evidence traced a narrow path for the theoretical model we developed to understand its determinants: the prominent excitable strip pointed towards a distribution of local excitability; the variability of wave modes and their lack of memory hinted at an important role of noise. In fact, we found that a multimodular large-scale in silico model of the slice that reproduced all these features must primarily incorporate a distribution of excitability in the form of heterogeneous synaptic self-excitation of the single modules. Furthermore, we also found that local self-coupling must be optimally balanced with the mutual excitation between modules. These conditions ensured that our in silico slice acted as a structured subexcitable medium in which each module displayed a marginally stable Down state from which, thanks to maximal sensitivity to endogenous finite-size noise and intermodular input, relaxation oscillations could be generated. As a result, the most excitable strip overlapped with the sites with maximal duration of the Up states, and with maximal firing during Up states, as confirmed by our experiments and further highlighting the pivotal role of L4/L5 in leading SW activity. All findings support the hypothesis that this default mode spontaneously expressed by the cerebral cortex (Sanchez-Vives and Mattia 2014; Sanchez-Vives et al. 2017) is due to a privileged configuration of the network that allows both dynamical richness and sensitivity to the underlying structure. Comparisons With Previous Studies Horizontal intracortical connections play an important role in determining the speed of activity propagation (Wadman and Gutnick 1993; Golomb and Amitai 1997; Sanchez-Vives and McCormick 2000; Wester and Contreras 2012), together with the excitability level modulated by the balance between excitation and inhibition (Sanchez-Vives et al. 2010), as has also been described in theoretical models of SW activity (Bazhenov et al. 2002; Compte et al. 2003). In the specific framework of in vitro spontaneous SO, the horizontal speed of SWs that we measured (5.5 ± 1.5 mm/s) matches well those from previous studies found to range from 7.2 (Sanchez-Vives et al. 2010) to 10.9 mm/s (Sanchez-Vives and McCormick 2000) in ferret cortical slices (same preparation adopted here), up to the 20 mm/s measured in rat brain slices (Wester and Contreras 2012). Note that the high speed reported in the latter work may be due not only to the different species, but also to the use of voltage-sensitive dye (VSD) imaging, known to probe subthreshold membrane potential rather than the suprathreshold MUA as we recorded here. An even higher horizontal velocity of 30 mm/s was previously measured with VSD imaging in rat neocortical slices submerged in low (0.1 mM) magnesium (Wu et al. 1999). The oscillatory frequency under these hyperexcitable conditions was 7–10 Hz, which is higher than that measured in our experiments (0.31 ± 0.12 Hz). Under higher excitability conditions our theoretical argument would support less variable wave propagation modes (reduced sensitivity to the underlying structure) and more correlated successive propagation modes. Indeed, Wu et al. (1999) also observed that waves usually propagate in the same direction in low magnesium, while Sancristóbal et al. (2016) reported higher spatiotemporal regularity and more stable ignition sites in higher potassium. Population bursts of activity are also known to spontaneously occur in cultures of dissociated neurons (Marom and Shahaf 2002). Such network activations develop like traveling waves (Orlandi et al. 2013). Interestingly, several of the initiation and propagation features we measured in this study are comparable to such unstructured networks. More specifically, in neuronal cultures activity fluctuations are a determining factor in the population burst onset. Ignition sites are randomly activated in time, therefore sharing the same memoryless properties we observed in cortical slices. Finally, the speed of the activation waves in cultured neuronal networks was about 10 mm/s (Orlandi et al. 2013), mildly faster than in slices. All of this further strengthens the hypothesis of a common organizing principle of cortical networks always leading to a default activity pattern, namely the SW activity (Sanchez-Vives and Mattia 2014; Sanchez-Vives et al. 2017). Relationship With in Vivo SWs In anesthetized rodents SW activity is well documented, and has been reported to be similar to SWs during deep sleep (Alkire et al. 2008; Cirelli and Tononi 2008; Siegel 2008), thus raising the question regarding their possibly similar dynamic origin as in vitro SWs. A first difference is that in vivo Down-to-Up wavefronts propagate faster, with a horizontal speed tangential to the cortical surface ranging from 25 to 40 mm/s (Fucke et al. 2011; Ruiz-Mejias et al. 2011; Stroh et al. 2013; Sheroziya and Timofeev 2014). This is irrespective of the recording method, ranging from intracellular and extracellular recordings to 2-photon calcium imaging. The columnar or vertical speed of the Up state onset has been studied in vivo by directly probing cortical columns during SO (Sakata and Harris 2009; Chauvette et al. 2010). Columnar activation in vivo was more asymmetric as deep layers (L5 and L6) activated simultaneously, and the activity spread towards superficial layers (L2/3) with a speed of about 30 mm/s, given that a depth of 700 μm was covered in about 25 ms (Sakata and Harris 2009). Such quantitative differences in asymmetry and speeds may be plausibly related to the known lack of long-range corticocortical connections in vitro (Stepanyants et al. 2009; Schnepel et al. 2015), which is expected to be more pronounced along the direction tangential to the cortical surface. This would explain why in the cortical slice the vertical speed is twice the horizontal one, while these 2 components of the speed have almost the same magnitude in the intact brain. The lack of long-range horizontal connections may also affect the optimal balance between local excitability and intermodular connectivity, which we find to be critical for the kinematics of SWs. Finally, the absence of thalamus in our brain slices may of course contribute to such discrepancies. Limitations of This Study Activation and silencing of wavefronts were reconstructed by relying on an interpolation procedure between the electrodes of the MEA, and by clustering together similar waves to improve the signal-to-noise ratio of the state transition time lags. Under the assumption of reasonably smooth wavefronts, this method allowed recovering information about the wavefronts in slice regions where no electrodes were available. A failure of this assumption might imply a wrong identification of the depth of the strip leading wave propagation. Experimental testing with an additional electrode in different positions of the MEA validated such assumption directly proving the reliability of the interpolation estimates (see Figs 2 and 4, and Supplementary Material). Having denser arrays of electrodes would allow capturing additional details about the shape of the activation wavefronts, and further confirm the reliability of our findings, although a rather accurate measure is granted by the multiscale geometry of the MEA we designed and implemented. The interpolation approach that we have described here is highly adaptable to MEA with different numbers of channels and spatial distributions of electrodes. We then expect that the approach could be applied also to other experimental frameworks such as the study of in vivo quasi-planar waves. Another possible limitation of this study is that we did not sort single units (SU) from the MUA. Yet, we do not expect any qualitative change in the characterization of in vitro SWs from relying on SU activity instead of MUA. Indeed, using rather different methodologies with different spatial and temporal resolutions ranging from intracellular and extracellular recordings, to calcium and VSD imaging, the observed propagating modes shared in all cases similar features both in vivo (Fucke et al. 2011; Ruiz-Mejias et al. 2011; Stroh et al. 2013; Sheroziya and Timofeev 2014) and in vitro (Wu et al. 1999; Sanchez-Vives and McCormick 2000; Sanchez-Vives et al. 2010; Wester and Contreras 2012). Functional and Dysfunctional Implications An appealing result of the modeling study was finding that, for the model to match the experiment, a proper balance of excitation between local and mutual modules was critical. This, together with the need of a weakly stable Down state, establishes a “sweet spot” of the network parameters where structure is influential but at the same time allows for a wide variety of propagation modes. It has been proposed that optimal sensitivity to structure can be a smart way to implement “self-organized instability” (Solé et al. 2002). However, whenever new stimuli are sensed, the cortical network moves to a different working point, likely no longer optimal. To bring the system back to an optimal working regime, a mechanism is needed to recover the proper weak stability of the Down state. A simple unspecific modulation of the network excitability would serve this purpose. As a result, the variability of stimulus-evoked waves would not be purely random, but rather sensitive to the network structure (Rabinovich et al. 2008). On the other hand, pushing the cortical tissue outside this sweet spot likely distorts the spontaneous activity, leading the network to display pathological dynamics, which can persist if optimality is not recovered. For instance, a wrong balance between excitation and inhibition is now known to influence the speed of the SWs traveling across the frontal cortex of animal models of Down’s syndrome (Ruiz-Mejias et al. 2016). An underexpressed synaptic inhibition seems to underlie the break of long-range coherence of SW activity in mouse models of the Alzheimer’s disease (Busche et al. 2015). In these cases our analytical approach, combined with properly tuned in silico models, could help to elucidate which combination of local excitability and intermodular synaptic coupling underlie SW activity changes, consequently unraveling the mechanistic roots of pathological dynamics. Authors’ Contributions M.M., C.C., M.V.S.V. designed research; B.R. performed and M.V.S.V. supervised experimental research; C.C. performed and M.M. and P.D.G. supervised simulations and analyses; A.M. and B.R. performed histology; X.I., M.M., and M.V.S.V. co-designed arrays and X.I. fabricated the new electrode array; M.M., P.D.G., C.C., and M.V.S.V. wrote the article. Funding CORTICONIC (EC FP7 grant 600806), the EU Horizon 2020 Research and Innovation Program under HBP SGA1 (grant no. 720270), BFU2014-52467-R, FLAGERA-PCIN-2015-162-C02-01, MINECO (Spain), and CERCA Program/Generalitat de Catalunya. We thank C. Gonzalez-Liencres for editing assistance. Notes Conflict of Interest: None declared. References Alkire MT , Hudetz AG , Tononi G . 2008 . Consciousness and anesthesia . Science . 322 : 876 – 880 . Google Scholar Crossref Search ADS PubMed Amigó JM , Monetti R , Tort-Colet N , Sanchez-Vives MV . 2015 . Layer 5 leads information flow during slow oscillations according to information directionality indicators . J Comput Neurosci . 39 : 53 – 62 . Google Scholar Crossref Search ADS PubMed Amit DJ , Brunel N . 1997 . Model of global spontaneous activity and local structured activity during delay periods in the cerebral cortex . Cereb Cortex . 7 : 237 – 252 . Google Scholar Crossref Search ADS PubMed Bazhenov M , Timofeev I , Steriade M , Sejnowski TJ . 2002 . Model of thalamocortical slow-wave sleep oscillations and transitions to activated states . J Neurosci . 22 : 8691 – 8704 . Google Scholar Crossref Search ADS PubMed Benucci A , Frazor RA , Carandini M . 2007 . Standing waves and traveling waves distinguish two circuits in visual cortex . Neuron . 55 : 103 – 117 . Google Scholar Crossref Search ADS PubMed Bringuier V , Chavane F , Glaeser L , Frégnac Y . 1999 . Horizontal propagation of visual activity in the synaptic integration field of area 17 neurons . Science . 283 : 695 – 699 . Google Scholar Crossref Search ADS PubMed Busche MA , Kekuš M , Adelsberger H , Noda T , Förstl H , Nelken I , Konnerth A . 2015 . Rescue of long-range circuit dysfunction in Alzheimer’s disease models . Nat Neurosci . 18 : 1623 – 1630 . Google Scholar Crossref Search ADS PubMed Capone C , Mattia M . 2017 . Speed hysteresis and noise shaping of traveling fronts in neural fields: role of local circuitry and nonlocal connectivity . Sci Rep . 7 : 39611 . Google Scholar Crossref Search ADS PubMed Chauvette S , Volgushev M , Timofeev I . 2010 . Origin of active states in local neocortical networks during slow sleep oscillation . Cereb Cortex . 20 : 2660 – 2674 . Google Scholar Crossref Search ADS PubMed Cirelli C , Tononi G . 2008 . Is sleep essential? PLoS Biol . 6 : 1605 – 1611 . Google Scholar Crossref Search ADS Compte A , Sanchez-Vives MV , McCormick DA , Wang X-J . 2003 . Cellular and network mechanisms of slow oscillatory activity (<1 Hz) and wave propagations in a cortical network model . J Neurophysiol . 89 : 2707 – 2725 . Google Scholar Crossref Search ADS PubMed Cross MC , Hohenberg PC . 1993 . Pattern formation outside of equilibrium . Rev Mod Phys . 65 : 851 – 1112 . Google Scholar Crossref Search ADS Ermentrout GB , Kleinfeld D . 2001 . Traveling electrical waves in cortex: insights from phase dynamics and speculation on a computational role . Neuron . 29 : 33 – 44 . Google Scholar Crossref Search ADS PubMed Ferezou I , Haiss F , Gentet LJ , Aronoff R , Weber B , Petersen CCH . 2007 . Spatiotemporal dynamics of cortical sensorimotor integration in behaving mice . Neuron . 56 : 907 – 923 . Google Scholar Crossref Search ADS PubMed Fucke T , Suchanek D , Nawrot MP , Seamari Y , Heck DH , Aertsen A , Boucsein C . 2011 . Stereotypical spatiotemporal activity patterns during slow-wave activity in the neocortex . J Neurophysiol . 106 : 3035 – 3044 . Google Scholar Crossref Search ADS PubMed Golomb D , Amitai Y . 1997 . Propagating neuronal discharges in neocortical slices: computational and experimental study . J Neurophysiol . 78 : 1199 – 1211 . Google Scholar Crossref Search ADS PubMed Han F , Caporale N , Dan Y . 2008 . Reverberation of recent visual experience in spontaneous cortical waves . Neuron . 60 : 321 – 327 . Google Scholar Crossref Search ADS PubMed Holcman D , Tsodyks M . 2006 . The emergence of Up and Down states in cortical networks . PLoS Comput Biol . 2 : e23 . Google Scholar Crossref Search ADS PubMed Illa X , Rebollo B , Gabriel G , Sánchez-Vives MV , Villa R . 2015 . A SU-8-based flexible microprobe for close and distal recordings from the cortical network. In: SPIE microtechnologies . International Society for Optics and Photonics . p. 951803 . doi:10.1117/12.2180454. Jancke D , Chavane F , Naaman S , Grinvald A . 2004 . Imaging cortical correlates of illusion in early visual cortex . Nature . 428 : 423 – 426 . Google Scholar Crossref Search ADS PubMed Koch C . 1999 . Biophysics of computation: information processing in single neurons . New York : Oxford University Press . Krause BM , Raz A , Uhlrich DJ , Smith PH , Banks MI . 2014 . Spiking in auditory cortex following thalamic stimulation is dominated by cortical network activity . Front Syst Neurosci . 8 : 1 – 24 . Google Scholar PubMed Lee S-H , Blake R , Heeger DJ . 2007 . Hierarchy of cortical responses underlying binocular rivalry . Nat Neurosci . 10 : 1048 – 1054 . Google Scholar Crossref Search ADS PubMed Luczak A , Barthó P , Harris KD . 2009 . Spontaneous events outline the realm of possible sensory responses in neocortical populations . Neuron . 62 : 413 – 425 . Google Scholar Crossref Search ADS PubMed Marom S , Shahaf G . 2002 . Development, learning and memory in large random networks of cortical neurons: lessons beyond anatomy . Q Rev Biophys . 35 : 63 – 87 . Google Scholar Crossref Search ADS PubMed Mascaro M , Amit DJ . 1999 . Effective neural response function for collective population states . Network . 10 : 351 – 373 . Google Scholar Crossref Search ADS PubMed Massimini M , Huber R , Ferrarelli F , Hill SL , Tononi G . 2004 . The sleep slow oscillation as a traveling wave . J Neurosci . 24 : 6862 – 6870 . Google Scholar Crossref Search ADS PubMed Mattia M , Del Giudice P . 2000 . Efficient event-driven simulation of large networks of spiking neurons and dynamical synapses . Neural Comput . 12 : 2305 – 2329 . Google Scholar Crossref Search ADS PubMed Mattia M , Pani P , Mirabella G , Costa S , Del Giudice P , Ferraina S . 2013 . Heterogeneous attractor cell assemblies for motor planning in premotor cortex . J Neurosci . 33 : 11155 – 11168 . Google Scholar Crossref Search ADS PubMed Mattia M , Sanchez-Vives MV . 2012 . Exploring the spectrum of dynamical regimes and timescales in spontaneous cortical activity . Cogn Neurodyn . 6 : 239 – 250 . Google Scholar Crossref Search ADS PubMed Mohajerani MH , McVea D a , Fingas M , Murphy TH . 2010 . Mirrored bilateral slow-wave cortical activity within local circuits revealed by fast bihemispheric voltage-sensitive dye imaging in anesthetized and awake mice . J Neurosci . 30 : 3745 – 3751 . Google Scholar Crossref Search ADS PubMed Muller L , Reynaud A , Chavane F , Destexhe A . 2014 . The stimulus-evoked population response in visual cortex of awake monkey is a propagating wave . Nat Commun . 5 : 3675 . Google Scholar Crossref Search ADS PubMed Orlandi JG , Soriano J , Alvarez-Lacalle E , Teller S , Casademunt J . 2013 . Noise focusing and the emergence of coherent activity in neuronal cultures . Nat Phys . 7 : 1 – 9 . Rabinovich MI , Huerta R , Laurent G . 2008 . Transient dynamics for neural processing . Science . 321 : 48 – 50 . Google Scholar Crossref Search ADS PubMed Reig R , Mattia M , Compte A , Belmonte C , Sanchez-Vives MV . 2010 . Temperature modulation of slow and fast cortical rhythms . J Neurophysiol . 103 : 1253 – 1261 . Google Scholar Crossref Search ADS PubMed Rubino D , Robbins KA , Hatsopoulos NG . 2006 . Propagating waves mediate information transfer in the motor cortex . Nat Neurosci . 9 : 1549 – 1557 . Google Scholar Crossref Search ADS PubMed Ruiz-Mejias M , Ciria-Suarez L , Mattia M , Sanchez-Vives MV . 2011 . Slow and fast rhythms generated in the cerebral cortex of the anesthetized mouse . J Neurophysiol . 106 : 2910 – 2921 . Google Scholar Crossref Search ADS PubMed Ruiz-Mejias M , Martinez de Lagran M , Mattia M , Castano-Prat P , Perez-Mendez L , Ciria-Suarez L , Gener T , Sancristóbal B , García-Ojalvo J , Gruart A , et al. 2016 . Overexpression of Dyrk1A, a down syndrome candidate, decreases excitability and impairs gamma oscillations in the prefrontal cortex . J Neurosci . 36 : 3648 – 3659 . Google Scholar Crossref Search ADS PubMed Sagués F , Sancho JM , García-Ojalvo J . 2007 . Spatiotemporal order out of noise . Rev Mod Phys . 79 : 829 – 882 . Google Scholar Crossref Search ADS Sakata S , Harris KD . 2009 . Laminar structure of spontaneous and sensory-evoked population activity in auditory cortex . Neuron . 64 : 404 – 418 . Google Scholar Crossref Search ADS PubMed Sanchez-Vives MV , Massimini M , Mattia M . 2017 . Shaping the default activity pattern of the cortical network . Neuron . 94 : 993 – 1001 . Google Scholar Crossref Search ADS PubMed Sanchez-Vives MV , Mattia M . 2014 . Slow wave activity as the default mode of the cerebral cortex . Arch Ital Biol . 152 : 147 – 155 . Google Scholar PubMed Sanchez-Vives MV , Mattia M , Compte A , Perez-Zabalza M , Winograd M , Descalzo VF , Reig R . 2010 . Inhibitory modulation of cortical up states . J Neurophysiol . 104 : 1314 – 1324 . Google Scholar Crossref Search ADS PubMed Sanchez-Vives MV , McCormick DA . 2000 . Cellular and network mechanisms of rhythmic recurrent activity in neocortex . Nat Neurosci . 3 : 1027 – 1034 . Google Scholar Crossref Search ADS PubMed Sancristóbal B , Rebollo B , Boada P , Sanchez-Vives MV , Garcia-Ojalvo J . 2016 . Collective stochastic coherence in recurrent neuronal networks . Nat Phys . 12 : 1 – 8 . Google Scholar Crossref Search ADS Sato TK , Nauhaus I , Carandini M . 2012 . Traveling waves in visual cortex . Neuron . 75 : 218 – 229 . Google Scholar Crossref Search ADS PubMed Schnepel P , Kumar A , Zohar M , Aertsen A , Boucsein C . 2015 . Physiology and impact of horizontal connections in rat neocortex . Cereb Cortex . 25 : 3818 – 3835 . Google Scholar Crossref Search ADS PubMed Sheroziya M , Timofeev I . 2014 . Global intracellular slow-wave dynamics of the thalamocortical system . J Neurosci . 34 : 8875 – 8893 . Google Scholar Crossref Search ADS PubMed Siegel JM . 2008 . Do all animals sleep? Trends Neurosci . 31 : 208 – 213 . Google Scholar Crossref Search ADS PubMed Solé RV , Alonso D , McKane AJ . 2002 . Self-organized instability in complex ecosystems . Philos Trans R Soc Lond B . 357 : 667 – 671 . Google Scholar Crossref Search ADS Stepanyants A , Martinez LM , Ferecskó AS , Kisvárday ZF , Ferecsko AS . 2009 . The fractions of short- and long-range connections in the visual cortex . Proc Natl Acad Sci USA . 106 : 3555 – 3560 . Google Scholar Crossref Search ADS PubMed Stroh A , Adelsberger H , Groh A , Rühlmann C , Fischer S , Schierloh A , Deisseroth K , Konnerth A . 2013 . Making waves: initiation and propagation of corticothalamic Ca2+ waves in vivo . Neuron . 77 : 1136 – 1150 . Google Scholar Crossref Search ADS PubMed Wadman WJ , Gutnick MJ . 1993 . Non-uniform propagation of epileptiform discharge in brain slices of rat neocortex . Neuroscience . 52 : 255 – 262 . Google Scholar Crossref Search ADS PubMed Wester JC , Contreras D . 2012 . Columnar interactions determine horizontal propagation of recurrent network activity in neocortex . J Neurosci . 32 : 5454 – 5471 . Google Scholar Crossref Search ADS PubMed Wilson HR , Cowan JD . 1972 . Excitatory and inhibitory interactions in localized populations of model neurons . Biophys J . 12 : 1 – 24 . Google Scholar Crossref Search ADS PubMed Wu J-Y , Guan L , Tsau Y . 1999 . Propagating activation during oscillations and evoked responses in neocortical slices . J Neurosci . 19 : 5005 – 5015 . Google Scholar Crossref Search ADS PubMed Xu W , Huang X , Takagaki K , Wu J . 2007 . Compression and reflection of visually evoked cortical waves . Neuron . 55 : 119 – 129 . Google Scholar Crossref Search ADS PubMed © The Author 2017. Published by Oxford University Press. All rights reserved. For Permissions, please e-mail: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - Slow Waves in Cortical Slices: How Spontaneous Activity is Shaped by Laminar Structure
JF - Cerebral Cortex
DO - 10.1093/cercor/bhx326
DA - 2019-01-01
UR - https://www.deepdyve.com/lp/oxford-university-press/slow-waves-in-cortical-slices-how-spontaneous-activity-is-shaped-by-xYDHPYomj4
SP - 319
VL - 29
IS - 1
DP - DeepDyve
ER -