A generative model of the connectome with dynamic axon growth

The network of axonal connections comprising a nervous system is known as the connectome (Hagmann, 2005; Sporns et al., 2005). Connectomes display nonrandom topological characteristics, such as small-worldness (SW) and modularity (Bassett & Bullmore, 2017; Sporns & Betzel, 2016) as well as rich diversity in the strength of connections and regions (Buzsáki & Mizuseki, 2014). While connectome topological properties are well characterized, the underlying wiring principles that give rise to these properties are poorly understood.

Generative models offer one avenue to investigate principles governing connectome development. Connectome-like networks have been generated in silico to model the micro-, meso-, and macroscale neural connectivity of many organisms, including C. Elegans, Drosophila, nonhuman mammals, and humans, through a variety of spatial, topological, and physiological wiring rules (Akarca et al., 2023; Betzel et al., 2016; Beul et al., 2018; Ercsey-Ravasz et al., 2013; Faskowitz et al., 2018; Henriksen et al., 2016; Kaiser & Hilgetag, 2004; Klimm et al., 2014; Oldham et al., 2022; Pavlovic et al., 2014; Priebe et al., 2017; Simpson et al., 2011; Vértes et al., 2012). In each of these models, the extent to which the generative process is guided by each wiring rule is determined by a set of tunable parameters. Typically, connections are more likely to be generated between regions that are close in spatial proximity to each other (Ercsey-Ravasz et al., 2013; Kaiser & Hilgetag, 2004) and/or for which inclusion of the proposed connection would enhance a desired topological criterion (Simpson et al., 2011; Vértes et al., 2012).

Axonal growth and guidance are important mechanisms that shape brain wiring. Since Ramon y Cajal’s discovery of growth cones and Sperry’s pioneering chemoaffinity hypothesis (Cajal, 1890; Chilton, 2006; Sperry, 1963; Zang et al., 2021), a variety of guidance molecules, such as netrins and slits, were found to contribute to axon pathfinding (Brose et al., 1999; Kennedy et al., 1994; Kidd et al., 1999; Serafini et al., 1994). Miswired connectomes are evident in model organisms deficient in guidance molecules and receptors; for example, abnormal optic chiasm development has been found in slit-deficient mice (Dickson, 2002). Given the importance of axonal guidance in brain network development and wiring, the modeling of pathfinding mechanisms may lead to improved connectome generative models that reflect multiple spatial phenomena, compared with preferential generation of connections between pairs of regions in close spatial proximity. Connectome generative models that consider axonal guidance may provide insight into connectome development, complementing the insight provided by current models and shedding light on the mechanisms that generate the characteristic geometry and spatial architecture of axonal fiber bundles.

Explicitly simulating axonal growth also provides an opportunity to generate weighted brain networks. Most established connectome generative models are unweighted—a connection is either present or absent between a pair of regions. As such, information about diverse connectivity strengths is overlooked and not modeled. Recent studies have proposed various methods to address this issue, such as through connectome community (Faskowitz et al., 2018) and communicability redundancy (Akarca et al., 2023). Despite the unique strengths of these approaches, axon counts remain a natural and straightforward representation for the strengths of physical neural connectivity. In this direction, researchers have established connectome generative models that simulate networks by growing axons in predetermined directions (Song et al., 2014). In contrast to tuning weights of connections themselves, this work parcellates a continuous space into discrete regions that are connected by multiple simulated axons. As a result, connection weights naturally arise from the axon counts between pairs of regions. Although a nodal correspondence between generated and empirical connectomes is missing, this approach has the advantage of being biologically tractable. The generated networks are found to replicate many topological properties of empirical connectomes, including degree, clustering, and triad distributions (Song et al., 2014). However, it remains unclear whether these attributes persist or new topological characteristics arise in the presence of dynamic axon guidance.

In this study, we establish a new spatially embedded generative model for weighted connectomes. Our model significantly builds on the seminal models of Kaiser et al. (2009) as well as Song et al. (2014), both of which feature axon outgrowth. Dynamic axon growth is a key novelty of our model, without which curved axons cannot form. Each brain region exerts a distance-dependent attractive force on an extending axon’s tip, steering the direction of axon growth. This emulates the process through which axon growth cones react to molecular guiding cues, whose concentration decays with the distance to chemical release sites. We find that our model can recapitulate a diverse array of topological features characteristic of nervous systems, at the edge, node, and network levels. We fit the two parameters of our generative model to individual connectomes, generating weighted networks that reflect interindividual variations in the brain network architecture. Overall, our work enables generation of connectomes in silico that are weighted, are spatially embedded, and feature axonal trajectories that appear biologically realistic.

We develop a model that generates weighted connectomes through dynamic axon outgrowth, an extension of the static generative model proposed by Song et al. (2014). Specifically, the direction of axon growth in a static model is governed in a one-shot manner, where fibers are generated in a direct, linear trajectory toward their targets. In contrast, axons in a dynamic model continuously assess the guidance gradient acting on them and respond accordingly. As such, dynamic models have an internal quasitemporal structure such that the position of the axon in the past influences its subsequent position over iterations. As a result, axons are generated with a richer, fascicle-like geometry.

For model simplicity and axon visualization purposes, we simplify the cerebral volume to a two-dimensional circular construct of radius R. The circumference and the internal space of the circle represent brain gray and white matter, respectively (Figure 1A). To discretize the circle into regions, we uniformly divide the circumference into N_n segments of equal length, where each segment represents a distinct brain region, otherwise known as a node. Given that brain regions differ in volume and surface area, geometric heterogeneity of nodes is introduced through a parameter ρ that perturbs node center coordinates (Figure 1B; detailed in the Methods section). Node regions are defined using a 1D Voronoi tessellation approach along the circumference, such that each point on the circle perimeter is assigned to the region of its nearest node center (Figure 1C). Note that perturbed node centers are not necessarily the geometric centers of the defined segments; however, we found that substituting perturbed node centers with geometric centers did not qualitatively impact our findings (Supporting Information Figure S1).

A distinguishing feature of the model is the dynamic axonal growth, echoing in vitro and in vivo evidence suggesting that axons actively modify growth pathways in response to local molecular and mechanical cues (Dickson, 2002; Oliveri & Goriely, 2022). The simulated axons grow progressively in a step-by-step manner, based on the attractive forces described by $F⇀s⇀$ (Figure 1D). For each step, an axon growth cone situated at position $s⇀$_i is extended in the direction of $F⇀s⇀i$ for a constant distance L_s. This step length L_s is the second tunable parameter of the model. We additionally impose weak regularity constraints on growth (see the Methods section and Supporting Information Figure S2) to avoid trajectories with biologically unrealistic sharp turns (Katz, 1985).

Axon propagation terminates as soon as its growth cone intersects with the circle, forming a connection between two points on the circle’s circumference (Figure 1E). Axons that successfully connect two points of the circumference are selected for network construction. Some axons fail to navigate to a point on the circumference, and they are excluded from subsequent analyses (Methods and Supporting Information Figure S3).

We construct networks by assigning the endpoints of successful axons to their nearest node centers (Figure 1F). The connectivity weight between a pair of nodes is given by the total number of axons linking them. Despite the intrinsic directionality of simulated axons, for simplicity, we focused on mapping weighted, undirected networks. Within-region connections are ignored. As shown in Supporting Information Figure S4, the number of excluded self-connections and unsuccessful axons is negligible for parameter combinations generating brain-like networks.

In this study, unless otherwise specified, we matched N_n to the number of brain regions in the Desikan-Killiany atlas (84 nodes) and simulated 2 × 10⁵ axons in each network. It should be noted that nodes in our generated and empirical connectomes do not have a one-to-one correspondence, an intrinsic limitation from using simplified brain geometry (Song et al., 2014). A complete summary of default parameters used can be found in Methods. We started by investigating how simulated networks behaved in response to variations in model parameters. Our model is governed by two key parameters (explained in the Parameter Specification section and Supporting Information Figure S5 in the Supporting Information): β—the distance-dependent decay of the attractive force, and L_s—the length of each extending step. Note that β determines the relative contribution of guiding cues exerted by each node, such that a larger β emphasizes the guidance from local, adjacent nodes, relative to distant nodes, whereas L_s governs the extent to which an axon can change its trajectory per unit length.

We first evaluated the generated networks in terms of the connectivity weight and nodal degree distributions. As spatially embedded networks, connectomes exhibit strong and abundant short-range connections and weak and rare long-range connections across a variety of scales and species. This property is typically modeled using an exponential distance rule (EDR) of connection weights (Betzel & Bassett, 2018; Ercsey-Ravasz et al., 2013; Gămănuţ et al., 2018; Horvát et al., 2016; Markov et al., 2014; Oh et al., 2014; Rubinov et al., 2015; Scannell et al., 1995). However, a small proportion of strong long-range connections that deviate from the expectations of the EDR model are consistently observed in empirical brain networks (Deco et al., 2021; Roberts et al., 2016). In addition, connectivity weights are typically lognormally distributed (Ercsey-Ravasz et al., 2013; Gămănuţ et al., 2018; Song et al., 2005; Wang et al., 2012), and nodal degrees are characterized by a scale-free distribution (Broido & Clauset, 2019; Eguíluz et al., 2005; Gastner & Ódor, 2016; Giacopelli et al., 2020; Sporns et al., 2004; van den Heuvel et al., 2008; Zucca et al., 2019). We examined whether our generative model recapitulates these properties. Figures 2 and 3 summarize the key findings in terms of variation in β and L_s, respectively.

Figure 2 (also see Supporting Information Figure S6) shows the effects of variations in β. We observed that small changes in β had marked effects in the topology of the generated networks, especially in terms of the formation of long-range bundles of axonal projections. Specifically, when β was either small (β = 0.98) or large (β = 1.02), generated networks lacked distant connections (Figure 2A–C). This was because local nodal guidance was too weak (strong) to allow axon terminations (outgrowths) of long-range connections, as detailed in Supporting Information Figure S3. Nevertheless, for moderate β values (β = 0.99, 1 and 1.01; Figure 2A–C), strong long-range connections emerged, decreasing in strength with greater β. Remarkably, these connections replicated the strong long-range connections observed in empirical connectomes that cannot be explained by a simple EDR (Deco et al., 2021; Roberts et al., 2016). That is, it is not only distance dependence that matters, but the outgrowth of axons depending on competing attractive factors also contributes to brain connectivity profiles (Cahalane et al., 2011), including long-distance connections.

Compared with null networks generated from a constrained random walk (see the Methods section), our model networks featured brain-like axonal projections resembling U-fibers and white matter tracts. As shown in Figure 2A (most evident for β = 1), axonal projections were organized into bundles at a distance from their origins and defasciculated before reaching their targets. This is a natural consequence of axon guidance in the presence of target-released diffusible chemoattractant (i.e., the node exerted a distance-dependent attractive force) that has been observed in past axon pathfinding studies (Hentschel & van Ooyen, 1999). Adjacent nodes were connected via U-shaped fibers that gradually steered according to the dynamically changing axon guidance. In contrast, axons generated by the null model failed to form organic fiber bundles. Additionally, the proportion of variance in the model generated fiber length explained by Euclidean distance is comparable with empirically observed values (Akarca et al., 2021), a characteristic missing in the random walk null model (Supporting Information Figure S8).

A negative association between connection weight and the Euclidean distance between nodes was evident across the range of β values considered (Figure 2C). Connection weights (normalized by the average nodal strength of connected regions; see the Methods section) spanned four orders of magnitude and were most parsimoniously modeled by lognormal distributions (Figure 2D, Supporting Information Figure S9). In contrast, connection weights for the random-walk null model were most accurately described by a gamma distribution (Supporting Information Figure S9) and exhibited less variability, distinguishing them from model networks and empirical connectomes (Figure 2B–D, Supporting Information Figure S10).

To investigate whether our generated networks showed scale-free degree distributions, we adopted the approach developed by Clauset et al. (2009), as detailed in the Methods section. In brief, the approach returned a p value describing the goodness of fit, and evidence of a scale-free distribution was deemed plausible if more than 50% of the networks in a population had a p > 0.1 (Broido & Clauset, 2019). For each representative parameter combination considered, we generated 1,000 networks and evaluated evidence for scale-free degree distributions. Figure 2E shows the networks with the median p value for each parameter combination (see Supporting Information Figure S11 for p value distributions). A scale-free behavior was evident for β = 0.98 and 1, yet it was missing in networks generated with other β values and our null networks. Crucially, β = 1 generated networks that were simultaneously characterized by all the connectomic properties were evaluated in this section.

We next characterized the impact of variation in step length, L_s, on properties of the generated connectomes (Figure 3; also see Supporting Information Figure S7). Greater L_s led to networks with higher connection density (ranging from 4% to 75% in Figure 3A). This was because a larger L_s resulted in fewer trajectory updates, and thus, the past guidance an axon received had a longer lasting impact on subsequent wiring. As a result, if two axons originated from the same node yet different coordinates, with a greater L_s, the difference in initial coordinates made their subsequent trajectories less likely to converge, which led to more diverse wiring and denser networks. We also found that the choice of L_s impacted connectivity weight distributions (Figure 3B–D). Certain values of L_s (L_s = 1 and 2; Supporting Information Figure S9) generated networks with lognormally distributed connectivity weights. Negative associations between connection weights and distances were consistently observed, and strong long-range connections that deviate from the EDR were found in generated networks except for L_s = 0.1 (Figure 3C). In addition, all representative L_s parameters generated networks with scale-free degree distributions (Figure 3E, Supporting Information Figure S11).

Collectively, these results suggested that combinations of β and L_s generate connectomes with realistic properties, including scale-free degree distributions, brain-like axonal bundles, negatively correlated connection weight and distance, lognormally distributed weights, and strong long-range connections that deviate from the EDR. Specifically, comparing the two model parameters β and L_s, variations in β had a stronger effect on degree distributions, whereas changes in L_s were more closely related to weight distributions. These observations remain robust if generated networks are controlled for network density via thresholding (Supporting Information Figures S12 and S13).

So far, we have focused on characterizing the distributions of connection weight and nodal degree generated by our model. In this section, we examine whether our generative model could give rise to connectomes exhibiting complex topological properties, including SW and modularity. We specifically investigated the network average clustering coefficient (CC), characteristic path length (CPL), SW, and modularity Q, benchmarked to weight and degree preserved null networks (see the Methods section). These topological features are hypothesized to relate to the functional segregation and integration of brain networks (Bassett & Bullmore, 2017; Fornito et al., 2015; Sporns & Betzel, 2016).

Using an exhaustive grid search (see the Methods section), we investigated how weighted topological properties of generated networks changed as a function of β and L_s. Figure 4A displays the variation of network topology among the parameter space that generated brain-like connection weight and nodal degree distributions. Topological properties were evaluated at a network density of 10%, yet the patterns of topological variations were insensitive to network densities (Supporting Information Figure S14). Across the investigated parameter space, generated networks consistently displayed small-world and modular structures. A higher CC and a longer CPL were evident relative to weight and degree preserved null networks. Variations in model parameters impacted the network topology in a continuous manner. Increasing β and decreasing L_s led to weighted networks with higher clustering, longer CPL, stronger SW, and weaker modularity (also see Supporting Information Figures S15 and S16).

To supplement this analysis, we visualized the axon organizations of example networks from different positions of the parameter space. As shown Figure 4B, a decrease in β and an increase in L_s improved the prevalence and strength of medium- to long-range connections, reducing segregation (measured by clustering) while promoting integration (measured by efficiency, i.e., the inverse of CPL).

In summary, variations in model parameters shifted network topology by adjusting the balance between short-range and long-range connectivity. While the generated networks consistently showed small-world and modular organizations, the degrees of these properties vary. Combined with the results from Figures 2 and 3 (i.e., analyses on bundle structure, connection weights, and degree distributions), certain combinations of model parameters (e.g., β = 1 and L_s = 1) were capable of generating networks that resembled all the evaluated connectomic features.

Finally, we focused on fitting the two parameters governing our generative model to individual human connectomes. Our motivations here were twofold. First, the results above suggested that networks generated by our model exhibited varying degrees of SW and modularity. It was unclear if certain parameters gave rise to these properties numerically close to empirical connectomes and if these parameters resided in the range that generated networks simultaneously manifested brain-like connection weights, nodal degrees, and axonal bundles. This question could be answered via a parameter optimization approach. Second, model optimization enables comparison of individuals in terms of their fitted parameters and can provide insight into interindividual variability in the processes guiding connectome development. For example, interindividual variation in parameters of existing generative models is associated with various traits, including age, sex, a schizophrenia diagnosis, and social-economic disadvantage (Akarca et al., 2021; Betzel et al., 2016; Carozza et al., 2023; Faskowitz et al., 2018; Simpson et al., 2011; Sinke et al., 2016; Siugzdaite et al., 2022; Zhang et al., 2021). Thus, a parameter fitting scheme is a prerequisite for future cohort studies to apply our model.

Parameters for the current state-of-the-art connectome generative models are typically optimized by minimizing an energy function that compares the degree, clustering, betweenness centrality, and Euclidean space edge lengths between generated and empirical connectomes (Betzel et al., 2016). However, our model does not generate networks with nodes that correspond one-to-one to empirical connectomes; as a result, the discrepancy in Euclidean space edge lengths cannot be evaluated and the classical energy cost is not applicable (despite this, when evaluated using the degree, clustering, and betweenness centrality, our model achieved a better fit relative to the established geometric model; see Supporting Information Figure S17). Therefore, we developed a new method to fit the parameters of our generative model to individual connectomes (detailed in the Methods section) and applied it to estimate parameters for 1,064 participants in the Human Connectome Project (HCP) (Van Essen et al., 2013).

The parameter estimation method aimed to minimize the dissimilarity between the generated network and an individual’s connectome, in terms of CC, CPL, and modularity Q. Figure 5A displays the best-fit model parameters for the HCP population, and networks generated by the fitted parameters numerically resembled topological features for which they were optimized (Supporting Information Figure S18). Individual parameters were found significantly different between males and females (p = 2e⁻¹⁵; Supporting Information Figure S19), confirming the model’s capability of capturing interindividual variations in connectomes. Using the group-averaged model parameters (β = 0.9968, L_s = 0.9531), we generated networks and investigated their weight, degree, and axonal bundle structures. As shown in Figure 5B, generated networks simultaneously manifested brain-like axonal bundles, negatively correlated connection weight and distance, lognormally distributed weights, scale-free degree distributions, and strong long-range connections that deviate from the EDR. These results suggested that a “sweet spot” of parameters can be identified to generate human brain-like connectomes possessing a host of empirically observed properties.

Generative models can provide insight into the wiring principles governing brain network organization. Most existing models generate unweighted connectomes in which connections are either present or absent (Betzel et al., 2016; Oldham et al., 2022; Simpson et al., 2011; Sinke et al., 2016; Vértes et al., 2012). However, it is well-known that connectivity weights are diverse and span multiple orders of magnitude (Fornito et al., 2016). In this work, we developed a novel model that generates connectomes weighted by axon counts and demonstrated that our model can recapitulate lognormal connectivity weight distributions. Our model is spatially embedded and characterizes the dynamics of axonal outgrowth. We demonstrated that our model manifests key topological properties of the connectome and yields axonal fiber bundle structures that resemble white matter fascicles. We were also able to fit our generative model to individual connectomes, enabling future cohort studies to apply the model.

Our work substantially builds on the seminal generative model developed by Song et al. (2014). Whereas this earlier model considered static axon propagation in a fixed direction, we established a dynamic axon pathfinding model in which the attractive forces guiding axonal outgrowth are continuously updated. With an appropriate selection of β and L_s, we observed that our model can generate networks with properties that are consistent with empirically mapped connectomes. In particular, negative associations between connectivity weights and the Euclidean distances between nodes were evident—a property that has been found in numerous studies (Betzel & Bassett, 2018; Ercsey-Ravasz et al., 2013; Roberts et al., 2016). We also found that connectivity weights were lognormally distributed and showed scale-free degree distributions. This suggests that a simple axon growth mechanism can generate key properties of a connectome’s topological architecture.

The choice of β and L_s determines whether the generated networks show complex topologies. As β is increased, distant nodes exert less influence on an axon’s growth and, thus, axons are attracted by neighboring nodes, forming short-range connections. This leads to the formation of clusters between spatially adjacent nodes and weakens the long-range connections that are critical to network efficiency. In contrast, a larger step length parameter L_s enables axons to propagate beyond a local nodal sphere of influence, increasing the prevalence and strength of long-range connectivity and improving network integration. While β and L_s impact generated axons through mechanisms that are both similar and distinct, their combined effects lead to generated networks with SW and modularity, akin to empirical connectomes.

Biological plausibility is a key characteristic of our model, where axon outgrowth is determined by the summed attractive forces exerted by each node. This is inspired by the observation that axons grow by responding to complex and combined effects of multiple guiding cues (Wadsworth, 2015). The two governing parameters also build on empirical observations in neural systems. The force decay parameter β determines the distance-dependent decay of attractive forces exerted by nodes, modeling the concentration decay of guiding cues with distance from releasing sites (Kaiser et al., 2009; Murray, 2002). The step length parameter L_s governs the extent to which an axon can change its trajectory per unit length, and it can represent the combined effects of multiple factors such as a growth cone’s growing speed and its sensitivity to molecular guidance (Alberts, 2017). While step length is seldom considered as a key parameter in applications such as tractography (Tournier et al., 2002), in vitro evidence suggests that it is a vital factor that models axonal growth. For example, variations in growth step length have been observed between frog and chick neurons, as well as between normal and regenerating frog neurons (Katz et al., 1984). Relating the parameter L_s to empirical data on neuronal growth rates of different species can help validate the biological interpretability of our model.

Our model generates connection weights by counting axons between nodes, a method distinct from other recently proposed models. Different weight inference approaches all have their unique strengths. Using a weighted stochastic block model, Faskowitz et al. (2018) inferred connection weights from network blocks, highlighting the community architecture of connectomes. Based on an unweighted connectome generative model, Akarca et al. (2023) introduced connection weights via minimizing the redundancy in network communicability, capturing dynamics in the strengthening and weakening of connections. In contrast, axon counts used in our model are intrinsically akin to streamline counts synonymous with structural connectomes, emphasizing the physical nature of connections as neural pathways.

Elucidating the mechanisms governing the formation of long-range connections remains a pivotal yet unsolved question in connectome generative model research. Early work suggested that, in addition to the distance rule, a topological homophily rule is required to promote the formation of long-range connections (Betzel et al., 2016; Vértes et al., 2012). Recently, the biological plausibility of topological rules was questioned, and homophily in gene expression and cytoarchitecture was hypothesized to contribute to long-range wiring (Kerstjens et al., 2022; Oldham et al., 2022). Nevertheless, existing frameworks failed to explain the specificity of long-range connectivity (Betzel & Bassett, 2018). Moreover, it is unclear how brain elements can perceive distant pairs without prior global knowledge of topology. Due to the spatial embedding of brain networks, a distance component might be required for brain elements to search for their wiring pairs. By including a pathfinding component, our model simulated long-range connections, including those that deviate from the EDR. Growth cones were sequentially guided by the strong local cues of a series of intermediate nodes before reaching their distant destinations (despite that weak attractive forces exerted by distant nodes still contribute). This mechanism is consistent with the hypothesis of intermediate targets in axon guidance, whose suggestive evidence has been observed in model organisms such as Drosophila and mice (Canty & Murphy, 2008; Dickson, 2002).

We conclude by acknowledging the limitations of our work and providing guidance for future improvement. Firstly, as the first attempt to generate connectomes from dynamic axon guidance, the model simplifies the brain as a two-dimensional circle and ignores complex brain structures such as sulci, gyri, deep gray matter, and cerebrospinal fluid. While this approach contributed to model simplicity and axon visualization, it also introduced limitations, such as the loss of nodal correspondence between generated and empirical connectomes. Using a realistic brain mesh/volume to incorporate three-dimensional neuroanatomical constraints in axonal outgrowth would not only naturally address these limitations but also entail higher computational demands. Additionally, neurodevelopment is patterned in space and time. Neurogenesis and axon growth initially occurs in the neural tube (Kintner, 2002); the development of long-range projection and commissural fibers occurs concurrently with the morphogenesis of primary foldings; the accelerated growth of short cortico-cortical association fibers coincides with the formation of secondary and tertiary foldings (Chavoshnejad et al., 2021). Future studies could investigate the evolution of connectivity in a spatially embedded developing brain. Secondly, we assume that all brain regions have the same distance-dependent attractiveness and that all axons are equally sensitive to guidance from brain regions. These assumptions are likely breached in the brain given the diversity in regional properties (e.g., cortical thickness, curvature of folds, laminar structure, cellular composition, and neuronal density), neuron types, and guiding cues (attractive and repulsive, chemical and mechanical). For example, we do not consider the impact of axonal repellents in our current model. Repellents are critical to the formation of long-range commissural fibers (Dickson, 2002), and modeling the forces resulting from repellents is an important avenue for future investigation. Recent efforts in generating high-resolution brain maps such as molecular and cytoarchitectural profiles (Amunts et al., 2013; Arnatkevičiūtė et al., 2019; Hansen et al., 2022; Markello et al., 2022) might provide an opportunity to refine the assumptions and improve the model’s capacity. Thirdly, while our model generates axon organization that is visually akin to axon bundles and white matter fascicles, factors that contribute to axon bundling are not considered. Incorporating fasciculation mechanisms such as the contact attraction between axons and axon-released guiding cues (Hentschel & van Ooyen, 2000) might help to build a more nuanced white matter and connectome organization. In addition, the model implements a deterministic axonal guidance rule, and as such, stochasticity, which is also fundamental to neural development (Carozza et al., 2023; Hassan & Hiesinger, 2015), was not taken into account. Future work could evaluate the robustness of the model with the presence of stochasticity, such as random noise in guiding cues and axon growth. Axonal branching that results in multiple synaptic connections is also ignored in our study. This property allows for one-to-many connections, promoting the formation of sophisticated connectome topology, whose impact could be evaluated in future work. Finally, our model generates macroscale connectomes, yet this is achieved by simulating axons that are microscale anatomical concepts. Future studies could investigate our model’s application in generating microscale connectomes.

An overview of our model is described in the Results section. Here, we provide finer details of the model, elaborating on aspects including node heterogeneity, path constraints, axon termination, and parameter specifications.

To parcellate the hypothetical gray matter, N_n node centers were evenly positioned along the circle perimeter, such that the angular distance between adjacent nodes equals 2π/N_n. Next, nodal heterogeneity was introduced by randomly perturbing node center coordinates. This was accomplished by applying a uniformly distributed angular displacement, ε ∼ ρ * U(−π/N_n, π/N_n), to each node center. Specifically, ρ = 1 was used in this study to maximize nodal heterogeneity while preserving the sequential arrangement of nodes along the perimeter.

Axons were simulated based on the distance rule in Equation 1. To encourage axons to traverse relatively noncurved trajectories, regularity constraints were applied to each axon from the second extending step onward. The regularity constraints stipulate that the angle formed between the direction of two consecutive steps cannot exceed the angle θ. In other words, if the angle between two consecutive steps exceeds θ, the second step is adjusted such that the angular difference is forced to θ (Supporting Information Figure S2).

Ideally, axons would terminate on the circle circumference, connecting two points of the hypothetical gray matter. However, not all simulated axons can successfully reach the circle perimeter. When the value of β was small, a “black hole” region emerged within the circle, as shown in Supporting Information Figure S3. Axons entering the “black hole” cannot escape, forming a circular trajectory of infinite loops. To address the problem, a parameter S_max was introduced to stipulate the maximum number of growing steps allowed. Axons failing to reach the circle circumference within S_max steps were considered unsuccessful and were excluded from network construction and analyses.

Eight parameters were defined in the model. Unless otherwise specified, default values of parameters (Table 1) were used. A comprehensive justification for the parameter choice was included in the Supporting Information.

We investigated the associations between edge weights and distances, and the weight distributions in generated networks. The edge weights were normalized by the average nodal strengths of the two regions they connect and were evaluated using a common logarithm scale.

The weight-distance associations were evaluated by calculating the Pearson’s correlation coefficient between edge lengths (i.e., the Euclidean space distance between two nodes connected by an edge) and the edge weights. Weight distributions were evaluated against fitted lognormal, gamma, normal, exponential, and Weibull distributions using one-sample Kolmogorov-Smirnov (KS) test.

Results of weight and degree analyses were visualized for representative parameters (L_s = 1, β = 0.98, 0.99, 1, 1.01, and 1.02; β = 1, L_s = 0.1, 0.5, 1, 2, and 5). These parameters were selected to generate diverse network properties while delineating the isolated effects of each parameter. Compared with an exhaustive grid search (used in a later section to evaluate global topology), this approach enabled us to uncover details (Figures 2 and 3) that were obscured in summary metrics (i.e., Pearson r, KS statistics, and p values).

Null networks generated from a constrained random walk were used to benchmark model networks. Specifically, axon growth directions were randomly sampled from U(−θ, θ) rather than being calculated from the distance rule in Equation 1. The step length parameter of L_s = 1 was used. All other parameters remained consistent with the model.

To characterize the global topology of generated networks, model parameters were drawn from a grid combination of β and L_s (0.99 ≤ β ≤ 1.01, 0.1 ≤ L_s ≤ 2.1; 101-by-101 grid). This parameter space was determined from preliminary experiments and was found to generate networks that replicated connectomic features. To account for the stochastic variability arising from node and axon sampling, 50 networks were generated for each parameter combination, forming 50 network landscapes. The network topology corresponding to each parameter combination was described by the average topological metrics over 50 landscapes.

We considered the weighted CC, characteristic pathlength, SW, and modularity Q of generated networks. Because topological measures are fundamentally related to network density and connectivity strengths, all generated networks were thresholded and normalized to have the same network density (10%) and total connectivity (2e⁵). Parameters whose generated networks have a density smaller than 10% were ignored. Topological measures were evaluated using the Brain Connectivity Toolbox (BCT), benchmarked to weight and degree preserved null networks constructed using the null_model_und_sign() function in BCT.

This study utilized the HCP Young Adults (1,064 subjects) datasets (Glasser et al., 2013; Uğurbil et al., 2013). A comprehensive description of data acquisition and connectome construction has been detailed elsewhere (Mansour L. et al., 2021). The HCP connectomes were mapped to the Desikan-Killiany atlas, comprising 68 cortical and 16 subcortical brain regions. Networks were thresholded to a density of 10%.

Parameters were optimized using a Monte Carlo method through an exhaustive grid search (see the Global Topology of Generated Networks section). To account for the stochasticity-dependent inaccuracy and unreliability, and to improve the computational tractability, we employed the fast landscape generation and the multilandscape method (generating 50 landscapes) developed by Liu et al. (2023). For each landscape, the best-fit parameters (with the smallest RMSE; values shown in Supporting Information Figure S18) were selected, and the average across 50 landscapes was deemed the optimal parameters.

Data were provided by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research and by the McDonnell Center for Systems Neuroscience at Washington University. The computation is supported by The University of Melbourne’s Research Computing Services and by the Petascale Campus Initiative. Y.L. is funded by the Melbourne Research Scholarship. C.S. is supported by the Australian Research Council (grant number DP170101815). R.F.B. is supported by the National Science Foundation (award ID: 2023985). D.A. is supported by the James S. McDonnell Foundation Opportunity Award and the Templeton World Charity Foundation, Inc. (funder DOI 501100011730) under the grant TWCF-2022-30510. M.A.D.B. was supported by an NHMRC Investigator Grant (1175754). A.Z. is supported by research fellowships from the NHMRC (APP1118153).

Supporting information for this article is available at https://doi.org/10.1162/netn_a_00397.

Yuanzhe Liu: Conceptualization; Formal analysis; Methodology; Software; Writing – original draft; Writing – review & editing. Caio Seguin: Conceptualization; Methodology; Supervision; Writing – original draft; Writing – review & editing. Richard F. Betzel: Methodology; Writing – original draft; Writing – review & editing. Daniel Han: Formal analysis; Writing – review & editing. Danyal Akarca: Methodology; Writing – original draft; Writing – review & editing. Maria A. Di Biase: Conceptualization; Methodology; Supervision; Writing – original draft; Writing – review & editing. Andrew Zalesky: Conceptualization; Methodology; Supervision; Writing – original draft; Writing – review & editing.

Yuanzhe Liu, Melbourne Research, Award ID: 720358. Caio Seguin, Australian Research Council (https://dx.doi.org/10.13039/501100000923), Award ID: DP170101815. Richard F. Betzel, National Science Foundation (https://dx.doi.org/10.13039/100000001), Award ID: 2023985. Danyal Akarca, Templeton World Charity Foundation (https://dx.doi.org/10.13039/501100011730), Award ID: TWCF-2022-30510. Maria A. Di Biase, National Health and Medical Research Council (https://dx.doi.org/10.13039/501100000925), Award ID: 1175754. Andrew Zalesky, National Health and Medical Research Council (https://dx.doi.org/10.13039/501100000925), Award ID: APP1118153.

The data of HCP are publicly available from https://www.humanconnectome.org/. The code for the generative model is available from https://github.com/yuanzhel94/connectome_from_pathfinding.

Small-worldness:

Brain regions are densely connected locally and can efficiently communicate globally.

Modularity:

Brain regions are grouped into modules such that connections are dense within modules and sparse between modules.

Generative model:

A model that can generate new samples of brain networks.

Chemoaffinity:

Neurons make connections based on the chemical signals that guide them to their targets.

Degree:

The number of brain regions that a region connects to.

Clustering:

Three interconnected brain regions form a closed triangle.

Scale-free:

The degree distribution follows a power-law function such that a few brain regions have many connections.

Characteristic path length:

The average smallest number of connections required to travel between two brain regions.