Self-organization of modular network architecture by activity-dependent neuronal migration and outgrowth

The spatial distribution of neurons and activity-dependent neurite outgrowth shape long-range interaction, recurrent local connectivity and the modularity in neuronal networks. We investigated how this mesoscale architecture develops by interaction of neurite outgrowth, cell migration and activity in cultured networks of rat cortical neurons and show that simple rules can explain variations of network modularity. In contrast to theoretical studies on activity-dependent outgrowth but consistent with predictions for modular networks, spontaneous activity and the rate of synchronized bursts increased with clustering, whereas peak firing rates in bursts increased in highly interconnected homogeneous networks. As Ca2+ influx increased exponentially with increasing network recruitment during bursts, its modulation was highly correlated to peak firing rates. During network maturation, long-term estimates of Ca2+ influx showed convergence, even for highly different mesoscale architectures, neurite extent, connectivity, modularity and average activity levels, indicating homeostatic regulation towards a common set-point of Ca2+ influx.


Introduction
Modularity is a fundamental design principle of neuronal systems and exists at the scale of cellular compartments, local circuits or interconnected brain areas. From a structural perspective, modularity can arise from inhomogeneities in the physical substrate that facilitate connectivity within a group of functional entities versus connectivity between such groups.
At the mesoscale level of local circuits, the cerebral cortex is organized in local clusters of tightly interconnected neurons (Feldman and Peters, 1974;Skoglund et al., 2004) that share common inputs and targets (Bosking et al., 1997;Voges et al., 2010), have similar functional properties (Ringach et al., 2016) and are thought to constitute a basic computational module (Buxhoeveden and Casanova, 2002;Casanova and Casanova, 2019;Mountcastle, 1997).
Although cortical architecture is largely genetically predefined at this level, blocking electrical activity during development disturbed the characteristic clustering of connections, suggesting that activity-dependent self-organization influences network modularity (Durack and Katz, 1996;Ruthazer and Stryker, 1996;Thompson, 1997). Intriguingly, computational models predict that modular connectivity, in turn, promotes spontaneous activity (Kaiser and Hilgetag, 2010;Klinshov et al., 2014;Mazzucato et al., 2015). Modularization and spontaneous activity may thus co-evolve in a self-enhancing process.
In early postnatal development, neuronal migration and neurite outgrowth are regulated by activity-dependent changes of the intracellular Ca 2+ concentration [Ca 2+ ] i (Kater and Mills, 1991;Komuro and Kumada, 2005;Spitzer, 2006;Zheng and Poo, 2007), suggesting that morphodevelopmental processes contribute to cellular Ca 2+ homeostasis (Zündorf and Reiser, 2011). Put simply, neurons would grow to increase neurite field overlap and the corresponding synaptic connectivity (Kossio et al., 2018;Shepherd et al., 2005;Stepanyants et al., 2002;Tetzlaff et al., 2010;van Ooyen et al., 1995) to establish the level of spike activity necessary to achieve some target value of [Ca 2+ ] i . As inter-neuron distance strongly affects the overlap of neurite fields and thus connectivity (Barral and D Reyes, 2016;Schnepel et al., 2015;Seeman et al., 2018), spatial clustering of neurons may play an important role in shaping modularity (Hernández-Navarro et al., 2017).
In the current study, we focus on the developmental self-organization that leads to different network architectures. In a simple computational model, varying the ratio of activity-dependent homeostatic growth versus migration was sufficient to modify neuronal clustering, mesoscale network organization, and the degree of modularity. Since controlled manipulation of network architecture and simultaneous activity monitoring is impractical in vivo, we tested this developmental interaction by modifying growth and migration in networks of cortical neurons in cell culture. These networks recapitulate major developmental processes such as cell migration and neurite outgrowth (Guan et al., 2007;van Huizen et al., 1987;van Pelt et al., 2004), develop varying degrees of clustering (Kriegstein and Dichter, 1983;Okujeni et al., 2017;Soriano et al., 2008;Teller et al., 2014) and produce a rich repertoire of spontaneous bursting events (SBE) (Kamioka et al., 1996;Okujeni et al., 2017;Wagenaar et al., 2006), similar to the developing cortex (Golshani et al., 2009;Minlebaev et al., 2007).
Increasing PKC activity in cultured networks amplified cell body clustering and local neurite entanglement at the expense of long-range connections, promoting local burst initiation and average firing rate (AFR) but reducing network recruitment during SBEs (Okujeni et al., 2017;Okujeni and Egert, 2019). This supports the theoretical predictions for modular networks mentioned above and is consistent with results from clustered networks created by mechanical constraints or modified growth substrates (Bisio et al., 2014;Tibau Martorell et al., 2018;Yamamoto et al., 2018).
Irrespective of network architecture, activity stabilized after approximately 21 days in vitro (DIV), suggesting that the target of homeostatic network development had been achieved. Different AFRs at this stage, however, conflict with previous studies assuming that AFR development reflects the homeostatic regulation of [Ca 2+ ] i (Abbott and Rohrkemper, 2007;Kossio et al., 2018;van Ooyen et al., 1995). Ca 2+ -influx, however, exponentially increases with membrane depolarization (Mazzanti et al., 1992) and thus depends on the temporal structure of spike activity. Our findings suggest that because of this non-linearity and specific differences in network-wide peak firing rates (PFR), long-term average Ca 2+ influx converges despite different AFRs and connectivity. Migration and neurite growth thus interact in a homeostatic process that defines the mesoscale architecture of neuronal networks.

Results
The connectivity between neurons depends on the overlap of their neurite fields and on their spatial distribution in the network. Like neurite growth, however, this distribution is dynamic because neurons migrate even in postnatal development. In a recurrent network, the input a neuron receives then depends on its embedding as well as the network's overall connectivity and activity structure. Here, we investigated how activity-dependent neurite growth and migration interact to establish connectivity and activity in neuronal networks.

Simulating activity-dependent neurite growth and migration
To gain insights into interdependencies between neurite growth and neuronal migration during the activity-dependent network self-organization, we extended a network growth model introduced by van Ooyen et al. (1995) that reproduces the outgrowth and subsequent pruning of neurites reported for developing neuronal networks (van Huizen et al., 1987;van Pelt et al., 2004). Following this, neurons were initially randomly seeded on a torus and their interconnectivity was modeled as degree of overlap between their circular neurite fields (no distinction was made between axons and dendrites). Input to neurons was calculated as the product of presynaptic firing rates and respective connectivity. A sigmoidal transfer function governed the relation between input-dependent membrane potential depolarization and firing rate ( Figure 1A). A growth process superimposed onto this framework allowed neurons to adjust their input by growing or shrinking their neurite fields, and thus the overlap with other fields, to establish a defined target firing rate ( Figure 1B,C). In addition to neurite growth, the final phase of neuronal migration observed in postnatal development is modulated by network activity and thus interacts with the formation of neurite fields and the regulation of connectivity. We therefore extended the original framework of the model by adding activity-dependent migration, where neuron somata migrated in the direction of the strongest input and gradually slowed down as their firing rates converged to the target level ( Figure 1B,D). In contrast to the bidirectional modulation of neurite fields, neurons were not repelled, however, if the activity level was above target. Prior to the formation of first contacts, migration was determined by erratic movements only. Neurons could thus increase their input by extending neurites and by migration to increase the overlap of neurite fields. The relative contribution of migration in network formation herein depended on its rate in relation to the net rate of neurite extension or pruning.

Migration and neurite outgrowth shape network architecture
Initially, neurite outgrowth ( Figure 2A) and migration ( Figure 2B) did not depend on activity. Once neurite fields began to overlap, directed migration towards areas that provided more input amplified statistical variations in the local cell density and led to clustering, indicated by decreasing clustering index (CI, Figure 2C). CI was calculated as the ratio between the average nearest neighbor distance in a network and the expected average nearest neighbor distance for random networks. CI above one indicates grid-like cell body arrangements and CI below one indicates clustering. Increasing clustering promoted connectivity buildup ( Figure 2D) and thus input to a neuron ( Figure 2E), which advanced the onset of spontaneous network activity ( Figure 2F). Migration and clustering of neurons ceased with the steep onset of network activity ( Figure 2B,C,F). In homogeneous networks, neurite fields had to grow larger than in clustered networks to establish the same degree of overlap and thus connectivity (Figure 2A,D). As a result, the size of neurites in mature networks correlated negatively with the degree of neuronal clustering (Figure 2-figure supplement 1). Connectivity, input activity and firing rates eventually converged to the same levels for different migration conditions ( Figure 2D-F).
Varying the rate of migration crucially impacted on the overall architecture of developing networks ( Figure 2I, Figure 2-videos 1-4). Without migration, networks developed the most homogeneous neurite field diameters and neurite coverage ( Figure 2G). Clustering led to more variable neurite field diameters as more isolated neurons required large fields to receive sufficient input, whereas within dense clusters, strongly overlapping neurite fields remained small.
The evolution of the largest connected subnetwork, that is the giant component, suggested that full network connectivity was established along the same developmental time line, irrespective of the degree of clustering ( Figure 2H, inset). In clustered networks, however, individual neurons played an important role in bridging subnetworks ( Figure 2I, arrows in the bottom panel). To quantify the tendency for modularity with different architectures, we calculated the giant component remaining after removing increasing subsets of randomly selected neurons in mature networks ( Figure 2H). In clustered networks, the giant component shrunk faster with an increasing fraction of neurons removed, demonstrating that individual neurons became critical bottlenecks in connectivity. Increasing activity-dependent migration relative to neurite growth thus increased the modularity (Q, Figure 2I) of the network.

Mesoscale network architecture in vitro
The growth model suggested that spatial clustering of neurons during development could play a crucial role in the formation of network connectivity by influencing the probability of neurites to overlap during outgrowth. We assessed this dependence experimentally by chronic activation or inhibition of PKC (PKC + and PKC À respectively), a regulator of neuronal migration, in developing networks of cortical neurons in cell culture. As described previously (Okujeni et al., 2017), PKC Figure 1. Model of activity-dependent network development. Neuronal wiring strategies may involve expansion of neurite fields and migration towards other neurons to increase connectivity modeled as neurite field overlap. (A) Transfer function of membrane depolarization between resting and maximal potential to firing rates. Dotted line: target firing rate. (B) Neurite growth (orange) and migration (green) were modulated as a function of [Ca 2+ ] i that corresponded to average firing rates. Neurites grew while the firing rate (corresponding to long-term average Ca 2+ influx) was below target and were pruned when above it. Migration rate decreased as neurons approached the target firing rate (dotted line). (C) The area of neurite field overlap, corresponding to connectivity in the model, can be increased by neurite outgrowth and neuronal migration towards neighboring neurons (D). DOI: https://doi.org/10.7554/eLife.47996.002 Figure 2. Model of activity-dependent growth and migration. (A) Activity-dependent growth produced a characteristic overshoot and subsequent pruning of neurite fields. The overall size of developing neurites decreased with increasing migration rates and clustering. (B) Mean migration distance of neurons after seeding (smoothed by 1 hr sliding average). (C) Migration promoted clustering of neurons, which saturated with the onset of network activity and neurite pruning (curves smoothed by 1 hr sliding average). All networks were initialized with the same spatial cell body distribution with CI Figure 2 continued on next page manipulation significantly altered the mesoscale architecture of networks with 600-800 neurons/ mm 2 ( Figure 3A), with striking similarity to mature networks generated with the growth model. Under control conditions (PKC N networks), networks appeared as inhomogeneous density landscapes with both, clustered and sparse regions ( Figure 3A, center panel). In particular in clustered areas, neurites formed tangles, which would increase the probability of local connections. Axons spanning several millimeters indicated monosynaptic connections between distant network regions. In comparison, PKC À networks with diminished migration had a more homogeneous distribution of cell bodies and coverage with dendrites and axons ( Figure 3A, left panel). Reduced fasciculation of neurites and a high density of long-range axons suggested a more isotropic embedding of neurons and more random-like connectivity. In turn, PKC + networks with enhanced migration had well delineated clusters of about 30-60 neurons with dense tangles of neurites that rarely reached into neighboring clusters ( Figure 3A, right panel), indicating high local connectivity and reduced inter-cluster connectivity.

Cell migration promotes neuronal clustering
To quantify the structural development, we seeded networks at lower densities of about 300 neurons per mm 2 that were more suitable for morphometric analyses (Figure 3-figure supplement 1). Within the first day of random seeding of neurons, rapid neurite outgrowth resulted in overlapping neurite fields between neighboring neurons. Simultaneously, neuronal cell bodies migrated across the substrate. Neuronal migration with concurrent outgrowth of neurites gradually increased neuron clustering within about three weeks in vitro ( Figure 3B). Chronic manipulation of PKC activity differentially modulated neuronal clustering during development ( Table 1). At 22 DIV, clustering was moderate in PKC N networks (CI = 0.75 ± 0.03) but significantly increased in the PKC + networks (CI = 0.67 ± 0.02, p=3.3*10 À2 ) and significantly reduced in the PKC À networks (CI = 0.88 ± 0.01, p=4.4*10 À4 ). CI did not change significantly after 22 DIV, indicating cessation of neuronal migration.
Note that the spatial patterning of somata depended on neuron density. Clusters in dense networks (~700 neurons/mm 2 at >22 DIV) typically contained 30-60 neurons (Okujeni et al., 2017), close to 1. Note that the fluctuations for zero migration results from the random jittering of neuron positions by half the cell body radius (6 mm). (D) Average connectivity increased more rapidly with stronger migration and clustering. (E) Input increased faster with increasing migration rate because clustering initially promoted connectivity. Input levels eventually converged. (F) Firing rates increased sharply once critical input levels were attained. Migration and clustering accelerated the onset of activity. With increasing migration, steps arise because of incremental integration and activation of clusters within the larger network. Note that clustering reduced the developmental overshoot of firing rates. (G) Moderate migration and clustering produced the highest variability of neurite field size across neurons in mature networks. (H) High migration rates increased modularity in mature networks. With increasing migration rate, the giant component more rapidly decreased in clustered networks when a certain fraction of neurons was randomly deleted, indicating that these networks break into disconnected subnets. Inset: the fraction of neurons in the giant cluster, that is the largest connected subnetwork, evolved similarly in different migration conditions. (I) Migration rates crucially determined the mesoscale architecture and modularity (increasing Q indicates stronger modularity) of developing networks. While average neurite fields were small in clustered networks, more isolated neurons generated larger fields (arrows) and formed bottlenecks for activity propagation by connecting otherwise unconnected or weakly connected subnetworks. DOI: https://doi.org/10.7554/eLife.47996.003 The following video and figure supplements are available for figure 2:

Clustering diminishes dendrite outgrowth
To address the interaction of neurite field extension, migration and clustering, we analyzed the average size of dendrites at several time points during development ( Figure 3C). Dendrite size was quantified as the ratio between the total length of detected dendrite stretches and the number of neurons within regions of interest ( Table 1). The measure estimates the average contribution of each neuron to the dendritic mesh. Chronic manipulation of PKC activity had little impact on dendrite size up to 8 DIV but significantly modulated dendrite outgrowth during subsequent development. At 22 DIV, dendrite size was significantly increased in the more homogeneous PKC À networks but significantly reduced in the more strongly clustered PKC + networks (PKC N : 1021 ± 41 mm; PKC À : 1413 ± 64 mm, p=7.9*10 À5 ; PKC + : 816 ± 24 mm, p=3.6*10 À4 ). In all conditions, dendrite size did not change significantly between 22 and 29 DIV, indicating stabilization of the dendritic network after the third week in vitro. As in the model (Figure 2-figure supplement 1), dendrite size in mature networks was negatively correlated with the degree of cell body clustering and, thus, the distance between neurons ( Figure 3D).

Dendrite outgrowth promotes synaptic connectivity
Network connectivity requires neurite overlap but further depends on the probability by which synapses are realized at axo-dendritic intersections. To assess how synaptic connectivity evolved in the different PKC conditions, we stained and detected presynaptic boutons (Figure 3-figure supplement 2) and determined the synapse density as the average number of presynaptic boutons per neuron ( Figure 3E, Table 1) and the dendritic occupancy as the number of synapses per unit dendrite length ( Figure 3F, Table 1). Manipulating PKC activity had no significant influence on early synaptogenesis up to 8 DIV, consistent with the comparable dendrite density in different PKC conditions at this stage. Paralleling dendritic outgrowth, synapse density increased significantly with increasing dendritic occupancy between 8 and 22 DIV in all conditions. Synapse densities and dendritic occupancy subsequently decreased between 22 and 29 DIV. This reduction was not significant in PKC À networks, however. Developmental manipulation of PKC activity profoundly affected mature synapse densities (PKC N : 1114 ± 56; PKC À : 2019 ± 110, p=4.5*10 À8 ; PKC + : 669 ± 21, p=5.7*10 À8 ) and dendritic occupancy (PKC N : 1.19 ± 0.04 mm À1 ; PKC À : 1.47 ± 0.04, p=3.0*10 À5 mm À1 ; PKC + : 0.9 ± 0.04 mm À1 , p=2.3*10 À5 ) at 29 DIV, both of which were significantly increased in the PKC À and reduced in the PKC + condition. Similar to dendrite densities, synapse densities were thus negatively correlated with the degree of clustering. Across PKC conditions and developmental stages, synapse In late development, dendrite size scaled inversely with the degree of clustering. For visualization the CI axis was inversed, so the degree of clustering increases from left to right. (E) The synapse density increased concurrently with dendrite growth. After 22 DIV synapse densities decreased in PKC N and PKC + networks, indicating synaptic pruning. (F) Dendritic occupancy with synapses differed slightly between conditions and decreased after 22 DIV. (G) The number of synapse per neurons increased with the dendrite size. Gray lines connect networks of the same age. The blue line illustrates a proposed quadratic scaling rule between dendrite size and synapse densities. (H) Neuron density declined with DIV to about one third of the seeding density. (I) Estimated upper bounds for connectivity based on the synapse density and the total number of neurons (on 113 mm 2 cover slips). PKC À at least doubled average connectivity. (J) In mature networks, maximum connectivity scaled inversely with clustering. All parameters are presented as mean ± SEM. Data from 4 to 24 images (Table 1, area 3.5 mm 2 ) taken in each of 2 networks per condition and age. Asterisks indicate p-values 0.05 (*), 0.01 (**) and 0.001 (***) tested against PKC N . DOI: https://doi.org/10.7554/eLife.47996.010 The following source data and figure supplements are available for figure 3: Source data 1. Source data and Matlab script for Figure 3B Table 1. Morphometric analysis of network development under different PKC conditions. Results are presented as mean ± standard error of mean (SEM). Significance was determined against PKC N , or between specified developmental time windows, using independent Student's t-test. N specifies the number of analyzed images taken from two networks per PKC condition and age. Neuron density 8 255 ± 6 (9.6*10 À7 ) 185 ± 11 168 ± 7 (2.0*10 À1 ) #/mm 2 15 214 ± 9 (5.9*10 À4 ) 131 ± 17 158 ± 12 (2.0*10 À1 ) #/mm 2 22 107 ± 8 (1.9*10 À1 ) 123 ± 6 85 ± 6 (5.7*10 À4 ) #/mm 2 29 87 ± 5 (3.0*10 À1 ) 96 ± 7 77 ± 4 (2.6*10 À2 ) #/mm 2 occupancy scaled approximately quadratic with the dendrite size ( Figure 3G), which could result from similarly modulated axonal densities (Okujeni et al., 2017) and the corresponding multiplicative increase in intersection probability.

Clustering reduces maximum global connectivity
Network connectivity is limited by the number of synapses per neuron and the overall number of neurons in a network since neurons obviously cannot have more partners than they have synapses.
The ratio between the number of synapses per neuron and the total number of neurons in the network defines an upper bound of connectivity for a network (maximum connectivity). The degree of connectivity realized, however, could be lower because of multiple structural synapses between neuron pairs. Although the density of neurons decreased during early development ( Figure 3H, Table 1), maximum connectivity increased significantly in all conditions between 8 and 22 DIV ( Figure 3I) and saturated between 22-29 DIV. At the same time, maximum connectivity almost doubled in PKC À networks compared to PKC N networks but was significantly reduced in PKC + networks (PKC N : 0.12 ± 0.01; PKC À : 0.23 ± 0.03, p=9.2*10 À4 ; PKC + : 0.08 ± 0.01, p=3.8*10 À2 ) and thus was negatively correlated with the degree of clustering ( Figure 3J).

Clustering decreases PFR and depolarization during SBEs
The hypothetical set-point of the homeostatic process, however, is not the firing rate per se but the associated [Ca 2+ ] i (Mattson and Kater, 1987), which is linked to molecular processes involved in growth and migration. Ca 2+ influx increases supra-linearly with increasing membrane depolarization (Mazzanti et al., 1992). This suggests that the long-term Ca 2+ gain is not a linear function of AFR but depends the depolarization of the membrane potential and thus on the temporal structure of activity. Depolarization depends on the number and synchronization of excitatory synaptic input, which becomes maximal during the peak phase of SBEs. Simultaneous intracellular and extracellular recording showed that higher SBE strength was indeed associated with stronger depolarization during SBEs (Okujeni et al., 2017). In PKC À networks, membrane depolarization high above spiking threshold frequently led to a depolarization block that outlasted the spike burst ( Figure 4D top trace). The fraction of time spent above threshold (À40 mV, Figure 4E) was significantly larger in neurons of PKC À networks (5.2 ± 0.7%, p=1.7*10 À4 , N = 30 neurons; mean ± SEM, independent Student's t-test) than in PKC N (1.7 ± 0.5%, N = 24) and PKC + (1.2 ± 0.7%, p=1.2*10 À3 , N = 24) networks (14-23 DIV). Depolarization was therefore not necessarily correlated with the individual firing rate of a neuron and the AFR in the network but rather reflected the network PFR during SBEs.

Homeostatic regulation of growth by long-term Ca 2+ influx
To assess how Ca 2+ influx depends on PFR, we determined the amplitude of Ca 2+ transients in excitatory neurons expressing GCaMP under the CAMKII promotor while simultaneously recording SBEs with MEAs ( Figure 5A). Most neurons indeed showed an exponential relation between PFR and the amplitude of Ca 2+ transients ( Figure 5B). PKC À networks realized much higher PFRs and had somewhat smaller exponents than PKC N (PKC N 0.12 ± 0.02, PKC À 0.11 ± 0.01, p=3.2*10 À18 ; Figure 5C,D, E).
In all network types, PFR increased steeply in early development and later declined concurrently with SBE strength. Throughout development, however, PFRs were highest in homogeneous networks and lowest in clustered networks ( Figure 5F, Table 2). Networks with low AFR thus had high PFR.  The following source data is available for figure 5: Source data 1. Source data and Matlab script for Figure 5B- Knowing the relationship between PFR and Ca 2+ influx allowed us to estimate Ca 2+ levels during development based on MEA recordings. We approximated the development of the average Ca 2+ influx per SBE ( Figure 5G) from their respective PFRs and the exponential Ca 2+ gain function with the average exponent of 0.11. Because higher PFRs, Ca 2+ influx per SBE was highest in the more homogeneous PKC À networks and lowest in clustered PKC + networks. Yet, in combination with the systematic increase of SBE rate with clustering, long-term Ca 2+ influx converged during late development for different PKC conditions, network architectures and AFR ( Figure 5H).

Differences in PFR reflect variations of network recruitment during SBEs
The predominately short-range connectivity observed in clustered PKC + networks could impair network-wide recruitment (Okujeni et al., 2017) and synchronization of activity. This would decorrelate inputs, explaining lower PFR and weaker membrane depolarization during SBEs. To test this, we determined network synchrony as the average spike correlations between all electrode pairs ( Figure 6A). Consistent with the rapid buildup of connectivity, network synchrony increased steeply between 3-15 DIV and reached stable levels already between 15-21 DIV, even though activity levels, connectivity and inhibition continued to develop. In line with connectivity estimates, synchronization was highest in PKC À networks (0.53 ± 0.04, p=4.6*10 À5 compared to PKC N ), intermediate in PKC N networks (0.35 ± 0.02) and lowest in PKC + networks (0.16 ± 0.03, p=1.7*10 À5 compared to PKC N ), that is network synchrony indeed decreased with the degree of clustering.

Discussion
Neuronal network architecture is not based on a genetic blueprint alone but is shaped by predefined rules of activity-dependent self-organization (Spitzer, 2006). Herein, neuronal migration (Komuro and Kumada, 2005;Zheng and Poo, 2007) and neurite outgrowth  are regulated by activity-related changes of [Ca 2+ ] i . Indeed, cell motility and growth is optimal within a narrow [Ca 2+ ] range and diminished otherwise, which led to the hypothesis that network connectivity and activity evolve under homeostatic control with the [Ca 2+ ] i as set-point parameter (Kater and Mills, 1991). However, basal cytosolic [Ca 2+ ] is very low due to efficient Ca 2+ -buffering and extrusion (Kater and Mills, 1991;Zündorf and Reiser, 2011) and remains relatively constant during development (Maravall et al., 2000). Free Ca 2+ for the regulation of growth is thus essentially determined by transient [Ca 2+ ] i elevations induced by synaptic input and spike activity. Accordingly, the developmentally attained spike rate was proposed to reflect the Ca 2+ set-point of growth (van Ooyen et al., 1995).
The overall capacity for neurite growth ultimately relies on gene expression for cytoskeletal building blocks, which crucially depends on nuclear [Ca 2+ ] (Berridge et al., 2000). Somatic membrane depolarization increases Ca 2+ influx close to the nucleus (Greer and Greenberg, 2008). In this context, intracellular stores like the endoplasmic reticulum can accumulate Ca 2+ over longer periods of time and then considerably amplify Ca 2+ signals by additional Ca 2+ -triggered release (Berridge et al., 2000;Pivovarova et al., 2002). This effectively acts as a low-pass filter and amplifier for Ca 2+ -signaling to the nucleus -modulating the expression of cytoskeletal proteins. Supporting this link, neurite tree morphology and size in different neuron types appear to depend on the expression of specific Ca 2+ -binding proteins that determine nuclear Ca 2+ buffering capacity (Mauceri et al., 2015). In contrast to nuclear Ca 2+ levels, local Ca 2+ transients in neurites direct migration and growth towards target neurons (Guan et al., 2007;Henley and Poo, 2004;Hutchins and Kalil, 2008), which promotes neurite overlap and synaptic connectivity (Shepherd et al., 2005;Stepanyants et al., 2002). Though local Ca 2+ influx activating PKC modulates cytoskeletal turnover involved in guided outgrowth and migration (Fogh et al., 2014;Kabir et al., 2001;Larsson, 2006), PKC may not be essential for constitutive neurite outgrowth (Flynn, 2013;Letourneau et al., 1987). We therefore speculate that local Ca 2+ transients and PKC activity regulate cytoskeletal motility to direct growth processes, whereas long-term accumulation of Ca 2+ in intracellular stores modulates signaling to the nucleus, transcription levels and thus the overall availability of cytoskeletal building blocks. This predicts a cessation of growth at a target longterm average Ca 2+ influx that is independent from PKC activity.

Migration contributes to homeostatic network development
Extending on growth models for homeostatic network formation based on activity-dependent neurite outgrowth, neuronal migration could likewise contribute to the regulation of connectivity and activity in developing networks. Eglen et al. (2000) already added migration implemented as repulsion between neurons to the neurite growth model by van Ooyen et al. (1995) to generate regular neuronal arrangements as observed in dense retinal cell mosaics. We showed that activitydependent attraction between migrating neurons leads to different degrees of modularity by the interaction of clustering with homeostatic regulation of neurite growth. While it is plausible to assume that neurons with small neurite fields and little connectivity may move, this seems less realistic once they are enmeshed in the network. In line with this, cell migration relies on localized Ca 2+ transients in leading neurites and the resulting Ca 2+ gradients across the cell (Guan et al., 2007) but ceases with increasing neuronal activity and frequency of Ca 2+ transients (Bando et al., 2016). We approximated this in the model by allowing attraction while input was below the setpoint but omitted repulsion with input above the set-point. In consequence, cell migration ended during the rapid increase of activity during development, similar to peaks in PFR and Ca 2+ influx, and cessation of clustering around 10-15 DIV ( Figure 3B, Figure 5F,H). Moreover, with rapid transitions to high network activity once neurite fields in the network overlapped sufficiently, the model showed a transient overshoot of connectivity. A more gradual build-up of activity diminished the average overshoot and pruning when the slope of the sigmoid mapping input to firing rate was reduced (Figure 2-figure supplement 2), in agreement with reports of varying degrees of growth overshoot or even saturating growth during development in vitro (Ito et al., 2013;Kondo et al., 2017;van Pelt et al., 2004). Neurons that connected to the network early in development, however, still showed an overshoot of connectivity, in agreement with Kossio et al. (2018).

Average Ca 2+ influx converges for different network architectures
Homeostatic regulation of growth processes by Ca 2+ was proposed to guide network development towards target firing rates (van Ooyen et al., 1995), which implies a quasi-linear relationship between Ca 2+ influx and AFR. In our model, connectivity, input activity and firing rates eventually converged to the same levels for different migration conditions and network architectures. In apparent conflict with the simulation, we found that different network architectures stabilized in vitro after about 3 weeks but with different AFR. Consistent with theoretical studies predicting that network modularity promotes spontaneous activity (Kaiser and Hilgetag, 2010;Klinshov et al., 2014;Mazzucato et al., 2015), SBE rates and AFRs increased with the level of clustering. Clustering, however, reduced network synchronization, lowered PFRs and weakened depolarization during SBEs. This strongly affected Ca 2+ -transients: Ca 2+ peak amplitude increased exponentially with PFR during SBEs, in agreement with reports of Ca 2+ currents through voltage-gated Ca 2+ channels increasing exponentially with depolarization (Mayer et al., 1987;Mazzanti et al., 1992). Because of the opposite modulation of SBE rates and PFRs with clustering, however, the estimated long-term Ca 2+ -gain converged for different network architectures during development, despite different AFR. The low spike rates during inter-burst intervals had negligible influence on Ca 2+ influx.
To account for the supra-linear increase of Ca 2+ with PFR we would need to use spiking neurons in our model. In addition, Ca 2+ influx would need to depend on the membrane potential, instead of on the average spike rate of a neuron as in extensions of the growth model with spiking dynamics (Abbott and Rohrkemper, 2007;Kossio et al., 2018). To accelerate the simulation of several weeks of network development, these studies initially increased the neurite growth rate and thus effectively decreased the temporal resolution until the networks approached the equilibrium state. The mesoscale structures forming in our networks, however, crucially depended on the continuous feedback between migration and neurite growth and activity. Low temporal resolution in the simulation would amount to a large decrease of the feedback speed, which leads to a random walk of neurons and more homogeneous network structures without clustering.

Interaction between growth and migration shapes network modularity
Increasing the rate of activity-dependent migration in the model promoted clustering, decreased neurite fields and accelerated the development of spontaneous activity by more rapidly increasing neurite overlap and connectivity. This resulted in network architectures covering a continuous gradient from homogeneous via partially clustered with scattered neurons to fully clustered networks with corresponding degrees of modularity. This was remarkably similar to the development in vitro, where PKC activity promoted clustering and SBE rates, and decreased neurite density. The model suggests that different network architectures can arise spontaneously based on simple rules regulating connectivity to achieve a target level of [Ca 2+ ] i .
Among the grand average developmental time courses of the most relevant aspects across all conditions, long-term Ca 2+ influx was the first property to peak while the impact of inhibition on network activity only started to increase when Ca 2+ influx stabilized ( Figure 6C).

Growth and migration shape the framework for synaptic connectivity
In our networks, synapse densities scaled approximately quadratically with the average dendrite size and thus negatively with the degree of clustering. This could be explained by the co-modulation of axonal and dendritic densities in the same direction (Okujeni et al., 2017), which multiplicatively increases the number of axo-dendritic contact sites, rather than their modulation in opposite directions as used in Tetzlaff et al. (2010). Such potential synapses realize into functional synapses with approximately constant probability in vivo (Stepanyants et al., 2002). The consistent relation of synapse density and dendrite size across developmental stages and PKC conditions ( Figure 3G) suggests that PKC manipulation did not critically impair synaptogenesis. Our estimates of maximum connectivity suggest a saturation of connectivity towards 10% in clustered and 20% in homogeneous networks, in the range of values reported for cultured (Marom and Shahaf, 2002) and native cortical networks (Feldmeyer, 2012).
The mesoscale network architecture formed early thus appears to determine the probabilistic framework for connectivity. PKC activity additionally influences synaptic plasticity, yet without general directionality towards LTP or LTD (Chung et al., 2000;Ferreira et al., 2011;Lan et al., 2001;Boehm et al., 2006;;Scott et al., 2007). Our model indirectly accommodates this influence. For example, synaptic depression, corresponding to reducing the synaptic weight factor s, would extend the outgrowth phase to increase connectivity and input necessary to reach the target level of [Ca 2+ ] i . Conceptually, this would be the inverse of the homeostatic scaling of synaptic weights with the level of connectivity (Barral and D Reyes, 2016;Okujeni et al., 2017;Wilson et al., 2007). This contribution of synaptic plasticity to the activity-dependent fine-tuning of connectivity likely gains importance with increasing developmental age and structural complexity of a network.

Conclusion
Based on our findings, we propose that interactions between neurite growth and neuronal migration affect the balance between local and global connectivity, thereby shaping network modularity. Cell migration defects were also proposed as a pathogenic mechanism involved in several neurological conditions associated with altered size and spacing of mini-columns in the cortex, aberrant neurite growth and hyper-or hypo-connectivity (Catts et al., 2013;Courchesne and Pierce, 2005;Di Rosa et al., 2009;Donovan and Basson, 2017;Fan et al., 2013;McKavanagh et al., 2015), suggesting that the mesoscale network organization could be a critical factor. The associated degree of modularity thus appears to have crucial impact on activity generation, propagation and perpetuation, neural synchronization as well as network function and dysfunction.

Network growth model
We adopted and modified the model of activity-dependent network growth introduced by van Ooyen et al. (1995). All simulations were carried out with Matlab (version 2017a, Mathworks, Natick, MA, USA; code available at doi 10.5281/zenodo.3459678). Networks were initialized by randomly seeding 500 neurons onto a torus surface of 1 mm 2 to avoid boundary effects. Newly introduced neurons conflicting with the minimal neuron distance of 12 mm, approximately the size of cell bodies, were discarded and the procedure continued until the required neuron density was obtained.
Neurite fields were modeled as circular fields, centered at cell bodies and were initiated with a radius of 12 mm. Connectivity between neurons W was nonsymmetrical and defined as the area A of neurite field overlap normalized by the area of the presynaptic neuron, which reflected the probability that dendrites of neuron i overlapped with the axons of presynaptic neuron k.
The gain s = 0.1 was chosen such that it produced networks with an intermediate degree of neurite field overlap (for s = 1, neurons would only connect to one or a few other neurons). Instead of simulating network growth with dimensionless equations (van Ooyen et al., 1995), we adjusted the time steps such that we could compare the dynamics to realistic developmental timescales. We estimated the loop-time across which activity is integrated based on the time constants for the accumulation of Ca 2+ in intracellular stores to be in the order of minutes (Pivovarova et al., 2002) and therefore set the temporal resolution of the simulation to 1 min.
Since inhibition is not explicitly relevant to the questions addressed here, we adapted the model for excitatory networks only. Long-term integration of activity in neurons was described by their state variable x i (ranging between 0 and 1), which increased with input from presynaptic neurons contributing with their firing rate f x k ð Þ times the synaptic strength W ki : where dt= t= 1 min was the time resolution of the simulation, corresponding to the time constant of long-term integration of activity. A sigmoidal transfer function for the depolarization state x i determined the firing rate f x ð Þ.
Þ=a where = 0.5 reflected the firing threshold and a = 0.12 determined the steepness of the function that crucially impacted on the developmental overshoot of connectivity and subsequent pruning of neurites. We chose a slightly shallower function than the original model by Van Ooyen (a = 0.1) to accommodate the degree of overshoot and pruning for cultured networks in recent reports (Ito et al., 2013;Kondo et al., 2017;van Pelt et al., 2004).
As in the original model by Van Ooyen, neurons were modeled to grow neurites and thereby increase input activity and firing rate to reach a target [Ca 2+ ] i . If this Ca 2+ level was surpassed, neurites were pruned, in turn. These bidirectional changes in the radius R of circular neurite fields, were determined by a sigmoidal function of the firing rate of a neuron multiplied with a fixed growth rate growth .
where " = 0.6 defined the target level for activity or [Ca 2+ ] i , b = 0.1 determined the steepness of the sigmoidal function and growth was the constant factor for the growth rate of neurite fields. We assumed that connectivity is mainly determined by the density of neurites rather than their maximal length. Given the homogeneous density of the neurite field used in the model, however, its radial expansion must be considerably slower than the average elongation rates of individual dendrites, which were reported to be 12 mm/day for isolated neurons in the first week in vitro . We therefore set growth = 4 mm per day.
In our model, neurons additionally migrated in the direction of presynaptic inputs, thus mimicking the guidance of migration by leading processes (Flynn, 2013;Guan et al., 2007) and consistent with the positive correlation between the rate of soma translocation and the amplitude and frequency of Ca 2+ transients (Komuro and Kumada, 2005;Zheng and Poo, 2007). We assumed synaptic activity in leading neurites as an important source of input, however, did not preclude contactmediated Ca 2+ signaling (Sheng et al., 2013), which may contribute in regulating migration early in development when activity levels are low. Changes in the spatial position of neuronal cell bodies S were caused by migration impulses that depended on [Ca 2+ ] i and, thus, on the firing rate f and a variable factor for the maximal migration rate migration .
where migration ranged 0-300 mm/day and m = -15 determined how strong migration impulses were diminished as neurons reached their target Ca 2+ level. We chose m to result in a negligible migration impulse at the target [Ca 2+ ] i . This mimicked a realistic migration process in which neurons are guided by local Ca 2+ transients in leading neurites and the resulting Ca 2+ gradients across the cell (Guan et al., 2007), but at the same time cease migrating when spiking-based Ca 2+ transients start to dominate (Bando et al., 2016). The migration speed of postnatal neurons in vitro indeed decays approximately exponentially during development from 0.7 mm/min (1008 mm/day) at 0 DIV tõ 0.05 mm/min (72 mm/day) at 12 DIV on Matrigel-coated substrates and with slower initial migration speeds of 0.1 mm/min (144 mm/day) on PEI coated substrates (Sun et al., 2011), as used in this study. In the model we varied migration rates within this range.
The direction of movement was determined involving a directed movement component and a random movement component to match erratic movements observed in time lapse videos. Movement direction of the directed component was determined by the vector sum v dir of direction vectors v ik that pointed to presynaptic neurons and were weighted by their input.
To obtain the final direction vector V, directed and the random component (updated every 10 min) were weighted (p = 0.9) and summed. The random directional component was necessary to mimic the erratic movement patterns observed in in vitro time lapse studies.
New neuronal cell body positions P were determined by multiplying the normalized final direction vector with the migration impulse.
In addition, neurons were set to jitter randomly around their current position by maximally their cell body radius to allow neurons to pass each other in the 2D simulation, which prevented unrealistic chains of neurons. This positional jitter decreased according to the exponential decay function modulating migration in dependence of [Ca 2+ ] i such that neurons stopped moving when reaching the target value. It was reset after each time step. Movements violating the minimal possible intersoma distance (12 mm) were discarded.
To assess the modularity of a network, we calculated the size of the largest subnetwork (the giant component) remaining after removing defined fractions of randomly selected neurons from the network as its fraction in the remaining total population. For each network, the results were averaged across 1000 repetitions of the procedure. We quantified the degree of modularity Q in the final networks based on the connectivity matrix using the Louvain method (Blondel et al., 2008) implemented for MATLAB by Mika Rubinov with gamma = 1 (Rubinov and Sporns, 2010). Q increases towards one with increasing modularity. Random networks yield Q = 0.

PKC modulation and disinhibition
PKC inhibitor Gö decke6976 (Gö 6976, 1 mM; Tocris Bioscience, Bristol, UK) and PKC agonist Phorbol-12-Myristate-13-Acetate (PMA, 1 mM; Sigma-Aldrich) 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%. GABAergic transmission was probed by acute application of the non-competitive GABA-A receptor antagonist Picrotoxin (PTX; 10 mM; Tocris Bioscience) during electrophysiological recordings. Recordings of spontaneous activity were started 10 min after application of PTX for 1 hr at different DIV. Changes of spike activity were calculated as mean burst strength across 1 hr with PTX vs. 1 hr baseline recording before application. Networks exposed to PTX were discarded.

Morphological analyses
The development of neuronal clustering, dendrite outgrowth and synapse densities was analyzed in sparse networks of~100 neurons/mm 2 that were more accessible for quantitative morphological analysis. Clustering of neuronal cell bodies was analyzed based on immunocytochemical staining of neuronal nuclei (NeuN; Rabbit-anti-NeuN, 1:500; Abcam, Cambridge, UK, RRID:AB_2744676) and of all cellular nuclei (DAPI; Sigma-Aldrich). Neuronal nuclei were detected based on NeuN and DAPI colocalization and evaluated for their degree of clustering using a modified Clark-Evans clustering index (CI) that accounts for cell body diameter as minimal possible inter-neuron distance (Clark and Evans, 1954;Galli-Resta et al., 1999;Okujeni et al., 2017). CI was calculated as the ratio between the average nearest neighbor distance in a network and the expected average nearest neighbor distance for random networks. Note that the degree of clustering increases with decreasing CIs below 1. CIs above one indicate grid-like cell body arrangements. Dendrite morphology was examined by immunocytochemical staining of microtubule-associated protein 2 (MAP2, Chicken-anti-MAP2; 1:500; Abcam, RRID:AB_2138147). To quantify the total length of dendrites, MAP2 images taken at 20-fold magnification (0.323 mm/pixel) were processed by median filtering (3 Â 3 kernel), background subtraction (lowest value in 7 Â 7 pixel field), contrast adjustment (saturation at highest and lowest 10%), thresholding and skeletonization of the resulting binary image, similarly to Pani et al. (2014). Synapses were detected based on an immunohistological staining of the presynaptic protein synapsin (Mouse-anti-Synapsin; 1:200; Synaptic Systems GmbH, Gö ttingen, Germany, RRID:AB_ 887805). Synaptic punctae were then determined by local maximum detection in high-pass filtered and contrast-enhanced images. We analyzed two networks per condition and age taken from images covering approximately 3.5 mm 2 . In each image, we typically analyzed 10-20 regions of interest with varying size (could overlap) and including dense and sparse network regions. The following measures were determined as the slope of the linear regression through data pairs from all regions of interest: Dendrite size, total length of dendrite stretches relative to the number of neurons; Synapse density: average number of synapses relative to the number of neurons; Dendritic occupancy: average number of synapses relative to the total length of dendrite stretches; Neuron density, average number of neurons per area; Maximum connectivity, ratio between the number of synapses per neuron and the total number of neurons in the network (extrapolated for the entire network area of~1.1 cm 2 given the image neuron density). All morphometric analyses were done with Matlab (versions 2014a -2017a). Results are presented as mean ± standard error of the mean (SEM) and significance was assessed with a two-tailed independent Student's t-test. Network architectures of dense networks (600-800 neurons/mm 2 ) were characterized qualitatively at 22 DIV with antibodies against MAP2 and phosphorylated neurofilament 200 kD (Rabbit-anti-neurofilament; 1:10; Abcam, RRID:AB_ 448148) to visualize dendritic and axonal compartments, respectively.

Extracellular recording and analyses
MEA recordings (MEA1060-BC and USB-MEA256-Systems; MCS, 25 kHz sampling frequency, 12 bit AD-conversion; MCRack software versions 3.3-4.5, RRID:SCR_014955) of multi-unit spike activity from individual networks were performed under culture conditions (37˚C, 5% CO 2 ) and lasted at least 1 hr. Action potentials were detected with a threshold set to À5 standard deviations of the high-pass filtered baseline signal (Butterworth 2 nd order high pass filter, 200 Hz cut-off; detection dead time 2 ms).
Raw data from MEA recordings was imported into Matlab using MEA-Tools (Egert et al., 2002) and the FIND toolbox (Meier et al., 2008). Spontaneous SBEs were detected as follows: Series of spikes with consecutive inter-spike intervals smaller than a threshold value (100 ms) were detected as bursts. SBEs were defined from periods in which a predefined fraction of electrodes showed simultaneous bursts (10% of all sites detecting spikes but minimally 3 and maximally 20 sites to keep criteria comparable between small and large MEAs). To account for buildup and fading phases of SBEs, spikes within a time windows of 25 ms prior to and following this SBE core were included into the SBE. Network activity was characterized by the following parameters: SBE rate in the recording period, SBE strength as the average number of APs per SBE divided by the number of electrodes with spikes at any time during the recording session (active sites); AFR as the grand average firing rate per active site during the recording session. PFR was calculated per SBE as the peak of the network-wide firing rate profile (box car filter applied to the global spike train; 0.2 s kernel width) divided by the number of active sites. Network synchrony was determined as average spike train correlation (30 ms bin width) between pairs of active sites.
For the developmental analysis of network activity, recordings from many networks were pooled within time windows of increasing width to account for the slowing development of activity dynamics as networks matured ( Table 2). Numerical results are presented as mean ± SEM and significance was assessed with a two-tailed independent Student's t-test.
For acute experiments with PTX, we defined as control period the last 1 hr section before application of PTX and excluded the first 10 min after application from the analysis to avoid transients due to handling. To determine the time course of the maturation of inhibition, changes in SBE strength following PTX application were quantified relative to the control period for different DIV. For visualization, trend lines were calculated with a sliding average (±7 DIV).
Data sets of at least 20 min were analyzed with Matlab. Membrane potential distributions for neurons with resting potentials between À64 ± 4 mV were determined for the entire recording period and averaged across neurons of the same PKC condition.

Calcium measurements and analyses
To assess neuronal Ca 2+ dynamics, cultures were transfected with AAV (=Adeno Associated Virus) vectors coding for GCaMP6s (AAV9.CAG.GCaMP6s.WPRE.SV40, titer:~10 11 ; Penn Vector Core, School of Medicine Gene Therapy Program, University of Pennsylvania) under control of the CAG promotor after 10-14 days in vitro. Ca 2+ dynamics were imaged at 20x magnification and 25 Hz frame rate (Examiner Z1 microscope, Zen software 2015, Carl Zeiss, Jena, Germany). Somatic regions were delineated by threshold detection in maximum projections of the Ca 2+ -movie with ImageJ (Schneider et al., 2012). The resulting regions of interest were corrected manually. Changes in the Ca 2+ signal DF/F were calculated as relative change to baseline following (Jia et al., 2011). For each SBE, the peak of the Ca 2+ signal (DF/F) within 200 ms after onset was related to the PFR determined from simultaneous MEA recordings. The exponential scaling between DF/F and PFR was assessed by fitting with the function DF=F ¼ e kÃPFR À 1 using the Matlab function fminsearch. Ca 2+ data were derived from five PKC N and four PKC À networks at 19-20 DIV in recordings of~30 min and analyzed with Matlab. Ca 2+ influx during SBEs was estimated as e 0:11ÃPFR À 1 to match the scaling found experimentally. Long-term Ca 2+ influx was approximated as the Ca 2+ influx integrated over all SBEs per hour. All results are presented as mean ± SEM. Significance was tested with a two-tailed independent Student's t-test.
German Research Foundation (DFG) and the University of Freiburg in the funding programme Open Access Publishing. The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Ethics
Animal experimentation: 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).

Data availability
Matlab code and source data files are provided for Figures 3-6. Data preprocessing is described in the methods. As the unprocessed data is considerably heavy (over 1TB), the raw data and analysis tools will be provided upon request. Code for the model simulation is available at https://dx.doi. org/10.5281/zenodo.3459678.
The following dataset was generated: Author (