Criticality and universality in neuronal cultures during “up” and “down” states

Abstract

<p>The brain can be seen as a self-organized dynamical system that optimizes information processing and storage capabilities. This is supported by studies across scales, from small neuronal assemblies to the whole brain, where neuronal activity exhibits features typically associated with phase transitions in statistical physics. Such a critical state is characterized by the emergence of scale-free statistics as captured, for example, by the sizes and durations of activity avalanches corresponding to a cascading process of information flow. Another phenomenon observed during sleep, under anesthesia, and in <italic>in vitro</italic> cultures, is that cortical and hippocampal neuronal networks alternate between “up” and “down” states characterized by very distinct firing rates. Previous theoretical work has been able to relate these two concepts and proposed that only up states are critical whereas down states are subcritical, also indicating that the brain spontaneously transitions between the two. Using high-speed high-resolution calcium imaging recordings of neuronal cultures, we test this hypothesis here by analyzing the neuronal avalanche statistics in populations of thousands of neurons during “up” and “down” states separately. We find that both “up” and “down” states can exhibit scale-free behavior when taking into account their intrinsic time scales. In particular, the statistical signature of “down” states is indistinguishable from those observed previously in cultures without “up” states. We show that such behavior can not be explained by network models of non-conservative leaky integrate-and-fire neurons with short-term synaptic depression, even when realistic noise levels, spatial network embeddings, and heterogeneous populations are taken into account, which instead exhibits behavior consistent with previous theoretical models. Similar differences were also observed when taking into consideration finite-size scaling effects, suggesting that the intrinsic dynamics and self-organization mechanisms of these cultures might be more complex than previously thought. In particular, our findings point to the existence of different mechanisms of neuronal communication, with different time scales, acting during either high-activity or low-activity states, potentially requiring different plasticity mechanisms.</p>

1 Introduction

The description of neuronal dynamics within the framework of critical phenomena has become commonplace amongst physicists in the last decade or so (Chialvo, 2010; Massobrio et al., 2015; O'Byrne and Jerbi, 2022). The first signatures of criticality in the brain were depicted as neuronal avalanches in organotypic cortical cultures (Beggs and Plenz, 2003, 2004), where sequences of high-frequency neuronal activations across a small population could be described by scale-free statistics of their sizes and durations distributions. Nowadays, similar critical signatures have been observed across many systems and preparations: from power-law statistics of correlations in whole-brain recordings (Tagliazucchi et al., 2012; Haimovici et al., 2013; Ponce-Alvarez et al., 2018) to neuronal avalanches in slices (Gireesh and Plenz, 2008), dissociated cultures (Pasquale et al., 2008; Yaghoubi et al., 2018), and in vivo (Petermann et al., 2009; Bellay et al., 2015; Yu et al., 2017; Priesemann et al., 2013; Curic et al., 2021; Rabus et al., 2023). They indicate the emergence of complex spatiotemporal dynamics with statistics compatible with a system being in the neighborhood of a critical point, in particular, near a second-order phase transition, or of a critical region (Moretti and Muñoz, 2013).

From a statistical physics point of view, neuronal avalanches are interpreted as a branching process (Beggs and Plenz, 2003; Korchinski et al., 2021), where an active neuron has a finite probability of activating its neighbors. When each active neuron induces, on average, the firing of a single neighbor, the system is thought to be in a critical state, where the activity neither explodes nor always dies out quickly. This results in a scale-free distribution of activation sequences observables, namely the size and duration of the neuronal avalanches. A key assumption in this picture is that there exists a separation of time scales, where the spontaneous activations of neurons (those that initiate an avalanche) happen much more slowly than the spreading of that avalanche across the population, i.e., there only exists a single avalanche at any given time. However, in most neuronal systems, that is far from the truth (Orlandi et al., 2013; Williams-Garćıa et al., 2017), and many avalanches can coexist simultaneously, resulting in a more difficult interpretation of the observed statistics (Korchinski et al., 2021).

Until recently (Yaghoubi et al., 2018; Yu et al., 2017), neuronal avalanches had mostly been described in systems where the network switches between periods of high and low-frequency activity such that the periods of high activity dominated the neuronal avalanche statistics. In vivo and in some slice preparations, such a switching behavior corresponds to the so-called “up” and “down” states, i.e., slow cortical oscillations present during slow-wave sleep (Sanchez-Vives and McCormick, 2000) that are initiated by pyramidal neurons near layer V and propagate toward the other cortical layers. In other preparations, like dissociated cultures (Orlandi et al., 2013), the bistable behavior is more akin to hippocampal sharp wave-ripples (Levenstein et al., 2018), which also have a well-defined size and duration. Theoretically, this phenomenon has often been described in terms of synchronization (Penn et al., 2016) or as spontaneous switching around a bifurcation (Millman et al., 2010). However, both cortical up and down states and hippocampal short-wave ripples have a strong spatial component, initiating at specific sites and propagating throughout the tissue, reminiscent of a classical spatially-extended excitable system (Orlandi and Casademunt, 2017).

The specific hypothesis of spontaneous switching around a bifurcation implies that the dynamics in the high-activity state are critical whereas it is subcritical in the low-activity state (Millman et al., 2010). Here, we test this hypothesis explicitly in neuronal cultures that do show alternating activity behavior by analyzing the two different states separately. We find that both high-activity and low-activity states can exhibit similar critical signatures if appropriate time scales are chosen for defining neuronal avalanches. These results point to the existence of different mechanisms of neuronal communication, with different time scales, acting during either high-activity or low-activity states. We also show that detailed model simulations of dissociated cultures follow the aforementioned hypothesis and are, thus, incompatible with the experimental observations for the low-activity state. This suggests the existence of processes with long timescales (of the order of a few hundred ms) that play a significant role in shaping the dynamics during the low-activity state that are currently not captured by existing models.

2 Materials and methods

2.1 Hippocampal cultures

Cultures from dissociated hippocampal neurons and glial cells, prepared from newborn P0 Sprague-Dawley rats, were plated on Si chips of 1 mm thickness, and 1 cm² surface area, placed on individual 24-well plate wells; as described previously (Colicos et al., 2001; Girotto et al., 2013). Each chip was Matrigel-coated (Beckton Dickinson) and placed in Basal Medium Eagle (BME). Cells were initially plated at a density of 30,000 cells/ml. The culture medium was not changed during the first week and every 4 days thereafter.

2.2 Calcium imaging

Cultures were grown for up to 2–3 weeks before imaging. Prior to imaging, cultures were incubated with Fluo-4 calcium indicator for 20 min. Afterward, cultures were washed and placed in an individual well with an extracellular bath solution (EBS) containing 135 mM NaCl, 10 mM glucose, 3 mM CaCl₂, 5 mM KCl, 2 mM MgCl₂ and 5 mM Hepes, pH was adjusted to 7.3 with NaOH, and osmolarity to 310 mOsm with Sorbital. Calcium imaging was performed following References (Colicos et al., 2001; Yaghoubi et al., 2018). In brief, we used a high temporal resolution camera that allowed us to record neuronal activity at 200 fps (Hammamatsu Orca-Flash 4.0) on an upright microscope with low magnification (field of view of L≈500 μm). We recorded the spontaneous, non-stimulated fluorescence neuronal activity of up to 1,500 neurons for 20 min (pixel resolution of 0.65μm/pixel). With this setup we are able to record most neuronal activity within a large field of view with single-cell resolution and a temporal resolution comparable to the fastest time-scales of synaptic integration; avoiding several of the drawbacks caused by spatial and temporal sub-sampling and hidden neurons (Ribeiro et al., 2014; Levina and Priesemann, 2017).

2.3 Pharmacology

In a subset of experiments, we blocked inhibitory connections by performing bath application of picrotoxin (PTX), a noncompetitive GABA_A receptor antagonist, with a concentration of 50 μM during 15 min at room temperature after the first imaging session. Cultures were then typically imaged in a different field of view for an additional 20 min.

2.4 Data preprocessing

Data preprocessing of calcium imaging experiments was performed as previously described in Fernández-García et al. (2020) using the NETCAL software (Orlandi et al., 2017) platform (see Figure 1). In brief, cell ROIs were automatically detected using a simple thresholding procedure on the time-averaged image of the recording and posteriorly cleaned up with morphological opening operations. Time series for each ROI were extracted, detrended, and normalized to ΔF/F₀ units (where F₀ was computed prior to the detrending operation). Spike inference was performed using the OASIS algorithm (Friedrich et al., 2017). See Table 1 for a summary of the list of recording experiments and their properties.

Figure 1

Figure 1. Calcium imaging analysis workflow. (A) 1,000–1,500 cells are automatically segmented within the field of view based on their average fluorescence. Inset: Zoom in on the field of view with four characteristic cells highlighted. Right: Timeseries from the same 4 highlighted cells as well as their inferred spike trains (black bars) before (B) and after (C) PTX application. For this particular example, we imaged the same cells before and after PTX. Only a smaller time window out of the 20 min recordings is shown.

Table 1

Table 1. Summary of the properties of different recordings.

2.5 Mathematical model and simulations

The dynamics of cultures from dissociated neurons were modeled and simulated following (Orlandi et al., 2013; Orlandi and Casademunt, 2017). In brief, for the 'homogeneous' simulations (simulations 1 and 2), single neuron dynamics are modeled by a quadratic integrate and fire model with adaptation (Alvarez-Lacalle and Moses, 2009; Izhikevich, 2003, 2007), i.e.,

C.v=k(v−vr)(v−vt)−u+I+I0+η, (1)Cv.=k(v−vr)(v−vt)−u+I+I0+η, (1)

τa.u=b(v−vr)−u, (2)τau.=b(v−vr)−u, (2)

if v≥vp,then v←vc,u←u+d. (3)if v≥vp,then v←vc,u←u+d. (3)

Here, Equation 1 corresponds to the dynamics of the soma membrane potential v(t) and v_r = −60 mV and v_t = −40 mV are the resting and threshold potentials, respectively. C = 100 ms is the normalized membrane capacitance and u models an inhibitory current that represents internal slow currents generated by the activation of ion channels. I accounts for the synaptic inputs from other neurons. I₀ represents the spontaneous release of synaptic vesicles (minis). The time between successive minis is often modeled as a memoryless, exponential process (Fatt and Katz, 1952). Hence, the number of minis released in a given interval will follow a Poisson process. The rate of this process depends on the experimental conditions (Ivenshitz and Segal, 2010), but in our simulations, and following previous work (Orlandi et al., 2013), we selected a frequency of λ = 0.1 ms⁻¹, which is chosen to produce a spontaneous firing rate of the neurons of the order of 0.1 Hz. Each mini generates a current with amplitude g_m = 30 mV which decays exponentially with a time constant τ_m = 10 ms. η is a white noise term with autocorrelation ⟨η(t)η(t′)⟩=2gsδ(t−t′)⟨η(t)η(t′)⟩=2gsδ(t-t′) with g_s = 30 mV. Equation 2 is the evolution of the slow currents and τ_a = 33 ms, k = 0.7 mV⁻¹, b = −2 and d = 100 mV are parameters that control recovery and adaptation.

Every time a given neuron i fires, it produces an (excitatory or inhibitory) current on its output neighbors of the form

Ii(t,tm)=gDi(tm)exp(t−tmτ)Θ(t−tm) (4)Ii(t,tm)=gDi(tm)exp(t−tmτ)Θ(t−tm) (4)

where t_m is the spike time and g is the synaptic strength. Its value for excitatory synapses is used as a control parameter to obtain the desired time interval separating subsequent up states (which is achieved for g≈40 mV) while the inhibitory strength is fixed at −50 mV. τ is the characteristic time constant (10 ms for excitatory currents and 20 ms for inhibitory ones). Θ(t) is the Heaviside function and D short-term synaptic depression (STD). The evolution of STD is described by

Ḋ=1τD(1−D)−(1−β)Dδ(t−tm), (5)Ḋ=1τD(1−D)−(1−β)Dδ(t−tm), (5)

where τ_D (2 s for excitatory and 0.2 s for inhibitory cells), is the synaptic recovery time and β (0.8 for excitatory and 0.95 for inhibitory cells) controls the level of depression after each spike.

To mimic the experimental cultures and establish realistic connectivity patterns, neurons were placed randomly on a square region with 10 mm sidelength until a density of ρ = 800 neurons/mm² was reached (80,000 neurons total). For each neuron, an axon was grown as a biased random walk with total length given by a Rayleigh distribution with variance 900 μm². Starting from the cell body a starting angular direction was picked randomly and a segment of 10 μm was grown. At the end of the segment, a new direction was chosen centered on the previous one following a Gaussian distribution with a standard deviation of 15 °. This process was repeated until the desired total length was reached. For each neuron, a dendritic tree was modeled as an effective circular area of interaction with a radius obtained from a Gaussian distribution with mean 150 μm and standard deviation 40 μm. If an axon crossed the dendritic tree of another neuron, a connection was established with probability 0.13. Finally, 20 % of the neurons were randomly chosen to be inhibitory and the remaining ones excitatory.

For the “heterogeneous” simulations (simulations 3 and 4) we lowered the standard excitatory population to 70% and added 10% of excitatory bursty cells (accomplished by changing v_c = −40 and d = 50). To all excitatory cells, we changed the resting membrane potential to a base value of v_r = −62 mV and added an offset drawn from a Rayleigh distribution with standard deviation σ = 2 mV, to add variability to the firing rates consistent with experimental data.

Each simulation had a fixed run length of 1 h and was simulated with a first-order Euler algorithm with a time step of 0.1 ms. Random numbers were generated with the MTGP32 implementation of the Mersenne Twister for the GPU (Saito and Matsumoto, 2013) and initialized with a random seed for each simulation. In total, we ran simulations in 4 different conditions: simulations 1 and 2 corresponded to homogeneous networks, and 3 and 4 to heterogeneous ones. In simulations 2 and 4 we mimicked the effect of blocking inhibition by setting the strength of inhibitory connections to 0.

2.6 Detection of up and down states

The detection of up and down states is done based on thresholding the population firing rate. The firing rate for each frame is defined as the number of spikes in that frame normalized by the number of neurons. A Gaussian kernel with σ= 5 frames for experimentally recorded data and σ= 20 frames for simulated data is used to smooth the firing rate traces. To separate up and down states we used a Schmitt trigger (Taub and Schilling, 1977), thresholding the normalized activity with an upper threshold = 0.001 and a lower threshold = 0.0003. We kept the same criteria for all recordings and simulations.

2.7 Neuronal avalanches and scaling collapse procedure

Following the standard approach for spiking data (Friedman et al., 2012; Pasquale et al., 2008), a neuronal avalanche is defined as the largest sequence of consecutive time bins containing spikes in every single time bin, separated by time bins during which none of the neurons in the culture fire (see Supplementary Figure S1 for some examples). The avalanche duration, T, corresponds to the number of time bins and the avalanche size, S, is the total number of spikes over the duration of an avalanche. Based on this definition, for a fixed number of neurons, it is expected that the choice of the size of the time bin affects the avalanche statistics, specifically the size S and duration T of avalanches. To capture this, we follow here the standard finite size scaling approach used in the context of phase transitions (see, for example, Christensen and Moloney, 2005). In our context, its original formulation in terms of a varying number of neurons can be recast in terms of a varying size of the time bin. It is based on the hypothesis that for a fixed number of neurons, the effect of temporal bin size on scale-free avalanche statistics (e.g., avalanche size distribution, P(S)) can be taken into account as follows:

P(S)~S−τ×f(S/binβS). (6)P(S)~S−τ×f(S/binβS). (6)

This functional form indicates, for suitably chosen scaling exponents τ and β_S, the distributions for different bin sizes can be collapsed onto the scaling function f by plotting P(S) × S^τ vs S/binβSS/binβS . In the case of exponential-like avalanche statistics (τ = 0), one can achieve a similar data collapse with a scaling exponent β_S when plotting P(S)×binβSP(S)×binβS vs S/binβSS/binβS . In both cases, to numerically determine the values of the scaling exponent(s) the distinctive features of each graph need to align horizontally for the largest arguments aka the right-hand side of the graph (Christensen and Moloney, 2005), which is typically assessed by eye as we do here and allows us to estimate the uncertainties.

To replicate the experimental field of view (which covers only a small fraction of the whole culture) in the simulations, we only analyzed the avalanche statistics of neurons within a circular patch of radius r = 0.5 mm for the E simulations and r = 0.8 mm in the E + I simulations. This radius was chosen to ensure that the total activity rate in the up states was of the order of 1 event per frame or time step (0.1 ms). Each of these patches contained in total 400 to 1,000 neurons. To increase the number of avalanches and obtain reliable estimates of the avalanche statistics, for each simulation and each experiment, we randomly selected 100 different local patches consisting of 50% of the neurons and combined the avalanches from all patches into a single distribution. This approach also allows us to obtain avalanche statistics for different temporal bin sizes.

2.8 p-value estimation for power-law distributions

After identifying the optimal exponents α and β, using the scaling collapse procedure (Equation 6), our next step is to assess the validity of the power-law model as a hypothesis for the specific recording or simulation. To achieve this, we check whether we can identify an extended range for the size or duration in the avalanche distribution function that gives a reasonably high p-value (>0.1). A detailed description of the process of identifying the range is presented in the captions of Supplementary Figure S2. To find the p-value, we first used the Kolmogorov-Smirnov (KS) statistic. KS statistic quantifies the distance between two probability distributions (Press et al., 2007), which is defined as the maximum distance between the cumulative distribution functions (CDFs) of the two distributions (here data and fitted model):

de=max∣∣Se(x)−Pe(x)∣∣,Xlow≤x≤Xhigh (7)de=max|Se(x)−Pe(x)|,Xlow≤x≤Xhigh (7)

Here, S_e(x) is the CDF of the empirical data, and P_e(x) is the CDF corresponding to the fitted model. The fitted model is estimated using the power-law scaling procedure related to Equation 6. The reported p-value is the probability of observing a KS value bigger than d_e for synthetic data generated by the fitted model and provides a measure of whether it is likely that the empirical data do indeed follow the fitted model. One can show that its value can be calculated from the following theoretical expression (Deluca and Corral, 2013):

p−value=2∞∑i=1(−1)i−1exp[−2i2(de√n+0.12de+0.11de/√n)2], (8)p−value=2∑i=1∞(−1)i−1exp[−2i2(den+0.12de+0.11de/n)2], (8)

where n is the number of samples in the data set. As mentioned above, to enhance our statistics we identified avalanches over 100 iterations (for both experimental and simulation data). The reported p-value is calculated over each of those iterations. The reported p-value is the mean ± standard error of the mean (SEM), where the mean value and SEM are calculated over 100 iterations. See Supplementary Figure S2 for a representative example of the procedure and Supplementary Table S2 for fitting details.

3 Results

3.1 Neuronal cultures

The overall activity of two to three-week-old neuronal cultures typically switches between high activity periods or “up states” and low activity periods or ‘down states‘ as shown in Figure 2. The up states—often also referred to as network bursts—occur quasi-periodically and involve the vast majority of all neurons. As Figure 2 also shows, their duration is quite regular as well. The activity during up states can be understood as a set of causal cascades of induced firings across the observed population of neurons in the culture, often modeled as a branching process (Beggs and Plenz, 2003). The cascade of neuronal activity is studied in the framework of neuronal avalanches as described in the Section 2. Previous studies have found that neuronal avalanches exhibit statistics of a branching process at or close to its critical point (Beggs and Plenz, 2003; Friedman et al., 2012). In the often-observed case of a mean-field branching process, the activity propagates on a tree-like network without feed-back loops and the avalanches follow a scale-free behavior, i.e.,

P(S)~S−τ (9)P(S)~S−τ (9)

P(T)~T−α (10)P(T)~T−α (10)

⟨T⟩(S)~Sγ (11)⟨T⟩(S)~Sγ (11)

where P is the probability distribution function (PDF) of the associated variable and τ = 1.5, α = 2.0, and γ = 0.5 are the critical mean-field exponents (Friedman et al., 2012). These critical exponents—whether they take on mean-field values or not—are necessarily related through the scaling relation

(α−1)/(τ−1)=γ−1, (12)(α−1)/(τ−1)=γ−1, (12)

which provides an additional test for the presence of critical behavior.

Figure 2

Figure 2. (A) Example of a raster plot for an experimental recording (top panel) and its corresponding normalized firing rate (bottom panel) are depicted. The up and down states are identified by thresholding the firing rate, where bins of data with firing rate > threshold are identified as up state and bins of data with firing rate less threshold are identified as down state (see section 2 for details). We define Tup and Tdown to measure the lengths of up and down states as visualized in the bottom panel. (B) Probability Density Functions (PDFs) of Tup, Tdown, Tup/ < Tup>, and Tdown/ < Tdown> are shown, where each curve represents one of the experimental recordings or simulations (see legends). (C) The mean values of Tup and Tdown for all the recordings and simulations are shown. Error bars correspond to 75 percentiles and the same color coding as panel (B) is used.

Here, we study the avalanche statistics of the up state and the down state separately. Due to the vastly different activity levels, the corresponding ISIs per neuron are also vastly different. For our experimentally recorded data the single cell ISIs are: ISI_{up state} = 217ms ± 43 ms and ISI_{down state} = 1, 065ms ± 302ms (mean ± SEM), as follows from Table 1. Nevertheless, we find similar statistical features across up and down states if the different activity levels are taken into account by choosing the time bins to define neuronal avalanches appropriately, being significantly larger for the down state. Examples of distributions of avalanche sizes and durations as well as the relation between sizes and durations are plotted in Figure 3 (see also Supplementary Figures S3–S5). Indeed, for both the up and down state in Figure 3, we recover Equations 9, 10, 11 and the exponents are close to the mean-field values. As Table 2 shows (see also Figure 3), this is often even true for the up state when we block all inhibitory connections by the application of a saturating concentration (50 μM) of picrotoxin (see Section 2 for details). Overall, we find that the majority of experimental recordings have an up state that is consistent with power-law statistics in their sizes over an extended range (see Supplementary Table S1). Such power-law behavior is slightly less prevalent in the down state (Table 2, Supplementary Table S1), but still prominent. Moreover, in all cases the scaling relation (Equation 12) holds within the statistical uncertainties (Supplementary Figure S6), consistent with critical behavior.

Figure 3

Figure 3. Example distributions of neuronal avalanche sizes and durations for up and down states of experimental data and simulated data are plotted for different temporal bin sizes, where the bin size of 1 corresponds to the temporal resolution of the recordings (5 ms) and the simulations (0.1 ms), respectively. (A) The up state is from experimental recording 3, the down state is from recording 1. (B) Simulation 1 is shown but all simulations exhibit almost identical avalanche statistics.

Table 2

Table 2. Critical exponents for all of the recordings and simulations.

Our analysis of the experimental data shows in particular that the statistical behavior of the neuronal avalanches in the down state can be statistically indistinguishable from the up state, especially if one focuses on τ, see Figure 3 and Table 2. Note that the different ranges in duration and sizes (see also Supplementary Table S1) are due to the shorter relative duration of the down states with respect to their corresponding ISI, as a comparison with the case of experimental recordings with continuously low steady-state activity (which can be interpreted as a system without an up state and instead being exclusively in a down state) (Yaghoubi et al., 2018) confirms.

As outlined in the Section 2, in this study we take advantage of a scaling analysis that also yields critical exponents denoted as β_S and β_T that can capture finite size behavior. This provides us with an additional tool for characterizing the dynamics of neuronal systems within the framework of neuronal avalanche statistics. As depicted in Figures 4A, B (see also Supplementary Figure S7), we studied the scaling properties of two types of distributions: (i) Power-law scaling which gives us β_S (β_T) and τ (α), and (ii) Exponential-like scaling that gives us only β_S (β_T), for avalanche sizes (durations). The estimate of the exponents becomes more reliable when the scaling collapse is obtained for a wider range of varying bin sizes. The summary of all estimated exponents for (i), along with the corresponding p-values for the reported critical exponents τ and α, (see Section 2 for details) for both up and down states, are plotted in Figure 4C. While the variation in τ and α is rather small, especially in the up state, this is not the case for β_S and, to a lesser degree, for β_T. This higher variability is also visible in Figure 4D, which displays the corresponding mean values and uncertainties for all exponents for both the up and down states. Figure 4E displays the mean values of β_S and β_T for distributions exhibiting exponential-like scaling (ii), which exclusively occurs in the down state, see Table 2.

Figure 4

Figure 4. (A) The power-law scaling collapse procedure for finding the exponents is visualized for an example experimental dataset (recording 3, up state, avalanche duration). For details of this procedure, see Section 2. (B) The exponential-like scaling collapse procedure for finding the scaling exponent [beta] for a sample experimental data is visualized (recording 3, down state, avalanche sizes). (C) Summary of the exponents for size and duration of all recordings during up and down states are presented here. Each point represents an experimental or a simulated dataset. We have only included recordings that exhibit power-law behavior over extended ranges as indicated in Table 2. p-values of the fitted power laws with exponents [tau] and [alpha], respectively, are shown in the bottom panels (mean + std of the ensemble of subpopulations). For each of the two columns, the same color scheme as the top panel is preserved. (D) Average of all power-law exponents for the data sets in (C). (E) Average of [beta]S and [beta]T for all data sets that show exponential-like behavior as the example in (B).

3.2 Model simulations

The similarity between the up state and the down state in terms of the neuronal avalanches for our experimental recordings suggests that either (a) there are different mechanisms of information processing in this neuronal system with different associated time scales, or (b) the concept of neuronal avalanches is not specific enough to establish insight into the differences between up and down states. To investigate this further, we employed several computational models that try to mimic the behavior and dynamics of neuronal cultures as described in Section 2.

In the literature, various computational models exist to describe up and down state dynamics. From mostly theoretical (Millman et al., 2010), to those mimicking the conditions of anesthetize animals (Holcman and Tsodyks, 2006), sleep (Bazhenov et al., 2002), or acute slices (Jercog et al., 2017; Camassa et al., 2022) among others. For the case of cultures from dissociated neurons, the up and down dynamics are better described with their own model. This is because in dissociated cultures, a down-to-up state transition can only occur after enough time has passed for the system to recover from a previous up state (Opitz et al., 2002). Once the system has recovered, spontaneous activations can be quickly amplified in a feedback loop due to the presence of recurrent connections, as described by Orlandi et al. (2013). On the other hand, transitions from up to down state are mediated primarily by the depletion of neurotransmitters caused by the high-frequency firings during the up state (Staley et al., 1998; Opitz et al., 2002). The models used here are able to reproduce most of these macroscopic observables as Figure 2 shows. However, when applying the same methodology to characterize the avalanche statistics as in the experimental data, there are important differences. In the simulations, the avalanches observed during the up states followed power-law distributions with slopes of τ≃−1.5 and α≃−1.8 for sizes and durations, respectively, see Figure 3B left and Table 2. These avalanche statistics also satisfied the scaling relation (Equation 12) within the statistical uncertainties (Supplementary Figure S6). During the down state, however, the avalanche statistics were always far from a power-law distribution, suggesting a sub-critical or exponential-like behavior instead (see Figure 3B right). As a result, the down state in the simulated data is absent from Figure 4D.

Although the scaling exponents τ and α were largely consistent during the up states with those observed experimentally (see Table 2), the scaling exponents β_s and β_T differed substantially from the experimental ones as follows from Figures 4D, E. The measured exponents across the different simulation conditions were robust, with little variability. The presence or absence of inhibitory connections, as well as heterogeneous cell populations (see Section 2), produced no significant differences in the avalanche statistics. Changes in connectivity strength, synaptic depression parameters (depression strength and recovery time constant), and the presence of other synaptic currents (NMDA), produced no changes in the reported exponents either (not shown). Similarly, although the metric properties of the network structure have a large impact on the presence of the up down transition (Orlandi and Casademunt, 2021), changes in connectivity, i.e., between metric-embedded to a random graph, also resulted in no changes in any of the exponents.

4 Discussion

In this work, we present for the first time the simultaneous characterization of neuronal avalanches in dissociated neuronal cultures during both up and down states. In several cultures, the avalanche statistics during both up and down states show critical behavior, as indicated by avalanche sizes and durations. In addition, the associated critical exponents satisfy the corresponding scaling relation (Equation 12), supporting the presence of critical behavior. These exponents are largely within the range of those of a critical mean-field branching process, in line with those previously reported across many preparations (Beggs and Plenz, 2003; Friedman et al., 2012; Yaghoubi et al., 2018). The exponents associated with finite-size scaling show some differences between up and down states. Yet, the variations are large and the statistics is rather limited such that it is difficult to make any definite statements. More importantly, we find that the scaling relation β_T = γβ_S (which follows directly from Equation 11) holds within the statistical uncertainties in almost all cases (Supplementary Figure S6), providing further evidence of criticality. The mechanisms by which these cultures can self-organize to maintain critical avalanche statistics across very different activity regimes are still unknown. If differences in the exponents were indeed present, this would further indicate that the two activity regimes are associated with distinct universality classes.

Of particular note is the fact that the temporal bins or time scales necessary to establish critical behavior based on neuronal avalanches are significantly different between up and down states. The time scales are closely related to the average interspike interval per neuron, being significantly larger for the down state. As the example of Figure 3A shows, the relative difference between the temporal bin sizes for neuronal avalanches showing critical behavior between up and down states largely corresponds to the relative difference in the average interspike interval per neuron (see Table 1). While in non-biological systems adjusting the temporal bin sizes over large ranges to evaluate avalanching or information spreading is quite common (Notarmuzi et al., 2022), the range considered for neuronal avalanches is rather limited, typically covering 2–20 ms (see, e.g., Pasquale et al., 2008; Friedman et al., 2012), which is comparable with the time scale of glutamatergic synaptic transmission (Ivenshitz and Segal, 2010). This is consistent with our time bins in the up state (see, for example, Figure 3A, Supplementary Figures S3, S5). Yet, for the down state the time bins necessary to recover neuronal avalanches exhibiting critical behavior are of the order of hundreds of ms (see, for example, Figure 3A, Supplementary Figure S4), suggesting communication mechanisms other than fast glutamatergic synaptic transmission, like those mediated by AMPA receptors, could be at play.

This hypothesis is consistent with our model findings. There has been extensive theoretical and modeling work to try and describe the up and down states as an emergent property of neuronal systems. However, only a few models can simultaneously reproduce the experimentally observed avalanche statistics and switch between up and down states. These include integrate-and-fire networks with structured connectivity through learning (Scarpetta and Candia, 2014); networks with scale-free connectivity coupled to a global field for the up and down switching (Lombardi et al., 2014); self-organized criticality around a saddle-node bifurcation (Millman et al., 2010); and a mesoscopic model on the verge of a synchronization transition (Santo et al., 2018). However, only Millman et al. (2010) treated the avalanche statistics of up and down states separately. They introduced a self-organized critical model with noise-driven, spontaneous transitions between up and down states. The transitions occur around a bifurcation such that the down state is subcritical and the up state is critical. Such a framework, however, is not applicable to our experimental system since its predictions are inconsistent with several experimental observations. Namely, (i) up states possessing characteristic durations, with a well-defined mean and variance (see, for example, Figure 2); (ii) the down states presenting a well-defined duration (usually called the interburst interval, IBI) that is correlated with the slow time-scale of synaptic depression (Cohen and Segal, 2009; Opitz et al., 2002); and (iii) the activity during the down states being possibly critical, as reported here and in Yaghoubi et al. (2018). Despite trying to take all of this into account, our computational model was still unable to reproduce some of the (new) experimental observations. Although the model was able to successfully capture many of the features of spontaneous activity in neuronal cultures (Orlandi et al., 2013), e.g., up and down transition statistics, distribution of states durations and their temporal correlations, up state exponents, etc., it could neither reproduce critical behavior observed during down states nor the values of the finite-size scaling exponents associated with the up states. For the model to simultaneously reproduce critical avalanche statistics during up and down states, its phase diagram would need to have critical points at two very different levels of average activity. The mechanism by which the model would be able to capture this is still unknown.

Conceptually, the fact that these cultures can dynamically transition between different levels of activity and still remain critical (or critical-like), suggests that we might have to move away from the traditional picture of a critical mean-field branching process, and even one of self-organization around a single critical point (Zapperi et al., 1995). Since the observed avalanche statistics need to be defined using substantially different time bins between the up and down states—which is related to their significantly different interspike intervals per neuron quantifying their activity levels—a model that can dynamically adapt the time-scales of synaptic integration (to maintain a constant effective probability of inducing firings from a neuron to their neighbors) could be a good candidate. Glial cells are known to modulate synaptic plasticity across different time-scales (Sancho et al., 2021) and have recently been shown to be involved in neural communication across long time-scales (Mu et al., 2019), hence being a likely candidate to support multi-scale critical dynamics. Investigating this hypothesis remains an exciting challenge for the future.

Data availability statement

The datasets presented in this study can be found in the following online repository: https://github.com/ComplexNSlab/multistateAvalanches.

Ethics statement

All animal care protocols were approved by the University of Calgary's Animal Care and Use Committee in accordance with the Canadian Council on Animal Care guidelines. All experimental procedures were carried out in accordance with these approved guidelines.

Author contributions

MY: Conceptualization, Investigation, Methodology, Writing – original draft, Writing – review & editing, Data curation, Formal analysis, Software, Visualization. JO: Conceptualization, Data curation, Investigation, Methodology, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing. MC: Funding acquisition, Methodology, Resources, Writing – review & editing. JD: Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Supervision, Validation, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This project was financially supported by an NSERC Discovery Grant (RGPIN/05221-2020) to JD, an NSERC Discovery Grant to MC, and the Eyes High Initiative of the University of Calgary (JO, JD).

Acknowledgments

We thank Y. Hernandez and A. Kipp for their help with culture preparation and maintenance.

Conflict of interest

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

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

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

References

Alvarez-Lacalle, E., and Moses, E. (2009). Slow and fast pulses in 1-d cultures of excitatory neurons. J. Comput. Neurosci. 26, 475–493. doi: 10.1007/s10827-008-0123-5

Bazhenov, M., Timofeev, I., Steriade, M., and Sejnowski, T. J. (2002). Model of thalamocortical slow-wave sleep oscillations and transitions to activated states. J. Neurosci. 22, 8691–8704. doi: 10.1523/JNEUROSCI.22-19-08691.2002

Beggs, J. M., and Plenz, D. (2003). Neuronal avalanches in neocortical circuits. J. Neurosci. 23, 11167–11177. doi: 10.1523/JNEUROSCI.23-35-11167.2003

Beggs, J. M., and Plenz, D. (2004). Neuronal avalanches are diverse and precise activity patterns that are stable for many hours in cortical slice cultures. J. Neurosci. 24, 5216–5229. doi: 10.1523/JNEUROSCI.0540-04.2004

Bellay, T., Klaus, A., Seshadri, S., and Plenz, D. (2015). Irregular spiking of pyramidal neurons organizes as scale-invariant neuronal avalanches in the awake state. Elife 4:e07224. doi: 10.7554/eLife.07224.019

Camassa, A., Galluzzi, A., Mattia, M., and Sanchez-Vives, M. V. (2022). Deterministic and stochastic components of cortical down states: dynamics and modulation. J. Neurosci. 42, 9387–9400. doi: 10.1523/JNEUROSCI.0914-22.2022

Chialvo, D. R. (2010). Emergent complex neural dynamics. Nat. Phys. 6:nphys1803. doi: 10.1038/nphys1803

Christensen, K., and Moloney, N. R. (2005). Complexity and Criticality. London: Imperial College Press Advanced Physics Texts. doi: 10.1142/p365

Cohen, D., and Segal, M. (2009). Homeostatic presynaptic suppression of neuronal network bursts. J. Neurophysiol. 101, 2077–2088. doi: 10.1152/jn.91085.2008

Colicos, M. A., Collins, B. E., Sailor, M. J., and Goda, Y. (2001). Remodeling of synaptic actin induced by photoconductive stimulation. Cell 107, 605–616. doi: 10.1016/S0092-8674(01)00579-7

Curic, D., Ivan, V. E., Cuesta, D. T., Esteves, I. M., Mohajerani, M. H., Gruber, A. J., et al. (2021). Deconstructing scale-free neuronal avalanches: behavioral transitions and neuronal response. J. Phys. Complex. 2:045010. doi: 10.1088/2632-072X/ac35b4

Deluca, A., and Corral, A. (2013). Fitting and goodness-of-fit test of non-truncated and truncated power-law distributions. Acta Geophys. 61, 1351–1394. doi: 10.2478/s11600-013-0154-9

Fatt, P., and Katz, B. (1952). Spontaneous subthreshold activity at motor nerve endings. J. Physiol. 117, 109–28. doi: 10.1113/jphysiol.1952.sp004735

Fernández-García, S., Orlandi, J. G., García-Díaz Barriga, G. A., Rodríguez, M. J., Masana, M., Soriano, J., et al. (2020). Deficits in coordinated neuronal activity and network topology are striatal hallmarks in Huntington's disease. BMC Biol. 18:58. doi: 10.1186/s12915-020-00794-4

Friedman, N., Ito, S., Brinkman, B. A. W., Shimono, M., DeVille, R. E. L., Dahmen, K. A., et al. (2012). Universal critical dynamics in high resolution neuronal avalanche data. Phys. Rev. Lett. 108:208102. doi: 10.1103/PhysRevLett.108.208102

Friedrich, J., Zhou, P., and Paninski, L. (2017). Fast online deconvolution of calcium imaging data. PLoS Comput. Biol. 13:e1005423. doi: 10.1371/journal.pcbi.1005423

Gireesh, E. D., and Plenz, D. (2008). Neuronal avalanches organize as nested theta- and beta/gamma-oscillations during development of cortical layer 2/3. Proc. Nat. Acad. Sci. 105, 7576–7581. doi: 10.1073/pnas.0800537105

Girotto, F., Scott, L., Avchalumov, Y., Harris, J., Iannattone, S., Drummond-Main, C., et al. (2013). High dose folic acid supplementation of rats alters synaptic transmission and seizure susceptibility in offspring. Sci. Rep. 3:1465. doi: 10.1038/srep01465

Haimovici, A., Tagliazucchi, E., Balenzuela, P., and Chialvo, D. R. (2013). Brain organization into resting state networks emerges at criticality on a model of the human connectome. Phys. Rev. Lett. 110:178101. doi: 10.1103/PhysRevLett.110.178101

Holcman, D., and Tsodyks, M. (2006). The emergence of up and down states in cortical networks. PLoS Comput. Biol. 2:e23. doi: 10.1371/journal.pcbi.0020023

Ivenshitz, M., and Segal, M. (2010). Neuronal density determines network connectivity and spontaneous activity in cultured hippocampus. J. Neurophysiol. 104, 1052–1060. doi: 10.1152/jn.00914.2009

Izhikevich, E. M. (2003). Simple model of spiking neurons. IEEE Trans. Neural Netw. 14:1569. doi: 10.1109/TNN.2003.820440

Izhikevich, E. M. (2007). Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. Cambridge, MA: MIT Press. doi: 10.7551/mitpress/2526.001.0001

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

Korchinski, D. J., Orlandi, J. G., Son, S.-W., and Davidsen, J. (2021). Criticality in spreading processes without timescale separation and the critical brain hypothesis. Phys. Rev. X 11:021059. doi: 10.1103/PhysRevX.11.021059

Levenstein, D., Buzsaki, G., and Rinzel, J. (2018). Excitable dynamics of nrem sleep: a unifying model for neocortex and hippocampus. bioRxiv 312587. doi: 10.1101/312587

Levina, A., and Priesemann, V. (2017). Subsampling scaling. Nat. Commun. 8:ncomms15140. doi: 10.1038/ncomms15140

Lombardi, F., Herrmann, H. J., Plenz, D., and Arcangelis, L. D. (2014). On the temporal organization of neuronal avalanches. Front. Syst. Neurosci. 8:204. doi: 10.3389/fnsys.2014.00204

Massobrio, P., de Arcangelis, L., Pasquale, V., Jensen, H. J., and Plenz, D. (2015). Criticality as a signature of healthy neural systems. Front. Syst. Neurosci. 9:22. doi: 10.3389/fnsys.2015.00022

Millman, D., Mihalas, S., Kirkwood, A., and Niebur, E. (2010). Self-organized criticality occurs in non-conservative neuronal networks during up states. Nat. Phys. 6, 801–805. doi: 10.1038/nphys1757

Moretti, P., and Muñoz, M. A. (2013). Griffiths phases and the stretching of criticality in brain networks. Nat. Commun. 4:2521. doi: 10.1038/ncomms3521

Mu, Y., Bennett, D. V., Rubinov, M., Narayan, S., Yang, C.-T., Tanimoto, M., et al. (2019). Glia accumulate evidence that actions are futile and suppress unsuccessful behavior. Cell 178, 27–43.e19. doi: 10.1016/j.cell.2019.05.050

Notarmuzi, D., Castellano, C., Flammini, A., Mazzilli, D., and Radicchi, F. (2022). Universality, criticality and complexity of information propagation in social media. Nat. Commun. 13:1308. doi: 10.1038/s41467-022-28964-8

O'Byrne, J., and Jerbi, K. (2022). How critical is brain criticality? Trends Neurosci. 45, 820–837. doi: 10.1016/j.tins.2022.08.007

Opitz, T., Lima, A. D. D., and Voigt, T. (2002). Spontaneous development of synchronous oscillatory activity during maturation of cortical networks in vitro. J. Neurophysiol. 88, 2196–2206. doi: 10.1152/jn.00316.2002

Orlandi, J. G., and Casademunt, J. (2017). Noise focusing in neuronal tissues: symmetry breaking and localization in excitable networks with quenched disorder. Phys. Rev. E 95:052304. doi: 10.1103/PhysRevE.95.052304

Orlandi, J. G., and Casademunt, J. (2021). Stochastic quorum percolation and noise focusing in neuronal networks. Europhys. Lett. 133:48002. doi: 10.1209/0295-5075/133/48002

Orlandi, J. G., Fernández-García, S., Comella-Bolla, A., Masana, M., Barriga, G. G.-D., Yaghoubi, M., et al. (2017). NETCAL: an interactive platform for large-scale, NETwork and population dynamics analysis of calcium imaging recordings. ZENODO Data Repository. doi: 10.5281/zenodo.1119025

Orlandi, J. G., Soriano, J., Alvarez-Lacalle, E., Teller, S., and Casademunt, J. (2013). Noise focusing and the emergence of coherent activity in neuronal cultures. Nat. Phys. 9:nphys2686. doi: 10.1063/1.4776497

Pasquale, V., Massobrio, P., Bologna, L., Chiappalone, M., and Martinoia, S. (2008). Self-organization and neuronal avalanches in networks of dissociated cortical neurons. Neuroscience 153, 1354–1369. doi: 10.1016/j.neuroscience.2008.03.050

Penn, Y., Segal, M., and Moses, E. (2016). Network synchronization in hippocampal neurons. Proc. Nat. Acad. Sci. 113, 3341–3346. doi: 10.1073/pnas.1515105113

Petermann, T., Thiagarajan, T. C., Lebedev, M. A., Nicolelis, M. A. L., Chialvo, D. R., Plenz, D., et al. (2009). Spontaneous cortical activity in awake monkeys composed of neuronal avalanches. Proc. Nat. Acad. Sci. 106, 15921–15926. doi: 10.1073/pnas.0904089106

Ponce-Alvarez, A., Jouary, A., Privat, M., Deco, G., and Sumbre, G. (2018). Whole-brain neuronal activity displays crackling noise dynamics. Neuron 100, 1446–1459.e6. doi: 10.1016/j.neuron.2018.10.045

Press W. H. Teukolsky S. A. Vetterling W. T. Flannery B. P. (2007). Numerical Recipes 3rd Edition: The Art of Scientific Computing. Cambridge: Cambridge University Press.

Priesemann, V., Valderrama, M., Wibral, M., and Quyen, M. L. V. (2013). Neuronal avalanches differ from wakefulness to deep sleep-evidence from intracranial depth recordings in humans. PLoS Comput. Biol. 9:e1002985. doi: 10.1371/journal.pcbi.1002985

Rabus, A., Curic, D., Ivan, V. E., Esteves, I. M., Gruber, A. J., Davidsen, J., et al. (2023). Changes in functional connectivity preserve scale-free neuronal and behavioral dynamics. Phys. Rev. E 108:L052301. doi: 10.1103/PhysRevE.108.L052301

Ribeiro, T. L., Ribeiro, S., Belchior, H., Caixeta, F., and Copelli, M. (2014). Undersampled critical branching processes on small-world and random networks fail to reproduce the statistics of spike avalanches. PLoS ONE 9:e94992. doi: 10.1371/journal.pone.0094992

Saito, M., and Matsumoto, M. (2013). Variants of mersenne twister suitable for graphic processors. ACM Trans. Math. Softw 39:12. doi: 10.1145/2427023.2427029

Sanchez-Vives, M. V., and McCormick, D. A. (2000). Cellular and network mechanisms of rhythmic recurrent activity in neocortex. Nat. Neurosci. 3, 1000–1027. doi: 10.1038/79848

Sancho, L., Contreras, M., and Allen, N. J. (2021). Glia as sculptors of synaptic plasticity. Neurosci. Res. 167, 17–29. doi: 10.1016/j.neures.2020.11.005

Santo, S. D., Villegas, P., Burioni, R., and Muñoz, M. A. (2018). Landau-ginzburg theory of cortex dynamics: scale-free avalanches emerge at the edge of synchronization. Proc. Natl. Acad. Sci. 115:E1356. doi: 10.1073/pnas.1712989115

Scarpetta, S., and Candia, A. D. (2014). Alternation of up and down states at a dynamical phase-transition of a neural network with spatiotemporal attractors. Front. Syst. Neurosci. 8:88. doi: 10.3389/fnsys.2014.00088

Staley, K. J., Longacher, M., Bains, J. S., and Yee, A. (1998). Presynaptic modulation of ca3 network activity. Nat. Neurosci. 1, 201–209. doi: 10.1038/651

Tagliazucchi, E., Balenzuela, P., Fraiman, D., and Chialvo, D. R. (2012). Criticality in large-scale brain fmri dynamics unveiled by a novel point process analysis. Front. Physiol. 3:15. doi: 10.3389/fphys.2012.00015

Taub, H., and Schilling, D. L. (1977). Digital integrated electronics. (No Title).

Williams-García, R. V., Beggs, J. M., and Ortiz, G. (2017). Unveiling causal activity of complex networks. Europhys. Lett. 119:18003. doi: 10.1209/0295-5075/119/18003

Yaghoubi, M., de Graaf, T., Orlandi, J. G., Girotto, F., Colicos, M. A., and Davidsen, J. (2018). Neuronal avalanche dynamics indicates different universality classes in neuronal cultures. Sci. Rep. 8:3417. doi: 10.1038/s41598-018-21730-1

Yu, S., Ribeiro, T. L., Meisel, C., Chou, S., Mitz, A., Saunders, R., et al. (2017). Maintained avalanche dynamics during task-induced changes of neuronal activity in nonhuman primates. Elife 6:e27119. doi: 10.7554/eLife.27119.012

Zapperi, S., Lauritsen, K. B., and Stanley, H. E. (1995). Self-organized branching processes: mean-field theory for avalanches. Phys. Rev. Lett. 75, 4071–4074. doi: 10.1103/PhysRevLett.75.4071