Long-term stability of avalanche scaling and integrative network organization in prefrontal and premotor cortex

Ongoing neuronal activity in the brain establishes functional networks that reflect normal and pathological brain function. Most estimates of these functional networks suffer from low spatiotemporal resolution and indirect measures of neuronal population activity, limiting the accuracy and reliability in their reconstruction over time. Here, we studied the stability of neuronal avalanche dynamics and corresponding reconstructed functional networks in the adult brain. Using chronically implanted high-density microelectrode arrays, the local field potential (LFP) of resting-state activity was recorded in prefrontal and premotor cortex of awake nonhuman primates. Avalanche dynamics revealed stable scaling exhibiting an inverted parabolic profile and collapse exponent of 2 in line with a critical branching process over many days and weeks. Functional networks were based on a Bayesian-derived estimator and demonstrated stable integrative properties characterized by nontrivial high neighborhood overlap between strongly connected nodes and robustness to weak-link pruning. Entropy-based mixing analysis revealed significant changes in strong link weights over weeks. The long-term stability in avalanche scaling and integrative network organization in the face of individual link weight changes should support the development of noninvasive biomarkers to characterize normal and abnormal brain states in the adult brain.


INTRODUCTION
Ongoing neuronal activity in the mammalian brain gives rise to distinct functional networks that reflect normal and pathological brain function (Bullmore & Sporns, 2009;Fox & Raichle, 2007;Raichle, 2015). Empirical measures of these networks, however, are often limited in the precision with which they reflect the organization of neuronal population activity and a lack of longitudinal follow-up. Networks are reconstructed either from indirect nonneuronal activity at low temporal resolution using functional magnetic resonance imaging (fMRI), or via more direct neuronal activity measures such as magnetoencephalography (MEG) and electroencephalography (EEG), which suffer predominantly from low spatial resolution (Bassett & Sporns, 2017). Accordingly, there is limited knowledge about the stability and nature of the spatiotemporal dynamics and corresponding functional connectivity of local neuronal populations that give rise to fluctuations in ongoing cortical activity.
In contrast to macroscopic and indirect methods, neuronal activity measured at singlecell resolution could in principle allow for precise and robust spatiotemporal reconstruction of functional networks in cortex. For example, in macaque cortex single-unit activity can be recorded over many weeks (Nicolelis et al., 2003) and identifies firing statistics of neurons over days (Dickey et al., 2009;Jackson & Fetz, 2007). Advanced approaches may involve tetrodes, automated classifiers (Tolias et al., 2007), or second-order statistics as classifiers (Fraser & Schwartz, 2011). Even advanced single-unit approaches, however, limit robust tracking of activity to small groups of individual neurons. Chronic expression of genetically encoded calcium indicators can increase the neuronal count and identify orientationselective neurons over months (Li et al., 2017), yet at low temporal resolution. The local field potential (LFP), alternatively, provides access to an intermediate spatial resolution of cortical activity that is between the macroscopic scale of fMRI, MEG and EEG, and singlecell accuracy. It robustly measures local neuronal population activity with high temporal and reasonable spatial resolution of several hundred micrometers (e.g., Bellay et al., 2021;Katzner et al., 2009) over many weeks. In nonhuman primates (NHPs), the LFP has been shown to correlate well with local unit-activity (Donoghue et al., 1998;Petermann et al., 2009;Rasch et al., 2008) and offers robust brain-machine interface performance (Mehring et al., 2003).
It is now well established that the ongoing LFP in superficial layers of nonhuman primates organizes into spatiotemporal activity clusters whose sizes and durations distribute according to power laws (Miller et al., 2019;Shew et al., 2009;Yu et al., 2017), the hallmark of neuronal avalanches (Beggs & Plenz, 2003). The spatiotemporal correlations underlying neuronal avalanches give rise to functional "integrative" networks, in which the topology and weights of links organize such that end nodes of strong links are more likely to have many common neighbors, as quantified by a high correlation between the link weight and the link clustering coefficient (Pajevic & Plenz, 2012; cf. edge-clustering coefficient in Radicchi et al., 2004), while being robust to the loss of weak links, in terms of maintaining the advantages of a small-world topology (Watts & Strogatz, 1998). This specific arrangement describes many complex networks ranging from genetic networks to social, that is, "friendship," networks (Pajevic & Plenz, 2012). Importantly, integrative networks can form using local plasticity rules that operate on neuronal avalanches (Alstott et al., 2015) in which weak links provide a large basin for plastic changes as they can be removed or rewired, permitting flexibility, while preserving the small-world network property. The stability over days and weeks of cortical avalanche dynamics and corresponding integrative networks has not been explored for cortex. Size and duration distributions of avalanches do not consider the spatiotemporal organization of avalanches, thus multiple functional connectivities could give rise to similar power law Functional connectivity: The activity between nodes can be positively or negatively correlated, resulting in a corresponding functional connectivity link that is either positive or negative. Often noted as link strength or link weight.
Local field potential: A neuronal signal that typically reflects the summed activity within a local neighborhood. It allows for robust tracking of local neuronal population activity, yet the origins of the contributing neural elements are ambiguous.
Neuronal avalanches: Scale-invariant spatiotemporal activity in ongoing and evoked activity of cortex. Sizes and durations of avalanches follow a power law and the mean size scales with avalanche duration according to a power law with slope of 2.
Small-world topology: A network organization where the local clustering between nodes is high while the average number of links to reach between nodes remains low, typically maintained by a low number of long-range links. statistics. Likewise, integrative networks allow for the reorganization of individual links while maintaining global integrative properties (Pajevic & Plenz, 2012).
Here we recorded ongoing LFP activity in prefrontal and premotor cortex of NHPs over weeks using chronically implanted high-density microelectrode arrays. We demonstrate stable avalanche dynamics with power law statistics, the temporal avalanche profile of an inverted parabola and corresponding scaling exponent of 2, which is the expected value for critical branching process dynamics (di Santo et al., 2017;Miller et al., 2019;Sethna et al., 2001). We further demonstrate that these avalanches give rise to nontrivial integrative functional networks despite progressive reorganization of the link strengths as identified using normalized count (NC) estimators (Pajevic & Plenz, 2009) and the "entropy of mixing" analysis (see Material and Methods). Our results suggest that avalanche dynamics and corresponding integrative network organization identify a robust cortical state in the adult brain, which should inform dynamical models of brain function and might allow for the early identification of pathological brain states.

RESULTS
The spatiotemporal organization of the ongoing LFP (1-100 Hz) was monitored with chronically implanted high-density microelectrode arrays (~10 × 10 electrodes; interelectrode distance 400 μm) in premotor (PM, n = 2 arrays) and prefrontal cortex (PF, n = 4 arrays) of three macaque monkeys (Macaca mulatta; K, V, N; Figure 1A, 1B). The animals were awake and seated in a monkey chair during the recording sessions but were not given any specific stimulus or task.

Stability of Average Pairwise Correlation and Avalanche Dynamics in Prefrontal and Premotor Cortex
The pairwise Pearson's correlation, R, is a standard measure to quantify the correlation in neuronal activity between distant cortical sites. We found that for the continuous LFP, the distribution of R, obtained from~5,000 electrode comparisons per array, was broad yet similar across days ( Figure 1D, 1E). The average correlation, <R>, while different between arrays, was constant over many days and weeks (Figure 1F, 1G; linear regression fit of 0.4 ± 0.7 × 10 −2 for all arrays; see Table 1 for individual arrays).
We next explored the stability of neuronal avalanche dynamics by discretizing the LFP from ongoing activity in the NHPs Yu et al., 2017). We previously demonstrated that negative deflections in the continuous LFP signal (nLFPs) track the local neuronal population activity and correlate with extracellular unit firing (Bellay et al., 2021;Petermann et al., 2009). Accordingly, discrete nLFPs were defined by crossing a threshold of −2 SD and extracting their peak times ( Figure 2A). The nLFPs were then grouped into avalanches at temporal resolution Δt of the mean time interval of successive nLFPs on the full array, <IEI>, which ranged between 3 and 4 ms ( Figure 2B; see Material and Methods). We found that the probability density distribution of avalanche size, which typically ranges from size 1 to the number of electrodes for each array , exhibited power laws for each day of recording (Figure 2C, 2D; Δt = <IEI>, truncated power law vs. exponential, LLR 1, p < 0.005; or vs. log-normal, LLR 1, p < 0.005; constant across all recording days). The stability of the power law in size was associated with an average slope α = 1.32 ± 0.05 for all arrays (Table 1), which is close to the expectations for the size exponent of a critical branching High-density microelectrode array: An array of microelectrodes that samples neuronal activity at discrete, spatial sampling points in accordance to the layout of the array.
Critical branching process: A dynamical process that describes the evolution of descendants from ancestors. In its simplest form, it is expressed in discrete generations (that is, discrete time steps), with no memory (in other words, ancestors contribute only to the next generation), and is stationary (that is, the average number of descendants per generation is constant). In a critical branching process, one ancestor on average leads to exactly one descendant. The resulting branched chains of ancestors/descendants over generations exhibit scale-invariant properties as found for neuronal avalanches in the brain tissue.
Entropy of mixing: An entropy measure to quantify the mixing over time of members that were initially separated into different categories.

Integrative network organization:
A network with small-world topology in which nodes connected by strong links also exhibit a relatively large common neighbor node overlap. This organization of strong links is also found in so-called friendship networks, in which agents with a high percentage of common friends communicate often.
Ongoing activity: The brain is spontaneously active even in the absence of specific sensory inputs or motor outputs. The specific organization of spontaneous activity might arise from numerous internal processes, such as posture maintenance, as well as internal cognitive processes, such as attention. process (Figure 2E, 2F;Harris, 1963). These results demonstrate the presence of avalanche dynamics in all NHPs for up to 2 weeks in PM and up to 10 weeks in PF.

Long-Term Stability in the Temporal Profile of Avalanches and Corresponding Scaling Collapse
The temporal profile of a neuronal avalanche describes how avalanches unfold in time by initiating locally, growing in magnitude, and eventually contracting spatially before ending. Recent work (Miller et al., 2019) confirmed the prediction that the temporal profile of neuronal avalanches obeys the characteristic motif of an inverted parabola and that empirically observed mean avalanche profiles collapse onto this scale-invariant motif regardless of avalanche duration. For this analysis, the temporal resolutions were set to Δt = 1 ms and Δt = 30 ms respectively, which demonstrated stable power law size distributions and expected change in slope α (Table 1; 1 ms: all LLR > 10 4 , all p 0.005; 30 ms: all LLR > 10, all p 0.005) as reported previously (Beggs & Plenz, 2003;Petermann et al., 2009). These resolutions allow for the collapse of temporal avalanche profiles in the presence of transient -oscillations in ongoing PF and PM activity of NHPs (Miller et al., 2019).
We explored the stability of these scaling predictions over time, finding that the characteristic motif for mean profiles of avalanche lifetime, L, with L = 3Δt − 5Δt did trace inverted parabolas for fine (Δt = 1 ms, Figure 3A, 3B) and coarse (Δt = 30 ms, Figure 3C, 3D) temporal resolutions (for statistical comparison see Table 2 for individual arrays). Importantly, the scaling exponent χ producing these parabolic, collapsed motifs was close to the value of 2 expected for a critical branching process (di Santo et al., 2017;Miller et al., 2019;Sethna et al., 2001) and did not change over time ( Figure 3B, 3D; for Δt = 1 ms, χ = 1.93 ± 0.14, linear regression fit of 1 × 10 −2 ± 0.2; for Δt = 30 ms, χ = 2.01 ± 0.24, linear regression slope of 2 ×10 −2 ± 0.04; see Table 3 for individual arrays). In contrast, time-shuffled data reliably produced flattened motifs ( Figure 3E, 3F) with a scaling exponent close to 1, as found for noncritical dynamics, in which sizes grow linearly with duration (Villegas et al., 2019). Table 1. Summary of average estimates for the mean Pearson's pairwise correlation coefficient, <R>, the exponent α of avalanche size distributions for each nonhuman primate, and array at temporal resolutions of Δt = 1 ms and Δt = 30 ms, respectively. Exponents were derived from maximum likelihood estimates. Stability of α over all recording days was quantified by the linear regression slope within 95% confidence limits (CfL)
reconstruction (Pajevic & Plenz, 2009). This method uses a shortcut to Bayesian estimation of network connections, making it approximate, but also computationally efficient, and hence applicable to reconstructing relatively large networks from lengthy event train recordings ( Figure 4A). It also reduces the influence of indirect common input correlations and has been found in simulations to accurately and efficiently reconstruct the directed networks and weights from simulated point-process dynamics (Pajevic & Plenz, 2009). An example of the directed connectivity matrix for a single recording session derived from spontaneous avalanche activity is shown in Figure 4B. We found that for all arrays and recording sessions, the weights of the directed functional connections were not particularly heavy-tailed and distributed between an exponential and log-normal function over more than 3 orders of magnitude ( Figure 4C-4E). Having reconstructed the directed, weighted links in our avalanche-derived functional networks, we proceeded to calculate several graph-theoretical properties. Integrative networks are characterized by a positive correlation, R CL , between the link cluster coefficient, C L , and the corresponding link weight, w, between two nodes ( Figure 5A; Pajevic & Plenz, 2012). This positive correlation demonstrates that strongly connected nodes have more common neighbors compared with weakly connected nodes ( Figure 5B). We employed degree-sequence preserved randomization to subtract potential contributions to link clustering from the node degree distribution and report here the resulting excess link clustering, ΔC L (see Material and Table 2. Summary of scaling function estimate for the temporal profile of neuronal avalanches. A parabolic fit (Δ parab ) is compared with a semicircle fit (Δ semi ) at temporal resolutions Δt = 1 ms and Δt = 30 ms for each nonhuman primate and array and its corresponding phase-randomized profile (random) at Δt = 1 ms  Table 3. Summary of estimates for the scaling exponent χ to optimally collapse temporal avalanche profiles at temporal resolutions Δt = 1 ms and Δt = 30 ms for each nonhuman primate and array. The stability of χ across all recording days was quantified by the linear regression slope reported and 95% confidence limits. Random: Phase-randomized control of temporal avalanche profile Scaling exponent χ (Δt = 1 ms) Methods). The monotonic increase in ΔC L with lowest-to-highest weight rank ( Figure 5C) clearly identifies integrative functional network organization in both PM and PF that was stable over time ( Figure 5D; linear regression slope of 0.005 +/− 0.01; 95% confidence limit; see Table 4 for individual arrays).   (Pajevic & Plenz, 2012) showed that integrative clustering motif was robust to weak-link pruning (solid lines; error bars denote mean ± SD over all recordings for an array), while the clustering "backbone" of the network quickly degraded when strong links were pruned first (dotted lines: color coded as in D). C − E, Δt = 1 ms. Table 4. Functional network stability based on NC reconstruction algorithm at Δt = 1. Shown is the stability of R CL over time (see also Figure 5D). Pruning analysis is an efficient approach to study the weight organization of networks (Onnela et al., 2007;Pajevic & Plenz, 2009). We found that pruning from the bottom, that is, pruning weakest links progressively, had negligible effect on the network clustering properties. This was evidenced by a near constant value of the excess node clustering, ΔC, until the vast majority (80-90%, Figure 5E; solid lines) of links were removed. In contrast, pruning from the top, that is, removing the strongest links first, rapidly decreased ΔC, a result consistently found for all NHPs ( Figure 5E). We conclude that the integrative property arises from a "backbone" of strong links that is robust to weak-link pruning, yet quickly degrades when strong links are removed.

Entropy of Mixing Identifies Link Progression Over Many Weeks
Given the stability observed over weeks in the global network dynamics, that is, <R>, α, χ, and integrative network organization, that is, R CL , we next sought to analyze in more detail the stability at the local level of individual link weights. Specifically, we examined whether the stability identified at the global level simply arises from stability at the local level. To quantify changes at the local level, we categorize the links, such as weak, medium, and strong, and analyze whether links switch between categories over time, for example from weak to medium Pruning analysis: The systematic change in network properties when links or nodes are removed systematically. This is contrasted by a perturbation analysis where network components are removed randomly to study something such as network failures. Figure 6. Entropy of mixing analysis demonstrates progressive changes in strong links over time for some networks. Entropy of mixing (H M ) analysis of the link weights obtained with NC reconstructions from Δt = 1 ms raster recordings. H M is plotted for the "top" half (strongest 50% of links; black curves) and "bottom" half (weakest 50% of links; green curves) versus the days after the initial recording. The black and green dashed lines indicate corresponding H M obtained when the labels were fully randomized, while the red dashed line indicates the theoretical maximum of H M . The six panels represent six different arrays from the three monkeys shown in two rows. (A) Arrays for which statistically significant progression of mixing was obtained for "top" links. (B) Arrays for which no statistically significant progression was observed. The "top" half of the links consistently showed lower H M than the "bottom" half of the links, which is mainly due to the smaller reconstruction error. The Day 1 points are not plotted because, by default of our labeling, the initial value for H M is 0. or vice versa. As a measure of mixing, we calculate the entropy of mixing, H M (see Material and Methods). Two main factors contribute to potential changes in H M , which originate from changes in link categories in a network from one day to another. First, independent random changes on each day will result in increased mixing and H M , but not progressively over time. For example, errors in network reconstruction alone will result in a constant H M over many days if the underlying network is not changing. Second, link weights could gradually evolve over days and weeks, and potentially change their categories, which is now progressively mixed and therefore is expected to lead to a monotonic increase in H M until the maximum in mixing is reached. Because reconstruction errors are higher for weak links compared with strong links, we quantified potential changes in link weight from NC-reconstructed networks separately for the strongest (top 50%) and weakest (bottom 50%) links. This separation is also motivated by our results for top-versus bottom-pruning, which demonstrated that networks are less robust to link changes in the top half of the link weight distribution. For each half, we partitioned link weights further into strength quintiles, yielding five categories, and then studied whether members from different categories will mix over time. In Figure 6, we show the data for the top half and bottom half and compare those curves to H M obtained if the category labels were fully randomized at each step and also present the theoretical maximum for H M (see also Material and Methods). Because all successive distributions used the link labels derived from Day 1 (making H M = 0), the average entropy of mixing was plotted beginning from the second recording for each NHP.
For all arrays, we found that H M for the bottom half of links reaches close to the maximal entropy by the following day, demonstrating that weak links fluctuate randomly. On the other hand, we found progressive and statistically significant increases of H M over subsequent days for strong links in three out of six arrays ( Figure 6A, 6B; see also Table 4). These findings demonstrate stability in avalanche dynamics and integrative network organization even though the network undergoes significant weight reorganization at the individual link level.

DISCUSSION
Fluctuations in ongoing neuronal activity give rise to functional networks that allow for the categorization of healthy and pathological brain states (Bullmore & Sporns, 2009;Fox & Raichle, 2007;Raichle, 2015), yet the dynamical and graph-theoretical markers underlying these ongoing fluctuations have been challenging to identify. By recording the LFP in prefrontal and premotor cortex of NHPs over weeks at high spatiotemporal resolution, we established two robust markers of ongoing activity: the dynamical category of neuronal avalanches and the graphtheoretical category of integrative networks. Avalanche dynamics revealed power law statistics, parabolic temporal profiles and a scaling exponent of 2 in support of critical branching process dynamics (di Santo et al., 2017;Miller et al., 2019;Sethna et al., 2001). We further demonstrate that avalanches form integrative networks despite reorganization of individual link weights as identified using NC estimators and an entropy of mixing approach (Onnela et al., 2007;Pajevic & Plenz, 2009, 2012. Our results suggest that avalanche dynamics and corresponding integrative networks identify a robust state of the awake frontal cortex. This work should be impactful for computational studies of cortex function and for translational research into the identification of brain dynamics in healthy subjects and corresponding deviation in brain disorders. Recent analysis of evoked cortical activity in nonhuman primates during behavior demonstrated the preservation of neuronal avalanche dynamics in the face of transient activity changes, specifically the nLFP rate (Yu et al., 2017). Given the robustness of integrative network properties observed in combination with avalanche dynamics in vitro (Pajevic & Plenz, 2009, 2012 as well as in vivo (the present study), we would expect integrative network properties to also persist during behavioral tasks.

Stability of LFP-Based Avalanches and Relation to Neuronal Activity
Technical advances in identifying long-term robustness in neuronal activity using cellular resolution approaches (Dickey et al., 2009;Fraser & Schwartz, 2011;Jackson & Fetz, 2007;Li et al., 2017;Nicolelis et al., 2003;Tolias et al., 2007) typically suffer from spatial and temporal subsampling (Levina & Priesemann, 2017;Petermann et al., 2009;Ribeiro et al., 2010;Ribeiro et al., 2014) limiting proper identification of avalanche dynamics and corresponding reconstruction of functional connectivity. The LFP, on the other hand, is a robust neuronal population signal, less prone to spatial subsampling and at high temporal resolution contains local spike information particularly when recorded from superficial layers (Donoghue et al., 1998;Mehring et al., 2003;Petermann et al., 2009;Rasch et al., 2008). Indeed, we recently demonstrated in nonhuman primates and in vitro slices that single neurons selectively participate in expansive, repeated LFP avalanches (Bellay et al., 2021), which suggest that the stability observed in the present study for LFP-based avalanches might extend to the single neuron level.
Experimental manipulations support a homeostatic regulation of avalanche dynamics in cortex. Transient pharmacological perturbation in vitro (Plenz, 2012) or sensory deprivation in vivo (Ma et al., 2019) initially abolishes avalanche dynamics, which is followed by a recovery over several days. Developmental findings further support an autonomous regulation of neuronal avalanches in isolated cortex in vitro in the absence of sensory input (Pasquale et al., 2008;Stewart & Plenz, 2007;Tetzlaff et al., 2010). Identifying the homeostatic regulation of avalanche dynamics and integrative networks will require precise perturbation approaches (Chialvo et al., 2020) in combination with advanced cellular techniques suitable for long-term monitoring such as ratiometric genetically encoded calcium indicators, which naturally compensate for cell-intrinsic expression variability (Lutcke et al., 2013).

NC Reconstruction of Integrative Networks
Reconstruction of our directed, weighted functional networks utilized a Bayesian-derived NC approach, which is easily scalable to large networks (Pajevic & Plenz, 2009, 2012 and reduces correlations from common input. Network simulations demonstrated ground-truth reconstruction for various small-world topologies from subcritical, supercritical, and critical dynamics, the latter in line with reconstructions from observed avalanche activity (Pajevic & Plenz, 2009).
Here we used significant events in LFP fluctuations, that is, nLFPs, for the reconstruction of functional connectivity, which differs from the more common ansatz based on continuous time series. Accordingly, our approach, which selectively reconstructs synchronization dynamics in the form of avalanches, introduces a timescale Δt to which the NC reconstruction is sensitive. In order to identify the spatiotemporal spread of avalanches, the use of a microelectrode array introduces a spatial discretization, which consequently enforces a discretization Δt in time. It was found empirically that Δt should be linked to the average propagation velocity of neuronal activity in the network, allowing Δt to be approximated by <IEI> (Beggs & Plenz, 2003;Petermann et al., 2009). For Δt smaller than <IEI>, errors of prematurely terminating avalanches increase, whereas errors of concatenating successive avalanches increase for Δt larger than <IEI>. An approximate balance regarding these two errors is struck for Δt = IEI. With respect to NC reconstruction, increasing Δt will increase the count in nLFP occurrences on the array per time step, resulting in an increase in the computational load and increased error in estimating priors for potentially activated links. Accordingly, we chose a Δt slightly smaller than IEI for a conservative and balanced estimate of the functional connectivity based on avalanche propagation.
Integrative networks differ from networks in which the weight is positively correlated with the node degree (Barrat et al., 2004;Bianconi, 2005;Pajevic & Plenz, 2012). Our shuffling correction further demonstrates that the organization of weighted links in integrative networks cannot be simply explained by the link-degree distribution and in which strong links are part of a continuous, in our case near exponential, distribution (cf. Figure 4D). Integrated networks are also not a simple consequence of critical branching processes generating avalanches, which can be realized in many topologies and architectures (Pajevic & Plenz, 2009, 2012.

Robust Avalanche Scaling of χ = 2 Over Many Weeks
Heterogeneous activity in the form of neuronal avalanches has been the hallmark of ongoing (Beggs & Plenz, 2003;Bellay et al., 2015) and evoked (Yu et al., 2017) brain activity in superficial layers of the cerebral cortex. The scale-invariant hallmark of avalanches has been found to capture ongoing neuronal activity in nonhuman primates at the mesoscopic level in the LFP (Klaus et al., 2011;Petermann et al., 2009;Yu et al. 2014;Yu et al. 2017;Yu et al., 2011) and macroscopic activity measured with EEG (Meisel et al., 2013), MEG , and fMRI Tagliazucchi et al., 2012). The observed power law distributions in cluster size and duration support the idea that the cortex operates close to a critical state at which networks gain numerous advantages in information processing (Beggs & Plenz, 2003;Gautam et al., 2015;Kinouchi & Copelli, 2006;Shew et al., 2009, Shew et al., 2011Yang et al., 2012).
Our present demonstration of an inverted-parabolic avalanche shape, collapse exponent of 2, and size distribution slope >−2 is in line with predictions that the unfolding of avalanches in the adult brain is governed by a critical branching process (for details, see Miller et al., 2019). This profile demonstrates that the initial spatial unfolding of an avalanche in time is truly expansive and should not be viewed as being dominated by a few strong links. Recently, avalanches with a scaling factor of 2 have been described for whole-brain neuronal activity in the zebra fish larvae (Ponce-Alvarez et al., 2018); however, the corresponding size and duration exponent were significantly steeper compared with reports for the mammalian cortex, suggesting a different dynamical model for the zebra fish.
Connectome analyses on a timescale of seconds have identified transient fluctuations in functional connectivity reflecting distinct states (Hutchison et al., 2013;Liu & Duyn, 2013;Zalesky et al., 2014), which based on computational modeling were found to be intermittent, a signature of metastability (Deco et al., 2017;Hansen et al., 2015). Because our analysis was based on nLFPs and not on continuous time signals, we were unable to obtain NCreconstructed networks for shorter time periods than 15-30 min. In fact, our entropy of mixing analysis allowed for an assessment of link weight progression in the face of reconstruction errors due to limited data availability. Transient fluctuations in functional connectivity estimates are not in contrast to the reorganization of links reported in our study. Short-term fluctuations in the face of long-term stability are hallmarks of critical dynamics Tagliazucchi et al., 2012).

Potential Advantages of Avalanche Dynamics and Integrative Networks for Cortex Function
Anatomically, the mammalian cortex is a sheet of local, functional modules that dynamically combine to support complex brain functions (Braitenberg & Schüz, 1991;Honey et al., 2007).
Local connectivity ensures diverse local operations, whereas long-range connections support global coordination (e.g., Bullmore & Sporns, 2009;Kirst et al., 2016). These structural hallmarks are dynamically realized in neuronal avalanches, which support a scale-free and selective organization of neuronal synchronization over all distances. They are graph-theoretically embedded in integrative networks, which maintain a high level of node clustering until about 90% of the weakest links are pruned, that is, about 10% of strong links are used for long-range connections. Intriguingly, simulations identified a particular local learning rule that operates on activity cascades such as avalanches to build integrative networks from random networks (Pajevic & Plenz, 2012). The rule increases link weights at locations of recent cascade failure, thereby facilitating the unfolding of future avalanches (Alstott et al., 2015). This mechanism opens bottlenecks in the network and further supports activity propagation over long distances. We therefore hypothesize that avalanche dynamics in combination with integrative network organization are beneficial for local operational diversity while supporting longrange, selective propagation of information.

Animal Procedures
All animal procedures were conducted in accordance with National Institutes of Health guidelines and were approved by the Animal Care and Use Committee of the National Institute of Mental Health. Three adult NHPs (Macaca mulatta; 1 male, 2 females; 7-8 years old) received two chronic implantations of high-density 96-microelectrode arrays each (Blackrock Microsystems; 4 × 4 mm 2 ; 400 μm interelectrode distance; 10 × 10 grid with corner grounds). To direct recordings towards superficial cortical layers II/III, electrode shanks of 0.6-mm length were used in prefrontal cortex (PF; n = 4), and shanks of 1-mm length were used for premotor cortex (PM; n = 2). During recording sessions, monkeys sat head-fixed and alert in a monkey chair with no behavioral task given. Portions of this dataset have been analyzed previously (Bellay et al., 2021;Meisel et al., 2013;Miller et al., 2019;Yang et al., 2012;Yu et al., 2014).

Local Field Potential Recordings in Awake Monkeys
Simultaneous and continuous extracellular recordings were obtained for 12-60 min per recording session (2 kHz sampling frequency), band-pass filtered between 1 and 100 Hz (sixth-order Butterworth filter) to obtain the local field potential (LFP), and notch-filtered (60 Hz) to remove line noise. About 2 ± 1% of time periods were removed from functional electrodes because of artifacts introduced by vocalization, chewing, sudden movements, and the like. These artifacts were identified by threshold crossing (SD > 7) and excised (± 0.25 s). Arrays on average contained 86 ± 8 functional electrodes that exhibited 64 ± 50 μV of spontaneous LFP fluctuations (SD). Channels that had been removed from an array at any recording day were discounted for all recording days. Electrode LFPs were z-transformed and recording sessions for each array were analyzed individually. The current study represents a combined 26 hr of ongoing cortical LFP activity.

Neuronal Avalanche Definition and Temporal Resolution
For each electrode in the array, peak amplitude and time of negative LFP (nLFP) threshold crossings (−2 SD) were extracted at the temporal resolution of Δt = 0.5 ms given our sampling rate of 2 kHz. We note that the negative peak amplitude of the LFP correlates with the probability of extracellular unit firing and synchrony in nonhuman primates (Bellay et al., 2021;Petermann et al., 2009). The mean inter-event interval, <IEI>, defined as the average time between consecutive nLFPs among all functional electrodes on the array, was calculated. We then binned nLFP times in steps of Δt = <IEI> for each electrode, which ranged between 3 and 4 ms across NHPs and arrays. All nLFP events from all electrodes were combined into a matrix, that is, raster, with rows representing electrodes and columns representing time steps. A population time vector was obtained by summing nLFPs in the raster for each time step. Avalanches were defined as spatiotemporal continuous activity in the population vector bracketed on each side by at least one time bin of duration Δt with no nLFP. The size of an avalanche, S, was defined as the number of nLFPs participating. Multiple nLFPs at an electrode during an avalanche are rare (Yu et al., 2014) and were counted in size estimates. Scaleinvariance of S was visualized by plotting probability distributions P(S) in double-log coordinates. We previously showed that analyzing nLFP at Δt = <IEI> results in a distribution of avalanches with a slope close to −3/2 in line with expectations for a critical branching process (e.g., Beggs & Plenz, 2003).
The maximum log-likelihood ratio (LLR) was calculated to test potential power law distributions in avalanche size S against the alternatives of exponential or log-normal distribution models (Clauset et al., 2009;Klaus et al., 2011). When tested positive for power law, the LLR estimate for best slope α was reported. Here we used the range of S = 1 to 40, which excludes the distribution cutoff (Klaus et al., 2011) close to the total number of functional electrodes on the array (n > 70). We note that analyzing avalanche dynamics at different temporal resolutions, such as at shorter Δt or longer Δt compared with the <IEI>, increases or decreases the slope α respectively, but not the power law form itself (see also Table 1; Beggs & Plenz, 2003;Petermann et al., 2009).

Avalanche Temporal Profile Collapse and Scaling Exponent
Given the difficulty of obtaining robust lifetime distributions in the presence of ongoing oscillations, the stability analysis for avalanches in the present work focused mainly on the calculation of size distributions and temporal profile collapse (for details, see Miller et al., 2019). The temporal profile, S(t), was described by the number of participating electrodes at each time step t = 1 up to the avalanche lifetime, L, which was defined by the number of time steps Δt that the avalanche persisted. In line with our previous results, profile analysis was carried out at Δt = 1 ms and Δt = 30 ms to avoid profile modulation by ongoing -oscillations (Miller et al., 2019). Avalanches were grouped by L in multiples of Δt and averaged to obtain the mean temporal profile for a given lifetime, hSi(t, L). After normalizing to dimensionless time units, that is, the relative lifetime, t/L, amplitudes were then rescaled via Equation 1.
The profile collapse function, shown in Equation 1, relates the mean profile for each lifetime L, hSi(t/L), with a characteristic temporal motif F F F F (t/L), and scaling factor, L χ−1 , which is independent of L according to Equation 1. To perform a shape collapse, we plotted hSi(t, L) from L − 1 through L + 1 for L min = 4Δt (to reduce finite size effects in shape caused by too few data points). The collapse error, Δ F, was quantified via a normalized mean squared error of height-normalized individual profiles to the combined normalized average of all collapsed profiles, F F F F (t/L). Minimized collapse error was calculated by scanning through χ = 0.5 to 3 at resolution of 0.001 to find the collapse in avalanche waveform associated with the smallest Δ F via χ collapse . A value of Δ F > 1 was considered a failure in collapse.
We note that for fractal objects such as spatiotemporal avalanches, changing the temporal resolution also partitions avalanches differently. For a Δt = 1 ms, avalanche durations in multiples of Δt will exhibit a power law as well as for Δt = 30 ms (see, e.g., Miller et al., 2019). Accordingly, the temporal profile of avalanches must be examined in multiples of generations of duration Δt. A comparison between timescales in absolute times may not be possible. In the present study, an avalanche of duration 150 ms encompasses L = 5 generations at Δt = 30 but would require L = 150 generations at Δt = 1, which is highly unlikely to occur.

Fit of the Average Temporal Profile of Avalanches
For the parabolic fit we used the approach by Laurson et al. (2013) as follows: The parabolic fit error, Δ parab , was quantified via a normalized mean squared error of individual profiles to an amplitude-matched parabola that was coarse-grained to match L. Comparison to a semicircle fit was conducted in the same manner to obtain Δ semi using (3)

Correlation-Based Network Reconstruction
Networks were reconstructed for each recording session individually and the links with low pairwise correlations were removed before subsequent network analysis. Rather than applying a strict minimum correlation threshold for removing the links, which would leave our analysis vulnerable to inter-NHP differences in electrode impedances and so on, we instead sequentially removed the weakest links in 0.01 increments of the correlation threshold until a predefined sparsity of approximately 40% was achieved. Weight bins were divided into deciles for pruning analysis (Figure 4).

Normalized Count Reconstruction of Network
To remove the influence of common input to electrodes on the array, which inflates the Pearson's pairwise correlation, we conducted a network reconstruction. This nonparametric method, called the NC approach, is used for the network reconstruction because its computational efficiency is comparable to the correlation analysis yet with much lower reconstruction error, which are important considerations for analysis of long LFP time series collected in awake in vivo preparations. For this reconstruction, the time window within which active nodes can influence each other was examined for 0.5, 1, <IEI>, and 30 ms. For each node activation (spike), all nodes that were active within Δt, prior to the spike, were deemed to be potential causes. Results were very similar for Δt = 0.5 and 1 ms (data not shown). For longer Δt, the precision of the algorithm decreases rapidly because of (a) the increase in ambiguous potential associations between the number of active nodes in successive frames and (b) a decrease in the number of transition time steps to estimate link weights (Pajevic & Plenz, 2009). Here we show results for Δt = 1 ms, which matches our avalanche scaling analysis.

Integrative Networks and Pruning Analysis
Once the networks were defined, we studied their integrative properties (Pajevic & Plenz, 2012). The link clustering coefficient, C L , was defined for each link as the ratio of the number of common neighbors of its end nodes to the number of all neighbors (see Figure 4A). In an integrative network, there is a positive correlation R CL between C L and the weight of the link, indicating that strongly connected nodes are likely to share more common neighbors. To visualize this correlation, we divide links based on their weight into equally sized 10 blocks ranging from the weakest block (rank 1) to the highest rank. For each block, we obtain the average C L and subtract the equivalent value that is obtained from the degree-sequence preserved randomized network, yielding the excess link clustering, ΔC L . We also conducted a "pruning analysis" in which we studied the robustness of the excess node clustering ΔC upon the progressive removal of the weakest (bottom pruning) or strongest (top pruning) links. Similarly, the excess node clustering, ΔC, was defined as the difference of the mean clustering coefficient of the original network and the corresponding randomized network after the degree-sequence preserved randomization.

Entropy of Mixing Analysis
The estimated link weights fluctuated over time both in pairwise correlation and in directed networks obtained with the NC reconstruction algorithm. We sought to explore whether these fluctuations predominantly represent an error in the network reconstruction or arise from a genuine and progressive change in the underlying network weights. To answer this question, we developed a novel method that utilizes the entropy of mixing to quantify the progression of these fluctuations in individual link strengths from one day to another. The entropy of mixing is a concept from thermodynamics and describes the increase in total entropy when partitioned (pure and equilibrated) subsystems are allowed to mix. It is quantified using the Shannon entropy in which probabilities are replaced by the fractions of each of the original species found at different partitions. Hence, to conduct our entropy of mixing analysis, we need to partition the links into different categories. Labels are assigned according to the magnitude of their weights obtained from the Day 1 data, and partitioning them into a finite set of categories, N c , ranging from the weakest links to the strongest links. In our case, we used N c = 5, so that the label c = 1 marked the weakest 20% of the sorted links and incrementally, c = 5 marked the strongest 20% of the links.
We conducted the entropy of mixing analysis on two different and complementary subnetworks. The first is the "top" network, consisting of the 50% of the strongest links, and the second is the "bottom" network, consisting of the weakest 50% percent of links. This was done mainly to have some test of the performance of the entropy of mixing analysis, as we expect that the "bottom" half is more dominated by the reconstruction errors of the NC algorithm than the "top" half. Note that with this two-stage stratification of the links, the strongest 20% links in the "top" network are effectively representing the 10% of strongest links in the full network, while the weakest 20% of the "top" network links have still greater magnitude than the strongest quintile of the "bottom" networks, and, vice versa, the weakest 20% of the links in the "bottom" represent the 10% of weakest links in the full network. Once the chosen links were labeled, we then monitored the changes of their weights over subsequent days. Specifically, we sort them again based on their new weights, but now partition them over a potentially larger number of weight groups, N g = k N c , where k is an integer, in order to make all groups contain a single category label for the Day 1 stratification (we use k = 3 and N g = 15, meaning that the links in groups 1 through 3 will all initially carry the label c = 1, the groups g = 4-6 will be labeled c = 2, etc.). Because of (a) reconstruction errors and (b) potential actual network weight changes from one day to another, each of the groups will on the subsequent days contain a mixture of the original labels, although it is expected that the low group indices, g, will early on predominantly carry the labels of the weak links, and vice versa. For each group, we calculate a fraction, f c , of the links that carry label c, which by definition are normalized to 1, that is, P Nc c¼1 f c = 1. For each group g we then calculate the Shannon entropy of such mixture according to H g = − P Nc c¼1 f c log f c . Finally, we report the average and the standard error over all non-edge H g values, that is, excluding the extreme groups, g = 1 and g = 15, as those have only unidirectional mixing (n = N g − 2 = 13). Hence, the average entropy of mixing reported in our results is H M = 1 Ng−2 P Ng−1 g¼2 H g .
To evaluate whether we have seen a progression in the weight evolution, we test it against the null hypothesis that there was no significant linear trend across the subsequent days. To do so we use a Monte Carlo (MC) procedure, in which we generate one million normally distributed replicates, N(0,σ G ) for each point, where σ G are the standard deviations of the means obtained using the non-edge groups described in the previous paragraph. For each of the N MC = 1,000,000 replicates we perform a linear regression and count the number of times the MC slopes were greater than the slope we obtained using the mean values for H M . From this count, N pos , our estimate of the p values is p = (N pos + 1)/(N MC + 1) (see Table 4).