## Abstract

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.

## Author Summary

The brain is spontaneously active even in the absence of specific sensations or movements. This ongoing activity arises from the interactions among hundreds of thousands of neurons, which has been captured by a variety of distinct functional networks predictive of healthy and pathological brain functions. The global dynamical states and overarching network principles that underlie such ongoing brain activity are not well understood. Here we identify avalanche dynamics and “friendship” networks as two major features of ongoing activity in the frontal cortex of nonhuman primates. We demonstrate their stability over weeks in the face of local network changes. Deviation from avalanche dynamics and “friendship” organization might provide robust biomarkers to identify normal and pathological brain states.

## 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 single-cell 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 orientation-selective 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 single-cell 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 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.

On average about 4 ± 2 hr of LFP activity during 9 ± 7 recording sessions over the course of 5 ± 4 weeks were analyzed per monkey (85 ± 8 electrodes/array; Figure 1C).

### 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).

**Table 1.**

. | <R> Stability
. | Avalanche size distribution (Δt = 1 ms)
. | Slope α stability (Δt = 1 ms)
. | Avalanche size distribution (Δt =
<IEI>)
. | Slope α stability (Δt =
<IEI>)
. | Avalanche size distribution (Δt = 30 ms)
. | Slope α stability (Δt = 30 ms)
. |
---|---|---|---|---|---|---|---|

Regression slope (95% CfL) (×
10^{−3})
. | Exponent α (± SD)
. | Regression slope (95% CfL) (×
10^{−3})
. | Exponent α (± SD)
. | Regression slope (95% CfL) (×
10^{−3})
. | Exponent α (± SD)
. | Regression slope (95% CfL) (×
10^{−3})
. | |

V-PF | −1 (−3, 2) | 1.71 ± 0.00 | 0 (−1, 0) | 1.26 ± 0.01 | −2 (−3, 3) | 0.65 ± 0.03 | −3 (−9, 2) |

V-PM | 0 (−62, 62) | 1.64 ± 0.02 | 6 (−40, 51) | 1.29 ± 0.04 | 16 (−28, 60) | 1.12 ± 0.05 | 20 (−2, 43) |

N-PF | 1 (0, 2) | 1.71 ± 0.02 | −1 (−1, 0) | 1.41 ± 0.01 | 0 (−1, 0) | 0.97 ± 0.02 | 0 (−1, 1) |

N-PM | 1 (−13, 15) | 1.73 ± 0.02 | 2 (−8, 12) | 1.31 ± 0.02 | 0 (−1, 1) | 0.75 ± 0.02 | 3 (−6, 12) |

K-PF1 | 1 (0, 1) | 1.65 ± 0.02 | −1 (−1, 0) | 1.29 ± 0.01 | 0 (−1, 0) | 0.90 ± 0.02 | 0 (0, 1) |

K-PF2 | 0 (−20, 20) | 1.77 ± 0.03 | 5 (−33, 43) | 1.36 ± 0.02 | 1 (−27, 29) | 0.74 ± 0.03 | 10 (−16, 35) |

Mean ± SD | 0 ± 1 | 1.70 ± 0.05 | 2 ± 3 | 1.32 ± 0.05 | 3 ± 7 | 0.86 ± 0.17 | 5 ± 9 |

. | <R> Stability
. | Avalanche size distribution (Δt = 1 ms)
. | Slope α stability (Δt = 1 ms)
. | Avalanche size distribution (Δt =
<IEI>)
. | Slope α stability (Δt =
<IEI>)
. | Avalanche size distribution (Δt = 30 ms)
. | Slope α stability (Δt = 30 ms)
. |
---|---|---|---|---|---|---|---|

Regression slope (95% CfL) (×
10^{−3})
. | Exponent α (± SD)
. | Regression slope (95% CfL) (×
10^{−3})
. | Exponent α (± SD)
. | Regression slope (95% CfL) (×
10^{−3})
. | Exponent α (± SD)
. | Regression slope (95% CfL) (×
10^{−3})
. | |

V-PF | −1 (−3, 2) | 1.71 ± 0.00 | 0 (−1, 0) | 1.26 ± 0.01 | −2 (−3, 3) | 0.65 ± 0.03 | −3 (−9, 2) |

V-PM | 0 (−62, 62) | 1.64 ± 0.02 | 6 (−40, 51) | 1.29 ± 0.04 | 16 (−28, 60) | 1.12 ± 0.05 | 20 (−2, 43) |

N-PF | 1 (0, 2) | 1.71 ± 0.02 | −1 (−1, 0) | 1.41 ± 0.01 | 0 (−1, 0) | 0.97 ± 0.02 | 0 (−1, 1) |

N-PM | 1 (−13, 15) | 1.73 ± 0.02 | 2 (−8, 12) | 1.31 ± 0.02 | 0 (−1, 1) | 0.75 ± 0.02 | 3 (−6, 12) |

K-PF1 | 1 (0, 1) | 1.65 ± 0.02 | −1 (−1, 0) | 1.29 ± 0.01 | 0 (−1, 0) | 0.90 ± 0.02 | 0 (0, 1) |

K-PF2 | 0 (−20, 20) | 1.77 ± 0.03 | 5 (−33, 43) | 1.36 ± 0.02 | 1 (−27, 29) | 0.74 ± 0.03 | 10 (−16, 35) |

Mean ± SD | 0 ± 1 | 1.70 ± 0.05 | 2 ± 3 | 1.32 ± 0.05 | 3 ± 7 | 0.86 ± 0.17 | 5 ± 9 |

We next explored the stability of neuronal avalanche dynamics by discretizing the
LFP from ongoing activity in the NHPs
(Petermann et al., 2009; 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 (Yu et al., 2011), 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 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 2.**

. | Temporal profile fit: Parabola vs. semicircle RMSE
(Δt = 1 ms)
. | Temporal profile fit: Parabola vs. semicircle RMSE
(Δt = 30 ms)
. | Random profile fit: Parabola vs. semicircle RMSE
(Δt = 1 ms)
. | |||
---|---|---|---|---|---|---|

Paired t test
(p value)
. | ANOVA (F)
. | Paired t test
(p value)
. | ANOVA (F)
. | Paired t test
(p value)
. | ANOVA (F)
. | |

V-PF | 4.36E-05 | 110.44 | 2.00E-04 | 62.26 | 5.96E-10 | 4,832.54 |

V-PM | 3.87E-04 | 121.17 | 3.40E-05 | 417.07 | 7.82E-07 | 2,776.65 |

N-PF | 6.12E-11 | 460.49 | 2.66E-13 | 1156.15 | 3.74E-14 | 1,607.96 |

N-PM | 6.81E-08 | 991.66 | 3.15E-06 | 272.38 | 2.41E-10 | 6,533.12 |

K-PF1 | 2.76E-22 | 472.20 | 2.41E-25 | 715.09 | 1.98E-50 | 18,530.05 |

K-PF2 | 2.60E-06 | 290.68 | 9.64E-06 | 186.09 | 1.91E-09 | 3,277.87 |

Fit Error Δ | Δ_{parab} <
Δ_{semi} | Δ_{parab} <
Δ_{semi} | Δ_{parab} <
Δ_{semi} | Δ_{parab} <
Δ_{semi} | Δ_{parab} >
Δ_{semi} | Δ_{parab} >
Δ_{semi} |

. | Temporal profile fit: Parabola vs. semicircle RMSE
(Δt = 1 ms)
. | Temporal profile fit: Parabola vs. semicircle RMSE
(Δt = 30 ms)
. | Random profile fit: Parabola vs. semicircle RMSE
(Δt = 1 ms)
. | |||
---|---|---|---|---|---|---|

Paired t test
(p value)
. | ANOVA (F)
. | Paired t test
(p value)
. | ANOVA (F)
. | Paired t test
(p value)
. | ANOVA (F)
. | |

V-PF | 4.36E-05 | 110.44 | 2.00E-04 | 62.26 | 5.96E-10 | 4,832.54 |

V-PM | 3.87E-04 | 121.17 | 3.40E-05 | 417.07 | 7.82E-07 | 2,776.65 |

N-PF | 6.12E-11 | 460.49 | 2.66E-13 | 1156.15 | 3.74E-14 | 1,607.96 |

N-PM | 6.81E-08 | 991.66 | 3.15E-06 | 272.38 | 2.41E-10 | 6,533.12 |

K-PF1 | 2.76E-22 | 472.20 | 2.41E-25 | 715.09 | 1.98E-50 | 18,530.05 |

K-PF2 | 2.60E-06 | 290.68 | 9.64E-06 | 186.09 | 1.91E-09 | 3,277.87 |

Fit Error Δ | Δ_{parab} <
Δ_{semi} | Δ_{parab} <
Δ_{semi} | Δ_{parab} <
Δ_{semi} | Δ_{parab} <
Δ_{semi} | Δ_{parab} >
Δ_{semi} | Δ_{parab} >
Δ_{semi} |

**Table 3.**

. | Scaling exponent
. χ (Δt = 1
ms) |
. χ Stability (Δt = 1
ms) | Scaling exponent
. χ (Δt = 30
ms) |
. χ Stability (Δt = 30
ms) | Scaling exponent
. χ Random
(Δt = 1 ms) |
. χ Stability Random (Δt = 1
ms) |
---|---|---|---|---|---|---|

Mean χ (± SD) | Regression slope (95% CfL) | Mean χ (± SD) | Regression slope (95% CfL) | Mean χ (± SD) | Regression slope (95% CfL) | |

V-PF | 1.78 ± 0.05 | 0.002 (−0.019, 0.022) | 1.61 ± 0.09 | −0.01 (−0.02, 0.01) | 1.15 ± 0.01 | 0.001 (−0.001, 0.003) |

V-PM | 1.93 ± 0.15 | 0.051 (−0.3, 0.4) | 2.25 ± 0.20 | 0.06 (−0.6, 0.7) | 1.15 ± 0.00 | 0.000 (−0.001, 0.002) |

N-PF | 2.10 ± 0.06 | 0.002 (−0.000, 0.004) | 2.15 ± 0.15 | −0.00 (−0.01, 0.00) | 1.14 ± 0.03 | 0.000 (−0.001, 0.002) |

N-PM | 2.05 ± 0.13 | 0.011 (−0.053, 0.075) | 1.85 ± 0.05 | 0.005 (−0.016, 0.027) | 1.14 ± 0.01 | 0.002 (−0.005, 0.008) |

K-PF1 | 1.92 ± 0.03 | 0.000 (−0.001, 0.001) | 2.11 ± 0.10 | −0.000 (−0.003, 0.003) | 1.15 ± 0.01 | 0.000 (−0.000, 0.000) |

K-PF2 | 1.76 ± 0.02 | 0.007 (−0.004, 0.018) | 2.12 ± 0.20 | 0.086 (−0.029, 0.200) | 1.15 ± 0.00 | 0.000 (−0.004, 0.005) |

Mean ± SD | 1.93 ± 0.14 | 0.012 ± 0.19 | 2.01 ± 0.24 | 0.02 ± 0.04 | 1.15 ± 0.01 | 0.001 ± 0.001 |

. | Scaling exponent
. χ (Δt = 1
ms) |
. χ Stability (Δt = 1
ms) | Scaling exponent
. χ (Δt = 30
ms) |
. χ Stability (Δt = 30
ms) | Scaling exponent
. χ Random
(Δt = 1 ms) |
. χ Stability Random (Δt = 1
ms) |
---|---|---|---|---|---|---|

Mean χ (± SD) | Regression slope (95% CfL) | Mean χ (± SD) | Regression slope (95% CfL) | Mean χ (± SD) | Regression slope (95% CfL) | |

V-PF | 1.78 ± 0.05 | 0.002 (−0.019, 0.022) | 1.61 ± 0.09 | −0.01 (−0.02, 0.01) | 1.15 ± 0.01 | 0.001 (−0.001, 0.003) |

V-PM | 1.93 ± 0.15 | 0.051 (−0.3, 0.4) | 2.25 ± 0.20 | 0.06 (−0.6, 0.7) | 1.15 ± 0.00 | 0.000 (−0.001, 0.002) |

N-PF | 2.10 ± 0.06 | 0.002 (−0.000, 0.004) | 2.15 ± 0.15 | −0.00 (−0.01, 0.00) | 1.14 ± 0.03 | 0.000 (−0.001, 0.002) |

N-PM | 2.05 ± 0.13 | 0.011 (−0.053, 0.075) | 1.85 ± 0.05 | 0.005 (−0.016, 0.027) | 1.14 ± 0.01 | 0.002 (−0.005, 0.008) |

K-PF1 | 1.92 ± 0.03 | 0.000 (−0.001, 0.001) | 2.11 ± 0.10 | −0.000 (−0.003, 0.003) | 1.15 ± 0.01 | 0.000 (−0.000, 0.000) |

K-PF2 | 1.76 ± 0.02 | 0.007 (−0.004, 0.018) | 2.12 ± 0.20 | 0.086 (−0.029, 0.200) | 1.15 ± 0.00 | 0.000 (−0.004, 0.005) |

Mean ± SD | 1.93 ± 0.14 | 0.012 ± 0.19 | 2.01 ± 0.24 | 0.02 ± 0.04 | 1.15 ± 0.01 | 0.001 ± 0.001 |

### Persistence of Integrative Global Network Properties

In order to quantify in more detail the spatiotemporal correlations underlying persistent avalanche activity, we reconstructed functional networks from neuronal avalanches using the NC 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 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).

**Table 4.**

. |
. R_{CL} stability
(Δt = 1 ms) |
. H_{M} trend
(Δt = 1 ms) | |
---|---|---|---|

Regression slope (95% CfL) (×
10^{−3})
. | Top 50% p
. | Bottom 50% p
. | |

V-PF | −1 (−28, 18) | 0.46 | 0.34 |

V-PM | 0 (−62, 62) | 0.007 | 0.05 |

N-PF | 1 (0, 0.002) | <10^{−5} | 0.03 |

N-PM | 14 (−13, 15) | 0.08 | 0.54 |

K-PF1 | 1 (0, 1) | 0.003 | 0.2 |

K-PF2 | 0 (−20, 21) | 0.27 | 0.35 |

Mean ± SD | 5 ± 10 |

. |
. R_{CL} stability
(Δt = 1 ms) |
. H_{M} trend
(Δt = 1 ms) | |
---|---|---|---|

Regression slope (95% CfL) (×
10^{−3})
. | Top 50% p
. | Bottom 50% p
. | |

V-PF | −1 (−28, 18) | 0.46 | 0.34 |

V-PM | 0 (−62, 62) | 0.007 | 0.05 |

N-PF | 1 (0, 0.002) | <10^{−5} | 0.03 |

N-PM | 14 (−13, 15) | 0.08 | 0.54 |

K-PF1 | 1 (0, 1) | 0.003 | 0.2 |

K-PF2 | 0 (−20, 21) | 0.27 | 0.35 |

Mean ± SD | 5 ± 10 |

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 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 graph-theoretical 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 (Shriki et al., 2013), and fMRI (Fraiman & Chialvo, 2012; 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., 2011; Yang 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 NC-reconstructed 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 (Fraiman & Chialvo, 2012; 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 long-range, selective propagation of information.

## MATERIAL AND METHODS

### 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. Scale-invariance 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

*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, 〈

*S*〉(

*t*,

*L*). After normalizing to dimensionless time units, that is, the relative lifetime,

*t*/

*L*, amplitudes were then rescaled via Equation 1.

*L*, 〈

*S*〉(

*t*/

*L*), with a characteristic temporal motif 𝓕(

*t*/

*L*), and scaling factor,

*L*

^{χ−1}, which is independent of

*L*according to Equation 1. To perform a shape collapse, we plotted 〈

*S*〉(

*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, 𝓕(

*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

_{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

### 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, $\u2211c=1Nc$*f*_{c} = 1. For each group *g* we then calculate the Shannon entropy of such mixture
according to *H*_{g} =
−$\u2211c=1Nc$*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} = $1Ng\u22122\u2211g=2Ng\u22121$*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).

## ACKNOWLEDGMENTS

We thank members of the Plenz lab and attendees of the Critical Brain Dynamics: 5th International Workshop on Criticality and the Brain, 2016, NIH, USA, for their support and useful discussions. This research utilized the computational resources of the NIH HPC Biowulf cluster (https://hpc.nih.gov).

## AUTHOR CONTRIBUTIONS

Stephanie R. Miller: Investigation; Methodology; Software; Writing – original draft. Sinisa Pajevic: Conceptualization; Formal analysis; Investigation; Methodology; Software; Writing – original draft. Shan Yu: Conceptualization; Formal analysis; Investigation; Methodology; Software; Writing – original draft. Dietmar Plenz: Conceptualization; Investigation; Methodology; Supervision; Writing – original draft.

## FUNDING INFORMATION

Dietmar Plenz, National Institute of Mental Health, Division of the Intramural Research Program, Award ID: ZIAMH002797. Sinisa Pajevic, National Institute of Child Health and Development, Division of the Intramural Research Program.

## TECHNICAL TERMS

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

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

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

## REFERENCES

## Author notes

Competing Interests: The authors have declared that no competing interests exist.

Handling Editor: Gustavo Deco