Figures
Abstract
The analysis of the activity of neuronal cultures is considered to be a good proxy of the functional connectivity of in vivo neuronal tissues. Thus, the functional complex network inferred from activity patterns is a promising way to unravel the interplay between structure and functionality of neuronal systems. Here, we monitor the spontaneous self-sustained dynamics in neuronal cultures formed by interconnected aggregates of neurons (clusters). Dynamics is characterized by the fast activation of groups of clusters in sequences termed bursts. The analysis of the time delays between clusters' activations within the bursts allows the reconstruction of the directed functional connectivity of the network. We propose a method to statistically infer this connectivity and analyze the resulting properties of the associated complex networks. Surprisingly enough, in contrast to what has been reported for many biological networks, the clustered neuronal cultures present assortative mixing connectivity values, meaning that there is a preference for clusters to link to other clusters that share similar functional connectivity, as well as a rich-club core, which shapes a ‘connectivity backbone’ in the network. These results point out that the grouping of neurons and the assortative connectivity between clusters are intrinsic survival mechanisms of the culture.
Author Summary
The architecture of neuronal cultures is the result of an intricate self-organization process that balances structural and dynamical demands. We observe that when the motility of neurons is allowed, these neurons organize into compact clusters. These neuronal assemblies have an intrinsic synchronous activity that makes the whole cluster firing at unison. Clusters connect to their neighbors to form a network with rich spontaneous dynamics. This dynamics ultimately shapes a directed functional network whose properties are investigated using network descriptors. We find that the networks are formed such that preference in connectivity between clusters is based on the similarity between their activity, a property that is called assortative mixing in networks' language. This particular choice of connectivity correlations must be rooted to basic survival mechanisms for the neurons constituting the culture.
Citation: Teller S, Granell C, De Domenico M, Soriano J, Gómez S, Arenas A (2014) Emergence of Assortative Mixing between Clusters of Cultured Neurons. PLoS Comput Biol 10(9): e1003796. https://doi.org/10.1371/journal.pcbi.1003796
Editor: Boris S. Gutkin, École Normale Supérieure, College de France, CNRS, France
Received: January 30, 2014; Accepted: July 6, 2014; Published: September 4, 2014
Copyright: © 2014 Teller et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Funding: Research was supported by the Ministerio de Ciencia e Innovaci\'on (Spain, www.idi.mineco.gob.es) under projects FIS2010-21924-C02-02, FIS2011-28820-C02-01, FIS2012-38266-C02-01, and FIS2010-09832-E. We also acknowledge the Generalitat de Catalunya (www10.gencat.cat/agaur_web) under projects 2009-SGR-00014 and 2009-SGR-838, and the EU FET (MULTIPLEX 317532). A.A. also acknowledges partial financial support from the ICREA Academia and the James S.\ McDonnell Foundation. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
The theory of complex networks [1]–[5] has proven to be a useful framework for the study of the interplay between structure and functionality in social, technological, and biological systems. A complex network is no more than a specific representation of the interactions between the elements of the system in terms of nodes (elements) and links (interactions) in a graph. The analysis of such resulting abstraction of the system, the network, provides clues about regularities that can be connected with certain functionalities, or even be related to organization mechanisms that help to understand the rules behind the system's complexity. Particularly, in biological systems, the characterization of the emergent self-organization of their components is of utmost importance to comprehend the mechanisms of life [6]–[8].
One of the major challenges in biology and neuroscience is the ultimate understanding of the structure and function of neuronal systems, in particular the human brain, whose representation in terms of complex networks is especially appealing [9], [10]. In this case, structural connectivity corresponds to the anatomical description of brain circuits whereas the functional connectivity is related to the statistical dependence between neuronal activity.
Network theory and its mathematical framework have provided, through the analysis of the distribution of links, statistical measures that highlight key topological features of the network under study. These measures have facilitated the comprehension of processes as complex as brain development [11], learning [12] and dysfunction [13], [14]. Particularly, these measures have unfolded new relationships between brain dynamics and functionality. For instance, synchronization between neuronal assemblies in the developing hippocampus has been ascribed to the existence of super-connected nodes in a scale-free topology [15]; efficient information transfer has been associated to circuits with small-world features [16], such as in the the nematode worm C. elegans [17] or the brain cortex [18], [19]; and the coexistence of both segregated and integrated activity in the brain has been hypothesized to arise from a modular circuit architecture [20]–[22].
A network measure that has recently caught substantial attention is the assortativity coefficient, which quantifies the preference of a node to attach to another one with similar (assortative mixing) or dissimilar (disassortative mixing) number of connections [23], [24]. Assortative networks have been observed in both structural [20] and functional [25] human brain networks. It has been proposed that assortative networks exhibit a modular organization [26], display an efficient dynamics that is stable to noise [27], and manifest resilience to node deletion (either random or targeted) [23], [28]. Resilience is ascribed to the preferred interconnectivity of high-degree nodes, which shape a ‘connectivity backbone’ [29] that preserves network integrity. The existence of such a tightly interconnected community is generally known as the ‘rich-club’ phenomenon [19], [22], [30]. On the other hand, disassortative networks, such as the ones identified in the yeast's protein-protein interaction and the neuronal network of C. elegans [23], are more vulnerable to targeted attacks. However, in these disassortative networks, the tendency of high degree nodes to connect with low degree ones results in a star-like topology that favors information processing across the network.
The assortativty coefficient is usually calculated through the Pearson correlation coefficient between the unweighted degrees of each link in the network [23]. To account for effects associated to large networks, the Spearman assortativity measure was introduced [31] and, later, weighted assortativity measures were proposed to include the weight in degree-degree dependencies [32].
To better understand the importance of these network measures in describing neuronal networks, in vitro preparations in the form of neuronal cultures have been introduced given their accessibility and easy manipulation [6], [33]. Two major types of cultured neuronal networks are of particular interest. In a first type, neurons are plated on a substrate that contains a layer of adhesive proteins. Neurons firmly adhere to the substrate, leading to cultures with a homogeneous distribution of neurons [34]–[37]. In a second type, neurons are plated without any facilitation for adhesion. Neurons then spontaneously group into small, compact assemblies termed clusters that connect to one another [38]–[41].
The formation of a clustered architecture from an initially isotropic configuration is an intriguing self-organization process [40], [41]. This feature has made clustered networks attractive platforms to study the development of neuronal circuits as well as the interplay between structural and functional connectivity at intermediate, mesoscopic scales [39], [42]–[45]. Moreover, the existence of a two-level network, one within a cluster and another between clusters, has made clustered cultures appealing to study dynamical and topological features of hierarchical [40], [41], [45], [46] as well as modular networks [46]–[49].
In this work we use spontaneous activity measurements in clustered neuronal cultures to render the corresponding directed functional networks and study their topological properties. We introduce a novel theoretical framework that uses the propagation of activity between clusters as a measure of “causality”, although strictly speaking we should refer to as a sequence of delayed activations, giving rise to functional connections that are both directed and weighted. Based on this weighted nature of the network, we propose a new measure of assortativity that explicitly incorporates the weight of the links. We observed that all the studied functional networks derived from clustered cultures show a strong, positive assortative mixing that is maintained along different stages of development. On the contrary, homogeneous cultures tend to be weakly assortative, or neutral. Finally, in combination with experiments that measure the robustness of network activity to circuitry deterioration, we show that the strongly assortative, clustered networks are more resistant to damage compared to the weakly assortative, homogeneous ones. Our work provides a prominent example of the existence of assortativity in biological networks, and illustrates the utility of clustered neuronal cultures to investigate topological traits and the emergence of complex phenomena, such as self-organization and resilience, in living neuronal networks.
Results
Experiments
We used rat cortical neurons in all the experiments. As described in Methods, neurons were dissociated and seeded homogeneously on a glass substrate. Cultures were limited to circular areas in diameter for better control and full monitoring of network behavior. The lack of adhesive proteins in the substrate rapidly favored cell-to-cell attachment and aggregation, giving rise to clustered cultures that evolved quickly (Figure 1A). By day in vitro (DIV)
, cultures contained dozens of small aggregates, which coalesced and grew in size as the culture matured. Connections between clusters as well as initial traces of spontaneous activity were observed as early as DIV
. Cultures comprised of
interconnected clusters by DIV
, and were sufficiently stable and rich in activity for measurements. Although the strength of the connections in the network and its dynamics evolved further, we observed that the size and position of the clusters remained stable. We therefore measured dynamics already at DIV
, and studied cultures up to DIV
.
A Bright field image of a network at day in vitro . Dark circular objects are aggregates of neurons (clusters), and filaments are visible physical connections between them. B Corresponding fluorescence image, integrated over 50 frames (
). Bright clusters at the top-left corner are active ones. C Spontaneous activity in the network (see Movie S1 for the actual recording). The top plot shows the average fluorescence signal of the clustered network shown in B, and along
of recording. The sharp peaks in fluorescence correspond to the fast sequential ignition of a group of clusters (burst). The bottom raster plot shows the clusters that ignite along the recording. The yellow bars relate a fluorescence peak with the ignition of a group of clusters, and highlights the tendency for the clusters to activate in specific groups. D Example of a particular ignition sequence in a region of the network containing
clusters. From left to right, the progress of cluster's activation is revealed by the increase in fluorescence signal of the downstream connected clusters. E Order of activation (black arrows) according to the analysis of the fluorescence signal. The clusters marked in yellow are those that fire simultaneously within experimental resolution. The ones in grey are clusters that do not participate in the firing sequence, and either fire independently or remain silent. F Detail of the fluorescence traces for the
participating clusters along two consecutive bursts, illustrating the accuracy in resolving the time delay in the activation of the clusters. The two bursts contain the same clusters, but the activation sequences are slightly different. Blue dots mark the ignition time, and yellow dots signal the clusters that fired simultaneously. The bottom orange boxes depict the final activation sequences of each burst. In the construction of the directed functional network, the influence of a cluster on another is conditioned by the time span between their activations. Close activations result in strong couplings (green arrows); far activations in weak ones (blue). Any two clusters whose activations are above
ms are considered functionally uncoupled (red).
The example shown in Figure 1A corresponds to a culture at DIV . Clusters appear as circular objects with an average diameter of
and a typical separation of
. Connections between clusters are visible as straight filaments that contain several axons.
We monitored spontaneous activity in the clustered network through fluorescence calcium imaging (Figure 1B). Fluorescence images of the clustered network were acquired at a rate of frames per second, and with an image size and resolution that allowed the monitoring of all the clusters in the network with sufficient image quality (see Movie S1). Activity was recorded for typically
hour, which provided sufficient statistics in firing events while minimizing culture degradation due to photo-damage.
The analysis of the images at the end of the measurement provided the variations in fluorescence intensity for each cluster and the corresponding onset times of firing (see Methods). As shown in the top panel of Figure 1C, the average fluorescence signal of the network is characterized by peaks of intense cluster activity combined with silent intervals. The accompanying raster plot reveals that this activity actually corresponds to the collective ignition of a small group of clusters, which fire sequentially in a short time window on the order of few hundred milliseconds. We denote by bursts these fast sequences of clusters' activations. Young cultures (DIV ) exhibited an activity of about
, while maturer cultures (DIV
) displayed about
(see also Table 1).
We observed that the time spanned between two consecutively firing clusters typically ranged between and
(see Methods), as also observed by others [47],[48]. These times are fairly large compared to the eventual scale of signal integration-propagation between single neurons (
), and is related to the large time scales associated to integration of the intra-clusters information. No consecutive activations were observed above
, signaling the termination of a burst. We therefore use this value of
as a cut-off to separate a given burst form the preceding one. Then, two clusters that fired above 200 ms cannot be influenced by one another and therefore are not causally connected.
Bursts occurred every 30 s on average for the experiment shown in Figure 1C and, as illustrated by the yellow bands in this figure, each burst typically encompassed a subset of clusters rather than the entire network. In general, however, the number of participating clusters within a burst depended on the details of the culture. Although in a typical experiment the collective firing comprised between and
clusters (see Movie S1), in some experiments the entire cluster population lighted up in a single bursting episode.
The analysis of the onset times of firing provides the cluster's activation sequence within each burst. As an example, Figure 1D depicts a highly active region of the network shown in Figure 1B. This region contains clusters, and
of them form a subset that regularly fires together. The series of frames show the progress in clusters' activation, revealed by the changes in fluorescence. Activity starts at the top-left cluster and progresses downwards. The time-line of sequence activation after image analysis is shown in Figure 1E, and the actual fluorescence traces are shown in Figure 1F. With our
temporal resolution we could resolve well the propagation of activity from a cluster to its neighboring ones (black arrows in Figure 1E). However, and for about
of the cases, the time delay between clusters' activation was either too short for detection or activation occurred simultaneously. The clusters associated to these ‘simultaneous’ events are marked in yellow in Figure 1E, and their inter-relation was treated as a bi-directional link (yellow arrow), since no causality can be inferred.
A typical recording provided on the order of bursting episodes. Some of them included the same group of clusters, although the precise sequence of activation could vary. An illustrative example is shown in Figure 1F, which depicts the fluorescence traces of
clusters along two consecutive bursts. The first sequence corresponds to the sketch of Figure 1E. The orange box at the bottom of the plot indicates the relative activation time of each cluster within the window, with two clusters treated as simultaneous. To introduce the construction of the directed functional network that is described later, we note that, intuitively, the firing of cluster #9 is most likely caused by #8 and therefore both clusters are (functionally) strongly coupled. At the other extreme, cluster #1 most likely did not trigger #9, and therefore their mutual coupling is very weak. For the second burst, we note that the activation sequence is very similar, but the relative delay times differ, therefore modifying the cluster's coupling strengths. Indeed, cluster #1 and #9 are now functionally disconnected given their long temporal separation.
We carried out measurements in different clustered networks, and labeled them with capital letters as networks ‘A’
‘O’. In order to compare their properties with the ones from cultures with a distinct structure, we applied the same measuring protocols and data analysis to
cultures characterized with a homogeneous distribution of neurons (see Methods and Figure S1), and labeled them as networks ‘P’
‘U’.
Construction of the directed functional networks
The above sequences of clusters' activations, extended to all the clusters and bursting episodes of the monitored culture, convey information on the degree of causal influence between any pair of clusters in the network. For instance, cluster #5 in Figures 1E-F can fire because of the first order influence of clusters #3 and #4, but also because of the second and third order influences of clusters #2 and #1, respectively. Hence, a realistic functional network construction should take into account these possible influences from the upstream connected clusters to build a network whose links are not only directed, but also weighted by the time delays in activation. This weighted treatment of the interaction between clusters is the major novelty of our work and the backbone of our model.
More formally, the interaction between any two clusters follows the principle of causality, i.e. the firing of cluster immediately after cluster
eventually implies that cluster
has induced the activity of
at that particular time. The likelihood of this relation between clusters is weighted according to its frequency along the full observational time, allowing to a statistical validation. Indeed, cluster
could induce the activity of various clusters, if all of them activate in a physically plausible short time window after cluster
. Such a construction is illustrated in Figure 2.
A Schematic representation of the experimental data, with firings of four different clusters (or neurons in the case of homogeneous cultures). B Stages of the method to construct the directed functional network. (1) The first step consists in detecting the different bursts in the whole sequence. Firings that are separated by more than
for clustered cultures (
for homogeneous ones) are not considered part of the same burst. (2) Calculation of the time lags
between consecutive firings inside the bursts, for the whole sequence. The frequency distribution of time lags is next fitted to a Gaussian distribution
, finally providing the variance
that will be specific for each culture. (3) Weighting procedure example for the first burst. Cluster #1 can activate #2 and #3. The weight of the links 1
2 and 1
3 depends on the time differences between clusters' activations, and is given by the function
determined by the previous fitting. Hence, weight
, and
. Cluster #3 can be activated as well by cluster #2, with
. (4) Schematic representation of the resulting directed functional connectivity network. The width of the connections is proportional to the weight of the links for clarity.
To construct the directed functional networks for each studied culture we proceed as follows. First of all, we divide the entire firing sequence into the bursts of clusters' activity (Figure 2A) using the cut-off of introduced in the previous section. Once the bursts have been detected, we compute the frequency distribution of time lags between pairs of consecutive firings (Figure 2B). This frequency distribution informs about the characteristic times expected between two consecutive firings within the same burst, and hence it is a good proxy of the causal influence of a cluster on another. We will use this information to weight the causal influence of firing propagation. The frequency distribution
presents a good fit to a universal Gaussian decay (
) in all the analyzed cultures, although the variance
is specific for each culture. We indeed observed (Figure S2) that
decreases with the culture age in vitro (correlation coefficient
, significance
), and increases with the number of clusters present in the network (
,
).
The last step in the construction of the directed functional networks consists in linking the interactions within each burst, and weighting them according to the previous frequency distribution (Figure 2B). The rationale behind this process is as follows: we hypothesize that every cluster influences other clusters (posterior in time) within a burst and, the larger the time after a cluster has fired the lower the influence we expect in the activation of another cluster (simply because the signal fades out). Then, the weighting of the interaction by the expected frequency observed in the distribution conveys the functional influence between clusters. The weights are reinforced every time the same pair of clusters' sequence is observed. After processing the full sequence we obtain a peer-to-peer activation map that is our proxy of the functional network.
We proceeded identically to construct the directed functional networks for homogeneous cultures (see Methods), with the only difference that the cut-off time corresponds to . We tested for both clustered and homogeneous cultures that the obtained functional networks were stable upon variations of the cut-off times (see Figure S3 and Discussion).
Analysis of the functional networks
We computed the functional networks of the (‘A’ to ‘O’) realizations of clustered cultures, as well as the
(‘P’ to ‘U’) homogeneous ones, and analyzed some major topological traits. Firstly, for each culture we obtained the number of nodes, the number of edges, the average degree of the networks, and its average strength (see Methods). The investigated networks and their topological measures are summarized in Table 1. Although young cultures display a richer activity, in general all networks present a similar number of nodes and a comparable functional connectivity, which is described by the number of edges, the average degree and the average strength.
Representative examples of the investigated functional networks for the clustered configuration are shown in Figure 3 (see Figure S1 for an example of the homogeneous ones). The position of the nodes and their size are the same as the actual clusters for easier comparison. Edges in the directed network are both color and thickness coded to highlight their importance, with darker colors corresponding to the highest weights. This representation reveals those pairs of clusters that maintain a persistent causality relationship over time. Nodes are also color coded according to their strength, i.e. the total weight of the in- and out-edges.
Top: Bright field images of representative neuronal cultures at different days in vitro. Bottom: Corresponding functional networks obtained from the directed and weighed construction described in Figure 2. From left to right, the pictures correspond to the cultures labeled D, H and O in Table 1. Only active clusters are used in the construction of the functional network. The size of the nodes is similar to the ones observed in the cultures, and facilitates the comparison of the functional network with the real culture. In the functional networks, the edges are both color and thickness coded according to their weight, while the nodes are only color coded according to their strength. The darker the color, the higher the value.
The functional networks exhibit some interesting features. First, there are groups of nodes that form tightly connected communities. These communities actually reflect the most frequent bursting sequences. Second, nodes preferentially connect to neighboring ones with some long-range connectivity, and often following paths that are not the major physical connections. This indicates that the structural connectivity of the network cannot be assessed from just an examination of the most perceivable processes. And third, as shown in Figure S4A, we observed that there is no correlation between the width of the physical connections and their weight (,
), or the size of the nodes and their strength (
,
, Figure S4B), and indicates that the dynamical traits of the network cannot be inferred from its physical configuration, stressing the importance of the functional study.
We also observed that the size of the clusters did not correlate with their average activity (,
, Figure S4C), i.e. small and big clusters displayed similar firing frequencies, and of
on average. However, since some clusters are initiators of activity and others just followers, we also computed the relative contribution of a given cluster size to initiate activity in the network. We found no significant correlation between initiation and cluster size (
,
, Figure S4D). These results strengthen the conclusion that one cannot predict the clusters that will initiate activity, or the most persistent sequences, by just a visual inspection of cluster sizes and their distribution over the network.
These analyses are important in the context of the work by Shein-Idelson and coworkers [41], who studied the dynamics of isolated clusters similar to ours, and observed that their firing rate increased from to
as the clusters' radii escalated from
to
. This remarkable difference in the dynamics between ‘isolated’ and ‘networked’ clusters reflects the dominant role of the network circuitry in shaping its dynamics.
Finally, to crosscheck that the results found for the functional networks presented here are robust to the inference method, we have also performed a classical mutual information analysis to construct functional networks for the same cultures (see Methods). The results obtained with the mutual information analysis are totally in agreement with the constructed functional networks using time delays.
Assortativity and rich-club properties
We determined the values of the weighted formulation of assortativity, both for the Pearson and Spearman
correlations, with values in
(see Methods for the generalization of assortativity to directed weighted networks). Positive values of the weighted assortativity indicate that nodes with similar strength tend to connect to one another, while negative values mean the preferred interconnectivity of nodes with different strength. In Table 1 we can observe that all clustered networks (labeled ‘A’-‘O’) exhibit a positive weighted assortativity, in the range
for the Pearson construction and
for the Spearman one. Although the values fluctuate across different cultures, the two assortativity measures provide the same value within statistical error, and reflect that network size corrections provided by the Spearman's treatment have little influence in strongly assortative networks.
To assess the importance of the measured assortativity values, we have also computed the weighted rich-club [50]. The rich-club phenomenon refers to the tendency of nodes with high degree to form tightly interconnected communities, compared to the connections that these nodes would have in a null model that preserves the node's degree but otherwise is totally random. Given the positive assortativity found, we analyzed whether this finding is also reinforced by the existence of rich-club structures.
The weighted formulation for the rich-club takes into account the node's strength instead of the degree, and is particularly useful in situations in which the weights of the links can not be overlooked [51]. The evaluation of the rich-club is performed by computing the ratio between the connectivity strength of highly connected nodes and its randomized counterpart, and for gradually higher values of the strength threshold
. The detailed calculation is described in the Methods section, and the results of the analysis for representative networks is shown in Figure S5. Ratios larger than
indicate that higher strength nodes are more interconnected to each other than what one would expect in a random configuration. On the contrary, a ratio less than
reveals an opposite organizing principle that leads to a lack of interconnectivity among high-degree nodes. After the calculation of the ratios for all the studied clustered networks, we found a positive tendency towards the creation of rich-clubs in all of them (Figure S5), which is in good agreement with the observed values of assortativity.
The above network measures were also analyzed in experiments with a homogeneous distribution of neurons (labeled ‘P’‘U’). The results are summarized in Table 1. Interestingly, the assortativity values are much lower (by an order of magnitude on average) than the ones for clustered cultures, in the range
for Pearson's and
for Spearman's. Accordingly, the rich-club phenomenon for the homogeneous cultures vanishes (Figure S5).
Network resilience
Several studies highlight the importance of assortative features for network resilience to damage. Given the strong assortativity of our clustered cultures, we carried out a new set of experiments to investigate the concurrent presence of resilient traits. As described in Methods, we considered two major ‘damaging’ actions to the network. In a first one, we gradually weakened the excitatory network connectivity by means of the AMPA-glutamate antagonist CNQX, and measured the decay in spontaneous activity as connectivity failed. In a second one, we continuously exposed a culture to strong fluorescence light, therefore inducing photo-damage to the neurons. This action resulted in random neuronal death across the network and hence a progressive failure of its spontaneous dynamics. The rate of activity decay upon radiation damage provided an estimation of the resistance of the network to node deletion. These investigations were carried out at the same time in clustered cultures (strongly assortative) and in homogeneous ones (weakly assortative or neutral). Their comparison provided a first reference to relate assortativity, network topology and resistance to damage.
Figure 4A shows the results for the application of CNQX to clustered cultures. We first monitored each cluster individually in the unperturbed case, and measured its average firing activity along
. We then applied a given drug concentration, measured the firing activity
for another
, and computed the relative changes in activity respect to the unperturbed case, as
. The protocol was repeated until activity ceased. Two illustrative examples of the action of CNQX on network activity are provided in Figure 4. In a clustered cultured and for weak CNQX applications (
) the activity in some clusters increases, while in some other decreases, and on average the network firing rate remains stable (
). As [CNQX] is increased to
, we observe that most of the clusters have reduced their activity, although there are still some that maintain a high activity or even increase it. This different behavior from cluster to cluster suggests that clustered networks are highly flexible, and that they may have mechanisms to preserve activity even with strong weakening of the connectivity. Conversely, homogeneous cultures (Figure 4B) lose activity in a more regular and faster way. These networks are characterized by a highly coherent dynamics [36], [37], and therefore all neurons in the network reduce activity similarly as CNQX is applied. Interestingly, for
the shown homogeneous culture has almost completely silenced (
), while the clustered culture is still highly active. We repeated this study on
different realizations of each culture type and observed that, on average, the critical concentration
at which activity complete stopped was
for clustered and
for homogeneous networks (Figure 4C).
A-B Examples of the degradation of neuronal activity in clustered and homogeneous cultures due to the gradual weakening of excitatory connectivity. Both culture types were investigated at the same day in vitro and contained a similar density of neurons. The weakening of connections is achieved by gradually increasing the concentration of CNQX, an AMPA-glutamate receptor antagonist in excitatory neurons. Network response upon weakening is quantified through the relative change in activity
between a given CNQX application and the unperturbed state. Activity variations are indicated separately for each cluster, and shown according to the cluster labeling number. A Clustered cultures show a mixed response upon weakening, with some clusters increasing activity and others reducing it. Only for relatively high concentrations of CNQX (
) the activity systematically decays up to the full silencing of the network. B In homogeneous cultures, activity is analyzed in
regions that cover in a regular manner the entire network. Activity decays almost equally in all regions. Relatively small drug concentrations of [CNQX]
practically suffice to fully stop activity. C Average critical concentration [CNQX]
at which spontaneous activity completed ceases, about
for clustered networks and
for homogeneous ones. Data is averaged over
network realizations of each type of culture. D Photo-damage experiments. Spontaneous activity is measured in cultures that are continuously exposed to strong fluorescence light, causing gradual neuronal degradation and ultimately the death of the entire network. The total radiation received by the neurons is calculated as the duration of the exposure times the area covered by the neurons in the culture (
and
on average for clustered and homogeneous clusters, respectively). The spontaneous activity in homogeneous cultures decays at a much faster rate than in the clustered counterparts. Data is averaged over
network realization of each type. Error bars show standard deviation.
Figure 4D shows the results for the resistance of the networks to node deletion as a consequence of direct photo-damage to the neurons. As can be observed, homogeneous cultures decay in activity much faster than the clustered ones, pinpointing the general resistance of clustered cultures to structural failure.
Discussion
Clustered neuronal cultures have a unique self-organizing potential. An initially isotropic ensemble of individual neurons quickly group to one another to constitute a stable configuration of interconnected clusters of tightly packed neurons. The formation of the clustered network is primarily a passive process governed by the pulling forces exerted by the neurites. Interestingly, aggregation occurs even in the absence of glial cells and neuronal activity [40], and is maintained up to the degradation of the culture [40], [52], [53]. Our work shows that this self-organizing process drives the network towards specific dynamic states, which shape a topology of the functional network that is distinctively assortative. We note that the number of clusters and their distribution are initially random. Therefore, a wide spectrum of physical circuitries and functional topologies are in principle attainable. However, in all the studied cultures, the network drives itself towards markedly assortative topologies with a ‘rich-club’ core. The emergence of these distinct topological traits, concurrently with a stronger network's resilience to activity deterioration, pictures a self-organizing mechanism that enhances network survival by procuring a robust architecture and dynamic stability.
We remark that the link between assortativity and resilience is based on the comparison between the response of clustered and homogeneous cultures upon the same perturbation. To obtain conclusive evidences that assortativity confers resilience traits exclusively from topology, we would require an experimental protocol in which we could arbitrarily ‘rewire’ the connectivity between clusters, or shape in a control manner different circuitries while preserving the number of nodes in the network. Although these strategies are certainly enlightening, they are of difficult development and a major experimental challenge.
We infer the functional connectivity maps of the clustered networks from their spontaneous dynamics. We considered small-sized networks to simultaneously access the entire population ( clusters). The approach that we have used to characterize this functional connectivity is based on the analysis of the time delays between consecutive clusters' activations. The uniqueness of our approach is to use these time delays to provide a direct measure of causality, giving rise to a functional network that is both directed and weighted, with the weights given by a decaying function that follows the frequency of the delay between pairs of clusters. Our formulation is simple and naturally derives from the intrinsic dynamics of the network.
We used two main parameters to quantitatively construct the directed functional network, namely the cut-off time for causality, and the variance of the Gaussian-like weighting function. The cut-off time is set to
, two times the maximum measured time delay between consecutive activations. The importance of the cut-off is first to discriminate two successive bursting episodes, and second to exclude individual firing events from an actual cascade of activations. Although these individual firings account for less than 2% of the total activations, they may occur in regions of the culture that are physically distant -though temporary close- from an actual sequence, and therefore they would add spurious, long-range functional connections to the network.
On the other hand, the variance is obtained from a Gaussian fit of the distribution of consecutive activation delays within bursts. The value of
is specific for each culture to take into account particular differences in the dynamics of the network, specifically the culture days in vitro or the number of clusters (Figure S2), parameters that could affect the delay times of activation. Young cultures for instance exhibit longer time delays between pairs of clusters, leading to a distribution
shifted towards higher values and therefore a larger
.
We tested that the obtained functional networks were stable upon variation of the above parameters. In particular, to examine whether the choice of the cut-off does or does not substantially affect the features of the generated functional network, we performed a sensitivity analysis on this parameter. As the process of generating the network from the sets of bursts is deterministic, we analyzed the influence of the cut-off value on the formed groups of firings. To quantify the variation on the bursts generated for different values of the cut-off, we calculated the variation of information [54] between the grouping of bursts at a certain cut-off value and the previous one as a measure to assess their difference (Figure S3). In the case of clustered cultures, we found that for values of cut-off of the variation of information is, on average, on the order of
. In the homogeneous case, for cut-off values of
ms, this value is on the order of
. This means that varying the cut-off values in these regions does not substantially change the grouping of the bursts, and therefore the generated networks are equivalent.
To assess the goodness of our construction in inferring the functional connectivity of the clustered networks, we compared our connectivity maps with those procured by information theoretic measures, such as Mutual Information or Transfer Entropy, applied to the original fluorescence recordings. These approaches have been used to draw the topological properties of neuronal networks in vitro, both in electrode recordings [55], [56] and calcium fluorescence imaging [57], [58]. The comparison of our method with these theoretic measures showed that the identified functional links were fundamentally the same, with small quantitative differences associated to the particular weighting procedures.
Our functional networks consistently maintained high assortativity values, and along a wide range of days in vitro. We also observed that, by contrast, the assortativity analysis in homogeneous cultures procured neutral or low assortativity values, a result that is supported by other studies in homogeneous networks similar to ours [55]. In our study, we have seen that the clustered, assortative networks exhibit a higher resilience of the network to damage compared to the homogeneous, non-assortative ones. Different studies also highlighted the importance of assortativity and the ‘rich-club’ phenomenon on higher-order structures of the network, in particular resilience, hierarchical ordering and specialization [10], [30].
Several studies in brain networks advocate that the functional connectivity reflects the underlying structural organization [59]–[61]. To shed light on this interrelation in our cultures, we would need the identification of all the physical links between clusters. The top images of Figure 3 indeed suggest that some structural connections could be delineated by a simple visual inspection. However, we observed by green fluorescence protein (GFP) transfection that physical connections have long extensions and may easily link several clusters together, and not just in a first-neighbor manner as seen in the images. Since the images provide a very poor subset of the entire structural layout, a complete description of the physical circuitry must be carried out before comparing the structural and functional networks. Such a detailed identification is difficult, and requires the use of a number of connectivity-labeling techniques. Nevertheless, for the connections that we could visualize, we draw two major conclusions. First, that neither the width of the physical connections nor the size of the clusters were related to a particular trait of the functional links, such as the weight of the connections or the strength of the nodes (Figure S4). And, second, that our construction inferred strong functional links between clusters that were not directly connected in a physical manner, highlighting the importance of indirect paths as well as long-range coupling in the flow of activity.
The identification of the full set of structural connections would certainly provide invaluable information to investigate the interplay between structure and function in our networks. In this context, the recent work by Santos-Sierra et al [52] is enlightening. They analyzed some major structural connectivity traits in clustered networks similar to ours, and observed that the networks were strongly assortative as well. Assortativity emerged at early stages of development, and was maintained throughout the life of the culture. Hence, in clustered cultures, the combined evidences of this study and ours hints at the existence of assortative properties in both structure and function.
An interesting peculiarity of our experiments is that, in most of the studied clustered cultures, the spontaneous bursting episodes comprised of a small subset of clusters rather than the entire network. This activation in the form of groups or moduli is often referred as conditional activity. It contrasts with the coherent activity of homogeneous cultures, where the entire network lights up in a short time window during a bursting episode. Given the acute differences in assortativity between clustered and homogeneous cultures, we hypothesize that the modular dynamics by itself increases or reinforces assortative traits in the functional network.
We finally remark that our neuronal cultures are spatial, i.e. embedded in a physical substrate, which imposes constraints to the layout of connections and, in turn, their assortative characteristics [52], [62]. Spatial networks have caught substantial interest in the last years to understand the restrictions—or advantages—that metric correlations impose on the structure and dynamics of complex networks [63], in particular brain circuits [64]. Vértes et al showed that spatial constraints delineate several topological properties of functional brain networks [65], and Orlandi et al showed that the initiation mechanisms of spontaneous activity are governed by metric correlations inherited by the network during its formation [36]. Strong spatial constraints in clustered networks can be attained by anchoring the neuronal aggregates in specific locations, for instance through carbon nanotubes [39]. The comparison of the functional maps of such a forced organization with our free one is enlightening, and would shed light on the importance of structural constraints in shaping functional connectivity.
To conclude, we have presented a simple yet powerful construction to draw the directed functional connectivity in clustered neuronal cultures. The construed networks present assortativity and ‘rich-club’ features, which are present concurrently with resilience traits. Our analysis has been based on spontaneous activity data, and may certainly vary from evoked activity. Hence, the combined experimental setup and functional construction can be viewed as a model system for complex networks studies, specially to understand the interplay between structure and function, and the emergence of key topological traits from network dynamics. Also, the spatial nature of our experiments may also procure invaluable data to understanding the role of short- and long-range connections in network dynamics; or to investigate the targeted deletion of the high degree nodes that shape the backbone of the network. The latter is a powerful concept that may assist in a detailed exploration of resilience in neuronal circuits, for instance to model the circuitry-activity interrelation in neurological pathologies.
Materials and Methods
Ethics statement
All procedures were approved by the Ethical Committee for Animal Experimentation of the University of Barcelona, under order DMAH-5461.
Clustered neuronal cultures
In our experiments we used cortical neurons from day old Sprague-Dawley rat embryos. Following standard procedures [36], [66] dissection was carried out in ice-cold L-
medium enriched with
glucose and gentamycin (Sigma-Aldrich). Cortices were gently extracted and dissociated by repeated pipetting.
Cortical neurons were plated onto glass coverslips (Marienfeld-Superior) that incorporated a poly-dimethylsiloxane (PDMS) mold. The PDMS restricted neuronal growth to isolated, circular cavities
in diameter. Prior plating, glasses were washed in
nitric acid for 2 h, rinsed with double-distilled water (DDW), sonicated in ethanol and flamed. In parallel to glass cleaning, and following the procedure described by Orlandi et al. [36], several
diameter layers of PDMS
thick were prepared and subsequently pierced with
diameter biopsy punchers (Integra-Miltex). Each pierced PDMS mold typically contained
to
cavities. The PDMS molds were then attached to the glasses and the combined structure autoclaved at
, firmly adhering to one another. For each dissection we prepared
identical glass-PDMS structures, giving rise to about
cultures of
in diameter. Neurons were plated in the PDMS cavities with a nominal density of
, and incubated in plating medium at 37°C,
CO2 and
humidity. Plating medium consisted in
of foetal calf serum (FCS, Invitrogen),
of horse serum (HS, Inivtrogen), and 0.1% B27 (Sigma) in MEM Eagle's-L-glutamate (Invitrogen). MEM was enriched with gentamicin (Sigma), the neuronal activity promoter Glutamax (Sigma) and glucose.
Upon plating, the absence of adhesive proteins in the glass substrate rapidly favored cell-cell attachment and, gradually, the formation of islands of highly compact neuronal assemblies or clusters that minimized the surface contact with the substrate. Clustered cultures formed quickly. By day in vitro (DIV) 2 the culture encompasses dozens of small aggregates that coalesce and grow in size as the culture matures. Spontaneous activity and connections between clusters were observed by DIV . Clusters at this stage of development also anchored at the surface of the glass and, although they continued growing and developing connections, their number and position remained practically stable along the next
weeks. At the moment of measuring, each PDMS cavity contained an independent culture formed by
interconnected clusters.
Clustered cultures were maintained for about weeks, as follows. At DIV
the medium was switched from plating to changing medium (containing
FUDR,
Uridine, and
HS in enriched MEM) to limit glial cell division. Three days later, the medium was replaced to final medium (enriched MEM with
HS), which was then refreshed periodically every three days.
Homogeneous neuronal cultures
Overnight exposure of the glass coverslips to poly-l-lisine (PLL, Sigma) provided a layer of adhesive proteins for the neurons to quickly anchor upon seeding, leading to cultures with a homogeneous distribution of neurons over the substrate. The remaining steps in the preparation and maintenance of the cultures were identical as the clustered ones, i.e. we used the same nominal neuronal density for plating, we included PDMS pierced molds to confine neuronal growth in cavities in diameter, and we refreshed the culture mediums in the same manner.
Experimental setup and procedure
Standard experiments.
To measure the spontaneous activity in the clustered networks we used cultures at day in vitro (DIV) , i.e. covering about two weeks of development. Cultures started to degrade by DIV
, and therefore we did not use cultures older than
weeks in our experiments.
Activity in neuronal cultures was monitored through fluorescence calcium imaging [67], [68], which allows the detection of neuronal activity by the binding of ions to a fluorescence probe upon firing. Prior to recording, the cultures under study were incubated for
in External Medium (EM, consisting of
NaCl,
CaCl2,
MgCl2,
sucrose,
glucose, and
Hepes; pH 7.4) in the presence of Fluo-4-AM (Invtrogen). We used
Fluo4 in a volume of
EM. We incubated a glass coverslip containing
cultures within the PDMS cavities at once, allowing for the simultaneous recording of different cultures or the selection of cultures with specific traits. After incubation, the cultures were washed with fresh EM and placed in the observation chamber, consisting of a standard glass bottom culture dish, filled with
ml EM, and with its wall and cover screened from external light. To minimize accidental damage to the aggregates during the manipulation of the cultures, the PDMS pierced mold was left in contact with the glass during both incubation and the actual experiment.
The observation chamber was mounted on Zeiss Axiovert inverted microscope equipped with a high-speed CMOS camera (Hamamatsu Orca Flash 2.8). We used an objective of X combined with a
X optical zoom. These settings provided a final field of view of
(
)
that supported the recording of
or
PDMS-confined cultures simultaneously. Individual frames were acquired as 8-bit grey-scale images, a size of
pixels, and a spatial resolution of
. All experiments were carried out at room temperature.
The fluorescence signal of the clusters' spontaneous activity was recorded with the software Hokawo 2.5, provided by the camera vendor (Hamamatsu Photonics). We used acquisition speeds in the range frames per second (fps), corresponding to a time interval of
between consecutive frames. These acquisition speeds were selected to optimize the balance between image quality, sufficient time resolution, and minimum light intensity. The latter was particularly important to minimize photo-damage and photo-bleaching, and allowed neuronal cultures to be studied with optimal conditions for at least 3 h. However, the combination of high acquisition speeds and high resolution images resulted in large data files—of at least
GB per hour of recording—that had to be stored and analyzed. We therefore limited most of our experiments to about 1 h of recording, which was sufficient to reliably build the functional networks, as described in the ‘Control experiments’ section.
Measurements in homogeneous cultures were carried out in the same way, with the only difference that the recording speed was increased to –
fps to take into account the fast propagation of activity fronts in these preparations, as observed for instance in the study of Orlandi et al. [36].
Resilience experiments.
We considered two groups of resilience experiments. In a first group, we monitored the gradual degradation of network activity due to photo-damage. In a second group, we measured the decay in activity as a consequence of the gradual weakening of the excitatory connectivity.
In the first group of measurements, we first considered a clustered culture and measured its spontaneous activity uninterruptedly along 2 hours, with neurons continuously exposed to a light radiation times stronger than normal. We then divided the sequence in blocks of 30 min, and determined, for each block, the average network activity by counting the number of bursting episodes within the block. Next, we switched to a homogeneous culture from the same batch (i.e. identical nominal density and age) and carried out the same protocol. In total we carried out
measurements for each kind of culture, and finally analyzed the decay in activity as a function of time. Although the radiation over the culture was homogeneous, the actual area occupied by the neurons was different between homogeneous and clustered cultures. Neurons in homogeneous cultures formed a monolayer that covered
of the available area, corresponding to about
for the
diameter wells. In clustered cultures, neurons were tightly packed in slightly three-dimensional structures, giving rise to a lower spatial coverage of
, i.e. about
. Hence, homogeneous cultures were effectively more exposed to light than their clustered counterparts. To take this spatial variability into account, the ‘total radiation’ received by the neurons in a given experiment was quantified as the duration of the light exposure times the area occupied by the neurons in the studied culture. We then averaged the results over the
different culture realizations of each type, and binned nearby values for clarity. The comparison of the activity-radiation plots between the two networks (Figure 4D) indicated which topology exhibited higher resistance to degradation in neuronal activity.
In the second group of measurements, we compared the change in activity between a clustered and a homogeneous culture during gradual weakening of neuronal connectivity. The weakening was achieved by progressive application of CNQX [37], [69], an AMPA-glutamate receptor antagonist in excitatory neurons (see also ‘Pharmacology’). We first measured the clustered network and thereafter the homogeneous one. In both cases, we first recorded spontaneous activity at , and used the average firing rate
as reference for the subsequent steps. We then increased the concentration of the drug to a preset value, waited
min for the drug to take effect and measured again for
, computing the new average firing rate
. We switched to a second preset values, and repeated the procedure until activity ceased. The relative decay in activity
is used to illustrate the gradual fall of activity in Figure 4C. The critical concentration [CNQX]
at which activity is absent along the
of recording (
) hints at the robustness of network dynamics to a global failure of its connectivity.
Pharmacology
The pharmacological protocols described below were used identically in clustered and homogeneous cultures.
Inhibitory connections.
The in vitro networks contain both excitatory and inhibitory connections. However, for sake of simplicity in the comparison between experiments and model, in the experiments at DIV and above we completely blocked γ-aminobutyric acid (GABA) inhibitory synapses with
of the antagonist bicuculine methiodide (Sigma). The drug was applied
before the actual recordings for the drug to take effect. Spontaneous activity in our experiments is therefore solely driven by excitatory connections. Although the balance between excitation and inhibition shapes the major traits of spontaneous activity [36], [37], such as the average firing rate of the network, we verified that the presence of inhibition did not qualitatively modify the results presented here.
We left active inhibitory synapses for experiments at DIV since at this early stages of development GABA has a depolarizing effect and therefore an excitatory action [37], [69], [70]. Its blockade would effectively reduce excitation and silence the network.
Network connectivity weakening through CNQX.
In the studies of network resilience to the weakening of connectivity, we studied the decay in spontaneous activity as a result of the gradual application of 6-cyano-7-nitroquinoxaline-2,3-dione (CNQX, Sigma), an AMPA-glutamate receptor antagonists in excitatory neurons. For the connectivity strength between neurons is maximum. As [CNQX] is administered, the efficacy of excitatory connections steadily diminishes, which is accompanied by a reduction in spontaneous activity (see e.g. [37] for illustrative data). High CNQX concentrations lead to a complete halt in activity. In the measurements we used CNQX concentrations in the range
, in quasi-logarithmic steps. We left the culture unperturbed for 5 min after each CNQX application for the drug to reach steady-state effects.
Data analysis for clustered cultures
The acquired images (recorded at a typical speed of fps) were first analyzed with the Hokawo 2.5 software to extract the fluorescence intensity of each cluster as a function of time. The regions of interest (ROIs) were chosen manually and typically covered an area of
pixels, each ROI corresponding to a single cluster. As illustrated in Figure 1C and Figure 1F, activity is characterized by a stable baseline (resting state) interrupted by peaks of fluorescence that correspond to clusters' firings. At the onset of firing, the fluorescence signal increases abruptly due to the fast intake of
ions. Fluorescence then reaches a maximum, and slowly decays back to the baseline in
s.
The algorithm that we used to detect the onset of firing for each cluster was as follows. We first corrected the fluorescence signal from small drifts, and calculated the resting fluorescence level
by discarding the data points with an amplitude two times above the standard deviation (SD) of the signal. The corrected signal was then expressed as
. We next took
and computed its derivative
in order to detect fast changes in the fluorescence signal. Finally, the onset of ignition in cluster was defined as the time where a maximum in
was accompanied by values of
two times above the SD of the background signal, and for at least 5 frames.
Reliability in detecting the clusters' ignition times.
Three major tests were carried out to assess the reliability of our analysis. In a first one, we measured spontaneous activity at fps, i.e. twice the standard recording speed, but used stronger light to compensate for the lower exposure time. We next analyzed the data, re-sampled the image sequence down to
fps and compared the results with the original acquisition. We observed that the detection of the onset times improved only by about 15%, which did not justify the excess of light and the associated damage to the neurons.
In a second test, we measured spontaneous activity in a culture using identical light settings but considering different acquisition rates, namely ,
, and
fps. We then selected ignition sequences that were as similar as possible in all three measurements, and compared the results. We observed that only in the few cases where the clusters fired with strong amplitudes the increased speed enhanced detection, and again by
. For the rest of the cases, the higher speeds actually worsened the analysis due to the poorer signal-to-noise ratio.
Finally, in a third test, we used sub-frame resolution analysis tools to evaluate the importance of finer ignition times. Following Orlandi et al. [36], we considered the approach of fitting two straight lines at the vicinity of each initially detected firing. A first fit included the points of the background signal that preceded ignition, and a second one extended to the
points that correspond to the fast rise in fluorescence. The crossing value of the two lines provided an onset time that refined the initially measured value. The better accuracy effectively increased the discrimination of sequences that were initially identified as simultaneous events. However, since these events are rare (by
), the finer temporal resolution had practically no effect in the construction of the functional networks and the derived analysis.
Activity propagation times.
The time delay in the propagation of activity between two connected clusters was measured in control experiments with high acquisition rates. We concluded that
varied in the range
ms, with an average value
. Other studies in clustered networks provided similar results [47]. With the detection algorithm described above and standard experiments at
fps, we could appraise the activation sequence in
of the cases. The remaining
corresponded to clusters that ignited in the same frame or time bin, and were treated as simultaneous events.
Data analysis for homogeneous cultures
Recordings in homogeneous cultures provided the activity of neurons in an circular area
in diameter. Neurons were marked individually as regions of interest in the images and the corresponding fluorescence time traces extracted using custom-made software. Ignition times for each neuron were next obtained by using the sub-frame resolution method described above (detailed in Ref. [36]), and that consisted in fitting two straight lines to the fluorescence data, a first fit encompassing the
points in the background region prior to firing, and a second fit including the
points during the fast rise in fluorescence that follows ignition. The crossing point of the two lines provided the onset of firing.
The extension of this analysis to all the active neurons within a burst, and along all the bursts, finally provided the entire set of ignition sequences. The construction of the directed functional networks for the homogeneous cultures was then carried out identically as the clustered ones.
Additional control experiments
Recordings in clustered cultures typically lasted for h and contained between
bursts in the quietest networks and
bursts in the most active ones. To test whether
bursts sufficed to draw the functional networks, we carried out a control experiment in which we monitored spontaneous activity along
h in a standard clustered culture, measured at DIV
and containing
nodes (Figure S6). We then analyzed the data using two different procedures. In the first one we drew the functional connectivity using the data extracted from the entire recording, and determined its assortativity values. In the second procedure, we separated the recorded sequence in three blocks, each
long, built the functional connectivity for each block, and computed the respective assortative values. The studied culture fired in a sustained manner at a rate of
, and procured a total of
bursts. Thus, each block typically contained about
bursts.
The results (Figure S6) led to two major conclusions. First, that the functional connectivity is very similar among the blocks, and between any of the blocks and the entire recording, providing assortativity values that are compatible within statistical error. And second, that the first 40 min of recording (with 45 bursts only) sufficed to shape the major traits of the functional network, therefore validating our strategy of using h of acquisition to procure a reliable estimate of the functional connectivity of the network and its assortative traits.
Network assortativity and rich-club
Here we describe the calculation of the assortativity coefficients, assortativity errors and the rich-club distributions. In the process, we have to define the assortativity for directed weighted networks.
Generalization of assortativity to directed weighted networks.
Newman [23] defined assortativity as the Pearson correlation between the degrees of every pair
of linked nodes in the network. More precisely, in the case of directed networks, if
is the output degree from node
,
the input degree to node
, and
scans all the edges in the network, then
(1)where
(2)
(3)
(4)
After some algebra, assortativity may also be written as (5)
(6)
This definition of assortativity is applicable to all kind of networks, either undirected, directed, weighted or unweighted. For weighted networks as in our case, the strength of the nodes carries important information about the structure of the network, and thus it would be useful to know the correlation between the strengths instead of the degrees. Since in this case each edge carries a weight, it seems logical that edges with higher weight should have a larger contribution to the correlation. Therefore, we define the weighted assortativity as the Pearson weighted correlation between the strengths of the nodes. In mathematical terms, if
is the weight of the link from node
to node
(zero if there is no link), then
and
are the output and input strengths, then
(7)where
(8)
(9)with
the total strength of the network, i.e.
(10)
Litvak et al. [31] showed that in disassortative networks the magnitude of the standard assortativity decreases with network size, a problem that was solved by replacing the Pearson correlation with the Spearman correlation, thus obtaining a Spearman assortativity
. Spearman rank correlation is calculated in the same way that the Pearson correlation but substituting the values (in this case, the degrees of the nodes) by their respective ranks, i.e. their position when the values are sorted in ascending order. This leads us to define the Spearman weighted assortativity
as the Spearman weighted correlation between the strengths of the nodes at both ends of each edge in the network.
The estimation of the error in the assortativity value (for any of the previous four variants) can be computed in several ways, for instance through the jackknife method, the bootstrap algorithm, or by using the Fisher transformation [24], [71]. We used in our work the bootstrap algorithm, and considered random samples of the data.
Rich-club analysis.
The rich-club analysis computes the degree-degree (or weight-weight) correlation distributions, and respect to a null case of non-correlated degrees (or weights). It allows us to reinforce the assortativity analysis presented before. Assortative mixing, in principle, will induce a rich-club effect that should be clearly detectable for a wide range of degrees (or weights). We will first introduce the formulation for the calculation of rich-club in weighted networks as presented in [51], and afterwards we will extend it to the case of weighted directed networks.
The rich-club score is calculated as follows: (11)where
is the rich-club score relative to the uncorrelated null case,
is the sum of the weights of the links of the subgraph formed only by those nodes whose strengths are higher than
,
(12)and
is the corresponding value in the case of uncorrelated strengths,
(13)
The term designates the subset of nodes
such that
.
is the total number of nodes in the network, and
and
are the first and second moments of the strength distribution. The ratio
is calculated for all values of
, and ranges from the minimum value of strength in the network to the maximum. This ratio indicates the presence or absence of a rich-club in the network: a network shows a rich-club effect when the high values of
give a ratio above 1.
To calculate this ratio on our functional networks, we have to consider that the network is not only weighted but also directed. Therefore, we need to adapt the former formulation for the case of weighted directed networks. For this reason we will consider the in- and out-strength of each node, expressed as and
respectively. In the directed formulation, Eqs. (11) and (12) remain unchanged. However,
must be redefined as
(14)and the term
becomes
(15)where the averages are calculated as
(16)
This formulation allows us to calculate the rich-club coefficient for weighted directed networks. Note that for an undirected network (where ) the latter formulation reduces to the original one.
Alternative construction of the functional network based on mutual information
Mutual information [56], [72] is a particular case of the Kullback-Leibler divergence [73], an information-theoretic measure of the distance between two probability distributions. In fact, the mutual information between two stochastic variables and
provides an estimation of the amount of information gained about
when
is known.
Let us indicate by the time series corresponding to the
-th cluster, with
and
the total number of time frames involved in the observation process. The time series adopted for the successive analysis are obtained by mapping the observed train of cluster activations to another time series termed walk, defined by
(17)
In the specific case of our analysis, the mutual information between two time series and
, corresponding to two different clusters, is interpreted as the amount of correlation between the dynamics of cluster
and
. In general, the time scale of the correlation between two time series is not known a priori. Such a time scale corresponds to the time delay required to maximize the gain of information. Therefore, in the spirit of Fraser and Swinney [74], we define the time delayed mutual cross information between
and
by
(18)where
and
are indices running over some partition of the observed time series. In Eq. (18),
indicates the probability to find a value of time series
in the
-th interval,
is the probability to find a value of time series
in the
-th interval, whereas
denotes the joint probability to observe a firing from the
-th cluster falling in the
-th interval and a firing from the
-th cluster falling in the
-th interval exactly
time frames later.
For the sake of simplicity, in the following we will adopt the more concise notation to indicate the time delayed mutual cross information. Finally, in order to gain the highest amount of information about the dynamics of cluster
by observing cluster
, we consider only the maximum value
of
with respect to the time delay
.
We estimate the importance of the observed amount of correlation by performing the above analysis on surrogate data. Surrogates adopted in this study are time series generated by randomly reshuffling the temporal observations of the firing series , for each cluster separately. Such a procedure destroys any correlation between pairs of time series while preserving the empirical probability distribution, thus allowing to test the null hypothesis that the observed correlation is obtained by chance.
We indicate by the walk corresponding to the surrogate obtained from time series
and with
the time delayed mutual cross information between
and
. We perform 200 independent random realizations of surrogates for each pair
and we estimate the corresponding expected value
of the maximum mutual cross-information, as well as the root mean square
of the underlying distribution.
Hence, we fix a priori the significance of the hypothesis testing and we estimate the
-score corresponding to each pair
by
. Therefore, the observed correlation between cluster
and
is said to be statistically significant if
, where
is the standard error function. Finally, we obtain the functional network of clusters by building the weight matrix
whose elements are defined by
if
, and
if
.
Supporting Information
Figure S1.
Homogeneous cultures. A The seeding of neurons in cover glasses previously coated with poly-l-lysine gives rise to neuronal cultures with a quasi-homogeneous distribution of neurons. A typical circular culture in diameter contains about
neurons that can be well identified either as bright spots in the fluorescence recordings or as circular objects in bright-field images (small panel on the right). Spontaneous activity in these homogeneous cultures is typically recorded at
frames/s, which suffices to extract the time delays between consecutive neuronal activations. The particular experiment shown here corresponds to network ‘P’ of the main text, with a total of
neurons manually selected over the images and monitored along
. The analysis of their spontaneous activity traces is analyzed in the context of our model, finally procuring the functional connectivity network and its topological properties. B Given the large number of nodes analyzed and the high average degree of the resulting network (
functional connections per neuron), a representation of the complete functional network is unpractical. As an example of the obtained functional networks, we here show all the functional links within a small region placed in the center containing
neurons. Connections are both color and thickness coded according to their weight. Nodes are color coded according to their strength. C As an alternative representation, we show here a
of the population (
neurons randomly chosen), each neuron showing the
of its links. Nodes and links are color coded according to their strength and weight, respectively. D A ring graph of the same neurons shows that most of them display a similar connectivity, in contrast to the strong modularity and variability in connectivity exhibited by the clustered cultures.
https://doi.org/10.1371/journal.pcbi.1003796.s001
(PDF)
Figure S2.
Variance and culture properties. A The variance
is obtained from the Gaussian fit of the activation delays between pairs of clusters. The plot shows that
decreases with the culture age in vitro, indicating that young cultures display a slower dynamics (larger delay times and therefore larger variance) than mature cultures. B The variance
increases with the number of clusters in the culture, indicating a broader and richer distribution of time delays as more clusters participate in the dynamics of the network. Errors bars in A show standard deviation.
https://doi.org/10.1371/journal.pcbi.1003796.s002
(PDF)
Figure S3.
Sensitivity of the functional network construction to the cut-off times. The cut-off determines the end of a sequence and therefore its variation modifies the set of burst chosen. The cut-off is set to for clustered cultures and
for homogeneous ones. To assess the sensitivity of the grouping of bursts to the cut-off, we have computed the Variation of Information (VI) between the grouping of bursts at a certain cut-off value and the previous one. It is computed as
, where
and
are two partitions,
is the entropy and
is the mutual information. In our case, each partition is the set of burst found at a certain value of the cut-off. We have screened the cut-off values from
to
. The analysis shows that
for homogeneous and
for clustered cultures are values for which the variation of information is already stabilized. Thus, the modification of the cut-off within the stabilization region does not change the grouping of clusters in each burst, and therefore the derived functional networks, as well as the corresponding network measures, remain the same.
https://doi.org/10.1371/journal.pcbi.1003796.s003
(PDF)
Figure S4.
Relation between the structural traits of the network and the functional connectivity. A Dependence of the weight of the functional connections on the width of physical connections between directly connected clusters. The sketch conceptually shows the comparison between structural and functional links. The plot represents the analysis of pairs of clusters, with data binned for similar widths. No significant correlation is observed. B The dependence of the node strength on cluster size shows no correlation, indicating that the functional connectivity cannot be assessed from the size of the clusters. Data is based on the analysis of
clusters. C For the same clusters, this plot shows that the activity of a cluster is independent of its size. D Activity within a burst is always initiated by a particular cluster, which triggers the sequential activation of all the downstream clusters. To quantify the importance of these ‘initiators of activity’ in network dynamics we computed the number of times that a cluster of a given size initiates a sequence of activations. The plot shows that there is no a significant correlation between initiation and size. The analysis is based on the study of
bursts. All these results indicate that the functional connectivity cannot be drawn from a visual inspection of the neuronal culture. Errors bars show standard deviation.
https://doi.org/10.1371/journal.pcbi.1003796.s004
(PDF)
Figure S5.
‘Rich-club’ analysis. The evaluation of the rich-club is performed by computing the ratio between the connectivity strength of highly connected nodes and its randomized counterpart,
, and for gradually larger values of the strength threshold
. The figure shows the rich-club analysis for 3 representative clustered and homogeneous cultures. Clustered networks exhibit values of
systematically higher than
for large values of the strength threshold
, evidencing the existence of a rich-club core of highly connected clusters in the network. On the contrary, homogeneous cultures display a mixture of positive and negative values, and with an average around
, ruling out the existence of the rich-club property.
https://doi.org/10.1371/journal.pcbi.1003796.s005
(PDF)
Figure S6.
Control experiment. A Bright field image of a clustered network whose spontaneous activity has been recorded for . The average bursting rate of the network is
. B Corresponding functional network. The size of the nodes is proportional to the size of the actual clusters, and their color is proportional to their strength. The weights of the links are both color and thickness coded. The darker the color, the higher the value of the observable. C Analysis of the 2 h recording in three blocks,
in duration each, and containing
bursts. The blocks show very similar traits between them, as well as with the entire recording. The blocks exhibit similar assortativity values, and share both the most important links and nodes’ strengths.
indicates the assortativity value of the depicted network, averaged over the Pearson and Spearman formulations.
https://doi.org/10.1371/journal.pcbi.1003796.s006
(PDF)
Movie S1.
Representative recording in clustered networks. The movie shows the spontaneous activity in the network depicted in Figure 1, showing the first of recorded activity. The original sequence was acquired at
frames/s. It is reproduced here 15 times faster. The ignition of a cluster is revealed by a sharp increase in the fluorescence signal. Bursts of activity correspond to the ignition of
clusters in a short time window. It is worth remarking the rich activity, the strong modularity, and the repetition of some sequences. Activity was driven solely by excitatory connectivity, with inhibition blocked with
bicuculline.
https://doi.org/10.1371/journal.pcbi.1003796.s007
(AVI)
Acknowledgments
The authors acknowledge S. Rüdiger and J. G. Orlandi for fruitful discussions and insight, and to E. Tibau for technical assistance.
Author Contributions
Conceived and designed the experiments: JS ST. Performed the experiments: ST JS. Analyzed the data: ST CG SG. Contributed reagents/materials/analysis tools: CG ST MDD SG AA. Wrote the paper: JS AA ST CG SG. Conceived the model: AA SG.
References
- 1.
Dorogovtsev SN, Mendes JFF (2003) Evolution of Networks: From Biological Nets to the Internet and WWW. Oxford University Press.
- 2. Alon U (2007) Network motifs: theory and experimental approaches. Nat Rev Genet 8: 450–461.
- 3.
Barrat A, Barthélemy M, Vespignani A (2008) Dynamical processes on complex networks. New York, NY, USA: Cambridge University Press.
- 4.
Newman M (2010) Networks: an Introduction. Oxford: Oxford University Press.
- 5.
Estrada E (2011) The Structure of Complex Networks. Oxford: Oxford University Press.
- 6. Eckmann JP, Feinerman O, Gruendlinger L, Moses E, Soriano J, et al. (2007) The physics of living neural networks. Phys Rep 449: 54–76.
- 7. Sun Y, Huang Z, Yang K, Liu W, Xie Y, et al. (2011) Self-organizing circuit assembly through spatiotemporally coordinated neuronal migration within geometric constraints. PLoS ONE 6: e28156.
- 8. Gross GW, Kowalski JM (1999) Origins of activity patterns in self-organizing neuronal networks in vitro. Journal of Intelligent Material Systems and Structures 10: 558–564.
- 9. Bullmore E, Sporns O (2009) Complex brain networks: Graph theoretical analysis of structural and functional systems. Nature Reviews Neuroscience 10: 186–198.
- 10. Sporns O (2011) The human connectome: a complex network. Annals of the New York Academy of Sciences 1224: 109–125.
- 11. Sporns O, Chialvo DR, Kaiser M, Hilgetag CC (2004) Organization, development and function of complex brain networks. Trends in Cognitive Sciences 8: 418–425.
- 12. Bassett DS, Wymbs NF, Porter MA, Mucha PJ, Carlson JM, et al. (2011) Dynamic reconfiguration of human brain networks during learning. Proceedings of the National Academy of Sciences 108: 7641–7646.
- 13. Seeley WW, Crawford RK, Zhou J, Miller BL, Greicius MD (2009) Neurodegenerative diseases target large-scale human brain networks. Neuron 62: 42–52.
- 14. Zhou J, Gennatas ED, Kramer JH, Miller BL, Seeley WW (2012) Predicting regional neurodegeneration from the healthy brain functional connectome. Neuron 73: 1216–1227.
- 15. Bonifazi P, Goldin M, Picardo MA, Jorquera I, Cattani A, et al. (2009) Gabaergic hub neurons orchestrate synchrony in developing hippocampal networks. Science 326: 1419–1424.
- 16. Latora V, Marchiori M (2001) Efficient behavior of small-world networks. Phys Rev Lett 87: 198701.
- 17. Watts D, Strogatz S (1998) Collective dynamics of ‘small-world’ networks. Nature 393: 440.
- 18. Sporns O, Zwi J (2004) The small world of the cerebral cortex. Neuroinformatics 2: 145–162.
- 19. Harriger L, van den Heuvel MP, Sporns O (2012) Rich club organization of macaque cerebral cortex and its role in network communication. PLoS ONE 7: e46497.
- 20. Hagmann P, Cammoun L, Gigandet X, Meuli R, Honey CJ, et al. (2008) Mapping the structural core of human cerebral cortex. PLoS Biol 6: e159.
- 21. Meunier D, Lambiotte R, Bullmore ET (2010) Modular and hierarchically modular organization of brain networks. Frontiers in Neuroscience 4: 200.
- 22. Zamora-López G, Zhou C, Kurths J (2010) Cortical hubs form a module for multisensory integration on top of the hierarchy of cortical networks. Frontiers in Neuroinformatics 4: 1.
- 23. Newman M (2002) Assortative mixing in networks. Phys Rev Lett 89: 208701.
- 24. Newman M (2003) Mixing patterns in networks. Phys Rev E 67: 26126.
- 25. Eguíluz VM, Chialvo DR, Cecchi GA, Baliki M, Apkarian AV (2005) Scale-free brain functional networks. Phys Rev Lett 94: 018102.
- 26. Avalos-Gaytán V, Almendral JA, Papo D, Schaeffer SE, Boccaletti S (2012) Assortative and modular networks are shaped by adaptive synchronization processes. Phys Rev E 86: 015101.
- 27. de Franciscis S, Johnson S, Torres JJ (2011) Enhancing neural-network performance via assortativity. Phys Rev E 83: 036114.
- 28. Rubinov M, Sporns O (2010) Complex network measures of brain connectivity: Uses and interpretations. NeuroImage 52: 1059–1069.
- 29. Achard S, Salvador R, Whitcher B, Suckling J, Bullmore E (2006) A resilient, low-frequency, small-world human brain functional network with highly connected association cortical hubs. The Journal of Neuroscience 26: 63–72.
- 30. Colizza V, Flammini A, Serrano MA, Vespignani A (2006) Detecting rich-club ordering in complex networks. Nature physics 2: 110–115.
- 31. Litvak N, van der Hofstad R (2013) Uncovering disassortativity in large scale-free networks. Phys Rev E 87: 022801.
- 32. Leung C, Chau H (2007) Weighted assortative and disassortative networks model. Physica A: Statistical Mechanics and its Applications 378: 591–602.
- 33. Wheeler B, Brewer G (2010) Designing neural networks in culture. Proceedings of the IEEE 98: 398–406.
- 34. Wagenaar DA, Pine J, Potter SM (2006) An extremely rich repertoire of bursting patterns during the development of cortical cultures. BMC Neurosci 7: 11.
- 35. Cohen E, Ivenshitz M, Amor-Baroukh V, Greenberger V, Segal M (2008) Determinants of spontaneous activity in networks of cultured hippocampus. Brain Research 1235: 21–30.
- 36. Orlandi JG, Soriano J, Alvarez-Lacalle E, Teller S, Casademunt J (2013) Noise focusing and the emergence of coherent activity in neuronal cultures. Nature Physics 9: 582–590.
- 37. Tibau E, Valencia M, Soriano J (2013) Identification of neuronal network properties from the spectral analysis of calcium imaging signals in neuronal cultures. Frontiers in Neural Circuits 7: 199.
- 38. Teller S, Soriano J (2013) Experiments on clustered neuronal networks. AIP Conference Proceedings 1510: 244–246.
- 39. Gabay T, Jakobs E, Ben-Jacob E, Hanein Y (2005) Engineered self-organization of neural networks using carbon nanotube clusters. Physica A: Statistical Mechanics and its Applications 350: 611–621.
- 40. Segev R, Benveniste M, Shapira Y, Ben-Jacob E (2003) Formation of electrically active clusterized neural networks. Phys Rev Lett 90: 168101.
- 41. Shein Idelson M, Ben-Jacob E, Hanein Y (2010) Innate synchronous oscillations in freely-organized small neuronal circuits. PLoS ONE 5: e14443.
- 42. Soussou W, Yoon G, Brinton R, Berger T (2007) Neuronal network morphology and electrophysiologyof hippocampal neurons cultured on surface-treated multielectrode arrays. Biomedical Engineering, IEEE Transactions on 54: 1309–1320.
- 43. Sorkin R, Gabay T, Blinder P, Baranes D, Ben-Jacob E, et al. (2006) Compact self-wiring in cultured neural networks. Journal of Neural Engineering 3: 95.
- 44. Macis E, Tedesco M, Massobrio P, Raiteri R, Martinoia S (2007) An automated microdrop delivery system for neuronal network patterning on microelectrode arrays. Journal of Neuroscience Methods 161: 88–95.
- 45. Shein M, Greenbaum A, Gabay T, Sorkin R, David-Pur M, et al. (2009) Engineered neuronal circuits shaped and interfaced with carbon nanotube microelectrode arrays. Biomedical Microdevices 11: 495–501.
- 46. Shein-Idelson M, Ben-Jacob E, Hanein Y (2011) Engineered neuronal circuits: A new platform for studying the role of modular topology. Frontiers in Neuroengineering 4: 10.
- 47. Tsai CY, Chang MC, I L (2008) Robustness and Variability of Pathways in the Spontaneous Synchronous Bursting of Clusterized Cortical Neuronal Networks In vitro. Journal of the Physical Society of Japan 77: 084803.
- 48. Yvon C, Rubli R, Streit J (2005) Patterns of spontaneous activity in unstructured and minimally structured spinal networks in culture. Experimental Brain Research 165: 139–151.
- 49. Berdondini L, Chiappalone M, van der Wal P, Imfeld K, de Rooij N, et al. (2006) A microelectrode array (mea) integrated with clustering structures for investigating in vitro neurodynamics in confined interconnected sub-populations of neurons. Sensors and Actuators B: Chemical 114: 530–541.
- 50. Zlatic V, Bianconi G, Díaz-Guilera A, Garlaschelli D, Rao F, et al. (2009) On the rich-club effect in dense and weighted networks. The European Physical Journal B 67: 271–275.
- 51. Serrano MA (2008) Rich-club vs rich-multipolarization phenomena in weighted networks. Phys Rev E 78: 026101.
- 52. de Santos-Sierra D, Sendiña Nadal I, Leyva I, Almendral JA, Anava S, et al. (2014) Emergence of small-world anatomical networks in self-organizing clustered neuronal cultures. PLoS ONE 9: e85828.
- 53. Shefi O, Golding I, Segev R, Ben-Jacob E, Ayali A (2002) Morphological characterization of in vitro neuronal networks. Physical Review E 66: 021905.
- 54.
Meilă M (2003) Comparing clusterings by the variation of information. In: Schölkopf B, Warmuth MK, editors, Learning Theory and Kernel Machines, Springer Berlin Heidelberg, volume 2777 of Lecture Notes in Computer Science. pp. 173–187.
- 55. Bettencourt LMA, Stephens GJ, Ham MI, Gross GW (2007) Functional structure of cortical neuronal networks grown in vitro. Phys Rev E 75: 021915.
- 56. Garofalo M, Nieus T, Massobrio P, Martinoia S (2009) Evaluation of the performance of information theory-based methods and cross-correlation to estimate the functional connectivity in cortical networks. PLoS ONE 4: e6482.
- 57. Stetter O, Battaglia D, Soriano J, Geisel T (2012) Model-Free Reconstruction of Excitatory Neuronal Connectivity from Calcium Imaging Signals. PLoS Comput Biol 8: e1002653.
- 58. Orlandi JG, Stetter O, Soriano J, Geisel T, Battaglia D (2014) Transfer entropy reconstruction and labeling of neuronal connections from simulated calcium imaging. PLoS ONE 9: e98842.
- 59. Deco G, Jirsa VK, McIntosh AR (2010) Emerging concepts for the dynamical organization of resting-state activity in the brain. Nature Reviews Neuroscience 12: 43–56.
- 60. Honey CJ, Sporns O, Cammoun L, Gigandet X, Thiran JP, et al. (2009) Predicting human restingstate functional connectivity from structural connectivity. Proceedings of the National Academy of Sciences 106: 2035–2040.
- 61. Honey CJ, Kötter R, Breakspear M, Sporns O (2007) Network structure of cerebral cortex shapes functional connectivity on multiple time scales. Proceedings of the National Academy of Sciences 104: 10240–10245.
- 62. Schmeltzer C, Soriano J, Sokolov IM, Rüdiger S (2014) Percolation of spatially constrained erdösrényi networks with degree correlations. Phys Rev E 89: 012116.
- 63. Barthélemy M (2011) Spatial networks. Physics Reports 499: 1–101.
- 64. Bullmore E, Sporns O (2012) The economy of brain network organization. Nature Reviews Neuroscience 13: 336–349.
- 65. Vértes PE, Alexander-Bloch AF, Gogtay N, Giedd JN, Rapoport JL, et al. (2012) Simple models of human brain functional networks. Proceedings of the National Academy of Sciences 109: 5868–5873.
- 66. Segal M, Manor D (1992) Confocal microscopic imaging of [Ca2+]i in cultured rat hippocampal neurons following exposure to n-methyl-d-aspartate. The Journal of Physiology 448: 655–676.
- 67. Takahashi N, Takahara Y, Ishikawa D, Matsuki N, Ikegaya Y (2010) Functional multineuron calcium imaging for systems pharmacology. Analytical and Bioanalytical Chemistry 398: 211–218.
- 68. Grienberger C, Konnerth A (2012) Imaging calcium in neurons. Neuron 73: 862–885.
- 69. Soriano J, Rodríguez Martínez M, Tlusty T, Moses E (2008) Development of input connections in neural cultures. Proceedings of the National Academy of Sciences 105: 13758–13763.
- 70. Ganguly K, Schinder AF, Wong ST, ming Poo M (2001) GABA itself promotes the developmental switch of neuronal GABAergic responses from excitation to inhibition. Cell 105: 521–532.
- 71.
Efron B (1979) Computers and the theory of statistics: Thinking the unthinkable. SIAM Review 21 : pp. 460–480.
- 72. Singh A, Lesica NA (2010) Incremental mutual information: A new method for characterizing the strength and dynamics of connections in neuronal circuits. PLoS Comput Biol 6: e1001035.
- 73. Kullback S, Leibler RA (1951) On information and sufficiency. The Annals of Mathematical Statistics 22: 79–86.
- 74. Fraser AM, Swinney HL (1986) Independent coordinates for strange attractors from mutual information. Physical review A 33: 1134.