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

- mesoscale architecture
- modularity
- network connectivity
- neuronal avalanches
- neuronal clustering
- self-organized criticality

## 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 * 10^{5} cells per network (∼1 cm^{2}). Networks developed in 1–2 ml growth medium in a humidified incubator (5% CO_{2}, 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/mm^{2} 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

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

#### 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% CO_{2}). 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; PKC^{N}, 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 *p* and *s* axis and normalized such that the first third of the curve *p*_{norm} (*s*_{norm}) crossed the point [1 −1] involving interpolation and resampling at 256 regular intervals, as follows. First, we calculated the logarithmic size range

These ranges were used to normalize ASDs

The initial ASD slope α was calculated as:

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 _{1}) weighted by scoring factor (PC_{1} score). The median α of close to power law ASDs with a PC_{1} 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:

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 *Q _{func}* 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

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

_{func}To visualize the gradual change from subcritical to supercritical as a function of *Q _{func}* (see Fig. 5

*B*) we rendered the array of ASDs as a surface by interpolation in s/N and

*Q*axes and boxcar smoothing with a kernel of 5% in

_{func}*s/N*-range and 20% in

*Q*-range.

_{func}#### 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 mm^{2} of a cultured network with 500 neurons/mm^{2}). 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

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

##### 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:
*N* = 10^{6}) 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

*Q*from the electrophysiological recordings) with negative correlations set to zero. For a given relative cluster spread

_{func}*Q*) of

_{est}*in vitro*networks was estimated by mapping from the dependence of

*Q** on the corresponding

*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, PKC^{N}; 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. 1*A*) in networks on coverslips (PKC^{−}, *N* = 5; PKC^{N}, *N* = 6; PKC^{+}, *N* = 7 networks) that developed under different PKC activity conditions with 400–900 neurons/mm^{2} (PKC^{−}, 824 ± 43; PKC^{N}, 491 ± 45; PKC^{+}, 519 ± 60 neurons/mm^{2}, mean ± SD). Neuronal cell bodies clustered in all networks with an average cluster size of 15–32 neurons (PKC^{−}, 30 ± 2; PKC^{N}, 19 ± 3; PKC^{+}, 27 ± 3 neurons/cluster, mean ± SD across networks) and an average ICD between 220 and 290 µm (PKC^{−}, 226 ± 4 µm; PKC^{N}, 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 *B*). Clusters were significantly broader in PKC^{−} networks (*p* = 2.6*10^{−3}; Student's unpaired two-sample *t* test) compared with PKC^{N} networks (^{+} networks (*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 ^{−} networks (1.36 ± 0.07, mean ± SD of *p* = 9.3 * 10^{−4}, Student's unpaired two-sample *t* test) compared with PKC^{N} networks (0.94 ± 0.18) and was decreased in PKC^{+} networks (0.45 ± 0.09, *p* = 5.5 * 10^{−5}; Fig. 1*C*) and thus positively correlated with CI. Across PKC conditions, *D*).

### 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 *E*) across the CI range found *in vitro* (Fig. 1*F*). Structural modularity *Q** was mostly controlled by *G*). Because of the pronounced positive correlation between *in vitro* networks, estimated structural modularity *Q _{est}* increased with decreasing CI, suggesting that neuronal clustering promoted structural modularity

*in vitro*(Fig. 1

*H*).

### 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; PKC^{N}, *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. 2*A*; PKC^{−}, 206 ± 74 µm, *N* = 13; PKC^{N}, 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/mm^{2} (PKC^{−}, 1070 ± 329; PKC^{N}, 798 ± 198; PKC^{+}, 1007 ± 120; mean ± SD). All networks produced synchronous bursting, regardless of their degree of clustering (Fig. 2*B*), 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 10^{4} and 10^{6} 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. 2*C*). To capture their fundamental shape, we normalized ASDs with respect to the initial third of each curve (Fig. 2*D*) and performed a principal component analysis. The first principal component (PC_{1}; Fig. 2*D*, inset) captured overall ASD curvature, explaining 86.4% of the variance (Fig. 2*E*, inset), which becomes apparent in the reconstruction of normalized ASDs using PC_{1} (Fig. 2*E*). PC_{1} continuously mapped the shape of normalized ASDs from concave to convex as function of the PC_{1} score with straight ASDs (power law) at zero (Fig. 2*F*). Distributions of the PC_{1} score differed significantly across PKC conditions [PKC^{−}, 3.7 ± 0.3 | 2.5, *p* = 6.8 * 10^{−11} (vs PKC^{N}), *p* = 4.3 * 10^{−22} (vs PKC^{+}); PKC^{N}, −0.7 ± 0.5 | 3.9; PKC^{+}, −5.4 ± 0.7 | 3.7, *p* = 1.8*10^{−7} (vs PKC^{N}); mean ± standard error of the mean (SEM) | SD, Student's unpaired two-sample *t* test]. In the morphologically analyzed networks, CI increased with PC_{1} score (Fig. 2*G*), 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 (PC_{1} scores between −2 and 2; Fig. 2*H*).

### 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. 3*A*), 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 PKC^{N} and PKC^{+} networks (Fig. 3*B*). Average correlation strength decreased with increasing clustering in the network. In addition, the distance dependence of the correlation became steeper with increasing clustering (Fig. 3*C*). 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 (*Q _{func}*) was significantly different between PKC conditions and their respective degrees of neuronal clustering [Fig. 3

*D*; PKC

^{−}, 3.2 * 10

^{−2}± 1.9 * 10

^{−3}| 1.4 * 10

^{−2},

*p*= 7.7*10

^{−4}(vs PKC

^{N}),

*p*= 2.4

^{−11}(vs PKC

^{+}); PKC

^{N}, 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 PKC

^{N}); 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

*Q*and PC

_{func}_{1}scores (Fig. 3

*E*) as well as between

*Q*and CI (Fig. 3

_{func}*F*), 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

*Q*(6.4*10

_{func}^{−2}) and close to power law ASDs (PC

_{1}score = 0). Variations of neuron density did not account for differences in

*Q*across networks (Fig. 3

_{func}*G*).

### Moderate functional modularity correlates with critical avalanche dynamics

To illustrate how *Q _{func}* was reflected in avalanche dynamics, we assigned networks of the different PKC conditions and respective ASDs (Fig. 4

*A*) into the three categories low (<33.3rd percentile), moderate and high (>66.7th percentile)

*Q*(Fig. 4

_{func}*B*). 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

*Q*group, with convex (supercritical) ASDs, whereas strongly clustered PKC

_{func}^{+}networks formed most of the high

*Q*group with concave (subcritical) ASDs. Moderately clustered PKC

_{func}^{N}networks, in turn, contributed most of the moderate

*Q*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. 4

_{func}*C*).

### 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 *Q _{func}* (Fig. 5

*A*,

*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

*C*, 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. 5

*C*, 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. 5

*D*).

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

*E*, 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. 5

*E*, 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. 5

*F*).

## 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 1*A*, 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