Research ArticleNEUROSCIENCE

Spatiotemporal ontogeny of brain wiring

See allHide authors and affiliations

Science Advances  12 Jun 2019:
Vol. 5, no. 6, eaav9694
DOI: 10.1126/sciadv.aav9694

Abstract

The wiring of vertebrate and invertebrate brains provides the anatomical skeleton for cognition and behavior. Connections among brain regions are characterized by heterogeneous strength that is parsimoniously described by the wiring cost and homophily principles. Moreover, brains exhibit a characteristic global network topology, including modules and hubs. However, the mechanisms resulting in the observed interregional wiring principles and network topology of brains are unknown. Here, with the aid of computational modeling, we demonstrate that a mechanism based on heterochronous and spatially ordered neurodevelopmental gradients, without the involvement of activity-dependent plasticity or axonal guidance cues, can reconstruct a large part of the wiring principles (on average, 83%) and global network topology (on average, 80%) of diverse adult brain connectomes, including fly and human connectomes. In sum, space and time are key components of a parsimonious, plausible neurodevelopmental mechanism of brain wiring with a potential universal scope, encompassing vertebrate and invertebrate brains.

INTRODUCTION

Cognition and behavior rely on the communication among distinct parts of the brain which is mediated by structural connections (1). Empirical studies in vertebrate and invertebrate brains have demonstrated that a graded and highly heterogeneous strength of structural connections is a basic feature of brain connectomes (25). Moreover, brain networks exhibit a characteristic global topology (Fig. 1) (6, 7), including features such as modularity. These network topology features pertain to brains as diverse as the Drosophila and the human brain and presumably facilitate integration and segregation of the distinct brain regions (6).

Fig. 1 Brain connectomes are characterized by interregional connectivity strength heterogeneity and exhibit characteristic global network topology.

Heterogeneity of interregional connectivity strength is parsimoniously described by the wiring cost and connectional homophily principles (right). Schematic depiction illustrates characteristic global network topology features observed in diverse brain connectomes (left). For simplicity, directionality of connections is omitted.

A parsimonious description of the heterogeneous strength of interregional connections is offered by the principles of wiring cost and homophily of connections (Fig. 1). These principles pertain to a diversity of brain connectomes (815). The wiring cost principle dictates that brain regions, presumably due to material and metabolic cost constraints, are strongly connected with brain regions that are spatially close, whereas connections with spatially distant brain regions, if any, are very weak (8, 13). Physical distance explains, from a statistical standpoint, a large percentage of connectivity strength heterogeneity (13). A more recently illustrated principle explaining strength heterogeneity is the principle of homophily of connections. The principle of homophily of connections is based on observations in social networks (16). In social networks, two persons with many common close friends are likely to become close friends themselves. In an analogous fashion, two brain regions, with many common brain regions that they strongly connect to, are likely to be themselves strongly connected (9, 10, 14, 15). Connectional homophily also explains, from a statistical standpoint, a large percentage of connectivity strength heterogeneity (10). In sum, empirical evidence demonstrates that the wiring cost principle, as reflected in the physical distance between brain regions, and the homophily principle offer a parsimonious description of the heterogeneous strength of structural connections. Wiring principles, that is, the wiring cost and homophily principles currently examined, not only offer a parsimonious description of connectomes but also guide approaches for deciphering the mechanisms underlying these principles. With respect to the homophily principle, tentative mechanistic explanations rely on Hebbian-like plasticity (10, 14). Other views are rooted in functional teleology, that is, the need for functionally similar regions to be interconnected or postulate that the connectional homophily of the brain is a consequence of its spatial embedding (9). With respect to the wiring cost principle, despite a plethora of studies demonstrating that neural systems are configured in such a way to minimize their wiring cost, the mechanisms that result in such an economical configuration remain elusive. A parsimonious suggested mechanism is based on the spatial embedding of the brain and the mere probabilistic fact that it is more likely for a brain region to establish connections with proximal rather than distant regions (17, 18). With respect to the global network topology of brains, generative models adopting a phenomenological approach based on rules that are grounded on empirically deciphered topological regularities can lead to synthetic connectomes with a global network topology that is comparable to the topology of empirical connectomes (15). Other models emphasize the importance of the temporal evolution of networks in shaping such empirically observed network topology properties (19).

To gain neurobiologically plausible insights into the ontogeny of brain connectomes, there is a pressing need to move beyond teleological and phenomenological explanations of brain wiring and adopt parsimonious modeling approaches that are interpretable at the neurobiological and mechanistic level. Moreover, there is a need to use modeling approaches in a comparative context to uncover the potentially universal scope of key components of plausible developmental mechanisms of brain wiring. Here, we use such a modeling approach and directly compare our modeling results with empirical observations in vertebrate and invertebrate brains. First, we use empirical connectomes and offer a comprehensive demonstration of the universality of the connectional homophily and wiring cost principles. Second, we use a computational modeling approach based on stochastic axonal growth and the simulation of heterochronous and spatially ordered neurodevelopmental gradients. Our model reconstructs a large part of the wiring principles and global network topology observed in empirical adult brain connectomes. Our modeling work highlights space and time, that is, spatial embedding and heterochronicity, as key components of a plausible, neurobiologically realistic developmental mechanism that can sculpt the wiring of both vertebrate and invertebrate brains and thus has a potentially universal scope.

RESULTS

Homophily and wiring cost principles

We analyzed Drosophila, mouse, macaque monkey, and human brain connectomes for which quantitative data on the strength of structural connections was available (25, 15). We used the cosine similarity between the wiring of regions A and B with the rest of the brain as a measure of connectional homophily. Both incoming and outgoing connections were considered for the computation of the connectional homophily. For each pair of regions, A and B, we excluded connections between regions A and B in computing the homophily of regions A and B. Cosine similarity ranges from 0 to 1, with (0) 1 denoting connectionally (dis)similar regions, thus exhibiting (low) high homophily. We used the physical distance between brain regions (Euclidean or geodesic) as a wiring cost metric.

We first used the empirical vertebrate and invertebrate brain connectomes for offering a comprehensive view of the principles of wiring cost and homophily. We binned the homophily and distance values, and for each bin, we calculated the average strength of connections (Fig. 2). For all connectomes and pairs of brain regions, increased homophily between two regions was accompanied by increased connectivity strength in-between them. Distance had the opposite effect, that is, increased distance was accompanied by decreased connectivity strength (Fig. 2). For a formal statistical assessment of the contribution of homophily and distance to the strength of connections, we used a least squares regression (see Materials and Methods). For the Drosophila connectome, fitting the homophily and distance factors to the strength of connections resulted in significant contributions of both factors (univariate model: β = 0.61 and −0.41, P < 0.001, model fit R2 = 0.37 and 0.18; multivariate model: β = 0.54 and −0.14, P < 0.001, model fit R2 = 0.39, for homophily and distance, respectively). Significant results were obtained for the rest of the connectomes, namely, for the mouse (univariate model: β = 0.63 and −0.38, P < 0.001, R2 = 0.40 and 0.15; multivariate model: β = 0.58 and –0.09, P < 0.001, R2 = 0.41, for homophily and distance, respectively), macaque monkey (univariate model: β = 0.71 and −0.49, P < 0.001, R2 = 0.51 and 0.24; multivariate model: β = 0.64 and −0.13, P < 0.001, R2 = 0.52, for homophily and distance, respectively), and human (univariate model: β = 0.85 and −0.68, P < 0.001, R2 = 0.71 and 0.47; multivariate model: β = 0.71 and −0.20, P < 0.001, R2 = 0.74, for homophily and distance, respectively). We used the aforementioned empirical model fit values as a measure for the quality of the fit of the synthetic connectomes to the empirical connectomes (see the “Comparing wiring principles in synthetic and empirical connectomes” section).

Fig. 2 Homophily, physical distance, and strength heterogeneity of empirical brain connectomes.

Increased homophily entails increased connectivity strength. Increased physical distance entails decreased connectivity strength. Depicted R2 values are derived from univariate regression. See the “Homophily and wiring cost principles” section for results and R2 values concerning multivariate regression. Note that the correlation, and thus the reported R2 values, between the strength of connections, physical distance, and homophily was computed on the nonbinned data and that the binning was performed for aiding visualization. Drosophila brain drawing from (6).

In sum, in both vertebrate and invertebrate brains, homophily and wiring cost constitute wiring principles that capture a substantial part of the heterogeneity of interregional connectivity strength. Increased physical distance between brain regions dictates decreased strength of connections between the brain regions, and increased homophily dictates increased strength of connections.

Simulation of neurodevelopmental gradients and connectivity formation

Heterochronicity in development was suggested to play an important role in the morphology of animals (20). In addition, the heterochronous neurogenesis of the developing brain seems related to the preferential connectivity patterns between brain regions of the rat and mouse brain (11, 21). Inspired by such theoretical work and empirical evidence, we aimed at simulating neurodevelopmental gradients and connectivity formation in a model based on heterochronicity (see Materials and Methods). Here, we use the term “heterochronous” to denote developmental events (population of the distinct parts of the brain with neurons and their subsequent wiring) that occur at different time points during development. In the current modeling context, we focused on the importance of the heterochronous and spatially ordered gradients of release and accumulation of neurons in the brain (22) and how such phenomena relate to the wiring of the brain. Therefore, detailed neurogenetic events, such as cell cycles, and other neurodevelopmental events, such as apoptosis, were not the focus of the current study and therefore were not modeled and examined (see Material and Methods for the possibility of inclusion of these phenomena into the model). We use the term “spatially ordered” to denote a mode of development where the population of the distinct parts of the brain with neurons occurs in such a way that, if a part of the brain, for instance, region A, is likely to be populated with neurons at time point t, then parts of the brain that are physically proximal to region A are also likely to be populated with neurons at time point t + 1. A two-dimensional (2D) synthetic brain constituted the niche of the simulations, with each unit surface hosting up to N neurons (Fig. 3A). Heterochronous development emanated from root(s)/origin(s), forming orderly, spatially embedded neurodevelopmental gradients, as empirical evidence dictates (Fig. 3A and fig. S1) (22, 23). The origins/roots of the gradients denote the spatial location in the developing brain, where neurons accumulate in early stages of development (22, 23). Subsequent population of parts of the brain with neurons emanates from these origins/roots (Fig. 3A and fig. S1) (22, 23). Each unit surface of the synthetic brain was assigned a time window dictating the probabilities of a neuronal population of n neurons migrating at each unit surface at a given time point t of the simulation (Fig. 3A). This time window was proportional to the physical distance of each unit surface from the origin(s). This modeling choice was informed by empirical evidence, suggesting spatially ordered waves of migration of neurons from neurodevelopmental origins (22, 23). When a neuron migrated to a unit surface, it attempted to form a connection by extending an axon toward a direction that was randomly sampled from a uniform distribution. If the straight line of this axonal extension intersected a circle with a radius r centered around each neuron that has already been placed to the synthetic brain, then a connection was formed (see Materials and Methods) (Fig. 3B). A straight line was used, as in previous modeling approaches (17, 24), in the absence of data capturing the exact trajectory of developing axons in vertebrate and invertebrate brains.

Fig. 3 Developmental modeling approach.

(A) Neurodevelopmental gradients were simulated in a synthetic 2D brain. Each surface unit was characterized by a time window that indicates at a given time point the probability of a neuronal population migrating to each surface unit. Each time window is shaped by the distance of each surface unit from the root(s), that is, origin(s), of the neurodevelopmental gradients. For instance, the surface unit close to the neurodevelopmental origin (petrol green) is more probable to be populated earlier than the surface unit further from the neurodevelopmental origin (magenta). (B) Heterochronous and spatially ordered ontogeny of synthetic connectomes. (C) Heterochronous, spatially random, and (D) tautochronous ontogeny of synthetic connectomes. (E) Creation of the synthetic connectome. Note that time units depicted in (A) and (B) denote “developmental ticks” occurring in a time frame from the onset of the simulation of development and connectivity formation until the end of the simulation (for details, see the “Modeling neurodevelopmental gradients and connectivity formation” section).

For assessing the importance of heterochronicity and the spatial order of development, two additional scenarios were simulated. In one scenario, heterochronous population of the brain with neurons and subsequent connectivity formation was simulated as previously described, but the population of the brain was spatially random (Fig. 3C). In the other scenario, population of the brain with neurons was tautochronous, that is, all neurons were simultaneously placed in the synthetic brain and formed connections (Fig. 3D). In all scenarios, the final synthetic brain contained the same number of neurons; thus, potential differences in model fit could not be attributed to this factor.

Since the empirical connectomes are region-to-region connectivity matrices, the synthetic brain was parcellated into regions with a Voronoi tessellation, resulting in a synthetic region-to-region connectivity matrix (Fig. 3E). The number of regions and connectivity density, that is, the fraction of connections that exist over the total number of connections that can exist given a number of brain regions, of the synthetic connectomes was determined on the basis of the number of regions and connectivity density of the empirical connectomes. We calculated connection strength as the number of projections to region A from region B divided by A’s total number of incoming projections. This normalization renders the connectomes comparable across simulations and empirical data and allows the depiction of comparable ranges of connectivity strength values across the different brain connectomes.

Comparing wiring principles in synthetic and empirical connectomes

The synthetic connectomes obtained from the simulations corresponding to each scenario (heterochronous spatially ordered, heterochronous spatially random, and tautochronous) were subject to a least squares regression analysis to obtain the parameters (β values) that correspond to the statistical association of homophily and physical distance to the strength of connections. We subsequently used these parameters, obtained exclusively from the synthetic connectomes, to estimate the fit to the empirical connectomes by predicting the strength of the connections of the empirical connectomes from the homophily and distance predictors derived from the empirical connectomes.

Model selection

We used Akaike information criterion (AIC), controlling for the complexity and fit of models, for ranking the models. We found that the heterochronous models better relate to the homophily principle, and tautochronous models better relate to the wiring cost principle, when these principles were examined separately. The heterochronous and spatially ordered models were better when considering the two principles simultaneously. Specifically, for the homophily principle, the heterochronous models resulted in the lowest AIC values for all brain connectomes (figs. S2 and S3). The temporal overlap of neurodevelopmental gradients, controlled by parameter a, and the number of origins of the neurodevelopmental gradients did not systematically influence the AIC values (fig. S2); thus, average summaries for each scenario are provided (fig. S3). When considering simultaneously the homophily and wiring cost principles, since both co-characterize adult brain connectomes, the heterochronous and spatially ordered models exhibited the lowest AIC values (figs. S2 and S3).

Model fit and null models

AIC is useful for ranking candidate models, but it does not offer information on the quality of the fit of the models and how different, from a statistical standpoint, this fit is in relation to a null model. Therefore, we used the coefficient of determination (R2) that is based on the residual sum of squares (RSS) of connectivity strength, as predicted from the parameters obtained exclusively from the synthetic connectomes, and the actual, empirical connectivity strength. An example of this empirical and predicted strength of connections is provided in fig. S4. The R2 values obtained from the fit of the parameters from the synthetic connectomes to the empirical connectomes were compared with the R2 values obtained from the parameters and fit when considering only the empirical connectomes. Moreover, we performed statistical null hypothesis testing (F-test) of the fit of the parameters estimated from the synthetic models to the empirical connectome data.

With respect to the homophily principle, we observed the best fit to the empirical connectomes for all brain connectomes for the heterochronous models (median R2 = 0.37, 0.40, 0.50, and 0.71 for the Drosophila, mouse, macaque monkey, and human connectomes, respectively. Note that all the reported R2 values are rounded to the nearest hundredth.) (Fig. 4). The parameter denoting the importance of homophily in predicting the strength of connections obtained from the synthetic models resulted in predictions of empirical connectivity with a fit that was very close to that obtained from a model fitted exclusively to the empirical data (R2 = 0.37, 0.40, 0.51, and 0.71 for the Drosophila, mouse, macaque monkey, and human connectomes, respectively) (Fig. 4). Hence, heterochronously developed and wired brains resulted in 100% (median) of the empirical fit (range, 98 to 100%) that corresponds to the relation between connectivity and homophily in empirical connectomes (Fig. 4).

Fig. 4 Interregional wiring principles captured by the synthetic connectomes.

Summary of R2 values for the predictions of empirical connectivity from parameters calculated exclusively from the synthetic connectomes. Values for R2 are depicted for the (A) Drosophila, (B) mouse, (C) macaque monkey, and (D) human connectomes when considering the homophily or wiring cost principles separately or simultaneously. Note the better fit of the heterochronously generated connectomes when considering simultaneously the wiring cost and homophily principles. Note that the predictions from model parameters estimated exclusively from the synthetic connectomes result in R2 values very close to R2 values obtained when using exclusively the empirical data (gray bars), thus replicating a high percentage of the empirically observed relations between connectivity strength, homophily, and distance across vertebrate and invertebrate brains. Drosophila brain drawing from (6).

With respect to the wiring cost principle, we observed the best fit to the empirical connectomes for all brain connectomes for the tautochronous models (median R2 = 0.15, 0.12, 0.24, and 0.45 for the Drosophila, mouse, macaque monkey, and human connectomes, respectively) (Fig. 4). The parameter denoting the importance of distance in predicting the strength of connections obtained from the synthetic models resulted in a fit to the empirical connectomes that was very close to that obtained from a model fitted exclusively to the empirical connectomes (R2 = 0.18, 0.15, 0.24, and 0.47 for the Drosophila, mouse, macaque monkey, and human connectomes, respectively) (Fig. 4). Hence, tautochronously developed brains resulted in 92% (median) of the empirical fit (range, 80 to 100%) that corresponds to the empirical relation between connectivity and distance (Fig. 4).

When considering the homophily and wiring cost principles simultaneously, since both co-characterize empirical connectomes, the heterochronous and spatially ordered model resulted in the best fit across the different vertebrate and invertebrate brain connectomes (R2 = 0.37, 0.29, 0.41, and 0.65 for the Drosophila, mouse, macaque monkey, and human connectomes, respectively. Note that, for the macaque monkey, heterochronous and spatially random models led to a similar fit with the heterochronous and spatially ordered models, that is, R2 = 0.41) (Fig. 4). Synthetic brains that were developed in a heterochronous and spatially ordered manner resulted in parameters, denoting the importance of wiring cost and homophily principles in predicting the strength of connections, which led to a fit to the empirical connectomes that was very close to that obtained from a model applied exclusively to the empirical connectomes (R2 = 0.39, 0.41, 0.52, and 0.74 for the Drosophila, mouse, macaque monkey, and human connectomes, respectively) (Fig. 4). Therefore, for the simultaneous assessment of the homophily and wiring cost principles, synthetic brains that were developed and wired in a heterochronous and spatially ordered fashion were very close to the fit obtained directly from the empirical connectomes, reaching 83% (median) of the empirical fit (range, 71 to 95%) (Fig. 4). Last, all the predictions of the strength of the empirical connections from the parameters obtained from the synthetic connectomes were statistically significant (P < 0.001, F-test) (table S1).

In sum, despite the fact that our simulations were based solely on heterochronous neurodevelopmental gradients and stochastic connectivity formation, and thus, were completely agnostic to the wiring principles that characterize empirical connectomes, the synthetic connectomes replicated a high percentage of the empirical relation between homophily, distance, and connectivity strength. Specifically, the heterochronous and spatially ordered models exhibited superior fit across the different vertebrate and invertebrate brains when considering both wiring principles simultaneously. Explaining the residual part of the empirical fit not captured by our model (Fig. 4) may involve additional mechanisms, such as pruning, axonal guidance, or activity-dependent plasticity.

Ontogeny of global network topology

After the examination of interregional wiring principles, we next sought to investigate to what extent the synthetic connectomes generated under the different scenarios reflected differences in global network topology. We focused on widely and commonly used global network topology metrics, that is, clustering coefficient, characteristic path length, diffusion efficiency, physical distance spanned by the connections, maximum clique size (defining the core of the network), and modularity (6, 7, 13, 25). First, for each brain connectome separately, we constructed a network topology profile for each of the synthetic connectomes that were generated under the three different scenarios, the latter defining three different groups. We performed a canonical discriminant analysis (CDA) to investigate whether the synthetic connectomes were distinct in terms of their network topology profile. The components obtained from the CDA defined a network topology morphospace (26) that separated the synthetic connectomes generated from the three different scenarios for all brain connectomes (Wilks’ Λ = 0.26, 0.21, 0.13, and 0.46; all P < 0.001, permutation test; for the Drosophila, mouse, macaque monkey, and human connectomes, respectively) (Fig. 5A). Nearly all between-group variance that separated the heterochronously from the tautochronously generated synthetic connectomes was explained by the first component (82, 95, 98, and 91% for the Drosophila, mouse, macaque monkey, and human connectomes, respectively) (Fig. 5A). Overall, heterochronously generated connectomes exhibited network topology metrics that differentiated them from the tautochronously generated connectomes, with the maximum clique size (defining the core) and characteristic path length (average shortest path length) among the most prominent network metrics segregating the synthetic connectomes in the network topology morphospace (Fig. 5A).

Fig. 5 Network topology morphospace of synthetic connectomes and classification of empirical connectomes.

(A) Heterochronously and tautochronously generated connectomes are separated on the basis of their global network properties across the first component. Core size (max clique size) and characteristic path length (average shortest path length) are the major network metrics differentiating the groups. (B) Posterior probabilities indicate that the global network topology signature of all empirical connectomes is more reminiscent of connectomes with a heterochronous ontogeny. Drosophila brain drawing from (6).

Next, we examined whether the global network topology of the empirical connectomes mostly resembled the topology of heterochronously or tautochronously generated synthetic connectomes. To this end, we situated the empirical connectomes in the network topology morphospace defined by the heterochronously and tautochronously generated connectomes. This was achieved by performing the CDA 100 times, each time with 70% of the synthetic connectomes, and using the derived CDA parameters to calculate the posterior probabilities of the empirical connectome belonging to one of the three groups. This procedure was separately repeated for each brain connectome. We obtained the higher posterior probabilities for the connectomes that were developed in a heterochronous and spatially ordered fashion (Fig. 5B).

Empirical connectomes mostly resemble the heterochronous and spatially ordered synthetic connectomes and are maximally dissimilar to random connectomes (fig. S5). Moreover, the synthetic connectomes generated in a heterochronous and spatially ordered fashion exhibited the lower AIC values, indicating that the heterochronous spatially ordered model should be selected as the model best approximating the global network topology of adult brain connectomes when the different number of parameters of each model is taken into account (fig. S6). However, empirical connectomes are also distinct from the synthetic connectomes (fig. S5), indicating that certain empirical network metrics are not accurately captured by the synthetic connectomes. We assessed how faithfully the synthetic connectomes reconstructed the empirical global network topology by plotting each network metric separately for the synthetic and empirical connectomes. Overall, good correspondence with empirical values was observed, with certain deviations, either underestimating or overestimating the empirical network topology (median percentage of deviations from 100% corresponding to a perfect reconstruction; thus, lower values denote better correspondence with the empirical network metrics; 9, 34, 32, and 6% for the Drosophila, mouse, macaque monkey, and human connectomes, respectively) (Fig. 6). Therefore, although the global network topology of the empirical connectomes is more reminiscent of connectomes that were developed in a heterochronous and spatially ordered fashion, deviations indicate that a more accurate reconstruction of certain aspects of the empirical global network topology may involve additional mechanisms, such as pruning, axonal guidance, or activity-dependent plasticity.

Fig. 6 Global network topology metrics for the synthetic and empirical connectomes.

The percentages below each bar denote the percentage of the values obtained for each network metric in relation to the empirical values. Thus, values over or below 100% denote overestimation or underestimation of the network metric, respectively. Drosophila brain drawing from (6).

DISCUSSION

Heterochronicity in development has been suggested as a key underlying principle of the observed morphological diversity in animals (20) and seems to underlie the intricate connectional configuration of rodent brains (11, 21). Here, we brought the concept of heterochronicity to the connectional realm, within a comparative framework that encompasses a diverse spectrum of brains, ranging from fly to human brains. We have demonstrated that heterochronous and spatially ordered developmental gradients and connectivity formation constitute key components of a universal and neurobiologically realistic mechanism that may underlie the empirically observed interregional wiring principles and global network topology of adult brains.

Our study offers mechanistic insights into the development of brain connectomes and thus extends beyond phenomenological modeling approaches that aim to offer a parsimonious description of brain connectomes (27). Our empirical results demonstrate that the homophily and wiring cost principles characterize the adult connectional configuration of not only the mammalian brain (9, 10, 14) but also the insect brain. Our modeling work demonstrates that such a wiring configuration can be laid down by the presence of spatially ordered and heterochronous neurodevelopmental events and a purely stochastic axonal growth, without gradients of axonal guidance molecules or activity-dependent remodeling of connections. Therefore, our results offer a solution to the conundrum of how two regions with high homophily, that is, similar patterns of connections, and thus, functional similarity end up preferentially and strongly interconnected in the adult brain. Heterochronous and spatially ordered development renders the emergence of such wiring configuration feasible, since two regions that are populated by neurons at proximal time points in development have similar available connectional targets with other regions, also occupying similar time points in development, leading to high connectional homophily of the two regions and the presence of strong connections in-between the two regions. In other words, we show how a structural skeleton that allows connectionally and functionally similar regions to be interconnected can be laid down with a parsimonious neurodevelopmental mechanism. This connectional configuration may be further refined through subsequent activity-dependent plasticity.

Our modeling work also offers a parsimonious mechanism for the wiring cost principle that pertains to diverse connectomes (6). The spatial embedding of the brain, and the mere probabilistic fact that the establishment of connections is more likely to take place with proximal rather than distant targets (18), is sufficient to replicate a large part of the empirically observed relation between physical distance and strength of connections from the fly to the human brain connectome. Thus, our results extend previous modeling and theoretical work (17, 18) by demonstrating a basic, plausible, and universal mechanism, leading to economically wired brains.

Brain connectomes exhibit nonrandom global topological properties such as modularity (6). Connectomes that have a heterochronous ontogeny exhibit a global network topology different from connectomes with a tautochronous ontogeny and are distinguishable from each other in a global network topology morphospace. The global network configuration of Drosophila, mouse, macaque monkey, and human brains is mostly reminiscent of synthetic connectomes with a heterochronous and spatially ordered ontogeny. Thus, heterochronicity may shape the empirically observed global network topology of brains. Deviant connectome topologies are observed in a number of pathologies (28). At the molecular level, the Notch and bone morphogenetic protein signaling pathways are involved in the initiation of neurogenesis (29) and thus influence the timing of the subsequent migration of neurons to the cortical plate. Our results show that pathologies involving these pathways and, consequently, the normative temporal population of the brain with neurons might also affect the global large-scale connectional configuration of the brain.

Our model adopted a parsimonious approach to connectivity formation and replicated a large percentage of empirically observed interregional wiring principles and global network topology features of diverse brains. Our results emphasize that, apart from spatial constraints, the temporal unfolding of neurodevelopmental events is crucial for the normative wiring configuration of brains. However, these results do not entail that other aspects, such as principles of chemotaxis, that is, chemoattraction and chemorepulsion for guiding axons to their proper targets (30), are irrelevant. These principles, in addition to activity-dependent plasticity rules, can be incorporated in future extensions of the model to reconstruct the currently residual part of the empirically observed interregional wiring principles and global network topology features of brain connectomes. Moreover, while currently we examined the wiring configuration of brains at the macroscale (interregional) level due to the existence of comprehensive data on multiple species, the progressive mapping of the wiring of the brain at the microscale (cell-to-cell connectivity) level will provide the necessary data for assessing if the currently illustrated mechanism can also explain more fine-grained wiring configurations at the synaptic level. Last, while we have demonstrated common principles and mechanisms across brain connectomes of different species, future studies should also aim at examining species-specific aspects of brain wiring.

In sum, we have demonstrated that heterochronous and spatially ordered development is a plausible key universal phenomenon that can sculpt the intricate wiring configuration of vertebrate and invertebrate brains. Hence, proper understanding of the ontogeny of brain connectomes should consider not only the spatial but also the temporal dimension of brains.

MATERIALS AND METHODS

Connectome datasets

For consistency across brain connectomes, all the connectomes described below concern intrahemispheric connections. Data are freely available from the indicated online resources.

Drosophila. The connectivity matrix for the Drosophila brain was reconstructed from the FlyCircuit 1.1 database (www.flycircuit.tw), a repository of images of 12,995 projections in the female Drosophila brain (4). Neurons were labeled with green fluorescent protein (GFP) using genetic mosaic analysis with a repressible cell marker. GFP-labeled neurons were delineated from whole-brain 3D images and coregistered to a female template brain using a rigid linear transformation. Individual neurons were partitioned into 49 local processing units (“brain regions”) with distinct morphological and functional characteristics. The result was a weighted and directed connectome composed of 1950 connections among the 49 local processing units. The Euclidean distances between the local processing units were used as a proxy for wiring cost. This dataset has been previously analyzed in, e.g., (31, 32).

Mouse. The mouse connectome was reconstructed on the basis of the freely available tract-tracing data from the Allen Institute Mouse Brain Connectivity Atlas (http://connectivity.brain-map.org). Anterograde recombinant adeno-associated viral tracer was injected into target areas in the right hemisphere of mouse brains. The brain was extracted 3 weeks after injection, and viral tracer projection patterns were reconstructed. Reconstructions were then smoothed and aligned to the common coordinate space of the Allen Reference Atlas. Details for the dataset are provided in (2). A denoised connectivity matrix of the mouse brain was used (5). Briefly, in total, 56 areas of the right hemisphere of the mouse brain were used, and the connections in-between these areas were weighted on the basis of the normalized connection densities, that is, the number of connections from a unit volume of a source area to a unit volume of a target area. We used the geodesic distances between brain regions as a proxy for wiring cost. See (5) for details. Thus, the mouse connectome is a weighted and directed connectome with 1975 connections among 56 areas.

Macaque monkey. The macaque monkey connectome was derived from http://core-nets.org, which is a database compiled from results obtained from invasive tract-tracing studies in the macaque brain (3). The injections involved the retrograde tracers fast blue and diamidino yellow. Quantitative information on cortico-cortical connectivity, i.e., strength and physical length, for a set of 29 regions spanning the frontal, parietal, temporal, and occipital cortex, was represented as a 29 × 29 weighted and directed connectome. The connectivity strength between regions A and B after an injection in region A was expressed as the fraction of labeled neurons involving region B, i.e., the number of projection neurons located in region B, divided by the total number of neurons observed across the whole brain after injecting region A (excluding “intrinsic” neurons, i.e., projection neurons located within region A). Thus, the macaque monkey connectome is a weighted and directed connectome with 536 connections among 29 areas. Physical distance between areas was estimated as the geodesic distance between the 3D coordinates of the barycenters of the cortical regions, thus approximating the trajectory of axons linking two regions. See (3) for details.

Human. Human brain connectomes were reconstructed from diffusion-weighted magnetic resonance imaging and a deterministic tractography algorithm. The network that we analyzed was a group-representative composite of subject-level networks (30 subjects). This network construction process entailed acquiring diffusion spectrum and T1-weighted anatomical images for each individual. Diffusion spectrum imaging (DSI) scans sampled 257 directions using a Q5 half-shell acquisition with a maximum b value of 5000, an isotropic voxel size of 2.4 mm, and an axial acquisition with repetition time (TR) of 5 s, echo time (TE) of 138 ms, 52 slices, and field of view of [231, 231, 125] mm. All procedures were approved by the Institutional Review Board of the University of Pennsylvania, and all participants gave informed consent. DSI data were reconstructed in DSI Studio (http://dsi-studio.labsolver.org/) using q-space diffeomorphic reconstruction (QSDR) (33). QSDR reconstructs diffusion-weighted images in native space, computes the quantitative anisotropy (QA) of each voxel, warps the image to a template QA volume in Montreal Neurological Institute (MNI) space using the statistical parametric mapping nonlinear registration algorithm, and reconstructs spin-density functions with a mean diffusion distance of 1.25 mm with three fiber orientations per voxel. Fiber tracking was performed using a modified FACT algorithm with an angular cutoff of 55°, a step size of 1.0 mm, a minimum length of 10 mm, a maximum length of 400 mm, and a QA threshold determined by diffusion-weighted magnetic resonance imaging signal in the colony-stimulating factor. For each individual, the algorithm terminated when 1,000,000 streamlines were reconstructed.

In parallel, T1 anatomical scans were segmented using FreeSurfer and parcellated using the Connectome Mapping Toolkit (www.connectomics.org), according to an atlas with 34 cortical regions (34). Each parcellation was registered to the B0 volume of subjects’ DSI data, and a B0-to-MNI voxel mapping was used to map area labels from native space to MNI coordinates. Streamlines were aggregated by the regions in which their starting and termination points were located. The connection weight between any pair of regions was defined as their streamline count normalized by the geometric means of their volumes. Group-representative matrices were generated using a consistency-based thresholding procedure (35). This approach has been described elsewhere and shown to be superior to alternative thresholding procedures (35). This procedure was applied separately to inter- and intrahemispheric connections. The resulting group average human connectome contained 337 undirected and weighted connections between 34 brain regions. This corresponds to a density threshold of 0.6. Note that only intrahemispheric connections were further analyzed.

Physical distance between regions, serving as a proxy for wiring cost, was calculated as the Euclidean distance between the barycenters of the brain regions. This human connectome dataset has been previously analyzed in (32).

Modeling neurodevelopmental gradients and connectivity formation

Populating the synthetic brain with neurons. Pillars of the current computational modeling approach are the studies of (19, 36). The brain was modeled as a 2D square, reminiscent of previous modeling work (24, 36). The size of the synthetic brain was set to 50 × 50 units for all the different brain connectomes (Fig. 3A). Variations of the exact size of the synthetic brain did not substantially influence the results. Each unit surface in the synthetic brain represents a structure of the developing brain, for instance, the cortical plate of mammalian cortices. Each unit surface can host N neurons. The population of the synthetic brain with neurons is achieved as follows. In the scenarios realizing a heterochronous development, a distinct time window was assigned to each surface unit. This time window dictates the probability of n neurons migrating at the corresponding unit surface at time point t and was defined as P(i,t,a)=(t2λ(tλ1)2)1m(i)α with m(i)=1k+1, λ=log(2)log(m), i denotes a unique integer identifying a time window assigned to a unit surface, k denotes the total number of time windows considered in the simulation, and α constitutes a parameter controlling the width of the time window of each unit surface of the synthetic brain (α > 0). Values for the parameter α closer to 0 lead to more narrow time windows and, hence, less overlap between the set of time windows of the surface units. Larger values for α lead to broader time windows, hence leading to more overlap between the set of time windows of the surface units (19). Here, we explored a range of values for parameter α, ranging from 0.2 to 0.8 with an increment of 0.2.

For each time window, the function P(i,t,a) was scaled to max[P(i,t,a)] for t = 0, …, 1 so that the outcome of the function is bound to [0 1] and can, thus, serve as a probability. The temporal resolution of the simulations was set to 0.05, where one developmental “tick” marked a temporal increment of 0.05. Thus, the simulations unfolded across a time frame defined by t = 0, 0.05, 0.1, …, 1. For each unit surface, the distance from the neurodevelopmental root(s)/origin(s) was computed. The total number of time windows was defined on the basis of the unique number of minimum distances from the neurodevelopmental root(s). The unique minimum distances of each unit surface served as a criterion for assigning to them a time window i. Time windows were assigned in such a way so that surface units with closer distances to the neurodevelopmental root(s) were assigned to time windows with the lowest corresponding unique integer i (i = 1, 2, …, k). Thus, surface units close to the neurodevelopmental root(s) were more likely to be populated with neurons at early stages of the simulations, that is, at low t values (Fig. 3A). The total number of neurons to populate the synthetic brain was increased exponentially, as empirical evidence dictates (37). Thus, the number of neurons was increased on the basis of the following relation: nt = ninit(1 + r)t, where ninit is the initial number of neurons, here set to 100, r is the rate of growth, here set to 0.2, and t is the current time point of the simulation. We should note that the exact phenomena of neurogenesis, such as number of divisions and apoptosis, were not within the scope of the current investigation and were not modeled. Instead, neuronal populations were positioned in the synthetic brain as described above, corresponding to the release and accumulation of neuronal populations across the developing brain (22). Modeling phenomena such as cell proliferation and apoptosis is feasible (38, 39), and thus, this future extension of our modeling approach is plausible.

In total, three scenarios were simulated, that is, heterochronous and spatially ordered, heterochronous and spatially random, and tautochronous development. In the first scenario, the population of the brain with neurons proceeded in a spatially ordered fashion as co-centric rings from the root(s) (Fig. 3B). In the second scenario, the population of the brain took place in a heterochronous fashion but in a spatially random manner. This was achieved by randomizing the relations between the position of the unit surfaces from the root(s) and their assigned time window. This randomization destroyed the spatially ordered co-centric population of the brain (Fig. 3C). Last, a tautochronous scenario consisted of the simultaneous placement of the neurons in the synthetic brain and was followed by subsequent connectivity formation. Note that, in all scenarios, the total number of neurons of the final “adult” synthetic brain contained the exact same number of neurons.

Connectivity formation. Once a neuron was placed in the synthetic brain, it could form a connection. Neurons eligible for establishing connections, that is, neurons that have not already established a connection, were selected uniformly and at random. Connectivity formation was established in a stochastic manner. A direction across the 2D space was randomly sampled from a uniform distribution, and if a straight line from the neuron that was about to develop a connection toward that direction was intersecting a circle placed around each of the other neurons in the synthetic brain, then a connection was established. A radius equal to the surface unit was used to define the circle dictating if a connection would be established. Each neuron could form one connection, corresponding to the fact that projection neurons have one axon. Axon collaterals were not modeled. Each neuron had a maximum capacity of constituting the target of up to 100 connections. If a neuron did not form a connection during a developmental time point, then it could retry establishing a connection during the next time point.

Constructing the synthetic connectomes. The empirical connectomes are interregional connectomes. Therefore, for a proper comparison of the synthetic and empirical connectomes, the following procedure was followed. First, the synthetic brain was parcellated with a Voronoi tessellation to “brain regions.” The number of the regions was dictated by the number of regions of each empirical brain connectome (Fig. 3E). Ten different parcellations were applied to each synthetic brain to avoid parcellation-dependent results. Second, a synthetic, directed and weighted connectome was assembled by creating a region-to-region synthetic connectivity matrix. A connection from region A to region B was formed by summing the number of neurons contained in region A that form connections with those in region B. The connectivity density of the synthetic connectomes, that is, the number of connections present in a connectome out of all possible connections given a number of regions, was equated to the connectivity density of each empirical brain connectome. Hence, the synthetic connectomes had the same number of regions and connections as the empirical connectomes. For the synthetic human connectomes, since the human empirical connectome is undirected, connectomes were symmetrized by computing the average of the weight of the connection from regions A to B and the weight of the connection from regions B to A.

Comparing empirical and synthetic connectomes

Relating homophily, physical distance, and connectivity. For assessing the relation of the homophily and physical distance of brain areas with the strength of connections, we used the least squares regression. For rendering the connectomes comparable across simulations and empirical data and across the different species, the connection strengths were normalized by calculating the fraction of projections to region A from region B, by dividing the number of projections to A from B with the total number of projections to A. This normalization was performed for the empirical data and the synthetic connectomes.

Model selection. For the synthetic connectomes, the following procedure was performed for each of the connectomes separately. For each simulated brain (×10) and parcellation (×10), that is, in total of 100 times, the relation of strength of connections, homophily, and distance was assessed with linear regression. The obtained regression coefficients (β values), calculated exclusively from the synthetic connectomes, denoting the relation of each of the regressors (homophily and distance) to the strength of connections, were used to predict the empirical connectivity strengths. Therefore, the empirical regressors (homophily and distance) were multiplied with the corresponding coefficients obtained from the regression applied only to the synthetic connectomes. In this way, the RSS, given the model parameters (coefficients estimated solely from the synthetic connectomes), could be estimated.

Subsequently, for comparing the models (heterochronous spatially ordered, heterochronous spatially random, and tautochronous), the AIC was used. AIC is a metric for model comparison that considers the fit to the empirical data (in our case, connectivity strength) and complexity of the model. AIC was calculated as AIC = 2 k + nlog(RSS), where k is the number of model parameters, n is the number of observations, RSS is the residual sum of squares of the actual and predicted strength of connections, and log is the natural logarithm. For generating the heterochronously and spatially ordered developed brains, as opposed to the tautochronously developed brains, the given time point must be known. Therefore, for computing the AIC, the heterochronously spatially ordered model was assigned one additional parameter than the tautochronous model. For generating the heterochronously spatially ordered brains, as opposed to the heterochronously spatially random brains, the distance from the root of the neurogenetic gradients must be known. Therefore, for computing the AIC, the heterochronously spatially ordered model was assigned one additional parameter than the heterochronously spatially random model. The AIC values for the 100 simulations for combinations of the parameters corresponding to the number of the roots of the neurogenetic gradients and the parameter α controlling the overlap of time windows are depicted in fig. S2. No systematic minimization of the AIC values was observed. Therefore, summaries across the number of neurogenetic gradients (roots) and α values were presented.

Model fit and null hypothesis. AIC is useful for model selection, but it does not offer an indication of the quality of the fit of the models to the empirical data nor does it offer the possibility for a statistical comparison with a null model. Therefore, the following steps were taken. The coefficient of determination (R2) was used to assess the fit of the model based on the homophily, wiring cost principle, or both. The R2 is defined as R2 = 1 − RSS/TSS, with TSS denoting the total sum of squares. The R2 was first computed only with the empirical data. This R2 served as a measurement for the quality of the fit of the synthetic to empirical connectomes. The synthetic to empirical connectomes fit was assessed by the R2, but now computed on the basis of the parameters obtained exclusively from the synthetic connectomes and subsequently applied to the empirical connectomes for obtaining predictions of the strength of the empirical connections. These predictions of the strength of connections were used for calculating the R2. The synthetic to empirical connectome fit was tested for statistical significance against a null model (a model including only a constant term and no additional regressors) with the F-test. Thus, the synthetic to empirical connectome fit not only was tested against a null model, which answers the question whether the model obtained from the synthetic connectomes offers a better fit than null models, but also compared against the fit obtained exclusively from the empirical connectomes, hence offering a rigid comparison for the faithfulness of the reconstruction of the empirically observed relations between strength of connections, homophily, and physical distance by the models derived from the synthetic connectomes.

Global network topology

Global network topology metrics. We estimated global network topology metrics that have been extensively used (6, 7, 13, 25) in the synthetic and empirical connectomes. Specifically, we estimated six metrics, that is, the clustering coefficient, characteristic path length, diffusion efficiency, physical distance spanned by the connections, maximum clique size (defining the network core), and modularity. All metrics, apart from the maximum clique size, were estimated with the Brain Connectivity Toolbox (https://sites.google.com/site/bctnet/) (40). The maximum clique size was computed with a MATLAB implementation of the Bron-Kerbosch algorithm with pivoting and degeneracy ordering (https://de.mathworks.com/matlabcentral/fileexchange/47524-find-maximal-cliques-for-large—sparsenetwork). All the metrics, except the core and the physical distance spanned from the connections, were computed by considering the weights of the connectomes. The core was defined as the union of all cliques with the maximum size in the network, irrespective of the weights of the connections (13).

Global network topology morphospace. The aforementioned topological metrics were used to construct the global network topology profile of both the synthetic and empirical connectomes. Subsequently, for assessing the distinct global network properties of the synthetic connectomes constructed under the three simulation scenarios (heterochronous spatially ordered, heterochronous spatially random, and tautochronous), a global network topology morphospace was constructed. To this end, we used CDA as implemented in the MATLAB Fathom Toolbox (https://de.mathworks.com/matlabcentral/fileexchange/68518-fathom-toolbox-for-matlab).

CDA aims at projecting multivariate data belonging to K groups to a low-dimensional space. CDA estimates a linear combination of the features that maximizes between-group variation in a similar way that the principal components analysis aims at constructing a coordinate system that maximizes total variation across all observations. Here, each observation was a synthetic connectome and the features corresponded to the six global network metrics. The three different simulation scenarios were used to assign the synthetic connectomes to three groups. Statistical significance of the separability of the groups in the CDA space, defining the global network topology morphospace, was assessed with Wilks’ Λ. The statistical significance of Wilks’ Λ was assessed with permutation tests, that is, the group assignments were shuffled, and the CDA and Wilks’ Λ were estimated with the random group labels. The procedure was repeated 1000 times, and a null distribution of the Wilks’ Λ values was obtained, allowing the estimation of a P value.

For uncovering the group that the empirical connectomes mostly resemble to, the CDA space defining the global network topology morphospace was used for classifying the empirical network topology profiles. For each brain connectome, the CDA was performed as described before, but now, 70% of the synthetic connectomes of each group were used. Subsequently, posterior probabilities of the empirical connectome belonging to each group were obtained. The procedure was repeated 100 times, resulting in a distribution of posterior probabilities for each group (heterochronous spatially ordered, heterochronous spatially random, and tautochronous). To take into account the different parameters used to create the synthetic connectomes and the similarity of their global network topology to the topology of the empirical connectomes, AIC values were computed, but in this case, the sum of the Euclidean distances of the empirical and synthetic connectomes in the 2D space defined by the CDA were used as a similarity measure, with lower Euclidean distances in this 2D space, indicating higher similarity/fit, in an analogous fashion that low RSS values indicate better model predictions/fit.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/6/eaav9694/DC1

Fig. S1. Neurodevelopmental gradients in the developing empirical and synthetic brain.

Fig. S2. AIC values for the different simulation scenarios.

Fig. S3. AIC values for comparison of the models.

Fig. S4. Examples of actual and predicted strength of connections for the different models across brain connectomes.

Fig. S5. Examples of the projections of the empirical brain connectomes to the morphospace defined by the global network topology of the synthetic brain connectomes.

Fig. S6. AIC values for the similarity of the global network topology of empirical and synthetic connectomes.

Table S1. F-test values for the predictions of the strength of connections of the empirical connectomes from parameters estimated exclusively from the synthetic connectomes.

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

REFERENCES AND NOTES

Acknowledgments: We would like to thank P. Vértes for helpful comments on previous drafts of the manuscript. Funding: This work was supported by the Humboldt Research Fellowship from the Alexander von Humboldt Foundation (to A.G.), and grants from DGF SFB 936/A1, Z3, SPP 2041 Computational Connectomics/HI 1286/7-1, and TRR 169/A2 (to C.C.H.) are gratefully acknowledged. Author contributions: A.G. designed the research; A.G., R.F.B., and C.C.H. performed research; A.G. and R.F.B. contributed new reagents/analytic tools; A.G. analyzed the data and performed the simulations; A.G. drafted the manuscript; A.G., R.F.B., and C.C.H. revised and finalized the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All nonhuman connectome data are freely available from the indicated online resources. Human connectome data may be requested from the authors. Additional data related to this paper may be requested from the authors.
View Abstract

Navigate This Article