Theory of branching morphogenesis by local interactions and global guidance

Branching morphogenesis governs the formation of many organs such as lung, kidney, and the neurovascular system. Many studies have explored system-specific molecular and cellular regulatory mechanisms, as well as self-organizing rules underlying branching morphogenesis. However, in addition to local cues, branched tissue growth can also be influenced by global guidance. Here, we develop a theoretical framework for a stochastic self-organized branching process in the presence of external cues. Combining analytical theory with numerical simulations, we predict differential signatures of global vs. local regulatory mechanisms on the branching pattern, such as angle distributions, domain size, and space-filling efficiency. We find that branch alignment follows a generic scaling law determined by the strength of global guidance, while local interactions influence the tissue density but not its overall territory. Finally, using zebrafish innervation as a model system, we test these key features of the model experimentally. Our work thus provides quantitative predictions to disentangle the role of different types of cues in shaping branched structures across scales.

In this Supplementary Note, we provide additional details on the model derivation, on the numerical simulations of branching with repulsion and external guidance, as well as on the data collection, analysis and tting to the model.

Supplementary Note 1: Derivation of the continuum model
In the continuum description, we concentrate on incorporating external guidance into a model of branching and annihilating random walks, in order to make simple predictions for the orientation of branches/tips during branching morphogenesis.

Alignment angle as a random variable
The growth of a branching tissue with active tips in an external potential can be in principle modelled in the framework of a biased and persistent random walk in two dimensions. However, in such a framework, a full expression for the time evolution of the probability density p(r, t) for an active particle to be at location r at time t in general cannot be obtained in a closed form [1], and will involve several bias and reorientation terms due to dierent frequencies arising from elongation vs. branching events. We therefore sought to reformulate the problem in a simpler framework for a single random variable. After realizing that if we focus on the alignment of branch segments with an external eld, rather than their exact position and polarity, we could restate the problem as nding how much the local angles of the segments diverge from the local direction of the external eld. To avoid confusion, we note that by branch segments we refer to the elementary vectors of a xed size , corresponding to the discrete steps taken by active tips, that are linked together consecutively to dene the branches of a tree.
Angular alignment of a branch segment is quantied by dierent angles for the two cases of external elds discussed in the main text: (i) For a horizontally oriented, axial external gradient (directed towards the positive x direction), the local angle ϕ of a branch segment with respect to the horizontal axis determines its alignment along the eld. (ii) For a radial external eld emerging from a central point of origin (which is relevant for the experimental geometry of the zebrash n), however, the angular alignment will be given by the angle dierence ψ: where θ denotes the angle of the active tip coordinates with respect to the origin of the external eld.
Supplementary Fig. 1(A) illustrates the alignment angles for the two dierent external potentials. Note that, the axial potential is a limiting case of a radial eld where the origin of the external eld is located at x → −∞. This choice then corresponds to setting θ = 0 where the alignment angle ψ becomes equal to the local angle ϕ. Because of its generality, we will derive our model for the radial case using the angle dierence ψ as the alignment angle in the following. Expressions for the axial case can then be obtained simply by setting θ = 0.

Eect of an external eld on the angle dierence ψ
We now consider an external eld that will inuence the directionality of growing tips. At each step, the angle of a given active tip will be modied both stochastically via the randomness associated to elongation direction and branching, but also deterministically via guidance from the external eld. For a branching network invading a circular region (in two-dimensions) we can dene a radial potential that will displace the tips towards the point P with coordinates given by the vector r + f c p c , where f c determines the strength of the external eld, r is the distance vector of length r between the active dimensions (with x, y coordinates as shown) highlighting the main variables used in the analysis: In a radial eld (right), for every branch segment with an active tip (orange node) one can dene the local angle ϕ that it makes with the horizontal axis (dotted lines), its angle θ to origin (star symbol), and the angle dierence ψ ≡ ϕ−θ. Accordingly, the alignment of a branch segment with the eld is determined by the angle dierence ψ, whereas for an axial eld (lefty), the origin of the external eld is located at x → −∞ (corresponding to θ = 0), and thus the alignment will be determined by the local angle ϕ only (ψ = ϕ). Blue solid lines represent (static) branch segments of size generated by tip elongation (blue circles indicate inactive tips after an annihilation event due to proximity to other branches) (B) Displacement of an active tip due to a radial eld. The external eld acts along the radial distance r of the active tip to the origin of the potential (labelled by star) with a strength determined by the coupling parameter f c . The branch segment is then aligned towards the extended point P but preserves its length . The angle dierence after displacement ψ * = ϕ * −θ * is then given by Eqs.(S2-S4), and can be approximated as a function of the eld strength f c and of its value ψ before displacement only, see Eq. tip and the eld origin, and p c is a vector of length parallel to r. We note that we discretise the problem by having tip elongation by a small length = 1 (relative to the overall size of the network) without loss of generality. After such a displacement by the external eld, the angle dierence is given where ∆ϕ ≡ ϕ − ϕ * and ∆θ ≡ θ − θ * denote respectively the changes in the local angle and in the angle to origin values. To express the new angle ψ * in terms of ψ and the external eld strength f c , we deduce the following trigonometric relations from simple geometric arguments, see Supplementary and where r * represents the distance of the active tip to the eld origin after the displacement by the external eld, and b determines the distance b between the point P and the branch node Q, see Supplementary Fig. 1 which is a function of f c and ψ only. Supplementary Fig. 1(C) provides a comparison between this approximation and the exact relation determined by Eqs.(S2-S4). This dependency on the external eld f c reects a sinusoidal reorientation model for the angle dierence ψ which was previously demonstrated for the gyrotaxis of spherical microswimmers [2], for cyanobacterial circadian oscillators [3] and for dipoles in a constant external eld [4].

Fokker-Planck equation with non-local jumps
In simulations of branching and annihilating random walks (BARWs), elongation events lead to a small rotational diusion of the branch segments, whereas branching events lead to large, abrupt changes in the local angle values of the active tips, as illustrated in Supplementary Fig. 2(A). These two dierent sources of noises are described by an a priori probability distribution λ(ϕ − ϕ ) for the dierence in the local angles before ϕ and after ϕ a jump. The external eld then subsequently acts to modify these jump probabilities. This simulation setup can thus be described by separating the jumps arising from elongation and branching events from the bias arising through the external potential. Because of the large, non-local jumps in the local angle values ϕ, the standard method of reaching the continuum limit by expanding around small ∆ϕ cannot be performed. We will therefore follow a generalized method developed for Lévy ights in Refs. [5,6] to derive the Fokker-Planck equation.
We rst derive the central result of Ref. [6] for generic non-local jumps in one spatial dimension (represented by the continuous variable x) and then apply it to the angle dierence ψ. We start with the continuum version of the dierence equation for the probability distribution P (x, t + ∆t) after a jump from site x to x within the time interval ∆t: where the transfer kernel encodes both the distance between the jumps x − x via the corresponding probability density function λ(x − x ), and the spatial dependency of the transition probabilities A(x) and B(x). H(z) denotes the Heaviside step function dened as H(z) = 1 for z ≥ 0 and H(z) = 0 otherwise. We furthermore assume symmetric jump distance distributions of the form λ(x − x ) = λ(x − x), i.e., the probability for dierent jump sizes only depends on the size of the jump and not on the direction. The directional biases are governed by the transition probabilities A(x) and B(x). The normalization for the transfer kernel is dened by integrating Eq.(S7) over the jump dierences z ≡ x − x : where we used the symmetry of λ(z) together with its normalization ∞ −∞ dz λ(z) = 1. We now multiply Eq.(S6) by e ikx and integrate over dx to obtain Substituting z = x − x in the integrals over x and using the Euler identity we get: where we introduced the Fourier transform of P (x, t) and dene By switching the integration bounds for λ − (k) and using symmetry properties of sine and cosine functions, we can express Eq.(S10) as follows (S12) We now introduce the cosine and sine transforms Using the convolution theorem for the second term in Eq.(S12) and recalling the relation A(x)+B(x) = 1, see Eq.(S8), we obtain the simple expression where * denotes convolution in Fourier space, as derived previously in Ref. [6].

Fokker-Planck equation for the angle dierence ψ
We now use this framework for non-local jumps to obtain an equation describing the temporal evolution of the probability distribution for the angle ψ. In general, changes in ψ after each elongation or branching event depend on the changes in the local angle ϕ − ϕ as well as on the changes in the angle to origin θ − θ values. However, after a small number of steps taken from the origin, the θ-dependent terms become negligible compared to changes in the local angle ϕ. The probability distribution for the jumps sizes of the angle dierence ψ can therefore be approximated by that of the local angle ϕ.
In fact, as we will show later, see Supplementary Fig. 2(B) and (C), the mean-squared displacements for these two angles will attain the same linear form for suciently large times, indicating that the free diusion of these two angular variables follow the same dynamics in this regime. As illustrated in Supplementary Fig. 2(A), the jumps in the angle values are described by uniform distributions with maximal jump sizes determined by δϕ e ψ e = π/10 and δϕ b ψ b = π/2 for elongation and branching events, respectively. Denoting the branching and elongation probabilities by p b and p e , respectively, and the change in the angle dierence after each jump byψ = ψ − ψ , we can express the jump size distribution as where for clarity we only express the positive part of the symmetric jump distribution, i.e. forψ ≥ 0.
Fourier cosine and sine transforms of λ(ψ) are then given by (S16) Now we use the Taylor expansions up to O(k 3 ) for the cosine and sine functions to obtain: can be dened as the diusion and mobility/advection coecients based on microscopic variables.
Inserting these expressions into Eq.(S14) and taking the inverse Fourier transform leads to We must now dene A and B to model the external eld which acts to reduce the alignment angle |ψ| at each jump. Recalling our discussion on a radial external eld, we introduce a sinusoidal force as described by Eq.(S5). Because A(ψ) + B(ψ) = 1, see Eq.(S8), and we want the drift on the particle to be determined by A(ψ) − B(ψ) = −f c sin(ψ), we get Inserting these expressions into Eq.(S19), and assuming identical mean stepping times τ ≡ t/n in the limit of large step numbers n, we obtain the Fokker-Planck equation In the steady state we impose no-ux boundary conditions such that At steady state, using the ansatz P st (ψ) = C exp µfc D cos(ψ) with constant C determined by the normalization condition π −π P st (ψ)dψ = 1, as well as the symmetry of the cosine function, this predicts that P st (ψ) should follow a von Mises distribution where I 0 (ν) = π −1 π 0 exp(νcos(ψ))dψ (S24) is the modied Bessel function of the rst kind of order zero and takes the role of the concentration parameter of the von Mises distribution. This is a central result of our analytical model, which we confront to experimental data in the main text. Importantly, this predicts the distribution of angles up to a single rescaled parameter ν, which in analogy to the Péclet number quanties the respective contribution of advection (arising from an extrinsic guiding eld) to diusion (arising from the intrinsic stochasticity of branching/elongation events) in the Fokker-Planck equation.
Note that the variance of the angle distribution ψ at steady state is thus directly related to this rescaled parameter: whereas the circular standard deviation takes the form [7] σ circ. ≡ −2 ln(R(ν)) , which deviates from the linear expression 1 − R(ν). Note that, unlike the latter, Eq. (S27) diverges for ν → 0 and therefore does not provide a reliable metric for suciently spread out distributions. For the parameter regime we investigate here, however, it is an adequate estimator for the standard deviation of the underlying von Mises distribution. It is instructive to compare Eq. (S27) with the standard deviation of the equilibrium distribution for a particle in a harmonic potential, which corresponds to a normal distribution for ψ and depends on the external eld f c as For suciently large values of ν, the von Mises SD given by Eq. (S27) exhibits a power-law behavior with a slightly larger decay coecient, see Fig.1 (F) in the main text for a comparison of the two scaling laws.
2 Supplementary Note 2: Simulation details Next, we briey summarize the details for the numerical algorithm to simulate branching and annihilating random walks (BARWs) in the absence of external eld or self-avoidance.

Branching and annihilating random walks (BARWs)
Similar to Ref. [8], we dene active tips that will elongate and branch randomly, taking discrete steps of unit size = 1 per discrete time interval τ = 1 and leaving behing inactive branch segments that remain immobile at coordinate r ≡ (x, y) on a two-dimensional plane. When an active tip comes in close proximity of an inactive branch segment, i.e., when the inactive branch is within an annihilation radius R a of the active tip, the latter will become inactive and immobile. We choose a rather small annihilation radius of R a = 1.5 to mimic contact-dependent membrane recognition. Note that the specic choice of the annihilation radius in two dimensions does not inuence the overall topology (up to global rescaling) of the networks because intersection of two branches occurs with probability 1 for all suciently small radii.
Simulations start with a single active tip at an initial position r 0 ≡ r(τ = 0) = (0, y 0 ) with a pre-dened polarity p 0 ≡ (1, 0), i.e., which is directed horizontally towards the right (higher x values), and proceed until a certain maximal time τ max = 200τ is reached. Because the active tips take a single step at each time interval τ , the maximal time τ max of the simulation also determines the maximal distance R max from the origin r 0 that the last surviving active tips can attain. For the experimental data, this parameter is therefore strongly constrained by the geometrical properties (overall size and shape) of the n, which we determine by counting the maximal number of steps in the coarse-grained reconstruction of the laments, see Supplementary Note 3.1. The initial polarity of the starting active tip allows us to dene an initial local angle ϕ 0 = 0 with respect to the x-axis associated with the tip.
The key parameter in this simulation setup is the branching probability p b , which will determine the frequency of branching events and thus also the nal density of the network.
The (i) elongation and (ii) branching events for the active tips are implemented in the simulation as follows: At each time point t, we draw n a (pseudo-)random numbers r j ∈ (0, 1) with j = 1, . . . , n a , where n a is the number of active tips at t. For each active tip, the elongation or branching occurs for r j > p b or r j ≤ p b , respectively. (i) When an active tip at coordinate r with local angle ϕ undergoes an elongation event, it takes one step forward such that its new local angle ϕ will take a random value uniformly distributed between ϕ ± δϕ e with δϕ e = π/10. The latter rule leads to a small rotational diusion of the tip during elongation events. The new position of the active tip will then be given by r = r + ( cos(ϕ ), sin(ϕ )). Note that for δϕ e = 0 the rotational diusion vanishes and the random walk becomes innitely persistent. (ii) When an active tip with local angle ϕ undergoes a branching event, it will produce two new active tips at positions r 1 = r + ( cos(ϕ 1 ), sin(ϕ 1 )) and r 2 = r + ( cos(ϕ 2 ), sin(ϕ 2 )), respectively, and will become inactive itself. The local angles ϕ 1 and ϕ 2 of the two new active tips respectively take values uniformly distributed in [ϕ + δϕ e , ϕ + δϕ b ] and , the two new tips will be located on the two dierent sides of the line determined by the polarity vector p = ( cos(ϕ), sin(ϕ)) of their parent while preserving a minimal angle of ∆ϕ min = ϕ 1 − ϕ 2 = 2δϕ e = π/5 between each other. The latter rule enforces a minimal distance between the two new tips to reduce the frequency of immediate annihilation of the two new tips. Furthermore, we also found that 85% of bifurcation angles obtained from the experiments were also distributed within the range [π/5, π], which justied this choice of the simulation setup. Supplementary Fig. 2(A) illustrates the elementary steps implemented in the simulation for the branching, elongation and annihilation of the active tips.
To test our analytical predictions for the free diusion of branch segments, we rst ran simulations for single branches that perform elongation and branching jumps without generating additional new active tips. We found that the mean-squared displacements of the nal local angle ϕ of single branches at time t max closely follow the relation ∆ϕ(t max ) 2 = 2Dt max with the diusion constant D given by the microscopic expression Eq. (S18), see Supplementary Fig. 2(B). We also obtained the same relation for the mean-squared displacements of the angle dierences, i.e., ∆ψ(t max ) 2 = 2Dt max for large values of t max , see Supplementary Fig. 2(C). The latter results justied our hypothesis that for suciently large times, the instantaneous jumps in the angle dierence values ψ are dominated by the changes in the local angle ϕ, and thus do not depend on the angle θ to origin.
For the full BARW simulations, we could then set a branching probability p b , the initial position r 0 , and the initial polarity vector p 0 , and produced branching networks that radially grow in all directions over time with active tips forming a front, leaving inactive branches of constant density in the inner regions. When the network growth is constrained within a spatial region delineated by xed boundaries, we observe an apparent directionality as predicted previously [8], see Supplementary We concluded that one needs to dene an external eld to guide the tips for a directed (anisotropic)

BARWs with external guidance and/or self-avoidance
In this subsection, we describe dierent possible microscopic ways to implement external guidance in the simulations of branching and annihilating random walks, and show importantly that the two give rise to the same qualitative predictions in the continuum model in terms of overall branch alignment.

Microscopic mechanisms for implementing external guidance
Observations of the movement of single neuronal tips in the presence of chemical cues (such as Netrin) indicate that the neuronal tips can reorient themselves even without neighboring tips or other branches in close proximity [9]. This change in the directionality can be eectively modelled as a displacement of the tip towards the external eld, as we discussed in Supplementary Note 1. Eq. (S20). We will therefore discuss these two dierent ways of implementing the external eld, and show that they generically result in the same behavior and the angular alignments follow the same statistics up to a constant prefactor for the eld strength.
External eld via biased probabilities. To closely follow the continuum model, we rst implemented biased jump probabilities in the simulation as follows: For each elongation event, we weighted the elongation probability p e = 1 − p b by the bias terms such that a forward and backward step that respectively increases or decreases the angle dierence value ψ occurs with the probability p e A(ψ) External eld via tip displacement. We also wanted to explore the modelling of the external eld in a dierent way, to show the generality of our approach and the insensitivity to the details of the microscopic tip behavior, and thus implemented guidance via tip displacements. In this simulation setup, for a radial external eld, the active tips would now be displaced at each time step by a factor f c p c , where p c ≡ ( cos(θ), sin(θ)) is the unit vector pointing along the eld lines, and θ denotes the angle to origin of the tip before displacement. To explore an axial eld, e.g. oriented along the positive x axis, one can dene the unit vector for the eld polarity to be p x c ≡ ( , 0). In this implementation of the external guidance, the eld thus directly corrected the directionality of the tip migration, as opposed to the eective correction which arose from the biased transition probabilities of the previous approach. Note that, the branch segments undergo a small rotational diusion upon elongation with a noise determined by the elongation angle δϕ e , as illustrated in Supplementary Fig. 2A. The displacement by the external eld is performed subsequently on the tips such that branch orientations are not completely deterministic. Even though these two ways of implementing the external eld are quite dierent, we found that they generated both qualitatively and quantitatively very similar network structures, see Supplementary Note 2.5 for a comparison of the results. Crucially, we observed that one only needs to tune the external eld strength f c by a constant prefactor α to obtain similar results using the two dierent algorithms.

Implementing self-avoidance
An additional feature observed in some branching tissues is the presence of self-avoidance of growing tips, e.g., as clearly evidenced in morphogenesis of starburst amacrine cells [10], where active tips sense, and are repulsed by, existing branch segments in close proximity. For neuronal branching it was shown that this self-avoidance of the growing dendrites is based on isoform recognition of the membrane proteins [11]. We modelled this eect by considering that active tips sense an average density vector p s [8] depending on the number of branch segments within a certain radius R s of self-recognition: Here we denote the position of the active tip at time t by r(t) and that of the remaining particles within R s by r j (t). Branch segments (particles) that belong to the same branch do not contribute to the sum, i.e. the active tip does not sense inactive segments from its own branch. The density vector p s will thus dene a normalized vector pointing away from nearby particles. We can now dene a self-recognition force by −f s p s , where the dimensionless parameter f s determines the strength of interaction with f s > 0 and f s < 0 corresponding to attraction and repulsion, respectively. Note that this self-repulsion force can eectively correspond to sequential recognition vs. retraction events such as seen in neuronal dendrites and as modelled e.g. in Ref. [12]. In contrast with the latter study, here we do not focus on the statistics of the retraction events and could therefore describe the self-recognition via a more simplied eective modelling. The position of the active tip after displacement by this self-recognition force will then simply change to r (t) = r(t) − f s p s . However, in order to preserve the step length , we will also correct this position such that the distance between the nal position r * (t) after displacement and the position at the previous time step is equal to , i.e., |r * (t) − r(t − τ )| = .
Note that, the new orientation of the branch segment does not align with p s but with the vector r * (t) − r(t − τ ). Supplementary Fig.5(A) illustrates the quantitites and displacement rule used in the implementation of self-avoidance.

Estimating the average opening angle
As briey discussed in the main text, the area of invasion of a branching network in a radial external eld will be proportional to the average opening angleθ between the outermost branch segments of the network. To quantify this, we conjectured that the overall opening angle is very likely to be dependent on the local dynamics at the growing boundary of a branching network. Indeed, using simple geometric arguments for the changes in the position and angles of the active tip (sine rule for the triangles dened by the origin and tip positions before and after tip displacement), similar to our analysis in Supplementary Note 1.1.1, we infer where ∆ϕ ≡ ϕ − ϕ * is the change in the local angle after the displacement by the external eld f c at time t, r(t) and r * (t) denote respectively the radial distance of the active tip before and after the displacement by the external eld, and δθ ≡ θ For the dynamics orthogonal to the growth direction that we are interested here, the changes in the radial distance between consecutive time steps can be neglected and thus described eectively by a single parameter r replacing the temporal parameter t. After these simplifying assumptions (and assuming δθ/dr is a well-dened innitesimal), we can integrate Eq. (S32) to obtain where we used the symmetry for the two growing boundaries of the branched network (prefactor 2).
Arguing that ∆ϕ is likely to be dominated by a small constant for large times, we can approximate the prefactor sin(∆ϕ) χ with a free parameter χ. For large times, the second term in Eq. (S33) inversely proportional to r becomes smaller and thus we arrive at Eq. (3) of the main text.

Simulation results for BARWs with self-avoidance 2.4.1 Strong self-avoidance without external eld
To clarify whether directed growth can arise in a self-organized way in the absence of an external eld just by the self-avoidance mechanism, we performed simulations with a strong self-avoidance of f s = −0.3 but setting the external eld strength to f c = 0. In the density regime set by a branching probability of p b = 0.1, we again observed isotropic growth of the networks, see Supplementary   Fig. 3(C). However, due to the the rather strong self-avoidance potential, the nal network exhibited a qualitatively much ordered structure, in agreement with the predictions in [8]. Interestingly, the active tips formed a sharper propagating front in contrast to the case without self-avoidance. We nally analyzed the low-density case by reducing the branching probability to p b = 0.02 to see if the network could preserve its pre-dened directionality uninterrupted by the outwards push due to branching events. However, such a self-organized anisotropic growth also could not be obtained for the low-density case, see Supplementary Fig. 3(D), which indicated that an external eld is required to break the isotropy in the tissue growth. To estimate the opening angleθ of a branched network (left), we focus on the active tip at the boundary of the network (orange tip in the boxed region): Tip displacement by the external eld with strength f c changes the local angle ϕ of the tip by a factor ∆ϕ and its distance to origin changes from r(t) to r * (t) (dashed lines to the tip). The angle to origin θ at time t after displacement is then given by θ * (t) = θ(t − τ ) + δθ, where θ(t − τ ) denotes the angle to origin at the previos time step t − τ . For a tip suciently distant from the origin, we can approximate r = r(t) = r * (t) = r * (t − τ ). Using the sine rule for the triangles dened by the origin and the tip positions before and after displacement, we obtain the relations leading to Eq.(S32).

Strong self-avoidance with external eld
Here we provide further results for simulations with a high self-avoidance potential and in the presence of external guidance. For increasing values of f s , we observed that the density and the space-lling properties of the network markedly increased, as shown in Fig. 5 in the main text. We also qualitatively observed that the networks with high self-avoidance exhibit better alignment with the external eld.
This increased alignment is indeed reected in the angle distributions, where the distributions for the angle dierence become narrower compared with those from simulations without self-avoidance, see

Large radius of self-avoidance
To explore the alternative of self-interactions that are controlled at larger distances (for instance e.g. longer retraction events after contact-mediated self-recognition) we set a large radius of self-recognition R s = 8 as compared with the radius R s = 3 we used otherwise. Interestingly, even on a visual level, we could observe strongly aligned branches in the morphology diagram, see Supplementary Fig. 6(A), which became quite pronounced for a large value of self-repulsion (f s = −0.2, rightmost column).
Analysis of the angle distributions revealed that this case indeed leads to large deviations from the analytical predictions, and the standard deviations decay also rather slowly with increasing f c , see Supplementary Fig. 6(B) and (C). For large f c , these eects of self-repulsion become rather dominant and lead to a saturation of the mean opening angles θ of the networks as shown in Supplementary   Fig. 6(E). Together, these results indicate that such long-range self-interaction eects can indeed inuence the angular alignment and morphology of the nal networks, especially in the presence of strong external potentials (although our data indicates that this it not the experimentally relevant limit, see  The two parameters Σ and Ω will then determine the new position of the tip after a branching or elongation event: (i) For a tip undergoing an elongation, we rst choose a random value for Σ = Σ e uniformly distributed in [0, π/10]. This thus xes the radius for the circle of the sphere where the new tip position will be located. We then pick a random value for Ω uniformly distributed in [0, 2π] to determine the position of the new tip on the chosen circle. Note that, for Σ = 0 the tip elongates while preserving its previous orientation. (ii) For branching, we now pick a random value for Σ = Σ b uniformly distributed in [π/5, π/2] because we nd that circles with an opening angle Σ smaller than π/5 lead to frequent annihilation events of newly born progenies upon branching. After xing Σ b we again choose Ω = Ω 1 uniformly distributed in [0, 2π] to determine the position of one of the new tips on the chosen circle. The second tip is then placed opposite the rst one by xing an angle parameter Ω = Ω 2 = Ω 1 + π, see Supplementary Fig. 8 for an illustration of these rules.
In the absence of external guidance, these rules then lead to isotropically growing branched networks in 3D, providing a similar phenomenology as seen in 2D. To explore the eects of external guidance, we now implement additional tip displacements along the eld polarity at every time step, in analogy to our implementation described in Supplementary Note 2.2.1 for a 2D geometry. For the simplest case of axial (linear) guidance, e.g. for a eld oriented towards the positive x-axis, we can dene additional tip displacements determined by a factor f c p x c where p x c ≡ ( , 0, 0) is the unit polarity vector of the external eld. After the displacement we update the branch position such that its length with respect to its previous position is conserved as , and transform back to spherical coordinates (Φ, Θ) in the local frame of reference of the active branch to obtain the local orientation of the branch segment after tip displacement. Finally, self-avoidance can be implemented in complete analogy to 2D, with the only modication being that the tip can sense neighboring branch segments within a sphere of radius R s around it, which will then lead to an additional displacement of the tip as described in Supplementary Note 2.2.2.

Results for BARWs with external guidance in a 3D geometry
We rst wanted to see if our main results for BARWs with external guidance could be reproduced qualitatively in a 3D setting. Indeed, building a morphology diagram for dierent values of external guidance strength f c (axial eld along x) and self-avoidance f s led to a qualitatively similar picture to the ones shown in Figs.1 & 2 in the main text: Increasing f c systematically decreased the overall territory of the branched networks, whereas large self-avoidance seemed to only increase the local density of branch segments, minimally inuencing the territory size, see Supplementary Fig. 9(A).
Next, to quantify these features, we generalized the opening angleθ which we used as a proxy for territory size in 2D by dening a right circular cone with its apex positioned at the origin of the branched network. The altitude vector of this cone (vector connecting the apex to the center of its base) is chosen to be aligned with the average polarity of the branched network p r ≡ j r j /| j r j |, where r j denotes the position vector of the j-th branch segment and j runs through all branch segments of the network. Using this construction, we could then dene an opening angleθ 3D of the cone that contains the entire branched network, where we determineθ 3D ≡ max(arccos(p r · r j )) after evaluating the dot product for each branch segment j. These estimated opening angles θ 3D (averaged over n = 30 simulations for each parameter choice) as a function of the external eld strength f c for vanishing (f s = 0) and strong (f s = −0.3) self-avoidance showed that, in complete analogy to 2D, self-avoidance indeed has a minimal inuence on the overall territory of the branched networks in 3D, see Supplementary Fig. 9(B).
To test the generalization of our analytical predictions on the branch orientations and scaling, we  π/2 Self-avoidance radius R s 3 , 8 Max. elongation angle δϕ e in 2D π/10 Annihilation radius R a 1.5 Branching angle Σ b in 3D U (π/5, π/2) Fit parameter χ for avoidance on the space-lling features. Interestingly, even though the fractal dimensions increased with increasing self-avoidance f s , this eect was weaker compared to the 2D case. Fractal dimensions did not change markedly for dierent values of the external eld strength f c , see Supplementary Fig. 9(E). We thus conclude that space-lling properties thus depend more sensitively on branching dimensionality.

Parameter values used in the simulations
In Table 1 we list the values for the diusion D and mobility µ coecients as predicted by our continuum model, see Eq. (S18), tting parameters and remaining parameters used in the simulation setup.
3 Supplementary Note 3: Experimental model system and methods

Coarse-grained reconstruction of neuronal laments
To obtain quantitative information on the branch numbers, branch lengths and probabilities of branching events, we needed to convert the manually reconstructed lament images into skeletonized datasets to extract coordinates. To achieve this, we rst conducted the following analysis using custom scripts in Python: The manual reconstruction images of lament and n borders were loaded and separated according to the channel information. After binarization, images were skeletonized with Lee's algorithm [13] implemented in scikit-image (v. 0.17.1) [14]. pre-dened stepsize , which had local angles with a ner distribution of values.
For the coarse-graining algorithm, we rst label all branching points in the network by identifying the three-valent vertices. We then start the coarse-graining loop by dening an active tip which is the vector closest to the origin. This active tip will evolve by scanning an area of a certain radius R for the underlying skeletonized data points and taking discrete steps (of average size s determined mainly by the scanned radius R) according to the following rules: (ii) Elongation: If there is no branching point in the scanned neighborhood, the active tip will try to move forward with respect to its current polarity by dening a polarity cone that provides a radial slice of the scanned neighborhood in the direction of the tip, and selecting the furthest data point of the underlying skeletonized network within the polarity cone in the next time step.
(iii) Branching: If the active tip itself is a branching point (which will necessarily be the case for a tip that undergoes the transition (i) above), it will search for two data points as progenies and produce two active tips at these coordinates in the next time step. The triangle connecting the two progenies and the active tip is required to have a minimal angle value at the vertex of the active tip in order to prevent branching events into the same branch of the underlying network.
(iv) Annihilation: If the above conditions are not fullled, the active tip will become inactive in the next time step, i.e., it will not be iterated further in the loop.
In some cases, two active tips start invading the same branch due to noisy regions in the original raw images. In order to prevent such events, we implement the additional rule that when an active tip moves on an already processed data point, it will immediately terminate. In the nal data processing step, we then revise the data of the corresponding branch such that it acquires the generation label and orientation of the older active tip that has a longer ancestral lineage.
This set of rules closely resembles the simulation setup for BARWs and provided a hierarchical tree for each sample from the experiments. After the coarse-graining loop, each processed datapoint was assigned a generation label and a local angle value determined by the vector connecting its positions at time t and in the previous time step t − τ . The resulting coarse-grained coordinates for the n = 8 networks analyzed here are shown in Supplementary Fig. 11 in black, where the underlying skeletonized networks are presented in blue. Coarse-grained coordinates of individual laments can be found as source data in Ref. [16].

Distribution of branch lengths and subtree topologies
A measure that is intimately linked to the probability of branching is the average branch length L j of a lament with j = 1, . . . , N , where N denotes the total number of branches in the tree. For a network with a high a priori branching probability, one would expect to nd short branches on average due to the frequent branching events. To quantify this, we estimated the branch lengths L j by calculating the total length Σ k s k of all branch segments (discrete steps of size s k ) between the starting and end points of each branch. The normalized histograms of branch lengths for the individual laments are displayed in Supplementary Fig. 12(A). Importantly, the combined data from n = 8 samples showed an exponentially decaying tail (see Fig. 2f in the main text) with an average branch length of L j 46, see Supplementary Fig. 12(B). This is what is expected in a stochastic branching process, validating a key assumption of the framework of BARWs. We also calculated the average step sizes corresponding to single steps of the coarse-graining algorithm described above, and obtained a mean step size of s 3.2 from the combined dataset, see Supplementary Fig. 12(C). The mean step size can be used to determine the normalized branch lengths L j / s 14 to compare the lengths with the simulation data (as shown in Fig. 2f in the main text).
Another signature of the BARW framework, as opposed to alternative theories of branching morphogenesis, that the subtree size structure of the networks should exhibit a strong heterogeneity with both large and small subtrees. The size of a subtree is dened by the total number of branches in it [17]. In agreement with this framework, tree topologies obtained from the experiment were quite non-stereotypical and showed a similar heterogeneity as those from simulations, see Supplementary

Estimation of the branching probability
A key parameter to determine for the comparison of our theoretical results with the experimental data is the branching probability p b of the networks. One can estimate p b from the measured distributions of branch lengths because the average branch length depends in general inversely on the branching probability, i.e. L j ∝ 1/p b . However, due to the frequent annihilation events the distributions of branch lengths in fact underestimate the average branch length, and thus lead to a high branching probability p b that generate networks qualitatively dierent from the experimental observations. We therefore decided to estimate p b by directly counting the number of branching events n b for each lineage and taking its ratio to the total number of branch segments (steps) n s for that lineage. The branching probability p b of a network then corresponds to the average ratio n b ns from all lineages of branches ending at a leaf. Supplementary Fig. 13(D) displays the normalized histograms for the ratio n b ns obtained from the combined dataset of n = 8 laments (D1) and the estimates of p b for the individual samples (D2). The average branching probability obtained from the combined data of p b 0.05 was nally used to generate the branching networks of the BARW simulations and to compare the results on angular alignment by determining the diusion and mobility coecients given in Eq. (S18).

Angle distributions and fractal dimensions for individual samples
To obtain the angle distributions we rst needed to clarify the coordinate of the origin for each individual sample. Because the initial branch of some laments was located in the spinal cord outside of the n region, we did not take the starting point of the initial branches as the origin. Instead, an alternative choice which had the advantage of being systematic was to rst locate the boundary of the n tissue with respect to the anterior-posterior axis, and then x a central point that has the same distance to the boundary for all samples from the same sh. After this identication, angles to origin θ and angle dierences ψ could be calculated for each node in the coarse-grained networks. The normalized histograms for the local angles ϕ, angles to origin θ, and the angle dierences ψ for the individual laments are shown in Supplementary Fig. 13(A-C). The histograms for the local angle ϕ and the angle to the origin θ attained rather irregular shapes, whereas the histograms for the angle dierences ψ for Finally, in Supplementary Fig. 14 we display the fractal dimensions for the individual samples obtained by the box-counting method. For this purpose, we selected box sizes gradually increasing from min 3 s to max 3 L j to describe a well-dened scaling behavior (i.e. to omit boundary eects arising from too large or too small boxes). The fractal dimensions showed a rather robust correlation with the inverse of the mean branch lengths and with the branching probability, even though these two measures were not strictly proportional for all samples, see Supplementary Figs. 14(B) and (C).