Constructive connectomics: How neuronal axons get from here to there using gene-expression maps derived from their family trees

During brain development, billions of axons must navigate over multiple spatial scales to reach specific neuronal targets, and so build the processing circuits that generate the intelligent behavior of animals. However, the limited information capacity of the zygotic genome puts a strong constraint on how, and which, axonal routes can be encoded. We propose and validate a mechanism of development that can provide an efficient encoding of this global wiring task. The key principle, confirmed through simulation, is that basic constraints on mitoses of neural stem cells—that mitotic daughters have similar gene expression to their parent and do not stray far from one another—induce a global hierarchical map of nested regions, each marked by the expression profile of its common progenitor population. Thus, a traversal of the lineal hierarchy generates a systematic sequence of expression profiles that traces a staged route, which growth cones can follow to their remote targets. We have analyzed gene expression data of developing and adult mouse brains published by the Allen Institute for Brain Science, and found them consistent with our simulations: gene expression indeed partitions the brain into a global spatial hierarchy of nested contiguous regions that is stable at least from embryonic day 11.5 to postnatal day 56. We use this experimental data to demonstrate that our axonal guidance algorithm is able to robustly extend arbors over long distances to specific targets, and that these connections result in a qualitatively plausible connectome. We conclude that, paradoxically, cell division may be the key to uniting the neurons of the brain.

sequences of expression templates. Growth cones terminate (or become dormant) when they can neither 211 improve their template match, nor replace their template. The growth cone does distinguish between no 212 gradient due to failure, and no gradient due to successfully reaching a peak. The cone will transition to 213 the next possible state in both cases. If the cone continues to make progress with its search templates, it 214 will ultimately reach a leaf state. If however, the local signals do not offer suitable gradients the growth 215 cone will fail after exhausting its options. 216 In executing this search algorithm, the initial cone and its clones extend axons along all the routes 217 in brain-space that offer contiguity in brain-space of familial expression patterns encoded in the lineage 218 tree, and as a result create a particular repeatable axonal arborization from source to target leaf nodes.

219
Results 220 We sought evidence of expression address maps by performing our hierarchical decomposition of covari-221 ances on experimentally measured gene expression data. We used the gene expression data published by 222 the ABI in their Developing Mouse Brain Atlas [51]. Their data are provided as 3D grids of isotropic 223 voxels. The expression energies of the ⇠2000 genes were measured by in-situ hybridization and take any 224 non-negative value (see Methods). 225 Our prediction, leaning on the principles outlined above, was that the measured covariances should 226 reveal a dual set of systematic patterns: a hierarchical ordering in expression space, and a nesting of 227 region in brain space. Here, 'expression space' denotes the abstract multidimensional space in which 228 the gene expression profile of a cell can be represented by a point; and 'brain space' or 'space' denotes 229 physical 3D space. We will use 'profile' for a fixed relative expression of all genes within a cell or voxel, 230 'pattern' to denote an evident spatial regularity in gene expression, and 'organization' to denote the 231 systematic order underlying these patterns.

232
Global spatial expression hierarchy 233 We performed our hierarchical decomposition (described in the Concept and Methods sections) on the 234 experimentally measured spatial gene expression, and consider to what extent that hierarchy could 235 constitute an estimate of the mitotic lineage tree. Figure 5 shows the results for P28. Other time points 236 from E11.5-P56 are shown in Figures S13-S20. 237 The hierarchical decomposition is based exclusively on the gene expression covariance measured across 238 unordered sets of voxels: The physical locations of voxels within the brain were not taken into account. 239 However, when the voxels are now considered in their correct physical locations, we observe that the sets 240 selected during the unordered hierarchical decomposition form nested and spatially continuous regions 241 spanning brain space. Thus, the hierarchy of nested voxel sets found by our analysis in expression space 242 parallels a hierarchy of spatially continuous regions in brain space, which is not entailed by the analysis. 243 To confirm that these patterns are due to intrinsic spatial organization of gene expression rather than 244 being induced by the analysis, we performed the same decomposition on voxels with randomly shuffled 245 gene expression values (see Methods). After shuffling, all spatial structure vanishes (Figures S13-S20). 246 We quantified the degree to which the hierarchy of asymmetries derived from gene expression projects 247 into physical brain space by measuring the spatial spread of the voxel set. If the nested voxel sets 248 form continuous regions in space beyond the first few generations, then the spread of their constituents 249 should decrease over successive generations as the sub-populations become more resolved (and therefore smaller). Otherwise, if the constituent voxels are scattered across space, their spread would remain 251 roughly constant. We find that indeed the spread consistently reduces with generation ( Figure 5c). 252 In addition to spanning brain space, the patterns are consistent over time from E11.5 to P56. The 253 temporal consistency cannot be measured directly, because there is no clear map between individual 254 voxels of subsequent time points. To side-step this problem we instead project the voxels of one time 255 point onto the hierarchy measured at another time point. In Figure 6 the voxels of all available time 256 points are projected onto the hierarchy measured at P28. 257 To quantify the temporal consistency of the spatial patterns we measure how many voxels are sorted 258 into the same hierarchical bin when projected into the hierarchy measured at the original time point 259 versus the hierarchy measured at P28. We find that the spatial patterns at the various time points agree 260 significantly above chance (Figure 6b). 261 We also measure the correlation between the axes of covariance C ! ov among the same hierarchical 262 node at different time points, and among different hierarchical nodes at the same time point. We find 263 that the axes of the same node across time correlate more than the axes of various nodes within a single 264 time point ( Figure S12).

265
Hierarchy exists in small gene sets 266 As we have selected genes without bias, the spatio-temporal patterning may be a general property of 267 gene expression, rather than the specific property of a particular set of specialized genes. We explored 268 this possibility by randomly selecting sets of genes of various sizes, and then measuring how well the 269 spatial patterns were maintained.

270
From the spatial gene expression data we found that small sets (⇠50) of genes selected randomly from 271 the database can already achieve an accuracy (ratio of voxels sorted in the same hierarchical bin) with 272 the signal over all genes (Figure 7). So, the lineage identity of a cell may be encoded in the profile of 273 any small set of genes.

274
Simulated mitosis induces spatial hierarchy 275 We verified that our division model is able to explain the above results by simulating the mitotic devel-  The simulation embodies the constraints of the model. The basic element of the simulation is a space-279 occupying (not a point) cell that expresses a profile of genes. When a cell divides, the expression profiles 280 of its two daughters are normally distributed variants of the expression of their mother. The cells are 281 additionaly positioned in a spatial 3D grid. When a cell divides it pushes a nearby cells in a random 282 direction to make space to place the mitotic daughters adjacent to one another (see Methods). In this 283 way, we are able to efficiently simulate the growth of a volume of cells, analogous to the embryonic brain.

284
A particular lineage tree of cells is generated by recursive application of this division rule.

285
The simulation results in a mass of 500,000 cells, each expressing 500 genes, and distributed over 286 100 independent lineages. This mass was divided into 21140 voxels of 3x3x3 cells (excepting voxels 287 on the outer surface, which may contain fewer cells), to emulate the voxelation of the ABI data. We 288 then applied the same hierarchical decomposition as was applied to the experimental data, yielding 289 qualitatively similar results ( Figure 8). This simulation indicates that the spatial and genetic constraints 290 of the model are sufficient to explain the spatial patterns that we observed in the developing mouse 291 brain.

292
Simulated axons traverse the brain 293 We simulate the process of axonal growth to demonstrate that they can use the Familial Guidance Model 294 described in the Concept to navigate through the voxels of the ABI gene expression atlas (Figure 9). The 295 Guidance Model in the Rationale was described in terms of individual cells and a given lineage tree. In 296 applying that model to the experimentally observed ABI data, we have to consider that the expression 297 data is voxelated, with each voxel containing many cells; the lineage tree is only estimated, as described 298 in the previous section; and we cannot yet predict which exact set of genes participates in the address 299 space. So, we need to test whether robust long range guidance is possible at all using our proposed 300 mechanism, rather than predict specific projections.

301
The growing axons were simulated using a spatial-state graph approach. In this approach, the axon 302 traverses a graph where each nodes represents the growth cone at a spatial location (i.e. a voxel) with a 303 specific receptor configuration. Two nodes are connected by an edge if they are spatially adjacent and 304 represent the same receptor configuration-this represents a move of the growth cone in space; or if they 305 represent the same spatial location, and the receptor configurations are adjacent in the lineage tree-this 306 represents a transition of the growth cone in state.

307
The axon can only traverse an edge to an adjacent node if the gene expression in that node is a better 308 fit to its receptors than the current node (i.e. if the currently held ancestral state is more correlated to 309 the adjacent node). This 'biological' algorithm resembles Dijkstra's algorithm for finding the shortest 310 route between two nodes of a graph [10].

311
Simulated axons were found to be able to extend up to about 10 µm source-tip distance, with axonal 312 lengths of up to 16 µm (Figure 9b). The axon length is longer than the source-tip covered distance 313 because the axons do not move in a straight line. Nevertheless, the relation between path and straight 314 distance is close to linear for many axons, particularly those of shorter ranges.

315
The axonal guidance is robust against signal noise. We demonstrated this by adding noise to the 316 gene expression after the lineage tree was reconstructed, but before the axons began guidance, effectively 317 reducing the signal to noise ratio. The added noise does not significantly change the trajectories of 318 the axons (Figure 9e). However, when gene expression is shuffled completely, either before or after 319 reconstructing the lineage tree, axons fail to navigate beyond 1 or 2 voxels (not shown).

320
Similar axon length distributions are found from the fully simulated tissue (Figure 9d). This result 321 further supports the hypothesis that simple constraints on mitosis can induce the address space.

322
The control case for axonal arborization is a random walk axon of identical total length, whose growth 323 cones are still self-avoiding, but able to move in random directions at every step. We find that the random 324 walk axons have much lower specificity, robustness, and spatial reach, than the navigating axons (See 325 Figure 9e). Note that for the navigation algorithm the total length of the axon is an implicit result, rather 326 than a set parameter. Thus, the random walk axon already over-informed compared to the navigating 327 axons. 328 We generated a typical connectivity matrix by simulating 500 axons rooted in voxels sampled uni-329 formly from the available data, (Figure 9f). The matrix is sparse and block structured, with many 330 off-diagonal components (rather than narrow diagonal band), indicating a specific, regionalized, con-331 nection pattern. Remarkably, these typical connections conform to reasonable anatomical patterns, as 332 can be appreciated by comparing the block-structured axonal connections with anatomical regions taken 333 from the independent annotations of the Allen Brain Institute (Figure 9f. We emphasize that these 334 anatomical annotations are not used at any point during simulation or analysis.

336
The literature describing the progressive organization of vertebrate brains over embryonic and evolution- the formation of their many synaptic connections require a global orchestration of guidance cues [17,50] 341 at various spatial scales.

342
In this paper we have explored the hypothesis that the mitotic lineage tree, which is implied by the  We further propose that the address space is navigable by axonal growth cones, which are able to grow 352 to specific target addresses by matching local gene expression patterns to those of successive nodes of a 353 lineage tree traversal. We call this process Familial Guidance ( Figure 4). In other words, the expression 354 of the brain, and the growth cone's ability to exploit it, are dual consequences of the brains developmen-355 tal process which both creates the Familial Address Space as a consequence of cellular differentiation, 356 and then exploits that differentiation for active cellular organization including the formation of axonal can be encoded within the genetic budget. Our proposed mechanism resolves this issue by showing 363 that the lineage tree can efficiently install unique labels in target tissue, and that navigating axons can 364 recognize them due to their shared origin in the cellular GRN. It also extends the scope of comprehensive 365 molecular labels from the retino-tectal projection to the brain at large. provide voxelated (rather than tied to pre-conceived anatomical regions) spatial expression data of de-368 velopmentally relevant genes throughout brain development. Previous analyses of these and related 369 atlases have been largely concerned with identifying profiles of co-expression that support anatomical 370 organization [37, 39, 38, 34, 51]; and also whether these regional profiles can be explained in terms of the 371 known expression of specific cell types [18, 58,60]. Although there are also systematic transcriptional 372 similarities across cortical areas [34,19], no global map-like organization has been reported as yet.

373
Our results now indicate that systematic spatial patterns of gene expression covariance do exist and are 374 widespread in the embryonic and postnatal brain. These patterns involve non-specific groups of genes, 375 occur on multiple spatial scales across the entire brain and spinal cord, transcend neuroanatomical 376 boundaries, and are consistent at least from E11.5 to P56. Interestingly, we found that the primary 377 axis of variance corresponds spatially to the dorso-ventral axis of the embryonic brain, rather than the 378 antero-posterior axis that is expected on the naive assumption of greater variance along a longer axis.

379
This suggests that the patterns do not simply reflect the geometry of the developing embryo, but are 380 related to controlled regionalization in embryogenesis itself. 381 We explored the embryological origin of these patterns by analyzing the statistical structure of the 382 expression covariance [20], rather than the relationship between expression and anatomy [34] or to pheno-383 typic expression of cells [18]. The essence of this structure is that the differential gene expression between 384 arbitrary sibling branches of a lineage tree (the asymmetry) in expression space has a dual expression 385 as covariance across the region of brain space occupied by the leaf nodes of those sibling branches, as   The covariance patterns indicate only that common gradients of expression exist across sets of genes, 391 and seem to be hierarchically organized. Our results do not of themselves indicate which genes contribute 392 most strongly to the patterns, nor which, if any, are actually utilized for addressing. It remains to be 393 tested whether the spatial organization observed in the current data is restricted only to the ⇠2000 genes 394 that the ABI chose for assaying [51], and consequently to the subset of ⇠1, 240 that we have analyzed. 395 However, both the experimental data and our simulations indicate that the organization does not arise 396 from the expressions of a specific set of genes dedicated to encode spatial structure, but rather can be ure S12) of the covariance nesting is suggestive of an address space. This putative address space has an 402 interesting property: Because the nested regions are a projection of the lineage tree onto 3D brain space, 403 regions composed of cells that are closely related in their respective lineage trees are also close in space.

404
Thus, the address map is a systematic arrangement of cells in terms of their ancestral gene-expression, 405 and so provides an implicit encoding of cell lineages that could be used as a relative localization mecha-  Thus, the exploration paths of the growth cone can be seen as the lineage tree hung from a leaf (with some 418 pruning). So, growth cone routes are anti-differentiating up the tree to some ancestral node, followed by 419 re-differentiating toward the leaf states accessible from that ancestral node.

420
An appealing aspect of this lineage-induced address is that it greatly simplifies the evolution of complex 421 spatial organization of cells. The systematic spatial labeling of cells is given as a direct consequence of 422 mitotic specialization and cell proximity. Evolution needed only to discover how to exploit this labeling 423 for organizational migration of cells and their components (e.g. growth cones). It could opportunistically 424 select a set of gene products for axonal growth cone guidance, because most gene sets will encode a 425 similar spatial pattern. This generality of the address space could also help to explain the wide range of  The Familial Address Space model is entangled by two factors. Firstly, the mitotic root of the 429 developing brain is difficult to define exactly. It seems reasonable to consider as starting point a small signal that is of interest here, rather than the absolute expressions of particular genes. We may also allow 448 that the stochastic profiles be subject to cell-external or internal factors, provided that these influences 449 are reproducibly regulated as part of the developmental process (and thus not due to environmental 450 noise external to the embryo).

452
There has been substantial progress in understanding how axonal growth cones respond to local guidance 453 cues [5]. They are exquisitely sensitive to local gradients, able to detect gradients on the order of a 454 few molecules across their span [41]. However, physical noise, ligand binding, and other signal detection 455 considerations indicate that molecular gradients alone are insufficiently robust to explain axonal guidance, 456 particularly at longer spatial scales [17,5]. Over these longer distances the algorithmic rather than the 457 reactive aspects of guidance rather relevant. and their resultant guidance cues secreted into the extracellular space or exposed on cellular surfaces. 488 The model asserts that the navigation sequence can be generated if the axon inverts its developmental 489 program to anti-differentiate to precursor expression states, and thereby traces a route through the 490 lineage tree. There is evidence that cells and also neurons are able to de-differentiate as a whole [14,30,491 57]. However, we require only that de-differentiation occur on some subset of the genome relevant the 492 familial marks. While there is as yet no systematic work on this question, there is nevertheless clear and 493 growing evidence that growth-cones used elaborate local context-dependent mRNA processing during 494 their guidance behavior [55,25,9].  be resolved in order to improve the prediction. Obviously, the range and resolution of experimental data 519 must be improved: the ABI data offers only subset of genes, and even these data degraded by averaging 520 over voxels and brains. Furthermore, the address space was inferred from the set of all genes. However, 521 this set is only one of many possible sets of genes that support the true address space (Figure 7). It is 522 likely that evolution has selected a particular set of genes to establish a particular address space that 523 supports well evolution's preferred connectivity. Unfortunately, this particular set of genes is as yet un-  Overall though, the address space induced by mitosis, as well as the guidance process that it supports, 531 is consistent with reasonable axonal projection pattern. Full agreement between our simulations and 532 experimentally observed projections will depend on the agreement between our estimated differential gene expression model, and the true differential gene expression generated by the actual gene regulation 534 network of the mouse. These simulations also confirm at least a partial projection of the expression 535 space onto brain space. This conformity is not self-evident, because the high dimensional expression 536 space cannot be faithfully projected into the lower, three dimensional brain-space. Even in the best 537 embedding, the pairwise distances between the embedded nodes in 3D Euclidean space cannot exactly 538 match the pairwise path distances between nodes in the tree. The error can be made arbitrarily small in 539 the limit of many dimensions. This limit is not relevant for 3D physical space because its dimensionality 540 is fixed; but for gene expression space it is relevant because the dimension can be increased by recruiting 541 genes for the embedding. Fortunately, constraints on mitotic daughter migration will result in at least 542 some regions of continuity in the lineage tree embedding. Thus, although not all traversals through the 543 lineage tree will be matched by traceable paths in brain space, those traversals whose embedding in 544 brains space provides for continuity of expression signals will be successful. This property is reflected 545 in that our connectivity matrix is not fully connected or diagonally structured, but sparse and block-546 structured (Figure 9f). The manifold of traceable paths, and so connection probabilities, will no doubt 547 be influenced by the anatomical distortions imposed on the growing cell mass by factors such as relative 548 mitotic rates, cell size, asynchronous axonal outgrowth, ventricular volumes etc.

549
Note that our algorithmic approach differs from more usual methods for the generation connections, and then deciding whether a connection exists between them by evaluating a probability distribution.

556
This method will establish suitable entries in a table, but does not explain how these connections should 557 be grown in space. To satisfy the rule by a developmental algorithm the growth cones must perform a 558 random walk in space, governed by a fixed probability of extension; and connect to encountered neurons 559 also with some fixed probability. This simple connection algorithm closely reproduces the empirically 560 observed axonal length distributions. However, because the behavior of these parameterized stochastic 561 models depend on random data (rather than fixed data), they are also unable to explain the repeatability 562 (across individuals) patterns of axonal trajectories and connectivity observed in biology. Repeatability 563 would require that the 'random' walk be decided by a frozen random number generator, so simulating a 564 deterministic guidance mechanism, whose data is that frozen random number. cells seem to preferentially target their clonal siblings [8,27,59], rather than simply nearby targets. Our analysis of gene expression in embryonic and postnatal mouse brains reveals a hierarchy of spatial 610 patterns of expression covariance that extend over the entire brain, and are stable over the available 611 data. This organization is present across the 1240 genes analyzed. However, they are also present in the 612 expressions of random subsets of as few as 50 genes. The organization is consistent with a multi-scale 613 address space that could be exploited for cell migration or growth cone guidance. Our simulation studies 614 confirm that this organization can be generated by persistent asymmetries of gene expression introduced 615 by the successive mitoses of the lineage trees that give rise to the brain, provided that mitotic daughters 616 do not stray too far from one another after their birth.

617
Due to the generality of these mitotic constraints, it is likely that similar map-like structures exist 618 also in other tissues, and may provide a fundamental scaffold for cell migration and tissue organiza-619 tion. However, the Familial Address Map has particular relevance for neurons, whose many stereotyped 620 connections cover distances up to the scale of the whole brain. 621 We conclude that the fundamental wiring of brain can be compactly encoded and expressed through the key to uniting the neurons of the brain. The resolution of the paradox is that division in reverse is 625 unification.

626
Future work must establish: which specific sub-set of genes is used for axon navigation; how the 627 growth cone reverts its host's differentiation and how receptors are generated to recognize an ancestral 628 state; and how the address space, that is the geometry of the brain and spatial gene expression, are 629 tuned to realize a specific observed connectivity.

630
Contrary to the prevailing reductive approaches to understanding the wiring of the brain, this paper 631 has taken a more global synthetic view. While much more effort will be required to confirm the var-632 ious implications of our approach, the theory and available data are remarkably consistent; and offer 633 the prospect that the connectome and its functioning can be more readily understood in terms of the 634 global mechanisms that generate it, rather than from interpretation of the final wiring diagram, just as 635 inspecting source code is more revealing of principles of operation than inspecting the compiled program. Figure 1: The connectome is the result of a constructive process that starts ultimately with the zygote, and involves the two aspects of first generating a mass of cells with various types, and then routing axons through this mass to their proper targets. An observer's description of the resulting detailed mouse connection matrix (right bottom) takes at least 10TB to encode. However, as development occurs largely in isolation, all instructions to construct this connectome must fit into the 1GB of genetic material of the zygote. This implies that neural progenitors have efficient methods for expanding the highly compressed wiring instructions into axonal trajectories. To do this, they need to, as they proliferate and differentiate, install a space of molecular addresses that axons can exploit for navigation. and physical space (c), so that cells of related differentiation have similar expression profiles (similar colors) and are nearby one another physically. Consequently, a trajectory from one leaf node to another through the lineage tree (a: red arrow) often corresponds to an unbroken trajectory through both gene expression space (b: red arrow) and physical brain space (c: red arrow). An axon navigates by inverting its source neuron's instance of the global genetic differentiation program (a: top right). This inversion generates a sequence of expression profiles that correspond to ancestral states and so act as guidepost profiles. d The axonal branch configures its growth cone to match the sensed expression to the internally generated expression, and so moves to the direction that improves that match . When the match can no longer be improved by moving, the axon updates the its internal state to the next ancestor, and repeats. If the match between internal and external expression can be improved by moving into multiple different directions, or by transitioning to multiple different states, the single axonal branch is split into two new branches that continue to execute the same algorithm, but whose independent states may subsequently diverge. When an axonal branch arrives at a leaf state, both in expression and physical space, navigation of that branch is complete and local synapses are formed. Figure 3: a Cells are points positioned in high-dimensional expression space, where each axis represents the expression of one gene. Here, this high-dimensional space is reduced to 2D dimensions for plotting purposes, so that their 2D distance approximates their high-dimensional distance. In our division model, the differential expression between a parent cell c 1 and its daughters c 2 , c 3 is a normally distributed random vector representing the genetic state transition from parent to daughter, denoted 2 = c 2 c 1 . (Here we use the division of the root progenitor 1 as a running example for any division.) The differential expression between two siblings, which we call the parent's asymmetry, is denoted 1 = c 3 c 2 = 3 2 . As a result, the correlation in gene expression between two cells reflects their distance through the lineage tree. (See c for verification of this process by numerical simulation.) b The expression of a progenitor can be estimated as the mean expression over its leaf progeny; and the asymmetry of a progenitor can be measured as the main axis of variance across its progeny. The diagram shows only the leaves of the lineage tree show in a-they have identical positions in embedded expression space. Each nested contour encloses the progeny of a progenitor; lines within the countour indicate the main axis of variance across the enclosed progeny; and dotted circles the average expression across the progeny. The sets of progenies for individual progenitors can be obtained by iteratively splitting the progeny along their main axis of variance, so with a decision boundary (black line with arrow) orthogonal to this axis. c Numeric simulation of expression profiles induced by our division model, and subsequent reconstruction of expression profiles and mitotic asymmetries from the leaves of the simulated tree. The root expression c 1 is drawn from a normal distribution with zero mean and unit variance. The expression profiles of other cells are generated recursively by adding differential expression patterns i , which are also normally distributed. (All random number are drawn independently.) The determination (squared correlation) was measured between between the true and reconstructed asymmetries (blue), and true and reconstructed expressions (orange). d Progenies group naturally in brain space according to their ancestry. Shown is a 2D simulation of growing tissue, started from a single root, only constrained to not detach from one another and not pass through each other. Figure 4: Navigation of an axon (red branching arrow) through the familial address space. Throughout the figure, similarity in color denotes similarity in gene expression profile. a The axon traverses the brain by traversing a sequence of familial states of the lineage tree that is implicit in its genome. The growth cone uses the sequence of familial states as successive search templates in brain space, and so navigates from a source leaf node to a number of target leaves. Familial states (colored circles) correspond to nodes of the encoded lineage tree. For purpose of explanation, the tree is hung from the leaf state corresponding to the axon's source neuron, rather than from its root node as in 3b. Terminal states of leaf (existant) nodes have a solid circumference, while ancestral states in the interior of the tree have a dotted circumference. Transitions between states occur downward, along the arrowed arcs, beginning at the source leaf (red encircled) and ending at (some) other leaves. The original tree root can be recognized as the only state having two edges, rather than three (since the root progenitor has no mitotic parent). b Various decision scenarios that the axon encounters during traversal. Each familial state is characterized by a profile of gene expression, whose distribution across all cells peaks at one or more locations in brain space. The gradient of a state in the familial address space is the frequency of encountered cells that test positively for a familial state. By selecting a particular familial template, the growth cone tunes into the corresponding expression gradient and filters out the others. If the tuned gradient is in range, the growth cone follows it to arrive at one of that gradient's peaks (case indicated by [1]). If the tuned gradient is not in range [3], the axonal branch of that growth cones fails. When the axon arrives at a peak, its growth cone tunes to the next downstream familial state, and so on, until a leaf state is found. If multiple downstream states are in range, the axon branches [2], with each branch tuned to one of the possible downstream states. The axon also branches if the gradient is bifurcated by a valley, so that the axon can follow an upward gradient in multiple directions [4]. Each branch pursues a different direction, but in this case they are tipped with growth cones in the same state (unlike the branches in scenario [2].) When a growth cone reaches a leaf state, guidance terminates [5]. c Cells have composite genetic identities, with one component (small inner circle) inherited from each ancestor state. The overall state of a leaf cell is the aggregation of these components (3). A growth cone can test whether a cell possesses a component by selecting the familial state template corresponding to that component, and then matching the internally produced gene expression to that of the tested cell. d Various regions of the brain correspond to branches of the mitotic lineage tree. Consequently, the regions are nested and each marked by the component of the genetic identity code corresponding to the common progenitor of the region. The bins are colored to represent the hierarchy: the parent bin has the average hue of the child bins. This coloring is applied throughout the paper. (c) Although regions are nested by construction (hierarchical decomposition), we quantified the extent to which the regions are also continuous by measuring their spatial spread (average distance from the region centroid) as a function of their depth in the hierarchy. At the root of the hierarchy the spatial spread covers the entire brain, and we expect that as the depth increases the spatial spread (i.e. the mean distance from the region centroid to the constituent voxels) decreases. To make the different time points and simulation comparable we present the spreads as a fraction of the root spread. The solid line indicates the median spread over all regions at that depth, and the gray area the first (below) and third (above) quartiles. As expected, the mean distance from the centroid decreases as the regions become more resolved at with depth. Figure 6: The root asymmetry measured at P28 is projected to the other available embryonic and postnatal developmental time points, and compared to the root asymmetry measured at the respective time point. (a) First division of the hierarchy, but the direction of variance used to sort the bins is derived from P28, rather than from the data of the time point itself (except Original P56). This temporally projected pattern only has small differences with the patterns derived from the original data (compare Original P56 to Projected P56). When the expression data is shuffled over voxels and genes, maintaining pooled expression statistics but destroying covariance structure, all spatial patterning disappears. Images are proportional to their actual brain sizes. (b) Quantification of the agreement between the original and projected hierarchy, measured as the proportion of voxels in matching bins, at different levels in the hierarchy. (Although the images in a are 2D, quantification is done on the 3D voxels.) The number of possible bins grows exponentially with tree depth, and so chance level decreases inverse-proportionally (dashed line), quantitatively verified by the shuffled case (yellow line). P28 projects onto itself, and is hence in perfect agreement. The other time points show an agreement consistently above chance. Consider that a mismatch at a shallow depth cannot be corrected at a deeper depth, and so mismatch can only accumulate. Figure 7: Random sets of genes of various sizes from embryonic age E11.5 were selected, and the spatial hierarchy they exhibit was compared to the hierarchy exhibited by the grand set of all genes at hierarchy. To compare hierarchies all voxels are projected onto both hierarchies. For each matching choice the score is incremented proportionally to the depth of the bin. As such, 1 indicates that the all voxels are sorted into corresponding nodes of the hierarchies, and the dotted line indicates the score if all voxels were sorted into hierarchical bins randomly (as in the shuffled case). The hierarchy established from a set of 20 random genes already agrees largely above chance with the original pattern. proposed to explain the results. Simulated 'brain' sphere composed of voxelated leaf cells was generated by 300,000 mitoses distributed over 10 independent lineages. Cells express 500 genes. Asymmetrical mitoses induce differential changes in gene expression. Each voxel contains 3 ⇥ 3 ⇥ 3 = 27 adjacent cells. Similar to experimental results, the hierarchical decomposition of covariance in gene expression voxels independent of location (left), is mirrored by matching decomposition in space (right). Figure 9: Simulated axons use familial guidance to navigate through the voxels of the ABI Developing Mouse Brain gene expression atlas. a Arborizations of 50 example axons, show in a sagittal projection of the ABI atlas. Each arbor is the collection of all branches that an axon could potentially navigate using this gene expression space. Each axon is colored according to its source region. The colorings correspond to those of Figure 5b. b Straight-line distance between the beginning of a branch (soma) and end of that branch (top) versus the actual path length. Branches are points sorted in hexagonal 2D bins, whose color intensity indicates the number of branches in that bin. c Same as a but on a tissue grown in simulation (as in Figure 8). d As b, but for the simulated tissue of c. e The dissimilarity between axons beginning from the same voxel (measured as average minimum distance), under varying levels (10% or 30%) of expression noise. (Because the navigation algorithm is deterministic the 0% noise case produces identical neurons.) The familial guidance dissimilarity is compared against a random walk axon of the same path length. f Connectivity matrix corresponding to the connections made by the axons of a. The connections conform to reasonable anatomical patterns. The anatomical regions marked on the matrix are taken from the annotations of the Allen Brain Institute. They are not used for the analysis.

783
All experimental data analyzed in this paper were originally published by the Allen Institute for Brain 784 Science [51], and are available at https://developingmouse.brain-map.org.

785
Code availability 786 All custom code is available at …(will be added with paper publication)  The data are provided as 3D grids of isotropic voxels of various sizes. The expression energies of the 801 ⇠2000 genes were measured by in situ hybridization and take any non-negative value, while -1 indicates an 802 invalid measurement in that voxel. 'Expression energy' is a combined measure of density and intensity.

807
The atlas data was retrieved through the API provided by ABI. The ABI expression grids were used as 808 published, without performing any additional re-sampling or interpolation (see below for preprocessing).

809
Thus, the voxel sizes were maintained as published by the ABI.  In order to make the gene expression energy levels roughly comparable across genes, the expression 826 values were normalized to unit variance and zero mean over the voxels at that developmental time.

828
The voxels measured at one time point are sorted into the leaves of an estimated lineage tree through our 829 hierarchical decomposition procedure. The procedure starts at the root of the tree, to which all voxels 830 are initially assigned. Each iteration of the procedure evaluates the voxels assigned to a node of the tree, 831 and reassigns each voxel to one of the node's two daughters. The procedure stops once every leaf node 832 is assigned exactly a single voxel.

833
An iteration considers the gene expression of the voxels assigned to a node. The collected expression 834 can be expressed as a matrix X, where each row corresponds to a voxel and each column to a gene. The 835 voxels will be split over the daughter nodes along the axis of greatest variance.

836
The axis of greatest variance is the eigenvector corresponding to the greatest eigenvalue of the covari-837 ance matrix. The covariance matrix is computed by first centering the data by subtracting the empirical 838 mean from each column X 0 ij = X ij P k X kj /n, where X ij is expression of the jth gene (column) in the 839 ith voxel (row), and n is the total number of voxels (rows). The covariance matrix is then Q = X 0 T X 0 .

840
The main axis of covariance is the eigenvector C

845
Based on these coefficient we sort the voxels into two subsets, namely one set (arbitrarily denoted L for 846 left) with L = {i|w i < 0}, and R = {i|w i > 0}. These voxels of these sets are assigned to the left and 847 right daughter nodes, respectively. The decomposition procedure is then repeated recursively on these 848 two daughter nodes.

849
If a node is assigned only a single voxel, the process terminates for that branch. The process as a 850 whole terminates when all branches have terminated.

852
Controls were performed to ensure that the observed spatial patterns are due to the spatial distribution 853 of experimental gene expression rather than being due to any inherent properties of the analyses. The Although randomness is used to establish the expression profiles, the (frozen) random deviations are 884 used as deterministic process with statistics that are indistinguishable from a random process. This is 885 analogous to fixing the seed of a random number generator: when the seed is fixed, the exact sequences 886 of numbers is reproduced, but the statistics still seem random. To create multiple lineage trees, the algorithm is still performed only once, but starting not from one, 899 but from multiple root nodes. So, the total number of divisions is possibly divided unequally over the 900 lineage trees.

902
The cells of the lineage tree are positioned spatially by Algorithm 3, as illustrated in Figure 10. We 903 developed this algorithm because it is simpler and more tractable than a direct simulation of soft (or 904 solid) body physics.

905
Our Familial Guidance model is simulated virtually, either on the voxel grid of measured expression, or 907 on a grid of simulated cells.

908
The axon begins in some chosen leaf voxel, and takes as its initial template the expression state of 909 its leaf voxel. Then, at each step, the growth cone senses the expression of the voxel it occupies, as 910 well as expressions of the immediately adjacent voxels. The cone then moves into the adjacent voxel 911 whose expression is most similar to its current template, so extending a new axonal segment between  In executing this search algorithm, the initial cone and its clones extend axons along all the routes 921 in brain-space that offer contiguity in brain-space of familial expression patterns encoded in the lineage 922 tree ( Figure 4).

923
For our axonal simulations brain space is discretized: Each spatial position corresponds to a (measured 924 or simulated) voxel. Spatially adjacent voxels are connected by an edge. These nodes and edges form 925 a graph encoding the geometry of the brain. The 3D positions of the nodes are used to establish the 926 spatial graph, but ignored thereafter.

927
The adjacency of nodes is established through a Gabriel tessellation, which is a subgraph of the 928 Delaunay tessellation [15]. In a Gabriel tessellation a edge of the Delaunay tessellation is kept only if the 929 sphere of which the edge is the diameter contains no other points. This criterion ensures the spatial graph 930 is connected, but that there are no edges across large empty spaces, such as ventricles and contours. This 931 is an improvement over the vanilla Delaunay tessellation, which always contains the convex hull of the 932 points, and therefore connects, for example, the rostral tip of the olfactory bulb to the cortex.

933
To navigate, axons follow signals on the spatial graph. A signal on a graph is a scalar value associated 934 with each node, and the gradient of the signal is a value associated with an edge, that is the difference 935 between the values of the nodes. The gradient depends on the direction the edge is traversed, and swaps 936 sign if the edge is traversed in the reverse direction.

937
The signal for a growth cone depends on the current state of the growth cone, and the expression of 938 the nodes of the graph. The state of the growth cone is a profile of gene expression that corresponds to 939 a node in the lineage tree. The signal over the nodes of the spatial graph relative to a growth cone state 940 is the correlation between the growth cone state's expression profile, and each node's expression.

941
For efficiency, the signal is only allowed to exist in the progeny of the ancestor whose state the growth 942 cone has adopted. This constraint reduces the search space of the growth cone significantly, without 943 significantly changing the routes taken by the growth cones.

944
The prominent action of the growth cone is to climb this signal by spatially moving across the graph,

947
In addition to moving, the growth cone can also change state. The state machine governing the 948 transitions the growth cone can take is (isomorphic to) the lineage tree, which is estimated through our 949 hierarchical decomposition. So, the growth cone can only transition to the parent state, or either of the 950 daughter states, of its current state.

951
To simulate this process, a spatial-state graph is constructed. The nodes of the spatial-state graph  On the spatial-state graph there is only a single guidance signal, attributing to each node the corre-956 lation between the node state's expression profile and node's expression.

957
The navigation of an axon starting from a voxel is simulated by executing Dijkstra's algorithm [10] 958 from a source node in the spatial-state graph to all possible nodes containing leaf states, allowing only 959 movements along positive gradients. For graph implementations the igraph library with python bindings 960 was used (https://igraph.org).

961
Axons were visualized using threejs (https://threejs.org) Figure S11: Hierarchical decomposition, as in Figure 5, but for all available time points. Only depth 3-the lowest tile in Figure 5-is shown, but other depths can be inferred by grouping similar colors. Decompositions were performed independently of one another (unlike Figure 6, where established hierarchies are projected across time points). The spatial spread of hierarchical regions goes down with hierarchy depth at each measured time point. Figure S12: Asymmetry profiles identifying regions of the hierarchical decomposition are poorly correlated within the hierarchy, but correlated across time points. Upper triangle Pairwise correlation coefficient between the estimated asymmetries C ! ov measured at the root of the hierarchies at various time points. Altough the asymmetry measurement is done independently at each time point, the main direction of covariance across all voxels is correlated. Generally, nearby time points are more correlated than distant time points. This correlation is surprising a priori, because the absolute gene expression changes from E11.5 to P56. Lower triangle Pairwise correlation coefficient between the estimated asymmetries C ! ov at the root of a hierarchy and other asymmetries within the same hierarchy. (Each column represents a time point, and each row a depth of the hierarchy, with the root at zero depth.) In contrast to standard principal component analysis, orthogonality between components is not enforced by our hierarchical decomposition. Nevertheless, we find that many pairs of components are poorly correlated. This implies that the direction of strongest covariance is not along any single direction for all subsets of voxels, but is rotated in high-dimensional expression space at each iteration of the decomposition. The model assumes that differential gene expression vectors , and consequently the asymmetrieŝ are independent. This matches the observation in the experimental data that the successive C ! ov i are poorly correlated in expression space. The poor correlation is not by construction, because unlike PCA (Principal Component Analysis), orthogonality is not enforced by our decomposition.