Research ArticleNEUROSCIENCE

Impact of modular organization on dynamical richness in cortical networks

See allHide authors and affiliations

Science Advances  14 Nov 2018:
Vol. 4, no. 11, eaau4914
DOI: 10.1126/sciadv.aau4914


As in many naturally formed networks, the brain exhibits an inherent modular architecture that is the basis of its rich operability, robustness, and integration-segregation capacity. However, the mechanisms that allow spatially segregated neuronal assemblies to swiftly change from localized to global activity remain unclear. Here, we integrate microfabrication technology with in vitro cortical networks to investigate the dynamical repertoire and functional traits of four interconnected neuronal modules. We show that the coupling among modules is central. The highest dynamical richness of the network emerges at a critical connectivity at the verge of physical disconnection. Stronger coupling leads to a persistently coherent activity among the modules, while weaker coupling precipitates the activity to be localized solely within the modules. An in silico modeling of the experiments reveals that the advent of coherence is mediated by a trade-off between connectivity and subquorum firing, a mechanism flexible enough to allow for the coexistence of both segregated and integrated activities. Our results unveil a new functional advantage of modular organization in complex networks of nonlinear units.


Modular organization is ubiquitous in many real-world systems, including biological networks such as the brain (13). Recent connectome analyses have revealed that such an organizational principle is conserved evolutionarily in animal brains, across multiple species (3, 4) and spatial resolutions (3, 5, 6). Modular organization is considered to be advantageous, since it is naturally small world, which facilitates robustness, adaptiveness, and resilience to damage (3, 7, 8). Computational studies have pinpointed that modularity underlies the emergence of complex spatiotemporal dynamics that are essential for cortical processing (2, 9) as well as rapid switching between locally segregated processing and globally integrated information flow with minimum wiring costs (3, 10). Integration is associated with the capacity of a neuronal network to operate as a whole and efficiently exchange information, while segregation refers to the capacity to distribute information onto localized communities that perform specialized tasks.

The inherent impossibility to architect different wiring schemes in natural neuronal tissues, however, has hindered the investigation of the causal relationships between modular organization and functional traits and called for the design of well-controlled experiments. Recently, microfabrication has emerged as a novel technology to engineer the growth of cellular systems (11), including neurons derived from the mammalian brain (1217). Such a technology has paved the way toward the construction of in vitro modular networks to investigate the relationship between modular architecture and function (1517), making it possible to directly address the question of how multicellular activity is synchronized among modules to realize integration-segregation balance (10, 18).

Here, we used this bottom-up approach to construct modular cortical networks in vitro and studied their functional significance. We show that dynamically rich spontaneous activity, in which both segregated and integrated dynamics coexist, emerges only when the intermodular coupling is at the verge of disconnection. Through computational modeling, we show that the latter is achieved via an excitatory feedback loop that amplifies inputs to a module, thus effectively reducing the threshold for neuronal inputs, or quorum (19, 20).


Modularity and dynamical richness

We considered a simple modular design consisting of four spatially segregated areas, which are represented as 200 μm by 200 μm squares and are connected either by zero, one, or three lines (Fig. 1A). We refer to these configurations as the no-bond, single-bond, and triple-bond patterns, respectively. The width of the line was designed to be 5 μm, on which multiple neurites can grow in either direction. The structures were fabricated as micropatterns of cell-adhesive proteins, which were stamped on a glass coverslip by microcontact printing (Fig. 1A). The micropatterns support the growth of primary cortical neurons within a predefined area (Fig. 1B). Immunostaining for axonal and somatodendritic markers confirmed the random growth of axonal and dendritic processes within the square unit with an occasional outgrowth to the interconnecting lines (fig. S1A). As a control structure, we fabricated a pattern in which the four areas were merged into a 400 μm by 400 μm square (merged; Fig. 1, A and B). A graph theoretical analysis of the patterns showed that their modularity, i.e., the tendency to exhibit distinct communities or modules (1, 21), increased as connections between the four areas were reduced (fig. S2), and this occurred in the order of merged, triple-bond, single-bond, and no-bond patterns.

Fig. 1 Engineering the modularity of living cortical networks.

(A) Confocal imaging of the polydimethylsiloxane (PDMS) stamps used for microcontact printing of the proteins. From left to right: merged, triple-bond, single-bond, and no-bond patterns. (B) Phase-contrast micrographs of cortical networks at 10 DIV. Modularity increases in the order of merged, triple-bond, single-bond, and no-bond networks. (C) A cortical network on the single-bond micropattern loaded with the fluorescence Ca indicator, Cal-520 acetoxymethyl ester (AM). Scale bar, 200 μm. (D) Fluorescence signals from a 20-min recording of the spontaneous activity of a single-bond network. Traces of eight neurons are shown, and the colors distinguish the affiliated modules. Black dots mark detected spikes. In this recording, network bursts were observed twice, in addition to 10 bursts that involved only two modules. In the third, fifth, and seventh traces (from the top), small random events were also observed, which, most probably, are sporadic action potentials. (E) The average number of network bursts per 20-min session. Note that for the no-bond network, the activity of individual modules was isolated (see fig. S4B), and an event refers to synchronized firing within each 200 μm by 200 μm island. n = 22 (merged), 21 (triple bond), 23 (single bond), and 20 (no bond). Error bars represent the SEM. (F) Distributions and cumulative distributions (inset) of IBI in the merged, triple-bond, and single-bond networks. Cum. prob., cumulative probability.

We investigated the effect of spatial modularity on network dynamics using fluorescence calcium imaging (Fig. 1C). In all patterns, the neurons formed a sufficient number of synapses by 10 days in vitro (DIV) to spontaneously generate coherent activity in the form of network bursts (Fig. 1D and fig. S1B) (14, 19) [merged, 77% of recorded cultures (n = 22); triple bond, 90% (n = 21); single bond, 84% (n = 23)]. In no-bond networks, we observed bursts in 50% (n = 20) of the cultures, with no correlation among modules. Since the number of neurons in a network affects its dynamics (14), we maintained the neuronal population constant at approximately 100 cells per culture (fig. S3).

Although all of the patterns supported network bursts, their temporal structure varied. As shown in Fig. 1E, the number of network bursts gradually decreased with increasing modularity, and a comparison of the distribution of interburst intervals (IBIs)—a signature of the temporal structure of the events—revealed that the IBI distribution of the single-bond networks, IBIsingle, deviated from that of the merged (IBImerged) and triple-bond (IBItriple) networks (Fig. 1F). We quantified these differences using the Jensen-Shannon divergence, DJS, and obtained DJS(IBIsingle||IBImerged) = 0.032, DJS(IBIsingle||IBItriple) = 0.036, and DJS(IBImerged||IBItriple) = 0.017. The differences among the three distributions were statistically significant, with P < 0.001 for IBIsingle against both IBImerged and IBItriple and P < 0.01 for IBImerged against IBItriple (two-sample Kolmogorov-Smirnov test). The broader distribution of IBIs in the single-bond network indicates a higher aperiodicity and temporal complexity of the events. The duration of individual bursts in each neuron, in contrast, was affected only weakly (fig. S4A).

Increased variability in the number of participating neurons in bursting events was another aspect that marked the dynamics of the single-bond pattern. This behavior is illustrated in Fig. 1D, in which the coherent activity switched between two and four modules. To quantify this phenomenon, we evaluated the correlation coefficient (CC; see Materials and Methods) averaged over neuronal pairs affiliated to different modules, i.e., intermodular correlations (fig. S4B). A comparison showed that the mean intermodular CC was more broadly distributed for the single-bond pattern than for the merged or triple-bond patterns, whose CCs were largely clustered at ~1.

We introduce dynamical richness Θ to provide a direct measure for this spatiotemporal variability. This quantity, as well as similar ones termed “complexity” in the literature (22, 23), reflects the ability of a neuronal network to exhibit a broad range of dynamical states. Intuitively, this richness characterizes a state in which neurons coactivate in groups of varying size and temporal occurrence and with Θ = 0 corresponding to a scenario of random activity or persistent whole-network synchronization. We define Θ as a combination of two assays, namely, the variability among the activity patterns of firing neurons and the variability among the size of global network activations (GNAs; see Materials and Methods). For the first assay, we analyzed the pairwise CCs, rij, to obtain the distribution p(rij). A distribution clustered at rij = 1 indicates fully correlated neuronal activity, and that clustered at rij = 0 indicates uncorrelated neuronal dynamics. Similarly, for the second assay, we analyzed the time series of GNA occurrences, Γt, to obtain the distribution pt) of coherent network activations, with a distribution clustered at Γt = 1 indicating persistent coherent activations and that clustered at Γt = 0 indicating the absence of any coactivation. The difference between the measured distributions and these extreme cases reflects the rich dynamics in the network. Following Zamora-López et al. (24), Θ is computed asEmbedded Image(1)where |·| denotes the absolute value and m = 20 is the number of bins used for estimating the distributions. ΘCC and ΘGNA indicate the richness associated with each assay. Our measure Θ is robust to errors in the inference of neuronal spike trains and to the presence of random neuronal activations.

Figure 2A shows the results of this analysis. The merged and triple-bond configurations exhibit highly periodic coherent activations and nearly identical activity patterns. The associated distributions yield ΘCC ≅ 0.2 to 0.3 and ΘGNA ≅ 0.1 to 0.2, and a final dynamical richness Θ ≅ 0.06. The single-bond pattern, however, exhibits variable activations that range from a few neurons in a spatial module to the entire network. CC and GNA are more broadly distributed, with ΘCC ≅ 0.61 and ΘGNA ≅ 0.56, and a final Θ ≅ 0.34.

Fig. 2 Dynamical richness and modularity.

(A) Computation of dynamical complexity in representative networks with gradually higher modularity, from merged to single bond. Left: Raster plot of spontaneous activity (violet) and corresponding GNA (red). Strong coherence in the merged network is revealed by the repeated network activations of size 1. Center: Matrix of pairwise CCs between neurons (i, j). The strong network coherence provides CC values close to 1. Right: Probability density functions (pdfs) of CC values rij (violet) and GNA values Γt (red-yellow). The pdfs for the merged network peak toward 1 and provide small values of ΘCC and ΘGNA, with a final dynamical richness of Θ ≅ 0.07. For the triple-bond modular network, the strong connectivity among modules facilitates integration and Θ ≅ 0.06. The single-bond modular network exhibits a much richer variability in activity patterns, with switching activations encompassing two or four modules. This variability is reflected in the GNA values and the matrix of CC values and provides broad pdfs that lead to relatively high values for both ΘCC and ΘGNA, and with Θ ≅ 0.34. (B) Distribution of GNA values for all observed concurrent neuronal activations, comparing the merged (M; 13 cultures and 345 activations), triple-bond (3-b; 13 cultures and 336 activations), and single-bond (1-b; 18 cultures and 529 activations) patterns. The cartoons on the right illustrate the relative size of the coherent activations (yellow glow). (C) Average Θ values. The single-bond pattern displays significantly higher Θ values than the merged or triple-bond ones. Merged and triple bond are not significantly different. Data are shown as means ± SEM. ***P < 0.001; **P < 0.01; n.s., not significant (two-tailed Student’s t test).

The high dynamical richness of the single-bond pattern was a widespread feature, as shown in Fig. 2B, where the GNA sizes for all of the observed activity episodes are compared among patterns. While the merged and triple-bond patterns essentially display full network activations or random events, the single-bond one practically covers the entire range of coactivation sizes. On average, the single-bond dynamical richness was significantly higher than the merged or the triple-bond ones, as shown in Fig. 2C. As a final remark, we note that the results shown here are robust with respect to other definitions of complexity, such as the ones described in Marshall et al. (23) or Zamora-López et al. (24), as shown in fig. S5. Thus, independently of the definition of complexity used, the single-bond configuration is able to exhibit unique dynamical traits as compared to the other patterns.

Effective connectivity and integration-segregation balance

The analysis of network activity using CC and Θ provided an insight into the behavior of the networks as a whole, without providing information on the communication among neurons or the functional interplay between spatial modules. To render a complete picture, we investigated the functional organization of the networks to expose the flow of neuronal activity, identify neurons that have a central role, and quantify the integration-segregation capacity of the patterned cultures.

By functional organization, we refer to a set of graph-derived network features that include effective connectivity, functional communities, hubs, and the normalized global efficiency G*EFF (see Materials and Methods) (2528). The latter captures the capacity of the network to efficiently exchange information among all neurons and therefore provides a direct measure of its degree of integration. Figure 3A shows representative effective connectivity maps of the studied patterns. The nodes and connections are color coded according to their strength and weight, respectively. For the merged and triple-bond patterns, neurons form a highly integrated graph that shapes a unique functional community with G*EFF ≅ 1. Modules in the triple-bond pattern are so strongly bound that they operate as a single functional community. Some neurons play an important role in this cohesiveness in that they exhibit high centrality characteristics (see Materials and Methods) and are ascribed as provincial hubs (2, 25). By contrast, the single-bond pattern exhibits a more heterogeneous effective connectivity, with the spatial modules being connected to one another with fewer links. Two main functional communities are observed in the illustrative example shown in Fig. 3A. These communities are also composed of neurons from different modules, but their internal connectivity is very different, with the community highlighted in yellow exhibiting an abundance of provincial hubs. The two communities are linked together through a kinless node (see Materials and Methods) and weak links, shaping a system that, as a whole, exhibits a very low global efficiency of G*EFF ≅ 0.2, much smaller than the previous cases and that portrays its segregated behavior.

Fig. 3 Effective connectivity and global efficiency.

(A) Representative effective connectivity maps for patterns with gradually higher spatial modularity. The squared contours with a dotted outline indicate the location of the spatial modules. The merged and triple-bond patterns provide highly connected, highly integrated effective networks, with G*EFF ≅ 1 and only one functional community. Provincial hubs in the triple-bond pattern facilitate the permanent functional cohesion among modules. The single-bond network exhibits a much weaker connectivity and a much lower integration capacity, with G*EFF ≅ 0.18. Its effective network is characterized by two main communities (orange and yellow areas) that extend across different modules. Provincial hubs maintain cohesion within communities, which are, in turn, linked through a kinless node. (B) Raster plot of spontaneous activity and temporal evolution of G*EFF for the same single-bond network. The graphs below show illustrative effective networks at different stages of the recording to highlight the rich functionality of the system, which switches from high integration (stages 1 and 4) to high segregation (stages 2 and 3). The values below each network indicate the number of main communities. (C) Average G*EFF values for all studied networks (merged, 13 cultures; triple bond, 13; single bond, 18). The merged and triple-bond patterns exhibit G*EFF values significantly higher than those of the single-bond patterns. The merged cultures are always integrated, the triple-bond ones are mostly integrated with rare or sporadic segregation, and the single-bond ones switch between the two scenarios. Data are shown as means ± SEM. ***P < 0.001, **P < 0.01 (two-tailed Student’s t test).

Figure 3B reveals the concert between dynamical richness and shifting functional organization by depicting the temporal evolution of G*EFF for the same single-bond pattern. The effective networks continuously change in structure and number of communities and switch from strong functional segregation with G*EFF ≅ 0.2 to strong functional integration with G*EFF ≅ 0.6. Although coherent episodes are frequent, the operability of the single-bond pattern is markedly richer than the merged and triple-bond patterns, which reside in a permanent integrated state.

The lower G*EFF of the single-bond pattern was consistent over realizations. As shown in Fig. 3C, the average G*EFF value for the single-bond cultures was significantly smaller than the merged or triple-bond patterns. The fact that the values for the latter do not reach 1 indicates that, on average, the integration capacity of the studied networks is smaller than that in their equivalent random graphs. This is possibly due to the inherent spatial embedding of the cultures, which promotes local connectivity (20) and thus an overall weaker integration.

To gain a deeper insight into the impact of intermodular connectivity on functional organization, we also targeted the strength of excitatory connectivity in the single-bond pattern by application of the AMPA–glutamate receptor antagonist CNQX (6-cyano-7-nitroquinoxaline-2,3-dione). As shown in fig. S6, a network that was regularly firing in a synchronous manner was driven into a more segregated state by application of [CNQX] = 500 μM. The reduction of excitation enhanced the dynamical richness of the system and shaped distinct functional communities. We therefore hypothesize that small variations of the excitatory connectivity among modules are sufficient to control the degree of integration or segregation in the neuronal network.

Physical connectivity among modules

The above results highlight the pivotal role of intermodular coupling in network functionality. Strong coupling results in persistent integration with systematically coherent dynamics, while coupling just at the edge of intermodular disconnection leads to dynamics that are spatially and temporally more complex.

To shed light on the role of intermodular connectivity in facilitating coherent dynamics in the weakly coupled, single-bond networks, we analyzed the relationship between the thickness of the neurite bundle and the CC of the modules coupled by it (Fig. 4A). An unsupervised classification analysis based on the k-means algorithm revealed that, as a general trend, CC between two modules was higher when the bundle was thicker (group 1, blue) and that a population existed in which the correlation was low (approximately zero) independent of bundle thickness (group 2, gray). The plots classified as group 2 derived from samples in which either one module or both modules were not active during the 20-min recording. However, in another population, the modules coupled with a thin bundle fired coherently (group 3, red). Coherence in sparsely coupled modules was also apparent at the level of three- and four-module interactions (Fig. 4, B and C). The existence of such a thin-but-strong coupling is most likely the reason for the occurrence of integrated—yet not fully persistent—activity in the single-bond networks.

Fig. 4 Subquorum coherence in micropatterned cortical networks.

(A) Relationship between the mean CC between two neighboring modules and the width of the interconnecting neurite bundle in the single-bond patterns. Each dot represents a module pair and is colored according to the result of the k-means algorithm. (B) Three-body interaction of neuronal modules. The mean CC between neurons in modules α and γ, Embedded Image, is plotted against the bundle thickness of modules α-β and β-γ. The value of Embedded Image is indicated in color, and the size of the dots is changed in proportion to the value to aid visualization. (C) Average intermodular correlation [Embedded Image] plotted against average bundle thickness [Embedded Image]. Each dot represents an individual network and is colored according to the result of the k-means algorithm. The dashed line represents a linear fit against the blue plots.

Network mechanism of subquorum integration

To clarify the mechanisms that support functional integrability within spatially segregated modules, we constructed an in silico model for the modular neuronal network. In the model, the number of inputs required to excite a neuron, or quorum (19, 20), was approximately 10 (see Materials and Methods). Network models were generated by forming a designated number of interconnections between four random networks, each containing 25 nodes (Fig. 5A).

Fig. 5 Subquorum synchronization in in silico model networks.

(A) Schematic representation of a modular network model. Yellow circles indicate neurons, gray lines denote the connections, and arrows indicate the direction of the connection. (B) Occurrence of network bursts (events) is shown as the average number per 20 min. A total of 300 networks were sampled per condition. The probability of a network exhibiting more than one event is shown on the counter axis. Results obtained over a wider range are shown in the inset. The dotted line indicates data obtained from networks with completely randomized connectivity with identical average node degrees. (C) Mean CC between two modules versus the number of connections between the modules (solid line). Individual data are plotted in violet. To ease visualization, we randomly displace the plots horizontally by ±0.5. (D) Raster plots of two network bursts that occurred in the same network with 26 intermodular connections. The red bar indicates the timing of an action potential. The activity is initiated in module 3 (neurons 51 to 75) and propagates to modules 2 (neurons 26 to 50) and 4 (neurons 76 to 100) and lastly to module 1 (neurons 1 to 25). (E) The mechanism of subquorum synchronization. The feedback excitation from the hidden layer nodes (orange) to the input layer nodes (blue) enhances the sparse input to each module. (F) Random rewiring of the feedback (FB) excitation to noninput layer nodes diminishes the occurrence of network bursts. A rewiring ratio of 1 corresponds to a case in which no feedback exists from the hidden layer to the input layer. The original network (rewiring ratio, 0) is a network with 26 intermodular connections, whose activity is shown in (D). Error bars represent 95% confidence intervals.

As a general trend, integrability increased with the number of intermodular connections. The probability of observing more than one network burst was >97% when the number of intermodular connections was >76, i.e., about 20 connections per module pair (Fig. 5B). However, network bursts were also observed in more segregated networks. For example, 25% of the networks generated such activity even when the number of intermodular connections was 40, i.e., 10 connections per module pair (Fig. 5B). As the connections were directed, only five connections were added from one module to another, which is 50% of that expected from the quorum of individual neurons. In Fig. 5C, the results are reorganized to illustrate the dependence of intermodular correlation on the number of intermodular connections. Although the CC was positively correlated with the number of intermodular connections, highly correlated modules were occasionally observed with less than 10 connections.

We first hypothesized that this subquorum synchronization, i.e., synchronization between neurons that are coupled with fewer number of connections than the quorum, was caused by the dense local coupling within modules, which would increase the resting membrane potential of the neurons. However, the membrane potential distribution for the modular networks was similar to those of random networks with identical numbers of nodes and connections (fig. S7A). This points to the importance of the temporal structure, underpinned by network connectivity, in regulating correlated activity. As illustrated in Fig. 5E, when a module receives several coherent inputs from another module, the neurons that directly receive the input (“input layer” nodes) activate other neurons within the module (“hidden layer” nodes). A fraction of the hidden layer nodes then provides a recurrent feedback to the input layer. The dense intramodular connectivity in a modular network enriches these recurrent connections.

To analyze the role of these recurrent connections in subquorum integration, we sampled in the model a network containing 26 intermodular connections (6.5 connections per module pair) that generated network bursts. We observed that all network bursts were initiated in module 3, which then propagated to modules 2 and 4 and lastly to module 1 (Fig. 5D). Assuming this directionality, we identified all the recurrent connections from the hidden layer to the input layer in modules 1, 2, and 4 and randomly rewired them to noninput layer nodes within each module. The effect of this rewiring on the occurrence of network bursts is summarized in Fig. 5F. The integrability of the network decreases as the recurrent connections are rewired, and their total elimination resulted in a ~95% suppression of the network bursts. Increasing the synaptic strengths decreased the number of intermodular connections required for generating network bursts, but subquorum characteristics caused by the intramodular feedback excitation remained apparent (fig. S7B).


Micropatterned cultures provide a simple yet elegant platform to investigate the dynamic capacity of modular neuronal networks. Our design is based on spatially segregated units coupled to one another, a structure that approaches in vitro the organization of the mammalian brain (29). The uniqueness of our approach is its potential for systematically exploring the functional impact of intermodular connectivity. Our system is scalable, and we can therefore go beyond our current design and investigate more complex structures such as modular networks with weighted couplings (fig. S8A), directed connectivity (fig. S8B), or hierarchical organization (3, 24, 30, 31).

In our experiments, abundant connections across the spatial modules in the triple-bond configuration render a dynamically stable yet excessively rigid system in which the different modules fire in a coherent manner. The excess of connections drives the system to operate as a unique entity, i.e., making the spatial embedding of neurons into modules irrelevant. This observation may have strong implications in naturally formed modular networks such as the brain, since an excess of connectivity among brain areas may severely disrupt its inherent integration-segregation capacity and functional operation.

The reduction of intermodular connectivity in the single-bond configuration shapes a flexible and dynamically rich system, where fully synchronized regimes coexist with fractional ones with variations in participating neurons. We note that a rich distribution of coordinated activity sizes reveals the existence of neuronal interactions across broad spatiotemporal scales. Thus, the single-bond patterns may operate near a critical point, a state that optimizes information flow, maximizes communication, and facilitates swift response to perturbations (9, 23). Analogous to the brain, the flexibility of the single-bond pattern capacitates the system to plunge into more segregated or more integrated computations depending on the states of the neurons and connections. This flexibility may be a signature for the ability of the network to sustain damage. Teller et al. (27) observed that dynamically rich, self-organized modular networks formed by a few dozens of modules are much more resilient to the failure of connections as compared to homogeneous, nonmodular networks. The high dimensionality of the dynamics is also relevant from the viewpoint of cortical information processing.

Our observations are consistent with those of previous theoretical investigations (22, 32) and other experimental designs (15, 16, 27, 33). These studies, however, did not examine the network mechanisms behind such a flexibility and dynamical richness. The central finding of our study is that the rich dynamics and flexibility of the single-bond patterns are mediated by a relatively few connections among modules. This connectivity trait facilitates the emergence of functional communities of varying sizes and structures. These dynamically rich episodes are combined with regimes of full coherence that are rooted on subquorum firing promoted by connectivity-induced amplification. We conjecture that short-term plasticity and fluctuations are the channels for a flexible dynamical switching that renders the dynamical richness of the single-bond patterns. Furthermore, we hypothesize that coherence within such a system could be locally and transiently controlled through stimulation, as occurs in the brain upon sensory input.

The observation that weakly connected spatial modules can reach full coherence through subquorum amplification may have important implications in other fields such as social sciences and epidemics (1, 34, 35), where geographical segregation is often a natural property. In epidemics, for instance, weak intermodular connectivity could be ascribed to air transport routes, while subquorum amplification could be related to the increased probability for passengers to become infected in densely occupied planes. Thus, despite the low importance of interregional airplane mobility as compared to intraregional transportation, a coherent, interregional disease outbreak could occur when recurrent amplification is brought into play.


Fabrication of patterned substrates

Microcontact printing was used to pattern a mixture of poly-d-lysine and extracellular matrix (ECM) gel on an agarose-coated coverslip. PDMS stamps were fabricated according to the protocol described previously (13). The ratio of prepolymer and curing agent of PDMS (SYLGARD 184, Dow Corning) was optimized for each micropattern design, which was 10:1 for the no-bond and single-bond patterns and 7.5:1 for the triple-bond and merged patterns. The three-dimensional topography of the stamp surfaces was analyzed using confocal microscopy (VK-X260, KEYENCE).

To print proteins, a glass coverslip (12-545-84, Thermo Fisher Scientific) was first cleaned by sonication in 100% ethanol and subsequent treatment in air plasma for 60 s (PM-100, Yamato). Four dots of paraffin wax (P3558, Sigma-Aldrich) were then placed at the periphery of the coverslip (approximately 0.5 mm in height), and the surface was then coated with 0.2% agarose (A9918, Sigma-Aldrich). After drying the agarose overnight, the coverslips were sterilized by ultraviolet irradiation for 30 min, and protein ink [ECM gel (E1270, Sigma-Aldrich; 1:100 dilution) + poly-d-lysine (50 μg ml−1; P0899, Sigma-Aldrich)] was stamped using a single-axis micromanipulator. The coverslip was dried overnight in a fume hood and was subsequently immersed in the neuronal plating medium [minimum essential medium (MEM; 11095-080, Gibco) + 5% fetal bovine serum + 0.6% d-glucose] a day before cell plating.

Cell culture

Animal experiments were approved specifically for this study by the Center for Laboratory Animal Research, Tohoku University (approval number: 2014BeA-005-1). Timed-pregnant Sprague-Dawley rats were obtained from Charles River Laboratories, Japan, and were used immediately after receipt. Primary neurons were obtained from the cortices of embryonic day 18 pups, plated on the microfabricated coverslip, and cocultured with astrocyte feeder cells in N2 medium [MEM + N2 supplement + ovalbumin (0.5 mg ml−1) + 10 mM Hepes] (13, 14). Half the medium was changed at 4 and 8 DIV with Neurobasal medium [Neurobasal (21103-049, Gibco) + 2% B-27 supplement (17504-044, Gibco) + 1% GlutaMAX-I (35050-061, Gibco)].

Calcium imaging

At 10 DIV, networks bearing 70 to 130 neurons (20 to 30 neurons) were randomly selected from each coverslip for the single-bond, triple-bond, and merged (no-bond) patterns. Cultured neurons were loaded with a fluorescence calcium indicator Cal-520 AM (AAT Bioquest) by first rinsing the cells in Hepes-buffered saline (HBS) containing 128 mM NaCl, 4 mM KCl, 1 mM CaCl2, 1 mM MgCl2, 10 mM d-glucose, 10 mM Hepes, and 45 mM sucrose and subsequently incubating in HBS containing 2 μM Cal-520 AM and 0.01% Pluronic F-127 for 30 min at 37°C. The cells were then rinsed in fresh HBS and incubated for an additional 10 min to complete the deesterification of the loaded AM dyes. The samples were imaged on an inverted microscope (IX83, Olympus) equipped with a 20× objective lens (numerical aperture, 0.75), a light-emitting diode light source (Lambda HPX, Sutter Instrument), a scientific complementary metal-oxide semiconductor camera (Zyla 4.2, Andor), and an incubation chamber (Tokai Hit). All recordings were performed at 37°C. Three networks were selected from a coverslip, and for each network, fluorescence imaging was performed for 20 min at 10 frames/s on the Solis software (Andor). The thickness of the neurite bundle that interconnects the modules was estimated from phase-contrast micrographs based on the brightness profile.

Neuronal activity analysis

The recordings were analyzed offline using either ImageJ (National Institutes of Health) or NETCAL (36). For analyses on ImageJ, 60 neurons (15 per module) were randomly selected, and their somas were manually defined as regions of interest (ROIs). Relative fluorescence intensity of cell i, (Fi(t) − F0)/F0 = fi(t), with Fi and F0 being the fluorescence intensity and the background fluorescence, respectively, was then analyzed as described previously (14). Alternatively, the recordings were analyzed using the NETCAL software. For this method, the image of the average fluorescence of the recording was first created and the background (areas with no change in fluorescence) was removed. Neuronal somas were next automatically detected and ascribed as ROIs. An average of about 90 neurons was identified per network, which were uniformly distributed across modules. The average fluorescence in each ROI along the recording was then extracted, and fi(t) was finally obtained as described above.

Spike event detection and network bursts

The time series of fi(t) were analyzed with NETCAL (36) to infer the timing of neuronal activations using the Schmitt trigger method. This method scans the fluorescence traces for events that first pass a high threshold and then remain elevated above a second lower threshold for at least a certain minimum duration. Shorter events are discarded. In our analysis, we used +3 SDs of the mean of the baseline noise as the high threshold, +1.5 SD as the low threshold, and 200 ms as the minimum event length. Schmitt spike inference was contrasted with the time-derivative thresholding method (14). The overall results were consistent across methodologies.

Network bursts provided a quick insight into the collective dynamics in cultures and were defined as those activity episodes in which more than 75% of the neurons fired simultaneously for more than 500 ms. Experiments with less than 10 network bursts were discarded in the analysis of dynamical richness and effective connectivity.

Correlation coefficients

The functional correlation between neurons i and j was quantified using the Pearson CC, rijEmbedded Image(2)where xi(t) is either fi or the inferred spike trains Si, and Embedded Image is the corresponding time-averaged values. Si(t) are binary signals, representing 1 for the presence of a spike at time t and 0 otherwise. For each network, the mean correlation was evaluated as Embedded Image, where N is the number of neurons, and was used as an index of the degree of global correlation. Mean neural correlation within a single module and between two separated modules was evaluated as the intramodular and intermodular correlations, respectively, by selecting corresponding i-j pairs.

Jensen-Shannon divergence

This was used to compute the similarity between the distributions of IBIs. For two distributions P and Q, the Jensen-Shannon divergence DJS is given byEmbedded Image(3)where DKL is the Kullback-Leibler divergence Embedded Image and M = (P + Q)/2.

Global network activation

Raster plots of neuronal activity were scanned using a Δt = 1 s window to detect episodes of network coherent firing. The fraction of active neurons in each episode provided the global network activity at a given time of the recording, Γt. Coherent episodes comprising less than 3% of the network were discarded to filter out sporadic, random activity. In addition, the details of the onset times for each participating neuron within the episode were termed “sequence” and were later used to compute the effective connectivity.

Effective connectivity

Directed and weighted graphs were constructed from the inferred neuronal spike trains. Following Teller et al. (27, 28), we used a formalism grounded on time delays to determine the strength of functional couplings among neurons. For each detected sequence of coherent activity, the weight of the functional links among all pairs of firing neurons was established as a decaying function of their time delay in activation. The analysis was then extended to all observed sequences, and the final value of the functional connection wij between neurons (i, j) was established as the sum of all computed weights. The significance of functional links was assessed through a random model wSij in which the number of firings per neuron was preserved, and a z score was introduced in the form ofEmbedded Image(4)where 〈wSij is the mean weight for 500 surrogates of the random model and Embedded Image is the corresponding SD. Values with Wij < 0 correspond to links that are less connected than those in the surrogates and were set to 0. The final matrix W = {Wij} (weighed and directed) provided the normalized weights between all the functionally connected neurons.

We verified (fig. S9) that the effective connectivity matrices derived using time delays were similar to those obtained with other approaches such as transfer entropy (6, 37). Both approaches shared 80% of the strongest links, procured identical communities, and exhibited very similar measures. In addition, we also verified that important network measures such as the global efficiency were robust to more restrictive thresholds in the z score (fig. S10).

Community detection

We first performed hierarchical clustering analysis using the Jaccard metrics as the distance (28) to obtain the pairwise similarity matrix, which was then normalized using a z score with respect to a null model that consisted of a set of 500 surrogates. The average linkage criterion was applied in the final matrix to establish the arrangement of the grouping of neurons shaped as a hierarchical tree (Statistics and Machine Learning Toolbox, MATLAB). To mathematically identify the communities, a threshold had to be established in the hierarchical tree (28). This threshold was set by taking advantage of the variation of information measure (38) to provide the final set of participating neurons in each community.

Centrality measures and hubs

We used a combination of centrality measures to identify hubs in the networks. Following Karrer et al. (38), we ascribed as hubs those nodes that excel (above population average) in at least two of the following categories: nodal strength, nodal betweenness, and local efficiency. Hubs are crucial for functional integration and facilitate the flow of information in a network. For networks that present more than one community, hubs were further classified to pinpoint their role in integrating the communities. Classification was established through the participation coefficient (39, 40), a measure that evaluates the distribution of the functional links of a node across all communities in a network. Hubs that present, above average, a high participation coefficient are referred to as connector hubs and promote communication among communities. On the contrary, hubs that show a low participation coefficient are referred to as provincial hubs. Their role is to facilitate communication within a community. Last, those nodes that present a high participation coefficient but were not initially classified as hubs are referred to as kinless nodes. These nodes typically do not belong to a particular community but facilitate intercommunity integration.

Normalized global efficiency

This provides a measure for the integration capacity of a network. It is computed by first evaluating the matrix of connection lengths, L = {Lij}, which are inversely related to the link weights, Wij. The higher the weight in the connection among two neurons, the shorter is their connection length. This matrix was then analyzed to compute the shortest path length dij between neurons i and j and, lastly, the global efficiency, GEFF, as the average inverse of the shortest path length (25, 26)Embedded Image(5)Embedded Image(6)where Luv is the length between nodes u and v, gij is the shortest geodesic path between nodes i and j calculated using the Dijkstra algorithm, and N is the total number of neurons. We note that dij = ∞ for disconnected pairs (i, j). To compare the global efficiency among different experiments, we considered its normalized form G*EFF = GEFF/GRGEFF, where GRGEFF is the global efficiency for the random graph equivalent (null model) of the studied network, computed from wS = {〈wSij}, as introduced before. We verified that the differences in global efficiency among the three configurations (merged, single bond, and triple bond) were robust to thresholds in the significance level of the procured effective networks (fig. S10).

Temporal evolution of normalized global efficiency

The 20-min recordings were divided into 10 sections of 2 min duration each. The network activity within each section was then analyzed as an independent dataset, and the effective network W and normalized global efficiency G*EFF were computed accordingly. One must observe that G*EFF analyzed along the different sections may take values higher than the corresponding value for the full 20-min recording. This is caused by the large statistics of firing events in the full recording, which filters out weak connections or coactivations that occur rarely.


Cortical neurons were fixed with a 4% paraformaldehyde (Sigma-Aldrich)/4% sucrose (Wako) solution, permeabilized with 0.25% Triton X-100 (Nacalai), and blocked in a 0.5% solution of fish skin gelatin (Sigma-Aldrich). Axons, somatodendritic regions, and presynaptic sites were stained with SMI-312 [837904, BioLegend; mouse immunoglobulin G1(IgG1)/IgM], MAP2 (ab5392, Abcam; chicken IgG), and Syn1 (106 001, Synaptic Systems; mouse IgG1) antibodies, respectively. The secondary antibodies were Alexa Fluor 488–labeled goat anti-mouse IgG1 (A21121, Molecular Probes) and Alexa Fluor 568–labeled goat anti-chicken IgG (A11041, Molecular Probes). Fluorescence images were obtained using either an epifluorescence microscope (IX83, Olympus) or a confocal microscope (SP8, Leica).

In silico modeling

The model network consisted of four modules each containing 25 leaky integrate-and-fire neurons, which are described in (14). Within each module, the neurons were randomly interconnected with directed connections at an average degree of 8, and the network modularity was tuned by adding interconnections between the four modules. To best model the experiments, intermodular connections were only added to neighboring modules, and connections between diagonal modules were forbidden.

Each simulation was run for 1200 s with a time step of 0.1 ms. Only excitatory synapses (synaptic conductance, 5 nS) were considered, and each contributed to a synaptic potential of approximately 2 mV. The resting and threshold potentials were set at −74 and −54 mV, respectively. Spontaneous activity was triggered by delivering Poisson input to each neuron at 0.5 Hz. The simulated activity of neurons was binned at 10 ms, and a network burst was detected when more than 80% of the neurons were continuously activated for more than 100 ms.


Supplementary material for this article is available at

Fig. S1. Axon outgrowth and synapse formation of micropatterned cortical neurons.

Fig. S2. Effect of intermodular connections on the network structure.

Fig. S3. Average number of cells in the networks.

Fig. S4. Analysis of network bursts.

Fig. S5. Other complexity measures.

Fig. S6. Role of excitatory connectivity on dynamical richness and functional organization.

Fig. S7. Membrane potentials and synaptic conductance in model networks.

Fig. S8. Impact of nonuniform intermodular connectivity in simulations.

Fig. S9. Robustness of effective connectivity inference.

Fig. S10. Dependence of global efficiency GEFF on the amount of the retained significant links.

Movie S1. Spontaneous activity recording of a merged network.

Movie S2. Spontaneous activity recording of a triple-bond network.

Movie S3. Spontaneous activity recording of a single-bond network.

Reference (41)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank S. Kono, K. Ishihara, S. Fujimori, and R. Matsumura for technical assistance. Funding: This study was supported by the Cooperative Research Project Program of the Research Institute of Electrical Communication at Tohoku University, JSPS KAKENHI (nos. 15K17449, 26390035, and 18H03325), and the JST CREST Program (JPMJCR14F3). J.S. and S.T. acknowledge financial support from the Spanish Ministerio de Economia y Competitividad through project no. FIS2016-78507-C2-2-P and from the Generalitat de Catalunya through grant no. 2014-SGR-878. This research is part of MESO-BRAIN. The MESO-BRAIN Project has received funding from the European Union’s Horizon 2020 Research and Innovation Programme under grant agreement no. 713140. Author contributions: H.Y. and J.S. conceived and designed the project with inputs from S.S., S.K., T.T., M.N., and A.H.-I. H.Y., K.I., and T.H. performed the experiments. S.M., H.A., and S.K. coded and performed the simulations. H.Y., S.M., K.I., and S.T. analyzed the results. H.Y. and J.S. wrote the manuscript. All authors reviewed and edited the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article