A generative model of the connectome with dynamic axon growth

Connectome generative models, otherwise known as generative network models, provide insight into the wiring principles underpinning brain network organization. While these models can approximate numerous statistical properties of empirical networks, they typically fail to explicitly characterize an important contributor to brain organization – axonal growth. Emulating the chemoaffinity guided axonal growth, we provide a novel generative model in which axons dynamically steer the direction of propagation based on distance-dependent chemoattractive forces acting on their growth cones. This simple dynamic growth mechanism, despite being solely geometry-dependent, is shown to generate axonal fiber bundles with brain-like geometry and features of complex network architecture consistent with the human brain, including lognormally distributed connectivity weights, scale-free nodal degrees, small-worldness, and modularity. We demonstrate that our model parameters can be fitted to individual connectomes, enabling connectome dimensionality reduction and comparison of parameters between groups. Our work offers an opportunity to bridge studies of axon guidance and connectome development, providing new avenues for understanding neural development from a computational perspective.


Constraints on axon growth directions
The model specified the maximum angular disparity  an axon can make at each increment step.smaller than or equal to .In growth step  − 1, the axon grew from point P1 (blue dot) to point P2 (black dot), and its growth course in step  − 1 is  ⃑ !(blue arrow).The black arrow denoted the summed force  ⃑ that the axon received at posi1on P2. ∆ represented the angular difference between  ⃑ and  ⃑ !, and the angle  between black and red dashed lines was the maximum angular disparity allowed.When ∆≤ , the axon's growth direc1on in step  ( ⃑ " ) equals  ⃑ .b) Same as a) but for ∆> , the angular difference between  ⃑ ! and  ⃑ " (red arrow) is forced to .
In addition, two special cases were considered.
Case 1: if F ,,⃑ and  ⃑ ! are in opposite direction ( ), the clockwise/counterclockwise direction of  could not be determined.In this scenario, we forced  ⃑ ( to be in the direction of  ⃑ ! .
Case 2: if  ⃑ has a magnitude of 0, its direction cannot be determined.In this scenario, the model randomly sampled a non-zero  ⃑ , after which the angular constraint  was applied.

"Black hole" region and 𝜷 values
In Results and Methods, we argued that growing axons could fail to land on the circle circumference because they were trapped within a "black hole" region.This happened when the values of  were small, such that the guidance from local nodes was too weak to allow axon terminations.As shown in Fig. S2, as  increases, trapped axons suddenly escape to form middle-range to long-range connections.However, if  is too large, long-range connections cannot form because of the strong attractive force from local nodes.This explains why small variations in  leads to dramatic change in generated networks: only a small critical range of  allows long-range connections to form. ) (the number of nodes) is determined by the data of interest.In this study, we used 84 nodes to correspond to the number of brain regions in the Desikan-Killiany atlas.In addition, 300 nodes were used in scale-free analysis to reduce the bias from finite network size.
The choice of  (controlling the node heterogeneity) has been explained in Methods. = 1 was used in this study to maximize nodal heterogeneity while preserving the sequential arrangement of nodes along the perimeter.
* (the number of axons simulated) was fixed to 2 + .This value was chosen as a balance between the favor to large  * and the computational affordability.The axon origins were sampled in a random manner; thus, a large  * is favored to capture the connectivity distribution of generated networks.Consider two networks (connectivity matrices  ! and  ( ) that only differ in the number of axons simulated ( * and  * , respectively;  * → ∞,  is a constant).The connectivity expectation would be  ( =  ! .This further guided the trick (linear scaling) used in parameter optimizations (See Parameter optimization section in supplementary materials).
The effects of R (circle radius),  (angular constraint), and  ' (growth step length) on generated networks were closely related.Fig. S3 demonstrated a linear relationship between R and  ' .The association between  and  ' is complex and non-linear; however, the ratio / ' could be approximated as the angular changes that an axon can make per unit growth length.Therefore, we decided to fix R and  and investigate the variations in  ' , making  ' a key tunable parameter in the model.R was set to a random constant ( = 30 was used), and the choice of  = 15° was inspired by the angle values typically used between successive steps in tractography.
,*-(the maximum number of growing steps allowed) ensured that simulations could stop even if axons were trapped by the "black hole" zone. ,*-= 3/ ' was used such that axons were able to connect furthest points on the circle (Euclidean distance of 2R), while a margin of R was included to allow curved axon trajectories.
describes the distance decay of attractive force, controlling the relative contribution of guidance exerted by adjacent and distance nodes.As a result, values of  are closely related to topology of generated networks, and we make it a key tunable parameter of the model.

Variance in fiber lengths explained by Euclidean distances
Previous work found that the variance in fiber lengths explained by Euclidean distances is inconsistent between studies, ranging from 22-79%.Despite the wide range reported, the results from our model generated axons fell within the range (generally between 70-80%), whereas values from the random walk null model generated axons were much larger (Fig. S4).

Weight distributions (null model, and alternative distributions ks test results)
In the main text, using representative parameter combinations, we suggested that the model was able to simulate networks whose connection weights were best described as lognormal distributions.Here, we statistically compared the fitted KS statistics among candidate distributions (lognormal, normal, gamma, exponential, Weibull, Fig. S5).50 networks were generated for each parameter combination, and the KS statistic of fit were compared between distributions using paired t-tests.Lognormal distribution was found to be the best fit (significantly smaller KS statistic compared to other distributions, all  < 0.05) in networks generated with  = 0.98, 0.99, 1, 1.01, and 1.02 ( ' was fixed to 1), and  ' = 1 and 2 ( was fixed to 1).It was also the best fit in networks generated with the group average parameter for the HCP population, despite the KS statistic difference between lognormal and Weibull distributions were not significant.It should be noted that connection weights in null networks were not best described as a lognormal distribution, despite the KS statistic is small (Fig. 2 and Fig. S5).To better characterize the difference between model and null networks, distributions of normalized weights were shown in Fig. S6.Connection weights in null networks exhibited less variability (strong connections are rare) relative to in model networks.

p-value distributions of scale-free test
In the main text, we suggested that the model was able to generate networks with scale-free degree distribution.The conclusion was drawn from the observations that more than 50% of networks generated by a parameter combination showed a  > 0.1 in the scale-free test.Fig. S7 showed a histogram of -values in 1,000 generated networks, for each evaluated parameter combination.Null networks did not show a scale-free property.Scale-free was evident in  = 0.98 and 1 ( ' was fixed to 1), and all  ' values considered ( was fixed to 1).The group

Contour plot patterns of topological properties are insensitive to network density
In Fig. 4, we evaluated the weighted topological properties of generated networks at a network density of 10%.We found that the patterns of contour plots are insensitive to network density and show the results for network density of 5% in Fig. S8.

Raw values of topological metrics
In Fig. 4a, we showed the weighted topological properties of generated networks, normalized by weight and degree preserved null networks.Here in Fig. S9, we show the raw values of weighted topological properties of CC, CPL, and Q. SW is not included here because it is by definition normalized.It should be noted that modularity Q changed with parameters in a different way compared to Fig. 4.This demonstrated that when  is increasing (or  ' is decreasing), modularity increases faster in null networks compared to model networks.where the matching index model used both geometric and topological information.In conclusion, geometry-dependent dynamic axon guidance generates networks better replicate empirical connectome than a simple distance rule, whereas this improvement is smaller than considering complex topological information.A related concern is that this may disrupt the inter-individual variations in total connectivity.To test this, we randomly selected a sample of 100 participants from the HCP cohort.Model networks were linearly scaled to the total connectivity matching each individual before parameter optimization.We found that fitted parameters derived from the two approaches were comparable.
As per previous studies , we report the cost value of selected model networks in Fig. S11.
Selected networks showed an RMSE of around 0.5 standard deviations on average, with comparable values among clustering coefficient, CPL, and modularity.

Associate model parameters to individual traits
We optimized the model parameters for individual connectomes participating in the HCP dataset and associate them (,  ' ) to age and sex.No significant association was found to age, probably because HCP is a healthy young adult dataset in which effects of development and aging are less

Fig
Fig.S1illustrates how this angular constraint is applied.

Figure S1 .
Figure S1.Angular constraint .a) axons grow in the direc1on of the summed force if the computed angular difference is

Figure S2 .
Figure S2.Effect of  on axon length.a)-e) each displayed a subset of 1,000 simulated axons from networks shown in Fig. 2a. = 0.98, 0.99, 1, 1.01, and, 1.02 respec1vely (from a to e). a) A "black hole" region is evident with small values of .Axons were trapped and failed to form long-range connec1ons.b)-c) As  increases, long-range connec1ons start to appear.d)-e) If  is too large, long-range connec1ons disappear due to strong local aNrac1ve forces.

Figure
Figure S3.R and   are linearly related.a) A network generated with R = 30 and  $ = 1.b) A network generated with R = 60 and  $ = 2.All other parameters, including the angular coordinates of nodes and axon origins, were the same.The two parameter combina1ons generated the same network.

Figure S4 .
Figure S4.Variance in fiber lengths explained by Euclidean distances is consistent with empirical observa@ons in model but not in null.
Figure S5.KS sta@s@cs of fiDed candidate distribu@ons.Legends were ranked by ascending mean KS sta1cs.

Figure S6 .
Figure S6.Probability density func@on of weight distribu@ons (normalized by nodal strengths) in simulated networks.Weights in null networks were small, deno1ng a lack of strong connec1ons.In contrast, strong weights could be found in model networks.
HCP population were found to be able to generate scale-free networks, despite the p value marginally above the threshold of  = 0.1.

Figure S7 .
Figure S7. value distribu@ons of scale-free tests.Each histogram displayed  values of 1,000 networks.

Figure S10 .
Figure S10.Our model achieved a beDer energy compared to the state-of-the-art geometric model, and a worse energy compared to the state-of the-art matching index model.All energy func@ons consider degree, clustering, and betweenness centrality only.

Figure S11 .
Figure S11.RMSE, and cost in CC, CPL, and modularity Q, of networks selected to op@mize parameters for the HCP cohort.
and females significant differ in  ' (Fig.S12), suggesting model parameters can capture inter-individual variations in connectomes.