Cover ArticleResearch Articles, Systems/Circuits

Structural Modularity Tunes Mesoscale Criticality in Biological Neuronal Networks

Samora Okujeni and Ulrich Egert
Journal of Neuroscience 5 April 2023, 43 (14) 2515-2526; https://doi.org/10.1523/JNEUROSCI.1420-22.2023

Abstract

Numerous studies suggest that biological neuronal networks self-organize toward a critical state with stable recruitment dynamics. Individual neurons would then statistically activate exactly one further neuron during activity cascades termed neuronal avalanches. Yet, it is unclear if and how this can be reconciled with the explosive recruitment dynamics within neocortical minicolumns in vivo and within neuronal clusters in vitro, which indicates that neurons form supercritical local circuits. Theoretical studies propose that modular networks with a mix of regionally subcritical and supercritical dynamics would create apparently critical dynamics, resolving this inconsistency. Here, we provide experimental support by manipulating the structural self-organization process of networks of cultured rat cortical neurons (either sex). Consistent with the prediction, we show that increasing clustering in neuronal networks developing in vitro strongly correlates with avalanche size distributions transitioning from supercritical to subcritical activity dynamics. Avalanche size distributions approximated a power law in moderately clustered networks, indicating overall critical recruitment. We propose that activity-dependent self-organization can tune inherently supercritical networks toward mesoscale criticality by creating a modular structure in neuronal networks.

SIGNIFICANCE STATEMENT Critical recruitment dynamics in neuronal networks are considered optimal for information processing in the brain. However, it remains heavily debated how neuronal networks would self-organize criticality by detailed fine-tuning of connectivity, inhibition, and excitability. We provide experimental support for theoretical considerations that modularity tunes critical recruitment dynamics at the mesoscale level of interacting neuron clusters. This reconciles reports of supercritical recruitment dynamics in local neuron clusters with findings on criticality sampled at mesoscopic network scales. Intriguingly, altered mesoscale organization is a prominent aspect of various neuropathological diseases currently investigated in the framework of criticality. We therefore believe that our findings would also be of interest for clinical scientists searching to link the functional and anatomic signatures of such brain disorders.

Introduction

For optimal computational performance, the brain is considered to operate at a critical phase transition between fading of activity and explosive recruitment (Beggs and Plenz, 2003; Shew et al., 2009; Beggs and Timme, 2012; Hesse and Gross, 2014; Gautam et al., 2015; Li et al., 2019; Plenz et al., 2021). Critical propagation dynamics were reported for various neuronal systems in vivo and in vitro, indicating that these operate at or near criticality (Beggs and Plenz, 2003, 2004; Pasquale et al., 2008; Petermann et al., 2009; Hahn et al., 2010; Tetzlaff et al., 2010; Beggs and Timme, 2012; Plenz et al., 2021). Yet, despite decades of discussion, it remains heavily debated how neuronal networks would self-organize to achieve criticality. Moreover, reports on highly correlated activity in clusters of several tens of neurons in vivo (Maruoka et al., 2017; Hosoya, 2019) or in vitro (Cohen et al., 2008; Teller et al., 2014; Lonardoni et al., 2017) suggest that neurons instead self-organize into supercritical local circuits with explosive recruitment dynamics. An extended theoretical framework predicts system criticality to arise from a mixture of regionally subcritical and supercritical dynamics (Wang and Zhou, 2012; Moretti and Muñoz, 2013; Hilgetag and Hütt, 2014) in modular networks. With appropriate modularity, runaway dynamics are confined within modules, allowing critical mesoscale dynamics to emerge more robustly than in random networks (Gutjahr et al., 2021) and at minimal wiring cost (Liang et al., 2022). Here, we propose that closeness to criticality in neuronal networks in vitro mainly depends on the homeostatic regulation of the neuronal interaction range as a function of neuronal clustering (spatial aggregation of cell bodies), which shapes the degree of modularity.

Modular mesoscale architecture reflects a prevailing design principle in neuronal networks in vivo (Mountcastle, 1997; Buxhoeveden and Casanova, 2002; Rockland, 2010) and in vitro (Stetter et al., 2012; Okujeni et al., 2017; Okujeni and Egert, 2019b). Synchronized activity and similar functional tuning of neurons, for example, within cortical minicolumns (Maruoka et al., 2017; Hosoya, 2019) or modules in the entorhinal cortex (Naumann et al., 2018; Tukker et al., 2022), indicates a binary mode of module activation that effectively relocates the elementary computational unit and computational benefits of criticality from the single neuron to the module level. Considered as a branching process (Beggs and Plenz, 2003), a critical state would then be achieved if the modules on average activate exactly only one downstream module during activation cascades termed neuronal avalanches. Yet, how local neuronal interactions shape appropriate modular connectivity in developing networks remains elusive.

Because it is currently not feasible to modify structural modularity systematically in intact brain tissue, we pharmacologically manipulated the self-organization of neuronal network architecture developing in vitro (Okujeni et al., 2017; Okujeni and Egert, 2019b). Here, we analyzed avalanche size distributions (ASD) in spontaneous network activity probed with microelectrode arrays (MEAs) and found a gradual transition from supercritical to subcritical ASDs for increasingly clustered networks. Assuming supercritical recurrent dynamics within clusters (Cohen et al., 2008; Teller et al., 2014; Lonardoni et al., 2017), this suggests weakening connectivity between clusters. ASDs followed a power law only in moderately clustered networks, indicating critical branching dynamics on the mesoscopic scale. Simulations of networks with realistic assumptions for neuronal interaction range and clustering indicate a strong correlation between the degree of clustering, network modularity, and the apparent branching dynamics.

We thus provide experimental support for the theoretical prediction (Wang and Zhou, 2012; Moretti and Muñoz, 2013; Hilgetag and Hütt, 2014) that moderate modularity creates a configurational corridor facilitating critical dynamics on the mesoscale level. Criticality and its relation to brain modularity is currently being investigated at different anatomic scales in the context of neurologic disorders and clinical diagnostics (Tagliazucchi et al., 2012; Jiang et al., 2018; Wang et al., 2019; Zimmern, 2020; Hagemann et al., 2021; Heiney et al., 2021; Alamian et al., 2022). Our findings offer insights into how neuroanatomical reorganization may have an impact on pathologic brain dynamics and criticality.

Materials and Methods

Cell culture techniques

Primary cortical cell cultures were prepared on MEAs (Multi Channel Systems; electrode grid layout/pitch distance in µm, 16 × 16/200) and standard coverslips (12 mm diameter, Roth). All substrates were coated with polyethyleneimine (150 µl 0.2% aqueous solution; Sigma-Aldrich; CAS number 9002-98-6) for cell adhesion. Cortical tissue was excised from brains of neonatal Wistar rat pups of either sex, minced with a scalpel, and transferred into PBS (Invitrogen). Tissue pieces were incubated with trypsin (isozyme mixture, 0.05%, 37°C, 15 min; Invitrogen). Proteolysis was stopped with horse serum (20%; Invitrogen). DNase (type IV, 50 μg/ml; Sigma-Aldrich) was added to eliminate cell trapping in DNA strands. Cells were dissociated by trituration with a serologic pipette, centrifuged (5 min, 617 × g) and resuspended in growth medium (Minimal Essential Medium supplemented with 5% heat-inactivated horse serum, 0.5–1 mm l-glutamine, 20 mm glucose, and 20 µg/ml gentamycin, 1 ml/pup; all from Invitrogen). Cells were counted with an automated cell counter (CASY, Schärfe System) and seeded at ∼3 * 105 cells per network (∼1 cm2). Networks developed in 1–2 ml growth medium in a humidified incubator (5% CO2, 37°C). Protein kinase C (PKC) inhibitor Gö6976 (1 μm; Sigma-Aldrich; CAS number 136194-77-9) and PKC agonist phorbol-12-myristate-13-acetate (1 μm; Sigma-Aldrich; CAS number 16561-29-8) were dissolved in dimethyl sulfoxide (DMSO; Sigma-Aldrich) and added to the culture medium directly after cell preparation. The maximal concentration of DMSO in the growth medium was 0.1%. Animal handling and tissue preparation were done in accordance with the guidelines for animal research at the University of Freiburg and approved by the Regierungspräsidium Freiburg (permits X-12/08D, X-16/07A, X-15/01H, X-18/04K).

Experimental design and statistical analysis

The statistical tools used are given below for morphologic, electrophysiological, and computational analyses separately.

Morphologic analyses

Cluster analysis

To evaluate the degree of neuronal clustering, networks with 400–1500 neurons/mm2 growing on coverslips or MEAs were fixed between 20–37 d in vitro (DIV) with paraformaldehyde (Sigma-Aldrich) and methanol (Sigma-Aldrich), respectively. Neuronal somata were either stained immunohistochemically (NeuN, rabbit-anti-NeuN, 1:500; Abcam; RRID:AB2744676; MAP2, chicken-anti-MAP2, 1:500; Abcam; RRID:AB_2138147) or labeled by transfection (pAAV-CaMKIIa-hChR2(H134R)-mCherry, titer, 1011 ml−1; Addgene; RRID:Addgene_26975). The detection of neuronal somata in immunofluorescent images (Axio Observer microscope, ZEN Digital Imaging for Light Microscopy, Carl Zeiss; RRID:SCR_013672) was based on the colocalization of one of the neuron-specific stains above with a general stain for cellular nuclei (DAPI, Sigma-Aldrich) using custom-designed image processing with MATLAB 2021a (MathWorks; background subtraction, spatial filtering, normalization, thresholding using Otsu's method). Clustering of neuronal somata was analyzed using a modified Clark-Evans clustering index (CI) accounting for the average cell body diameter as minimal possible interneuron distance (Clark and Evans, 1954; Okujeni et al., 2017). Cluster affiliation was determined by the proximity of a neuron to cluster centers identified as significant peaks in the neuron density landscape. These were obtained by convolution of the coordinates of each neuron with a 2D-Gaussian kernel (SD, 35 µm), z-scoring, and detection of local maxima above 0.1. Clusters were analyzed for size (average number of neurons within clusters) and cluster spread (square root of the average squared distances between neurons and their respective cluster center). Intercluster distance (ICD) was calculated as the average edge length (excluding perimeter edges) in a mesh obtained by Delaunay triangulation of cluster centers. The relative cluster spread σc of each network was calculated by dividing the average cluster spread by the average ICD.

Neurite analysis

Neurites were detected in immunofluorescence micrographs for axons [rabbit-anti-neurofilament (NF), 1:10; Abcam; RRID:AB_448148] using MATLAB (spatial high-pass filtering, peak detection in line scans at different orientation, skeletonization) as described previously (Okujeni et al., 2017). Axonal projection range was approximated statistically based on the average length of axon per neuron calculated as the sum of axon pixels divided by the number of neurons in the image and multiplied with the pixel resolution and a correction factor of 1.12. The latter accounted for the orientation dependence of neurite length per pixel (i.e., √2 for diagonal vs 1 for horizontal or vertical neurite orientation; Smit et al., 1994). The relative axon length λa was calculated as the ratio between the average length of axon per neuron and the average ICD in the network.

Electrophysiological measurements and analyses

MEA recordings

Multiunit spike activity was recorded with MEAs (MEA1060-BC and USB-MEA256 systems, 25 kHz sampling frequency, 12 bit analog/digital conversion, Multi Channel Systems; MC Rack software versions 3.3–4.0; RRID:SCR_014955) under culture conditions (37°C, 5% CO2). Recordings lasted at least 1 h per network. Action potentials were detected with a threshold set to −5 SD of the high-pass-filtered electrode signal (Butterworth second-order high-pass filter, 200 Hz cutoff, detection dead time 2 ms). Raw data from MEA recordings was processed with MATLAB for further analysis.

Avalanche analysis

Avalanche detection was modified from Pasquale et al. (2008). In our study, we considered that bin width should be a function of the propagation velocity and distance between recording sites (Beggs and Plenz, 2003). Accordingly, the shortest possible delay in the subsampled branching process is the distance-dependent propagation delay between neighboring electrodes. Instead of using fixed bin widths, we therefore calculated bin width as the average distance between neighboring electrodes with spike activity (≥0.2 mm) obtained by Delaunay triangulation divided by the action potential velocity of ∼1 m/s reported for cortical networks in vitro (Bakkum et al., 2013). The resulting average bin widths were slightly larger in clustered networks with larger spacing between electrodes detecting spike activity (PKC−, 0.26 ± 0.05 ms; PKCN, 0.28 ± 0.10 ms; PKC+, 0.32 ± 0.12 ms). Following the bin size calculation for each network, global spike trains were binned and avalanches were determined as continuous series of bins with activity that were delimited by empty bins. Avalanche size was defined as the number of unique recording sites that participated in the avalanche. ASDs captured the relative frequency p of a given s, omitting avalanche size s=1 to reduce the impact of uncorrelated noise on the initial slope of the distribution. To compare ASDs across recordings and networks, raw ASDs were logarithmized in p and s axis and normalized such that the first third of the curve pnorm (snorm) crossed the point [1 −1] involving interpolation and resampling at 256 regular intervals, as follows. First, we calculated the logarithmic size range Δs of the first third of the ASD: Δs=13log(smaxsmin), where smin is the smallest and smax the largest observed avalanche. Second, we calculated the logarithmic probability range Δp of the same ASD section as follows: Δp=log(p(smin + eΔs)p(smin)).

These ranges were used to normalize ASDs p(s) → pnorm(snorm) as follows: snorm=1Δs*log(ssmin) pnorm=1|Δp|*log(p(s)p(smin)).

The initial ASD slope α was calculated as: α=ΔpΔs.

The procedure highlighted the overall shape of a distribution as either straight, convex, or concave on a double logarithmic scale. Principal component analysis (PCA; the set of ASDs were treated as observations and their pnorm(snorm) values as dimensions) of the normalized ASDs showed that this fundamental shape aspect was captured by the first principal component (PC1) weighted by scoring factor (PC1 score). The median α of close to power law ASDs with a PC1 score between −2 and 2 was −1.89. We used this exponent to rescale ASDs to a common s/N-axis with s normalized by the number of sites with spike activity N, adapting similar analyses (Klaus et al., 2011; Levina and Priesemann, 2017), as follows: s→s/N p(s)→presc(s/N)=(NmaxN)∝*p(s), where Nmax = 256 is the number of MEA electrodes.

The exponent α = −1.89 was in the range of critical exponents previously reported for neuronal avalanches in cultured networks (Tetzlaff et al., 2010; Levina and Priesemann, 2017). Functional modularity Qfunc was calculated based on the correlation of avalanche participation between pairs of recording sites with spike activity. The resulting correlation matrices with negative correlations set to zero were analyzed by the Louvain method (Blondel et al., 2008). Network classification into low, moderate, and high Qfunc was based on the 33.3rd and 66.7th percentiles. To test how MEA size affected the avalanche analysis, we defined subsets of electrodes by successively removing perimeter electrodes from the dataset to subsample the 16 × 16 arrays down to 8 × 8 electrodes.

To visualize the gradual change from subcritical to supercritical as a function of Qfunc (see Fig. 5B) we rendered the array of ASDs as a surface by interpolation in s/N and Qfunc axes and boxcar smoothing with a kernel of 5% in s/N-range and 20% in Qfunc-range.

Computational model

To study dependencies between the structural network architecture and modularity, we implemented a spatial network model in which we could systematically modulate the degree of clustering and the projection range of neurons. We then explored how different network architectures influenced probabilistic activity propagation and resulting ASDs with respect to criticality.

Modulating the degree of clustering

We defined networks of 1000 neurons forming 37 clusters on a hexagonal grid with an ICD of 227.9 arbitrary units (a.u.; with a.u. roughly reflecting µm in vitro, this would correspond to 2 mm2 of a cultured network with 500 neurons/mm2). Neuron coordinates were determined from the coordinates of a randomly selected cluster center plus a deviation in x and y directions drawn from a Gaussian distribution with a varying SD. Positions violating a minimal distance (10 a.u.) to other neurons were rejected to account for the diameter of a neuronal soma. Neuronal clusters were subsequently redefined by assigning neurons to their respective closest cluster center. The effective relative cluster spread σc* was subsequently calculated for each simulated network as the SD of neuron coordinates from their respective cluster centers divided by ICD, therefore σc*<0.5 (superscript asterisks indicate model parameters corresponding to experimentally derived parameters). We modulated σc* between 0.16 and 0.32 to cover the range determined experimentally. The iterative neuron seeding procedure yielded clusters with 27 neurons on average and neuron distributions resembling those in cultured networks with a homogeneous distribution of clusters of varying compactness.

Modulating projection range

In neuronal networks in vitro, the probability to find a connection between two neurons decays with distance following a Gaussian profile (Barral and Reyes, 2016). Hence, we modulated projection range based on the SD of a 2D-Gaussian profile defining the distance-dependent probability to form synaptic connections. The relative axonal length λa determined from experiments was mapped to 2 SDs accounting for 95% of connections. The simulated axonal projection range relative to the ICD λa* was then varied between 0.2 and 1.5, which covered the corresponding experimental λa range. Each neuron was connected to 100 presynaptic neurons drawn according to their relative distance-dependent probability. The procedure fixed the in-degree but entailed variable out-degrees to account for the homeostatic regulation of connectivity in real networks. However, the average out-degree equaled the fixed in-degree. Individual presynaptic neurons could be drawn repeatedly, creating multiple synapses with a postsynaptic neuron (multapses).

Simulating avalanches

Neuronal avalanches were simulated based on probabilistic activity propagation in the network with each synapse conveying activity with a certain probability. First, we balanced synapses to establish an average branching ratio β of 1, that is, a neuron would activate one postsynaptic neuron on average in the next time step as follows: wij=β■1Nsyn, where synaptic strength wij reflects the probability of a synapse of neuron i to activate the postsynaptic neuron j, and Nsyn=100 is the number of input synapses per neuron. Multapses were treated as multiple individual synapses, each with the same probability to activate the target neuron. Neuronal avalanches (N = 106) were initiated by activating a randomly selected neuron in the network and ended when no further neuron was activated. ASDs were calculated based on the number of recruited neurons during avalanches. In a second step, we then increased β to 1.05 in the same networks to establish supercritical recruitment dynamics. In all simulations, avalanches were terminated once all neurons had been recruited (maximal avalanche size) to save computation time, which was particularly necessary in case of persistent avalanches occurring with branching ratio >1. In real networks, synaptic depression would also enforce avalanche termination following explosive recruitment dynamics and recurrent activation of neurons and their synapses.

Structural modularity Q* was calculated based on the connectivity matrices using the Louvain method (Blondel et al., 2008) and accounted for multapses effectively forming stronger connections. Functional modularity Q*func was calculated based on the correlation matrices for the participation of neurons (or subsampled neurons) in avalanches (analogous to Qfunc from the electrophysiological recordings) with negative correlations set to zero. For a given relative cluster spread σc and relative axon length λa, structural modularity (Qest) of in vitro networks was estimated by mapping from the dependence of Q* on the corresponding σc* and λa* in simulated networks, assuming the same dependency in vitro.

Results

In the current study, we investigated how the mesoscale network structure influences avalanche dynamics in networks of rat cortical neurons. Toward this, we manipulated the structural self-organization of neuronal networks in vitro by pharmacological modulation of PKC activity (PKC inhibition, PKC−; normal condition, PKCN; PKC activation, PKC+) as described previously (Okujeni et al., 2017; Okujeni and Egert, 2019b). Experimental results and simulations show that criticality assessed at the mesoscopic scale is a function of network modularity.

Neuronal clustering and axonal projection range in vitro

In spatial networks, modularity would depend on the relationship between neuronal clustering and the range of axonal projections. To assess this, we first analyzed the spatial distribution of neurons (Fig. 1A) in networks on coverslips (PKC−, N = 5; PKCN, N = 6; PKC+, N = 7 networks) that developed under different PKC activity conditions with 400–900 neurons/mm2 (PKC−, 824 ± 43; PKCN, 491 ± 45; PKC+, 519 ± 60 neurons/mm2, mean ± SD). Neuronal cell bodies clustered in all networks with an average cluster size of 15–32 neurons (PKC−, 30 ± 2; PKCN, 19 ± 3; PKC+, 27 ± 3 neurons/cluster, mean ± SD across networks) and an average ICD between 220 and 290 µm (PKC−, 226 ± 4 µm; PKCN, 235 ± 9 µm; PKC+, 273 ± 10 µm, mean ± SD across networks). The overall contrast in the neuron density landscape was determined as the relative spatial cluster spread σc (average SD of neuron distances to their respective cluster center divided by average ICD in a network), which decreased with increasing clustering, that is, with decreasing CI (Fig. 1B). Clusters were significantly broader in PKC− networks (σc = 0.29 ± 0.01; mean ± SD of σc across networks; p = 2.6*10−3; Student's unpaired two-sample t test) compared with PKCN networks (σc = 0.26 ± 0.02) and compact clusters in PKC+ networks (σc = 0.20 ± 0.01, p = 9.6 * 10−6). Second, we assessed how the projection range scaled with the degree of clustering. The true distribution of projection ranges cannot be reconstructed histologically in high-density networks. We therefore used the length of axon per neuron as a statistical proxy for axonal projection range, calculated by dividing the total length of axon segments by the number of neurons detected in an image. In contrast to intercluster connections, axons entangled within clusters cannot be readily resolved microscopically. Hence, the detection procedure underestimated the absolute length of axons within clusters and predominantly reflected intercluster connectivity. We therefore calculated the ratio λa of the length of axon per neuron and ICD as an indicator of connectedness across clusters. λa was increased in PKC− networks (1.36 ± 0.07, mean ± SD of λa across networks; p = 9.3 * 10−4, Student's unpaired two-sample t test) compared with PKCN networks (0.94 ± 0.18) and was decreased in PKC+ networks (0.45 ± 0.09, p = 5.5 * 10−5; Fig. 1C) and thus positively correlated with CI. Across PKC conditions, λa and σc were strongly correlated (Fig. 1D).

Figure 1.
Figure 1.

Clustering and structural modularity in developing networks. A, Chronic pharmacological modulation of PKC activity altered the mesoscale architecture of self-organizing networks in vitro. Top, Fluorescence micrographs showing somata (red, NeuN), axons (green, NF) and cellular nuclei (blue, DAPI) in networks with modified PKC activity at 20 DIV (somata appear white because of colocalization of NeuN with NF and nuclear stain). Bottom, Cluster identification (marker color) and assignment of neurons to clusters based on peaks (white crosses) in the landscape of relative neuron density (background). Clusters become more compact and delimited from left to right. B, Cluster spread relative to ICD (σc) decreased with increasing clustering (decreasing CI). C, Average length of axon per neuron relative to ICD (λa) decreased with increasing clustering. D, Relative cluster spread σc and relative axon length λa were positively correlated. E, Networks resembling those observed in vitro were modeled by increasing relative cluster spread (σc*) and axonal projection range (λa*) relative to the ICD, respectively (bottom left, labels denote σc*; λa*), within experimentally approximated ranges. Perimeter color refers to networks representing PKC+, PKCN, and PKC− conditions. Colored lines are projection patterns of individual neurons. F, As in vitro, CI* in simulated networks increased toward 1 (random distribution) with increasing σc* (arrow colors refer to networks in E; error bars indicate SD across networks). G, Structural modularity Q* mainly depended on λa* and was weakly influenced by σc* only in combination with a short projection range. Gray values were interpolated from combinations of λa* and σc* in simulated networks (white markers; colored arrows refer to networks in E). λa and σc of in vitro networks were within the parameter range simulated (colored markers, data from D). H, Structural modularity of in vitro networks estimated from the map in G (Qest) was strongly correlated with CI. In B-D, F and H, magenta line and label show the linear regression with respective coefficient of determination (R2) .

Neuronal clustering and estimated structural modularity in vitro

Estimating modularity in networks with several tens of thousands of neurons is difficult because of the inaccessibility of the full neuronal connectivity in vitro. However, distance-dependent connectivity in conjunction with clustering can be expected to create modular networks. Using simulated networks to gain insight into dependencies between the degree of clustering, axonal connectivity range, and modularity, we varied σc* as well as λa* (Fig. 1E) across the CI range found in vitro (Fig. 1F). Structural modularity Q* was mostly controlled by λa* and was affected by σc* only if λa* ≪ ICD (Fig. 1G). Because of the pronounced positive correlation between λa and σc in in vitro networks, estimated structural modularity Qest increased with decreasing CI, suggesting that neuronal clustering promoted structural modularity in vitro (Fig. 1H).

Neuronal clustering correlates with altered avalanche statistics

To determine dependencies between avalanche dynamics and network architecture, we recorded spontaneous spike activity from networks grown on MEAs (PKC−, N = 55/24; PKCN, N = 69/36; PKC+, N = 29/15 recordings/networks). Although MEA recordings undersampled network activity on the neuron level, the electrode pitch of 200 µm matched the spatial scale of ICDs determined for a subset of recorded networks also analyzed morphologically (Fig. 2A; PKC−, 206 ± 74 µm, N = 13; PKCN, 221 ± 79 µm, N = 11; PKC+, 258 ± 85 µm, N = 2; mean ± average SD). Neuron densities in this subset of networks ranged between 400 and 1600 neurons/mm2 (PKC−, 1070 ± 329; PKCN, 798 ± 198; PKC+, 1007 ± 120; mean ± SD). All networks produced synchronous bursting, regardless of their degree of clustering (Fig. 2B), but stronger clustering correlated with higher burst rates and activity levels and weaker synchronization (Okujeni et al., 2017; Okujeni and Egert, 2019a, b). Neuronal avalanches were detected as periods of continuous network activity delimited by silent periods. Each network was typically recorded twice for 1 h between 18 and 57 DIV, yielding between 104 and 106 avalanches per recording session. Critical branching dynamics produce ASDs following a power law. Raw ASDs for the number of recruited electrodes differed considerably across networks with some following a power law and others showing under-representation (concave curve in log-log plots) or overrepresentation (convex curve) of large avalanches, suggesting subcritical and supercritical dynamics, respectively (Fig. 2C). To capture their fundamental shape, we normalized ASDs with respect to the initial third of each curve (Fig. 2D) and performed a principal component analysis. The first principal component (PC1; Fig. 2D, inset) captured overall ASD curvature, explaining 86.4% of the variance (Fig. 2E, inset), which becomes apparent in the reconstruction of normalized ASDs using PC1 (Fig. 2E). PC1 continuously mapped the shape of normalized ASDs from concave to convex as function of the PC1 score with straight ASDs (power law) at zero (Fig. 2F). Distributions of the PC1 score differed significantly across PKC conditions [PKC−, 3.7 ± 0.3 | 2.5, p = 6.8 * 10−11 (vs PKCN), p = 4.3 * 10−22 (vs PKC+); PKCN, −0.7 ± 0.5 | 3.9; PKC+, −5.4 ± 0.7 | 3.7, p = 1.8*10−7 (vs PKCN); mean ± standard error of the mean (SEM) | SD, Student's unpaired two-sample t test]. In the morphologically analyzed networks, CI increased with PC1 score (Fig. 2G), reflecting the transition from concave (subcritical) to convex (supercritical) ASDs with increasing degree of neuronal clustering. The rescaling exponent α = −1.89 was determined as the median initial slope of ASDs close to a power law (PC1 scores between −2 and 2; Fig. 2H).

Figure 2.
Figure 2.

Neuronal clustering and avalanche statistics. A, Spontaneous spike activity was recorded at 256 sites in neuronal networks grown on MEAs (white circles, 0.2 mm pitch). Staining and cluster assignment was performed at 37 DIV as in Figure 1A. B, In all PKC conditions, spike activity was organized in periods of synchronous bursting and little interburst activity. Burst rates increased, and peak firing rates decreased with the degree of somata clustering. C, Raw ASDs with probability p plotted against avalanche size s (number of recruited recording sites; s = 1 was omitted to diminish the influence of noise) for PKC−, PKCN, and PKC+ networks. Thick lines show averaged ASDs of PKC conditions. D, E, ASDs normalized to the first third of the logarithmic curve (crosshairs; normalization yields a slope αnorm = −1) were analyzed by PCA. The first principal component (PC1, inset) weighted by PC1 score (E) captured the continuous spectrum of convex to concave ASD shapes by accounting for 86.4% of the variance (inset). F, Surface derived by sliding average across the array of ASDs illustrates the gradual transition from concave to convex ASDs as a function of PC1 score with ASDs close to a power law at PC1 score = 0 (white line). Surface color indicates the pointwise difference between the normalized avalanche size probability and a normalized power law (αnorm = −1) shown as a tilted grid. G, PC1 score increased with increasing clustering (decreasing CI, asterisks indicate data derived from networks shown in Fig. 1A; magenta lines show linear regression with 95% confidence intervals). Power law ASDs predominantly appear in networks with moderate clustering of somata (CI = 0.82 at PC1 score = 0). H, Slopes of the initial third of ASD curves (no normalization) α as function of PC1 score. Magenta line shows the moving median (window ±2) of α. ASDs close to a power law around PC1 score = 0 had a median α = −1.89 (crossing of blue lines) and were dominated by PKCN networks.

Neuronal clustering and functional modularity in vitro

Estimating modularity based on multisite recordings relies on a meaningful description of the functional network connectivity. We derived a measure of functional modularity from the pairwise correlation of sites in their participation in avalanches. The substructuring of the correlation matrices increased with the degree of clustering in networks (Fig. 3A), indicating the appearance of functional submodules. The strongest correlations in the most homogeneous PKC− networks showed little dependence on the distance between sites. Correlations were increasingly dominated by nearby and neighboring sites in the more clustered PKCN and PKC+ networks (Fig. 3B). Average correlation strength decreased with increasing clustering in the network. In addition, the distance dependence of the correlation became steeper with increasing clustering (Fig. 3C). This suggested progressively more local connectivity and spatial modularization of the network with increasing clustering. Functional modularity calculated from the correlation matrices (negative correlations set to zero) using the Louvain method (Qfunc) was significantly different between PKC conditions and their respective degrees of neuronal clustering [Fig. 3D; PKC−, 3.2 * 10−2 ± 1.9 * 10−3 | 1.4 * 10−2, p = 7.7*10−4 (vs PKCN), p = 2.4−11 (vs PKC+); PKCN, 5.3 * 10−2 ± 5.0 * 10−3 | 4.2*10−2; PKC+, 1.7 * 10−1 ± 2.4 * 10−2 | 1.3*10−1, p = 1.5*10−9 (vs PKCN); mean ± SEM | SD, Student's unpaired two-sample t test]. Combining electrophysiological and morphologic analyses for the subset of recorded networks mentioned above showed a strong correlation between the Qfunc and PC1 scores (Fig. 3E) as well as between Qfunc and CI (Fig. 3F), indicating a crucial role of neuronal clustering and associated network modularity in the gradual shift from supercritical to subcritical avalanche dynamics. Networks with moderate CI (0.82) displayed moderate Qfunc (6.4*10−2) and close to power law ASDs (PC1 score = 0). Variations of neuron density did not account for differences in Qfunc across networks (Fig. 3G).

Figure 3.
Figure 3.

Neuronal clustering promotes functional modularity. A, Pairwise correlation matrices for the participation of recording sites with spike activity determined for all avalanches became increasingly structured with the degree of clustering in PKC−, PKCN, and PKC+ networks (examples for networks shown in Fig. 1A), indicating the emergence of functional submodules. Negative correlations were zeroed. Only for illustration, correlations were saturated at the lowest 1% and highest 99% of their range. B, Functional modularization illustrated for the highest 400 correlations (edges) between MEA recording sites with spike activity (black dots; gray dots, no spike activity). C, Pairwise correlation weakened with distance between recording sites and decreased overall with the degree of clustering in PKC−, PKCN, and PKC+ networks. Isoclines indicate the cumulative fraction of equal or lower correlation coefficients at a given distance between recording sites. D, Functional modularity Qfunc was calculated based on the pairwise correlation matrices using the Louvain method and significantly differed between PKC conditions (box plot shows median and 25th and 75th percentiles with whiskers extending to the most extreme datapoints not considered as outliers). Asterisks denote significant differences between distributions (p < 0.001, Student's unpaired two-sample t test). E, Qfunc decreased exponentially with PC1 score increasing from negative (concave ASD shapes) to positive (convex ASD shapes). ASDs close to power law (PC1 score = 0) were associated with moderate Qfunc (intersection of blue lines) prominent in PKCN networks. F, Qfunc increased exponentially with decreasing CI (increasing clustering) calculated from micrographs. The correlation suggested that moderate Qfunc appeared with a moderate CI ∼0.8. G, Qfunc was only weakly correlated with neuron density. Magenta lines (E–G) show linear regression with logarithmic Qfunc-axis and 95% confidence intervals. Asterisks refer to networks shown in Figure 1A.

Moderate functional modularity correlates with critical avalanche dynamics

To illustrate how Qfunc was reflected in avalanche dynamics, we assigned networks of the different PKC conditions and respective ASDs (Fig. 4A) into the three categories low (<33.3rd percentile), moderate and high (>66.7th percentile) Qfunc (Fig. 4B). Note that networks recorded several times in the course of development can contribute to more than one category. Homogeneous PKC− networks dominated in the low Qfunc group, with convex (supercritical) ASDs, whereas strongly clustered PKC+ networks formed most of the high Qfunc group with concave (subcritical) ASDs. Moderately clustered PKCN networks, in turn, contributed most of the moderate Qfunc group with ASDs close to a power law (approximately critical). ASDs showed little dependence on spatial extent of the sampling of activity as estimated with computationally reduced MEA sizes (Fig. 4C).

Figure 4.
Figure 4.

Moderate functional modularity promotes close to power law ASDs. A, ASDs from the respective PKC conditions with avalanche size s normalized by the number of recordings sites N using α = −1.89 (obtained for close to power law ASDs, Fig. 2H) for rescaling to a common axis. Individual networks were typically recorded twice in the course of development after 18 DIV. B, Networks were assigned into the three categories representing low (<33.3rd percentile), moderate, and high (>66.7th percentile) Qfunc. PKC− networks were mostly part of the low Qfunc group with convex ASDs, whereas PKC+ networks dominated in the high Qfunc group with concave ASDs. PKCN networks contributed most to the moderate Qfunc group with ASDs close to a power law but also constituted large fractions in the other groups. ASD line color refers to PKC conditions in A. Blue lines show the respective Qfunc group averages. C, Reducing the virtual MEA size by successively removing perimeter electrodes did not fundamentally change ASDs.

Mesoscale criticality in simulated networks

Avalanche statistics derived from mesoscale subsampling of spike activity in in vitro networks showed a gradual transition from subcritical to supercritical ASDs as a function of Qfunc (Fig. 5A,B). To understand how network modularity affected recruitment dynamics across different scales, we determined avalanche statistics in model networks. First, we simulated probabilistic activity propagation with an average branching ratio β = 1, corresponding to a critical branching process, by balancing synaptic weights relative to the number of synapses Nsyn. Any neuron would thus on average activate one postsynaptic neuron. Increasing modularity in spatial networks shifted full network ASDs from a power law distribution with α = −1.5 toward concave subcritical distributions with a shift toward smaller avalanche sizes (Fig. 5C, left; Note that the over-representation of network-size avalanches is a finite-size effect collapsing all larger avalanches that could appear in networks with infinite size.). This can be explained by an increasing likelihood of redundant neuron activation through multapses and convergent connections with increasing modularity (coalescence). As supercritical ASDs dominated in homogenous PKC− networks, we increased β to 1.05 in the simulation. As expected, random networks with such dynamics displayed supercritical ASDs. Increasing modularity tuned full network ASDs toward a power law with a critical exponent α = −1.5 and then toward subcritical ASDs (Fig. 5C, right). Recruitment statistics within clusters showed the opposite dependence during avalanche dynamics, with supercritical intracluster recruitment in networks with high modularity transitioning to subcritical recruitment dynamics in homogeneous and random networks. In simulated networks, Q*func increased exponentially with Q* toward convergence at 1 (fully isolated neurons; Fig. 5D). β=1.05 effectively increased functional connectivity and thereby decreased Q*func, leading to a steeper exponential dependence between Q*func and Q*. MEA recordings sample avalanche dynamics on the mesoscale and typically show multiunit activity. To understand how this influences ASDs, we calculated avalanches for the simulated networks from the combined time series of the two neurons closest to each cluster center. With β = 1, mesoscale subsampling in random connectivity produced steeper ASDs with α ∼ −1.9, close to the exponent determined from our recordings (Fig. 5E, left). Increasing modularity led to subcritical recruitment dynamics. However, with β = 1.05, increasing modularity led to a gradual transition from supercritical to subcritical dynamics (Fig. 5E, right). Moderate structural modularity (Q* ∼0.8) was associated with close to power law ASDs with α approximately −1.9. For mesoscale subsampling, Q*func showed a similar exponential decay with decreasing Q* as assessed from the full network, with a steeper decay for β = 1.05 (Fig. 5F).

Figure 5.
Figure 5.

Clustering and modularity determine mesoscale avalanche statistics. A, B, ASDs sorted as a function of Qfunc. ASDs gradually transition from concave (subcritical) to convex (supercritical) shapes with decreasing Qfunc. Close to power law ASDs occur at moderate Qfunc. Line or surface color in A and B indicate the pointwise ratio of avalanche size probability and a power law with α = −1.89 shown as a tilted plane or grid. B, Surface obtained by moving median smoothing of the array of ASDs in A to visualize the gradual transition. White lines delimit the moderate Qfunc group in Figure 4B. C, Full network ASDs derived from simulations of probabilistic activity propagations (N = 1000 neurons; networks in Fig. 1E–G). Left, ASDs become increasingly subcritical with increasing structural modularity Q*, although neurons on average activate one postsynaptic neuron (branching ratio β = 1). ASD line color denotes pointwise ratio with a power law, indicated by a plane with α = −1.5. The ASD with finite-size effects expected for random connectivity is shown at Q* = 0.6 for reference. Right, Increasing β to 1.05 resulted in supercritical ASDs in networks with low Q*. Increasing Q* dampened supercritical recruitment toward critical and subcritical ASDs (color code and reference lines as in left). D, Q*func derived from the recruitment statistics in simulated networks (as in Fig. 3) scaled exponentially with Q* and depended on β. Increasing β from 1.0 to 1.05 reduced the functional modularity Q*func for a given structural modularity Q*, with increasing reduction toward lower Q* (steeper slope; magenta lines, regression with logarithmic Q*func-axis). E, Avalanche dynamics subsampled at the mesoscale level (using only activity of the two neurons closest to each of the respective cluster centers; N = 37) mimics MEA recordings (same simulations as in C). s/N hence denotes the fraction of sampled sites (1 per cluster) recruited in an avalanche. Mesoscale ASDs show the same dependence on Q* as full network ASDs, however, with a steeper slope α = −1.9 and reduced finite size effect (shown in A and B). ASD line color denotes pointwise ratio with a power law (α = −1.9). ASD effects expected for random connectivity shown at Q* = 0.6 for reference. F, The reduction of Q*func with decreasing Q* was more pronounced with subsampling. For β = 1.05 and low Q*, Q*func decayed toward zero indicating full recruitment of the sampled sites during avalanches.

Discussion

Ample experimental evidence suggests that neuronal networks self-organize toward a critical point where activity on average neither tends to grow nor decay (Plenz et al., 2021). Although a tempting concept, given the associated computational benefits (Gautam et al., 2015), fine-tuning of complex biological networks toward a narrow range in a vast configurational space seems unlikely (Heiney et al., 2021). Theoretical studies proposed that network modularity could be an important ingredient for critical dynamics by providing more configurational leeway on the basis of regional mixtures of subcritical and supercritical dynamics (Rubinov et al., 2011; Wang and Zhou, 2012; Moretti and Muñoz, 2013; Hilgetag and Hütt, 2014). Whether and how neuronal network self-organization establishes adequate modularity levels for critical dynamics is not known.

Dependencies between clustering and modularity

The degree of modularity in a neuronal network may be quantified in terms of its structural connections or based on the functional correlations between neurons. Because of the inaccessibility of the full connectivity, given the large number of neurons and synapses and synaptic plasticity even at short timescales, it is virtually impossible to quantify the true structural modularity in biological networks. Nevertheless, in neuronal networks in vivo and in vitro, the vast majority of connections are formed locally (Levy and Reyes, 2012; Barral and Reyes, 2016), and simulations of abstract networks show that modularity is influenced by clustering if average axonal projection ranges are in the range of ICDs (Gilarranz, 2020). Homeostatic regulation of connectivity by adjusting neurite field growth to the local neuron density (Okujeni and Egert, 2019b) would amplify short-range connectivity at the expense of long-range intercluster connectivity. Consistent with this, we find a negative dependence of average axon length normalized to ICD on the degree of clustering in vitro. Because of this control loop, axons would become shorter in overall denser networks, which explains why functional modularity was largely independent of the average neuron density. Functional modularity, however, increased exponentially with the degree of neuron clustering and estimates of structural modularity. The same relations appeared in our simulations of modular networks. The exponential dependence may be explained by structural modularity reflecting direct connections, whereas functional modularity also comprises indirect links that contribute to activity correlations. The resulting increase of functional connectivity decreases modularity. The effective strength of indirect connections decays exponentially with the number of synaptic steps by multiplication of activation probabilities. Consequently, in clustered networks with shorter axons, the number of effective indirect pathways and their influence on functional modularity decreases.

Impact of modularity on avalanche dynamics

At criticality and branching ratio = 1 synaptic weights need to be balanced by the number of synapses formed by each neuron. This could be reflected in the continuous redistribution of synaptic resources (Minerbi et al., 2009; Hazan and Ziv, 2020) and synaptic scaling (Wilson et al., 2007; Turrigiano, 2008) in biological networks. However, in recurrent networks, multiple synapsing of individual presynaptic neurons and convergent synapsing of neurons with correlated firing can result in excess depolarization of postsynaptic neurons, which effectively absorbs synaptic input (coalescence) and reduces postsynaptic firing (Zierenberg et al., 2020), decreasing the branching ratio with increasing modularity. Previous work indicated multapses in highly clustered PKC+ networks and a large number of weak divergent synaptic connections in homogeneous PKC− networks (Okujeni et al., 2017). The latter would have fewer multapses and less convergence; and, hence, activity propagation would be affected less by coalescence. This is in agreement with the transition from supercritical to subcritical dynamics with increasing clustering and modularity reported here.

Criticality assessment with subsampled networks

In modular networks with supercritical local circuits but critical mesoscale branching dynamics, sampling ASDs from local populations would entail strong deviations from power laws (Levina and Priesemann, 2017). In fact, most experimental studies reporting criticality in neuronal systems used electrode spacing too large (>0.1 mm) to differentiate intracluster dynamics but corresponded well with mesoscopic activity propagation across neuron clusters (Pasquale et al., 2008; Petermann et al., 2009; Hahn et al., 2010; Tetzlaff et al., 2010; Plenz et al., 2021). The more negative critical exponents reported for cultured networks in vitro, around −2 rather than −1.5 associated with a critical branching process, could be readily explained by subsampling (Carvalho et al., 2020).

Mesoscale criticality in vitro and in vivo

Our experimental data indicates a strong correlation between the degree of clustering and modularity in neuronal networks in vitro and their closeness to criticality. We propose that neurons in homogeneous networks will exhibit supercritical recruitment dynamics. Reports indicating supercritical dynamics within local clusters in vitro (Cohen et al., 2008; Teller et al., 2014; Lonardoni et al., 2017) support this concept in that each cluster constitutes an intrinsically homogeneously connected network. Appropriate modularity with weaker intercluster connectivity would then dampen activity propagation across clusters toward criticality (Wang et al., 2017; Gutjahr et al., 2021). In such a scenario, criticality would emerge at the mesoscopic network level of interacting clusters.

Mesoscale criticality may apply to the neocortex with its minicolumnar architecture and to other modular brain regions where clustered network architecture (Rockland, 2010; Ji et al., 2015) combines with highly correlated intracluster activity (Maruoka et al., 2017; Hosoya, 2019). Conversely, neuronal migration deficits that have an impact on the minicolumnar organization of the neocortex have been found in several neurodevelopmental disorders (Di Rosa et al., 2009; McKavanagh et al., 2015; Casanova and Casanova, 2019). In line with reports of disrupted criticality in such conditions, typically sampled from mesoscale and macroscale structures, clinical diagnostics increasingly investigates indicators of criticality in brain activity (Arviv et al., 2016; Jiang et al., 2018; Zimmern, 2020; Hagemann et al., 2021; Heiney et al., 2021; Alamian et al., 2022).

In conclusion, network modularization could be a powerful mechanism contributing to self-organized mesoscale criticality in neuronal networks.

Footnotes

  • This work was supported by the German Research Foundation Grants EXC 1086 (BrainLinks-BrainTools) and INST 39/963-1 FUGG, State of Baden-Württemberg through bwHPC, Federal Ministry of Education and Research Grant FKZ 01GQ0420 (Bernstein Focus Neurotechnology Freiburg-Tuebingen), and a Carl Zeiss Foundation Grant (Computational Neuroscience of Brain Disease). We thank Ute Riede for technical assistance, Ruth Eneida Montaño Crespo for the PKC+ micrograph in Figure 1A, and Ad Aertsen for comments on the manuscript.

  • The authors declare no competing financial interests.

  • Correspondence should be addressed to Samora Okujeni at okujeni{at}bcf.uni-freiburg.de

References

  1. ↵
    2. Alamian G,
    3. Lajnef T,
    4. Pascarella A,
    5. Lina J-M,
    6. Knight L,
    7. Walters J,
    8. Singh KD,
    9. Jerbi K
    (2022) Altered brain criticality in schizophrenia: new insights from magnetoencephalography. Front Neural Circuits 16:630621.
    OpenUrl
  2. ↵
    2. Arviv O,
    3. Medvedovsky M,
    4. Sheintuch L,
    5. Goldstein A,
    6. Shriki O
    (2016) Deviations from critical dynamics in interictal epileptiform activity. J Neurosci 36:12276–12292. https://doi.org/10.1523/JNEUROSCI.0809-16.2016 pmid:27903734
    OpenUrlAbstract/FREE Full Text
  3. ↵
    2. Bakkum DJ,
    3. Frey U,
    4. Radivojevic M,
    5. Russell TL,
    6. Müller J,
    7. Fiscella M,
    8. Takahashi H,
    9. Hierlemann A
    (2013) Tracking axonal action potential propagation on a high-density microelectrode array across hundreds of sites. Nat Commun 4:2181. https://doi.org/10.1038/ncomms3181 pmid:23867868
    OpenUrlCrossRefPubMed
  4. ↵
    2. Barral J,
    3. Reyes AD
    (2016) Synaptic scaling rule preserves excitatory-inhibitory balance and salient neuronal network dynamics. Nat Neurosci 19:1690–1696. https://doi.org/10.1038/nn.4415 pmid:27749827
    OpenUrlCrossRefPubMed
  5. ↵
    2. Beggs JM,
    3. Plenz D
    (2003) Neuronal avalanches in neocortical circuits. J Neurosci 23:11167–11177. https://doi.org/10.1523/JNEUROSCI.23-35-11167.2003 pmid:14657176
    OpenUrlAbstract/FREE Full Text
  6. ↵
    2. Beggs JM,
    3. Plenz D
    (2004) Neuronal avalanches are diverse and precise activity patterns that are stable for many hours in cortical slice cultures. J Neurosci 24:5216–5229. https://doi.org/10.1523/JNEUROSCI.0540-04.2004 pmid:15175392
    OpenUrlAbstract/FREE Full Text
  7. ↵
    2. Beggs JM,
    3. Timme N
    (2012) Being critical of criticality in the brain. Front Physiol 3:163.
    OpenUrlCrossRefPubMed
  8. ↵
    2. Blondel VD,
    3. Guillaume J-L,
    4. Lambiotte R,
    5. Lefebvre E
    (2008) Fast unfolding of communities in large networks. J Stat Mech Theory Mech 2008:P10008. https://doi.org/10.1088/1742-5468/2008/10/P10008
    OpenUrlCrossRef
  9. ↵
    2. Buxhoeveden DP,
    3. Casanova MF
    (2002) The minicolumn and evolution of the brain. Brain Behav Evol 60:125–151. https://doi.org/10.1159/000065935 pmid:12417819
    OpenUrlCrossRefPubMed
  10. ↵
    2. Carvalho TTA,
    3. Fontenele AJ,
    4. Girardi-Schappo M,
    5. Feliciano T,
    6. Aguiar LAA,
    7. Silva TPL,
    8. de Vasconcelos NAP,
    9. Carelli PV,
    10. Copelli M
    (2020) Subsampled directed-percolation models explain scaling relations experimentally observed in the brain. Front Neural Circuits 14:576727. https://doi.org/10.3389/fncir.2020.576727 pmid:33519388
    OpenUrlPubMed
  11. ↵
    2. Casanova MF,
    3. Casanova EL
    (2019) The modular organization of the cerebral cortex: evolutionary significance and possible links to neurodevelopmental conditions. J Comp Neurol 527:1720–1730. https://doi.org/10.1002/cne.24554 pmid:30303529
    OpenUrlPubMed
  12. ↵
    2. Clark PJ,
    3. Evans FC
    (1954) Distance to nearest neighbor as a measure of spatial relationships in populations. Ecology 35:445–453. https://doi.org/10.2307/1931034
    OpenUrlCrossRef
  13. ↵
    2. Cohen E,
    3. Ivenshitz M,
    4. Amor-Baroukh V,
    5. Greenberger V,
    6. Segal M
    (2008) Determinants of spontaneous activity in networks of cultured hippocampus. Brain Res 1235:21–30. https://doi.org/10.1016/j.brainres.2008.06.022 pmid:18602907
    OpenUrlCrossRefPubMed
  14. ↵
    2. Di Rosa E,
    3. Crow TJ,
    4. Walker MA,
    5. Black G,
    6. Chance SA
    (2009) Reduced neuron density, enlarged minicolumn spacing and altered ageing effects in fusiform cortex in schizophrenia. Psychiatry Res 166:102–115. https://doi.org/10.1016/j.psychres.2008.04.007 pmid:19250686
    OpenUrlCrossRefPubMed
  15. ↵
    2. Gautam SH,
    3. Hoang TT,
    4. McClanahan K,
    5. Grady SK,
    6. Shew WL
    (2015) Maximizing sensory dynamic range by tuning the cortical state to criticality. PLoS Comput Biol 11:e1004576. https://doi.org/10.1371/journal.pcbi.1004576 pmid:26623645
    OpenUrlCrossRefPubMed
  16. ↵
    2. Gilarranz LJ
    (2020) Generic emergence of modularity in spatial networks. Sci Rep 10:8708. https://doi.org/10.1038/s41598-020-65669-8 pmid:32457480
    OpenUrlPubMed
  17. ↵
    2. Gutjahr N,
    3. Hövel P,
    4. Viol A
    (2021) Controlling extended criticality via modular connectivity. J Phys Complex 2:035023. https://doi.org/10.1088/2632-072X/ac202e
  18. ↵
    2. Hagemann A,
    3. Wilting J,
    4. Samimizad B,
    5. Mormann F,
    6. Priesemann V
    (2021) Assessing criticality in pre-seizure single-neuron activity of human epileptic cortex. PLoS Comput Biol 17:e1008773. https://doi.org/10.1371/journal.pcbi.1008773 pmid:33684101
    OpenUrlPubMed
  19. ↵
    2. Hahn G,
    3. Petermann T,
    4. Havenith MN,
    5. Yu S,
    6. Singer W,
    7. Plenz D,
    8. Nikolic D
    (2010) Neuronal avalanches in spontaneous activity in vivo. J Neurophysiol 104:3312–3322. https://doi.org/10.1152/jn.00953.2009 pmid:20631221
    OpenUrlCrossRefPubMed
  20. ↵
    2. Hazan L,
    3. Ziv NE
    (2020) Activity dependent and independent determinants of synaptic size diversity. J Neurosci 40:2828–2848. https://doi.org/10.1523/JNEUROSCI.2181-19.2020 pmid:32127494
    OpenUrlAbstract/FREE Full Text
  21. ↵
    2. Heiney K,
    3. Huse Ramstad O,
    4. Fiskum V,
    5. Christiansen N,
    6. Sandvig A,
    7. Nichele S,
    8. Sandvig I
    (2021) Criticality, connectivity, and neural disorder: a multifaceted approach to neural computation. Front Comput Neurosci 15:611183. https://doi.org/10.3389/fncom.2021.611183 pmid:33643017
    OpenUrlPubMed
  22. ↵
    2. Hesse J,
    3. Gross T
    (2014) Self-organized criticality as a fundamental property of neural systems. Front Syst Neurosci 8:166.
    OpenUrlCrossRefPubMed
  23. ↵
    2. Hilgetag CC,
    3. Hütt M-T
    (2014) Hierarchical modular brain connectivity is a stretch for criticality. Trends Cogn Sci 18:114–115. https://doi.org/10.1016/j.tics.2013.10.016
    OpenUrlCrossRefPubMed
  24. ↵
    2. Hosoya T
    (2019) The basic repeating modules of the cerebral cortical circuit. Proc Jpn Acad Ser B Phys Biol Sci 95:303–311. https://doi.org/10.2183/pjab.95.022
    OpenUrl
  25. ↵
    2. Ji W,
    3. Gămănuţ R,
    4. Bista P,
    5. D'Souza RD,
    6. Wang Q,
    7. Burkhalter A
    (2015) Modularity in the organization of mouse primary visual cortex. Neuron 87:632–643. https://doi.org/10.1016/j.neuron.2015.07.004
    OpenUrlCrossRefPubMed
  26. ↵
    2. Jiang L,
    3. Sui D,
    4. Qiao K,
    5. Dong H-M,
    6. Chen L,
    7. Han Y
    (2018) Impaired functional criticality of human brain during Alzheimer's disease progression. Sci Rep 8:1324. https://doi.org/10.1038/s41598-018-19674-7
    OpenUrl
  27. ↵
    2. Klaus A,
    3. Yu S,
    4. Plenz D
    (2011) Statistical analyses support power law distributions found in neuronal avalanches. PLoS One 6:e19779. https://doi.org/10.1371/journal.pone.0019779 pmid:21720544
    OpenUrlCrossRefPubMed
  28. ↵
    2. Levina A,
    3. Priesemann V
    (2017) Subsampling scaling. Nat Commun 8:15140. https://doi.org/10.1038/ncomms15140 pmid:28469176
    OpenUrlCrossRefPubMed
  29. ↵
    2. Levy RB,
    3. Reyes AD
    (2012) Spatial profile of excitatory and inhibitory synaptic connectivity in mouse primary auditory cortex. J Neurosci 32:5609–5619. https://doi.org/10.1523/JNEUROSCI.5158-11.2012 pmid:22514322
    OpenUrlAbstract/FREE Full Text
  30. ↵
    2. Li M,
    3. Han Y,
    4. Aburn MJ,
    5. Breakspear M,
    6. Poldrack RA,
    7. Shine JM,
    8. Lizier JT
    (2019) Transitions in information processing dynamics at the whole-brain network level are driven by alterations in neural gain. PLoS Comput Biol 15:e1006957. https://doi.org/10.1371/journal.pcbi.1006957 pmid:31613882
    OpenUrlCrossRefPubMed
  31. ↵
    2. Liang J,
    3. Wang S-J,
    4. Zhou C
    (2022) Less is more: wiring-economical modular networks support self-sustained firing-economical neural avalanches for efficient processing. Natl Sci Rev 9:nwab102.
    OpenUrl
  32. ↵
    2. Lonardoni D,
    3. Amin H,
    4. Di Marco S,
    5. Maccione A,
    6. Berdondini L,
    7. Nieus T
    (2017) Recurrently connected and localized neuronal communities initiate coordinated spontaneous activity in neuronal networks. PLoS Comput Biol 13:e1005672. https://doi.org/10.1371/journal.pcbi.1005672 pmid:28749937
    OpenUrlCrossRefPubMed
  33. ↵
    2. Maruoka H,
    3. Nakagawa N,
    4. Tsuruno S,
    5. Sakai S,
    6. Yoneda T,
    7. Hosoya T
    (2017) Lattice system of functionally distinct cell types in the neocortex. Science 358:610–615. https://doi.org/10.1126/science.aam6125
    OpenUrlAbstract/FREE Full Text
  34. ↵
    2. McKavanagh R,
    3. Buckley E,
    4. Chance SA
    (2015) Wider minicolumns in autism: a neural basis for altered processing? Brain 138:2034–2045. https://doi.org/10.1093/brain/awv110 pmid:25935724
    OpenUrlCrossRefPubMed
  35. ↵
    2. Minerbi A,
    3. Kahana R,
    4. Goldfeld L,
    5. Kaufman M,
    6. Marom S,
    7. Ziv NE
    (2009) Long-term relationships between synaptic tenacity, synaptic remodeling, and network activity. PLoS Biol 7:e1000136. https://doi.org/10.1371/journal.pbio.1000136 pmid:19554080
    OpenUrlCrossRefPubMed
  36. ↵
    2. Moretti P,
    3. Muñoz MA
    (2013) Griffiths phases and the stretching of criticality in brain networks. Nat Commun 4:2521. https://doi.org/10.1038/ncomms3521 pmid:24088740
    OpenUrlCrossRefPubMed
  37. ↵
    2. Mountcastle VB
    (1997) The columnar organization of the neocortex. Brain 120:701–722. https://doi.org/10.1093/brain/120.4.701
    OpenUrlCrossRefPubMed
  38. ↵
    2. Naumann RK,
    3. Preston-Ferrer P,
    4. Brecht M,
    5. Burgalossi A
    (2018) Structural modularity and grid activity in the medial entorhinal cortex. J Neurophysiol 119:2129–2144. https://doi.org/10.1152/jn.00574.2017 pmid:29513150
    OpenUrlPubMed
  39. ↵
    2. Okujeni S,
    3. Egert U
    (2019a) Inhomogeneities in network structure and excitability govern initiation and propagation of spontaneous burst activity. Front Neurosci 13:543.
    OpenUrl
  40. ↵
    2. Okujeni S,
    3. Egert U
    (2019b) Self-organization of modular network architecture by activity-dependent neuronal migration and outgrowth. Elife 8:e47996. https://doi.org/10.7554/eLife.47996
    OpenUrl
  41. ↵
    2. Okujeni S,
    3. Kandler S,
    4. Egert U
    (2017) Mesoscale architecture shapes initiation and richness of spontaneous network activity. J Neurosci 37:3972–3987. https://doi.org/10.1523/JNEUROSCI.2552-16.2017
    OpenUrlAbstract/FREE Full Text
  42. ↵
    2. Pasquale V,
    3. Massobrio P,
    4. Bologna LL,
    5. Chiappalone M,
    6. Martinoia S
    (2008) Self-organization and neuronal avalanches in networks of dissociated cortical neurons. Neuroscience 153:1354–1369. https://doi.org/10.1016/j.neuroscience.2008.03.050 pmid:18448256
    OpenUrlCrossRefPubMed
  43. ↵
    2. Petermann T,
    3. Thiagarajan TC,
    4. Lebedev MA,
    5. Nicolelis MAL,
    6. Chialvo DR,
    7. Plenz D
    (2009) Spontaneous cortical activity in awake monkeys composed of neuronal avalanches. Proc Natl Acad Sci U S A 106:15921–15926. https://doi.org/10.1073/pnas.0904089106 pmid:19717463
    OpenUrlAbstract/FREE Full Text
  44. ↵
    2. Plenz D,
    3. Ribeiro TL,
    4. Miller SR,
    5. Kells PA,
    6. Vakili A,
    7. Capek EL
    (2021) Self-organized criticality in the brain. Front Phys 9:1–23.
    OpenUrl
  45. ↵
    2. Rockland KS
    (2010) Five points on columns. Front Neuroanat 4:22.
    OpenUrlPubMed
  46. ↵
    2. Rubinov M,
    3. Sporns O,
    4. Thivierge J-P,
    5. Breakspear M
    (2011) Neurobiologically realistic determinants of self-organized criticality in networks of spiking neurons. PLoS Comput Biol 7:e1002038. https://doi.org/10.1371/journal.pcbi.1002038 pmid:21673863
    OpenUrlCrossRefPubMed
  47. ↵
    2. Shew WL,
    3. Yang H,
    4. Petermann T,
    5. Roy R,
    6. Plenz D
    (2009) Neuronal avalanches imply maximum dynamic range in cortical networks at criticality. J Neurosci 29:15595–15600. https://doi.org/10.1523/JNEUROSCI.3864-09.2009 pmid:20007483
    OpenUrlAbstract/FREE Full Text
  48. ↵
    2. Smit AL,
    3. Sprangers JFCM,
    4. Sablik PW,
    5. Groenwold J
    (1994) Automated measurement of root length with a three-dimensional high-resolution scanner and image analysis. Plant Soil 158:145–149. https://doi.org/10.1007/BF00007928
    OpenUrl
  49. ↵
    2. Stetter O,
    3. Battaglia D,
    4. Soriano J,
    5. Geisel T
    (2012) Model-free reconstruction of excitatory neuronal connectivity from calcium imaging signals. PLoS Comput Biol 8:e1002653. https://doi.org/10.1371/journal.pcbi.1002653 pmid:22927808
    OpenUrlCrossRefPubMed
  50. ↵
    2. Tagliazucchi E,
    3. Balenzuela P,
    4. Fraiman D,
    5. Chialvo DR
    (2012) Criticality in large-scale brain FMRI dynamics unveiled by a novel point process analysis. Front Physiol 3:15.
    OpenUrlCrossRefPubMed
  51. ↵
    2. Teller S,
    3. Granell C,
    4. De Domenico M,
    5. Soriano J,
    6. Gómez S,
    7. Arenas A
    (2014) Emergence of assortative mixing between clusters of cultured neurons. PLoS Comput Biol 10:e1003796. https://doi.org/10.1371/journal.pcbi.1003796 pmid:25188377
    OpenUrlCrossRefPubMed
  52. ↵
    2. Tetzlaff C,
    3. Okujeni S,
    4. Egert U,
    5. Wörgötter F,
    6. Butz M
    (2010) Self-organized criticality in developing neuronal networks. PLoS Comput Biol 6:e1001013. https://doi.org/10.1371/journal.pcbi.1001013 pmid:21152008
    OpenUrlCrossRefPubMed
  53. ↵
    2. Tukker JJ,
    3. Beed P,
    4. Brecht M,
    5. Kempter R,
    6. Moser EI,
    7. Schmitz D
    (2022) Microcircuits for spatial coding in the medial entorhinal cortex. Physiol Rev 102:653–688. https://doi.org/10.1152/physrev.00042.2020
    OpenUrl
  54. ↵
    2. Turrigiano GG
    (2008) The self-tuning neuron: synaptic scaling of excitatory synapses. Cell 135:422–435. https://doi.org/10.1016/j.cell.2008.10.008
    OpenUrlCrossRefPubMed
  55. ↵
    2. Wang R,
    3. Lin P,
    4. Liu M,
    5. Wu Y,
    6. Zhou T,
    7. Zhou C
    (2019) Hierarchical connectome modes and critical state jointly maximize human brain functional diversity. Phys Rev Lett 123:038301.
  56. ↵
    2. Wang SJ,
    3. Zhou C
    (2012) Hierarchical modular structure enhances the robustness of self-organized criticality in neural networks. New J Phys 14:023005.
  57. ↵
    2. Wang Z,
    3. Tian C,
    4. Dhamala M,
    5. Liu Z
    (2017) A small change in neuronal network topology can induce explosive synchronization transition and activity propagation in the entire network. Sci Rep 7:561. https://doi.org/10.1038/s41598-017-00697-5
    OpenUrl
  58. ↵
    2. Wilson NR,
    3. Ty MT,
    4. Ingber DE,
    5. Sur M,
    6. Liu G
    (2007) Synaptic reorganization in scaled networks of controlled size. J Neurosci 27:13581–13589. https://doi.org/10.1523/JNEUROSCI.3863-07.2007
    OpenUrlAbstract/FREE Full Text
  59. ↵
    2. Zierenberg J,
    3. Wilting J,
    4. Priesemann V,
    5. Levina A
    (2020) Description of spreading dynamics by microscopic network models and macroscopic branching processes can differ due to coalescence. Phys Rev E 101:022301. https://doi.org/10.1103/PhysRevE.101.022301 pmid:32168601
  60. ↵
    2. Zimmern V
    (2020) Why brain criticality is clinically relevant: a scoping review. Front Neural Circuits 14:54.
    OpenUrl
Back to top

In this issue

Email
Print
View Full Page PDF
Citation Tools
Respond to this article
Structural Modularity Tunes Mesoscale Criticality in Biological Neuronal Networks
Samora Okujeni, Ulrich Egert
Journal of Neuroscience 5 April 2023, 43 (14) 2515-2526; DOI: 10.1523/JNEUROSCI.1420-22.2023
Twitter logo Facebook logo Mendeley logo

Jump to section

Keywords

Responses to this article

Respond to this article

Jump to comment:

No eLetters have been published for this article.