{'fraction_non_alphanumeric': 0.046085189112760695, 'fraction_numerical': 0.013108283256745611, 'mean_word_length': 4.208023318248063, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A model of fasciculation and sorting in mixed populations of axons

We extend a recently proposed model (Chaudhuri et al., EPL 87, 20003 (2009)) aiming to describe the formation of fascicles of axons during neural development. The growing axons are represented as paths of interacting directed random walkers in two spatial dimensions. To mimic turnover of axons, whole paths are removed and new walkers are injected with specified rates. In the simplest version of the model, we use strongly adhesive short-range inter-axon interactions that are identical for all pairs of axons. We generalize the model to adhesive interactions of finite strengths and to multiple types of axons with type-specific interactions. The dynamic steady state is characterized by the position-dependent distribution of fascicle sizes. With distance in the direction of axon growth, the mean fascicle size and emergent time scales grow monotonically, while the degree of sorting of fascicles by axon type has a maximum at a finite distance. To understand the emergence of slow time scales, we develop an analytical framework to analyze the interaction between neighboring fascicles.


I. INTRODUCTION
Reaction-diffusion phenomena arise in diverse fields such as physical chemistry [1] or developmental biology [2]. In certain reaction-diffusion systems, the process of path aggregation occurs, in which preferred paths of the diffusing elements are established and evolve in time. The path aggregation process is found in diverse realm of nature, e.g., in formation of insect pheromone trails [3][4][5] and human walking trails [6,7], in aggregation of trails of liquid droplets moving down a window pane, in river basin formation [8,9], etc.
One class of the mathematical models in which path aggregation processes have been studied is the activewalker models [7,10] in which each walker while moving through the system changes the surrounding environment locally, which in turn influences the later walkers. The ant trail formation is an example of such a process [3,4]. An ant leaves a chemical trail of pheromones on its path which the other ants can sense and follow. Evaporation of pheromone leads to an aging of these trails. Similarly, the mechanism of human and animal trail formation is mediated by the deformation of vegetation that generates an interaction between earlier and later walkers [6,7]. This deformation, and therefore its impact, decays continuously with time [7].
In a recent Letter [11] we analyzed the dynamics of path aggregation using a simple model that belongs to the class of active walker systems discussed above. In contrast to the active-walker models, in our model the individual paths do not age gradually, but rather maintain their full identity until they are abruptly removed from the system. This particular rule for path aging was chosen to allow application of our model to the process of axon fasciculation (formation of axon bundles [12]), which we discuss below and more in detail in Sec. II.
In order to develop neuronal connections, sensory neurons born in peripheral tissues project their axons (long tubular part of the neuron cell that conducts electrical excitations) towards target regions in the brain. Frequently, multiple axons come together to form axon fascicles, and may sort according to the cell type of the neuron to which the axon belongs. This fascicle formation and sorting can be driven by inter-axon interactions leading to, e.g., a pre-target spatial map in the mammalian olfactory system [13][14][15].
In our model, the axons are represented as directed random walks in two spatial dimensions. In Ref. [11] we formulated and analyzed the simplest version of the model, in which all axons belong to the same type and have strong adhesive interactions, so that each newly growing axon encountering an existing fascicle will join the fascicle and never detach. In the presence of axon turnover (aging of the paths), a steady state characterized by a distribution of fascicle sizes [11] is eventually established. The focus of Ref. [11] was on the analysis of the surprisingly long time scales that emerge from this simple dynamics.
In the current paper, we significantly extend this theoretical analysis. We develop an analytical description of the dynamics of two neighboring fascicles, and show how the slowest mode of their interaction gives raise to the slow time scales observed in Ref. [11]. We also systemat-ically discuss the limited analogies that can be made between our 2-dimensional model and 1-dimensional models of particle coalescence [16], aggregation [17], and chipping [18,19]. These analogies are useful for the understanding of stationary quantities of our model such as the distribution of fasicle sizes and the distribution of inter-fascicle separations.
The main contribution of the current paper, however, is to generalize the previous model of Ref. [11] to attractive interactions of finite strength (so that detachment of axons from fascicles is possible) and to multiple axon types with type-specific interactions. Such a generalization is necessary to allow the biological application of the model.
In the following section we give a detailed biological motivation for the model we consider. In Sec. III we introduce the model and the Monte Carlo (MC) simulation scheme that we use to investigate its properties numerically. Followed by this, in Sec. IV, we give a brief overview of guiding concepts that will recur in the rest of the paper. In Sec. V, we present a detailed analysis of the system containing axons of a single type that follows the "always attach, never detach" rule. We extend the numerical results of Ref. [11] for the properties of the steady state and for the emerging time scales. We review the analytical framework of single-fascicle dynamics, developed in Ref. [11], and significantly extend it by deriving results for the interaction dynamics of two neighboring fascicles. In Sec. VI we numerically study the effects of non-vanishing detachment rates of axons from a fascicle. In Sec. VII we discuss some limited analogies of our model to one dimensional aggregation and coalescence processes. In Sec. VIII we discuss the sorting of fascicles by axon types in a system containing two types of axons, the simplest manifestation of a mixed population of axons. Finally, we provide a summary of main results in Sec. IX and conclude in Sec. X by discussing the outlook for biological applications of our model.

II. BIOLOGICAL MOTIVATION
Sensory neurons located in peripheral tissues connect to more central locations of the nervous system via axons [20]. During the development of an organism, axons of newly maturing sensory neurons must establish connections to the proper location. Axon growth is initiated at the soma (main cell body) of each neuron, and proceeds with a typical rate of extension 100µm/h [21]. The direction of growth is controlled by the dynamic growth cone structure at the tip of the axon. The growth cone probes the environment in its vicinity, and can detect gradients of spatially distributed chemical signals. In the absence of strong directional signals, the path of the growth cone is highly stochastic [22,23], while in the presence of appropriate guidance cues, the direction of motion becomes strongly biased. The overall direction of axon growth may be guided by spatial gradients of chem-ical cues generated by the target. A number of distinct molecular guidance cues that influence neuronal development have been identified in recent years [24,25], and the response of the growth cone to graded cues has been studied theoretically [26][27][28]. In this work, we do not directly model the axon guidance by graded chemical cues, but subsume their influence into the setup of our model by giving all axons a common preferred growth direction.
In this article, we study the collective effects that arise from direct local interactions among the growing axons. When such interactions are attractive, the growth cone of a newly growing axon tends to follow the tracks (i.e., the axon shafts) of older axons. The strength of this interaction is governed by the type and expression level of the relevant cell adhesion molecules [21,29]. The resulting dynamics can lead to selective formation of fascicles of axons [30][31][32], a common and essential phenomenon in the developing nervous systems.
An additional important aspect included in our model is that of neuronal turnover. During development, a significant portion of sensory neurons with fully grown axons may die, and be replaced by younger sensory neurons which attempt a new connection to the brain. For example, up to 80% of retinal ganglion cell axons are lost during the development of the visual system in the cat [33,34]. In the mammalian olfactory system, both neuronal birth and death persist throughout the life of the animal, leading to a dynamical steady state pattern of connectivity. In particular, the average lifetime of an olfactory sensory neuron in the mouse is of the order of 1-2 months [35], which is less than one tenth of the mouse's lifespan.
To motivate the introduction of multiple types of axons into our model, we now briefly discuss the intricate connectivity pattern of the mammalian olfactory system, which implements the sense of smell. In the mouse, the adult nasal epithelium contains approximately 10 6 olfactory sensory neurons, which send their axons through the olfactory tract to the olfactory bulb in the forebrain. Remarkably, the sensory neurons belong to approximately 1200 distinct types [14,36], and the axons of each type connect to a distinct neuronal structure, a glomerulus, on each olfactory bulb [36,37]. Such precise connectivity is fully established only after several turnover periods, while in newborn mice, split glomeruli and glomeruli that mix several axon types are often observed [35,38].
In olfactory sensory neurons, elegant genetic analysis shows that the axonal type is determined by the expression of a specific odorant receptor gene [36,[39][40][41]. Physiological experiments on mice show that the expression of specific types of cell adhesion molecules, that dictates the strength of adhesive forces between axons, is strongly correlated with this axonal type [40,42]. A wide range of strengths of interactions between axons may be generated through combinatorial expression of multiple types of cell adhesion molecules.
In Fig. 1 we show configurations of olfactory axons as observed in in vivo [41] and in vitro [43] experiments.   [41]). Only axons belonging to one type of olfactory sensory neurons (expressing the M71 receptor gene) are labeled; axons of other types are present but not visible. The axons progressively fasciculate and terminate in a glomerulus visible in the center of each half-image. Fig. 1(b) shows fasciculation of axons growing from an explant of the rat olfactory epithelium ( Fig.7(a) of Ref. [43]). In this case, the fluorescent labeling does not distinguish the axonal type, and (with a high probability) the visible axons belong to multiple types.
Our model aims to provide a quantitative framework for evaluating the contribution of axon-axon interactions to the formation of patterns described above. The presence of turnover and multiple axon types in our model distinguishes our work from previous theoretical studies of axon fasciculation [44,45]. Our implementation of the individual axon dynamics is particularly simple, to allow us to concentrate on collective effects arising from interactions within a population of axons.

A. Setup and interactions
In our model, each growing axon is represented by a directed random walk in two spatial dimensions ( Fig. 2(a)). The random walkers (representing the growth cones) are initiated at the epithelium (y = 0, random even x) with a birth rate α, and move towards the bulb (large y) with constant velocity v y = 1. In the case of multi-type systems, a type is assigned to each newly initiated random walker (specifically in the simulations of Sec. VIII, the type is decided randomly with equal probability for each of the two types). The trail generated by a random walker (growth cone) is regarded as an axon shaft. A forward moving directed random walker (growth cone) interacts with trails (axon shafts) of other walkers. In the numerical implementation on a tilted square lattice, at each time step the growth cone at (x, y) can move to (x − 1, y + 1) (left) or (x + 1, y + 1) (right). The probability p {L,R} to move left/right is evaluated based on the axon occupancies at the (x±1, y +1) sites and their nearest neighbors (x ± 3, y + 1) (see Fig. 2(a)). At a given y, two axons are considered to be part of the same fascicle if they are not separated by any unoccupied sites (i.e., they are not separated by more than two lattice spacings).
We assume a short-range attractive interaction between each growth cone and the close-by axon shafts. The range of interaction (two lattice spacings in our model) corresponds biologically to the range of extension of sensory filopodia from the growth cone (of the order of 10µm). The attractive interaction is mediated by cell adhesion between the growth cone and the axon shafts. We assume that the interactions are additive and typespecific. For a given growth cone, the model assumes a weak nearest neighbor attraction E o < 0 if the neighboring axon shaft is of a different type and a stronger attraction E h (< E o ) if the neighboring axon shaft is of the same type. In each time step, a growth cone at (x, y) attempts a Monte Carlo move to the left (x − 1, y + 1) and to the right (x + 1, y + 1) with probabilities 1/2. The moves are accepted with probabilities (δE r = −δE l ) and the occupancy number n h denotes the number of axons belonging to the same type as the growth cone, while n o is the number of axons of other types. Notice that in calculating the difference in energy, the occupancy of the positions (x ± 1, y + 1) does not appear, as their contributions to the energy cost mutually cancel. Periodic boundary conditions are used in the xdirection. In this kinetic MC scheme, we used parallel updates of all random walkers in each MC step.
We now clarify the relation of the general model described above to the model studied by us in Ref. [11]. In this simple version of the model, all axons belonged to the same type and the interaction between them was governed by the "always attach, never detach" rule, which is the "zero temperature" version of the above-mentioned general MC scheme, i.e., the dynamics was a pure energy minimization process: p L = 1 (p R = 1) when δE l < 0 (δE r < 0); p L = p R = 1/2 in all other cases. This is in contrast to the "finite temperature" dynamics of the general model, in which there is a non-zero rate for the detachment of growth cones from fascicles. Ref. [11] used sequential updates, in contrast to parallel updates used in this paper. An additional difference is that in Ref. [11], the interaction was not assumed to be additive, i.e., the strength of interaction of a growth cone with a fascicle did not depend on the number of axons in the fascicle.
Note that two fasciculated axons run parallel to each other with their separation in x-direction restricted within two lattice spacings, the interaction range. Thus the typical width of a fascicle containing multiple axons remains 2-3 lattice spacings in our model. We do not implement any on-site repulsion, i.e., axons are free to grow on top of each other.
In our model, we do not consider any relaxation dynamics of axon shafts. This corresponds to the assumption of strong adhesion of axons to the substrate, so that the line tension on axons can not straighten out the local curvatures.

B. Turnover
To capture the effect of neuronal turnover, each random walker is assigned a finite lifetime θ from an exponential distribution of lifetimes Π(θ)dθ = 1 T exp(−θ/T )dθ with mean θ = ∞ 0 θ Π(θ)dθ = T . When the lifetime expires, the random walker and its entire trail (i.e., the whole axon) is removed from the system. The mean number of axons in the system reaches the steady-state value N (y) = N 0 exp(−βy), where β = 1/T is the mean death rate per axon and the steadystate occupancy at y = 0 is N 0 = α/β. In the simulations, we use T = 10 5 time steps, and restrict ourselves to y ≤ T /10. The birth rate α is chosen so as to obtain the desired number of axons N 0 , or equivalently, the desired axon density ρ = N 0 /L (ρ = 1/2 implies an average occupancy of one axon per site), where L is the system size in x-direction. As we will show in detail, the time scales needed to achieve the steady state of fascicle size distribution can be very long compared to the time scale T needed to achieve the steady state value of the total number of axons. A typical late-time configuration for a system of axons involving two distinct axonal types with type-specific interactions is shown in Fig.2(b).

C. Parameters
In this section we briefly discuss the biological meaning and physical values of the parameters in our simulations. In our model the interaction range in x-direction is chosen to be 2 lattice units. Since we assume only contact interactions, this range of 2 lattice units corresponds to the length of a filopodium which typically is 10 µm. Thus the lattice spacing in x-direction ∆x = 5 µm.
The time step ∆t in our model needs to be large enough to allow the growth cone to integrate a signal and react to it. Ref. [23] suggests this time scale to be of the order of tens of seconds; we choose ∆t = 60 s. This corresponds to a diffusion constant (in x-direction) (∆x) 2 /2∆t = 12.5 µm 2 /minute which compares well with in vitro observation for short time scales up to tens of minutes [22]. Note that with this choice of ∆t the mean lifetime of an axon T = 10 5 ∆t used in our simulations corresponds to 69.4 days, quite typical of axons of mouse olfactory sensory neurons [35].
The lattice unit in y-direction is now chosen to give a reasonable growth velocity v y . Choosing ∆y = 1 µm, we have v y = 60 µm/hour, which is a typical value for growing axons of sensory neurons [21].
Note that with the above mentioned choices, Fig. 2(b) corresponds to a system size of 800 ∆x = 4 mm in xdirection and 10 4 ∆y = 10 mm in the y-direction. These dimensions are comparable to the size of the olfactory bulb in mice [36].
The effective interaction energies E h and E o should be chosen to match the observed rates with which axons detach from fascicles. We introduce the quantity π d = exp(E)/∆y that expresses the rate of detachment of one axon per unit length of a two-axon fascicle given that the two axons interact via E = E h or E o . Thus a growth cone interacting with a fascicle of n axons will follow the fascicle over the mean distance L y = π −n d before it detaches. It is not straigthforward to use published experimental images to deduce E as usually, the size of the fascicle is not known, and the location at which a growth cone first attached to the fascicle is not recorded. The observed typical distance L y varies widely depending on the specific neural system, ranging from tens of µm to centimeters. In our simulations, we use the range of homotypic interaction strength E h = −4 to -1, which corresponds to detachment rates π d = 0.02 to 0.37 µm −1 . In Table I we list the meaning and values of the parameters used in our model.

IV. OVERVIEW
We show in this paper that simple directed growth of axons that interact via a short-range attraction leads to reliable formation of axon fascicles, in absence of any external chemical guidance cue. Once a fascicle is formed its position does not move appreciably, however, the fascicle size (number of axons present in the fascicle) fluctu- ates. The turnover of individual axons generates a slow dynamics of reorganization of fascicles at a fixed y-level.
In the simplest case (Sec. V), the system contains only a single type of axons which grow and form fascicles via an energy-minimizing dynamics (strong inter-axon interaction), so that once attached the axons do not leave a fascicle. For this case one can uniquely assign a basin of each fascicle at any specified y-level (see Fig. 2(c)). The basin size D of a fascicle is the interval at the level y = 0 between the right-most and left-most axons belonging to the fascicle (Fig. 2(c)). Any axon growing from within this basin must contribute to the fascicle size unless it dies before reaching the specified y-level. Thus the average number of axons that survives in a fascicle at level y (in the steady state) isn = D ρ exp(−βy). The axons initiated at the opposite edges of the basin are expected to meet each other in y ≃ (D/2) 2 steps of random walk in x-direction. Therefore, one obtains the mean-field prediction for the mean fascicle sizen(y) ≃ 2y 1/2 ρ exp(−βy) up to y ≃ (L/2) 2 , where complete fasciculation occurs, i.e.,n = N (y).
In a system with finite detachment rates (Sec. VI), however, growing axons can leave one fascicle and attach to another. Thus the fascicle basins overlap and axons introduced in the basin of one fascicle can end up in a different fascicle. Still the mean-field estimate of the increase of mean fascicle size with y, shown above, turns out to remain approximately valid (Sec. VI A). For higher detachment rates (weaker interactions), the prefactor of the power-law growth is reduced, corresponding to smaller fascicles. In the limit of extremely weak interactions, each axon would grow independently of the others, and no fasciculation is possible. This overall picture remains intact even for systems having multiple axon types. The dynamic steady state is characterized by a position (y) dependent distribution of fascicle sizes which shows a scaling law. The peak of the distribution shifts towards larger fascicle sizes at higher y-levels (Sec. V B and Sec. VI B).
In a steady state configuration at fixed time t (such as in Fig. 2(c)) the axon fasciculation with increasing y may be formally viewed as the evolution of an one-dimensional diffusion-aggregation or diffusioncoalescence process [46]. These limited analogies can be used to approximately understand steady state properties like distribution of fascicle sizes and inter-fascicle separation (Sec. VII). However, the full dynamics of our model has no simple one-dimensional counter-part [11]. A projection of the dynamics onto one dimension would involve complicated long-time correlations between the random walkers.
The dynamics of fascicle reorganization can be characterized by the slow approach to steady state, or by the steady state auto-correlation time for the mean fascicle size. These time scales grow with y and may reach values orders of magnitude larger than T . In absence of detachment, the slowest mode of fascicle reorganization occurs via partial exchange of neighboring basins (Sec. V D).
In the presence of detachment, in addition, the basin of one fascicle can easily drain to another. Thus the time scales decrease with decreasing inter-axon attraction (Sec.VI A).
In systems containing multiple types of axons with type-specific interactions, we evaluate the degree of sorting S that quantifies the type-wise purity of the local environment of axons (Sec. VIII). At steady state, S(y) shows a non-monotonic variation, with a maximum at some intermediate y. The position and the value of the maximum depend on system parameters like the mean density of axons ρ and the interaction strengths. This non-monotonicity is due to the attractive heterotypic interaction which merges mid-sized, relatively pure fascicles to form large impure fascicles.

V. SINGLE TYPE OF AXONS, NO DETACHMENT
In this section we analyze the collective behavior of axons belonging to a single type following the energyminimizing "always attach, never detach" rule. This model has been investigated in detail in a previous publication [11]. In this section, we extend the numerical and analytical results of Ref. [11]. In the simulations, we use a modified implementation of the Monte Carlo update rules. In contrast to Ref. [11], the strength of interaction between a growth cone and a fascicle is assumed to be proportional to the number of axons present in the fascicle. Another difference is that we use parallel updates, instead of sequential updates used in Ref. [11]. We therefore include a comparison to the main results we reported in Ref. [11], to show that these are not altered. A. Approach to steady state: mean fascicle size and time scale A typical late-time configuration for a system with L = 800 and N 0 = 100 (density ρ = N 0 /L = 1/8 at y = 0) is shown in Fig. 2(c). With increasing y, the axons aggregate into a decreasing number of fascicles. The number of axons in a fascicle is referred to as the fascicle size n. At steady state, the mean fascicle sizen at level y may be estimated asn ≃ 2ρy 1/2 exp(−βy) up to y ≃ (L/2) 2 , where complete fasciculationn = N (y) is expected [11].
The measured mean fascicle size, obtained by averaging over all the existing fascicles at a given y (Fig. 3), where τ ap (y) defines the time scale of approach to the steady state value n ∞ (y). The same behavior was observed earlier in simulations reported in Ref. [11]. The semi-log plot in the inset of Fig. 3 shows clearly the slow exponential approach to the steady state mean fascicle size. Using the above-mentioned double-exponential fitting we extract the time scale τ ap and the steady-state mean fascicle size n ∞ at all the y-levels.
The approach-to-steady-state time scale τ ap increases with y. τ ap can exceed the mean axon lifetime T by orders of magnitude ( Fig. 3 and 4). Note that τ ap is longer in a system with larger density of axons ρ (Fig. 4). Ref. [11] discussed this point in detail. Further, asymptotically in y, we find n ∞ = c + 2ρy b exp(−βy), with b ≈ 1/2  Impact of interaction range: To test the impact of the range of inter-axon interaction, we have simulated a similar system with purely "contact" interaction, i.e., the interaction range is taken to be zero. With this reduction in the range of interaction, we find that the emerging time scales decrease. For instance, for a system of L = 800 and N 0 = 100 the approach-to-steady-state time becomes τ ap 10T . Thus an increase in the range of interaction increases the emerging time scales. The steady state distribution of fascicle sizes shows the same scaling behavior as in the case of nearest-neighbor interaction discussed in the following.

B. Steady state
The steady state is characterized by the stationary distribution of fascicle sizes P s (n, y), defined as the number of fascicles of size n at level y.

Steady state: scaling regime
For a system with L = 800 and N 0 = 100, P s (n, y) is shown at a series of y-levels in Fig. 6. Within the range y = 10 3 − 10 4 all data collapse onto a single curve after appropriate rescaling (Fig. 6). This data collapse implies the scaling law [11] P s (n, y) = n(y) −r φ(n/ n(y) ) with r = 2.1 and the scaling function φ(u) = N u exp(−νu − λu 2 ). Note that the steady state averaged fascicle size n(y) is a quantity equivalent to the asymptotic n ∞ (y) discussed in the previous subsection. The scaling law in Eq. 1 can be justified starting from the assumption of homogeneity of fascicle size distribution P s (n, λy) = λ −p P s (λ −q n, y), with the exponents p and q undetermined at this stage. Noting that the mean number of axons N (y) = dn n P s (n, y) and the mean number of fascicles B(y) = dnP s (n, y), the homogeneity condition leads to the relations N (λy) = λ −p+2q N (y), B(λy) = λ −p+q B(y). Since by definition the mean fascicle size n(y) = N (y)/B(y), n(λy) = λ q n(y) . Invoking the mean-field prediction n(λy) = λ b n(y) with b = 1/2 (Sec. IV and V.A), we find q = b. If N (y) were independent of y, we would have had p = 2b. However, in fact N (y) = N 0 exp(−βy). In the region βy < 1, we can write p = 2b + δ with δ ≈ βy/ ln y ≪ 2b. Note that the relation P s (n, λy) = λ −p P s (λ −b n, y) can be recast in the form P s (n, y) = n −r φ(n/ n ) where r = p/b = 2 + δ/b 2, in agreement with Eq. 1. As δ is y-dependent, the scaling of P s (n, y) is only approximate.

Steady state: crossover to complete fasciculation
The steady state distribution changes its shape drastically beyond the scaling regime. Near y = (L/2) 2 , on an average all the axons are expected to collapse onto a single fascicle, thereby generating a distribution sharply peaked at n = N (y). To demonstrate this fact we take a system of small size and high density ρ = 1/2 (L = 100 and N 0 = 50). Within the scaling regime (y 500) the distribution function maintains the scaling form u exp(−νu − λu 2 ). However, at higher y-levels the distribution becomes bimodal with a new maximum appearing, characteristic of the complete fasciculation. This shows a coexistence of two preferred fascicle sizes. Finally, at y ≃ (L/2) 2 the whole weight of the distribution shifts to this new maximum and the distribution becomes unimodal again (see Fig. 7). In systems with larger L, the regime of coexistence shifts towards higher y-levels.

Steady state: correlation time
The dynamics in the steady state is characterized by the auto-correlation function for the mean fascicle sizē n(t) at a fixed y-level: c(t) = n(t)n(0) which fits to the form p + q exp(−βt) + r exp(−t/τ c ) (as in Ref. [11]). The correlation time τ c increases with y and significantly exceeds the axon lifetime T . We show this behavior for systems at ρ = 1/2 (L = 100, N 0 = 50) in Fig. 4. This shows a regime of approximate power law growth of the time scale τ c ∼ y 2b with b ≈ 1/2.

C. Effective single-fascicle dynamics at fixed y
In this subsection, we review our analytical results from Ref. [11]. The following subsection presents new results for the time scales arising from the interaction of two neighboring fascicles.
The concept of effective single-fascicle dynamics has been introduced in Ref. [11]. The dynamics of a mean fascicle at level y with n(t) axons can be viewed as a stochastic process with gain rates u + (n) (for transitions n → n + 1) and loss rates u − (n) (for transitions n → n − 1). A fascicle loses axons only by the death of individuals, therefore, u − (n) = βn [11].
In absence of detachment events, any axon introduced within the basin (size D) of a fascicle (see Fig. 2(c)) can not escape the fascicle. Moreover, some of the axons born in the neighboring inter-basin gaps (size E) eventually join the fascicle under consideration. These two processes contribute to u + (n). As was shown in Ref. [11], the time series of D(t) and n(t) tend to co-vary. Thus treating the dynamics of D as slave to n, we get a form u + = a + bn [11]. Note that the basin size D can not exceed 2y or L, and D > 2y 1/2 occurs with low probability. Therefore a saturation of u + (n) is expected for large values of n. The measured average gain and loss rates u ± (n) obtained from our current simulations agree with the functional forms u + (n) = a + + b + n − c + n 2 and u − (n) = βn (data not shown). The quadratic correction to linear growth captures the saturation of u + (n) at large n.
The master equation of the growth-decay process for the effective single fascicle of size n at level y may be written aṡ P (n, t) = u + (n − 1)P (n − 1, t) + u − (n + 1)P (n + 1, t) − [u + (n) + u − (n)]P (n, t), for n > 1. For the boundary state (n = 1) where J + (y) represents the rate with which new single axons appear between existing fascicles at y. The solution of the master equation at steady state was derived in Ref. [11] and has the form, where γ = a + /β − 1 and ℓ = 1 − b + /β and κ = c + /2β.
From the master equation one can estimate the approach-to-steady-state time τ ap and the correlation time at steady state τ c [11]. The correlation time τ c for the fascicle size n, near the macroscopic stationary point n s [u + (n s ) = u − (n s )] can be expressed [47] as τ c = 1/(u ′ − (n s ) − u ′ + (n s )) = 1/(β − b + + 2c + n s ). Under the linear approximation of u + (n) = a + + b + n the approach-to-steady-state time scale for the average fascicle size n can be written as τ ap = 1/(β − b + ) [47].
Further, the mean lifetime of fascicles can be defined as τ f = [ ∞ 1 P s (n, y)dn]/J + (y) and is evaluated to obtain [11] Notice that the above derivations of the time scales already involved the numerical observation u + (n) = a + + b + n − c + n 2 . Using the y-dependence of b + and c + obtained from numerical simulations, we found power law growth of the time scales, τ c ∼ y b , τ ap ∼ y b and τ f ∼ y 2b with b ≈ 1/2 [11].
In the following subsection, using a purely analytical, deterministic treatment of the dynamics of two neighboring fascicles, we show that a time scale growing as y 2b emerges due to an exchange of basin size between the fascicles.

D. Effective dynamics of two interacting fascicles
In this section we analytically explore the effective dynamics of two neighboring fascicles. Under the "always attach, never detach" rule the basins of neighboring fascicles do not overlap. Three dynamical variables characterize the dynamics of the fascicles: (i) the number of axons n i present in the fascicle, (ii) the basin size D i and (iii) the separation E i between the i-th and (i + 1)-th basin (Fig. 2(c)). We explore the dynamics at a fixed y. We consider the thermodynamic limit of large N 0 and L (with a fixed density ρ = N 0 /L), and where δ = 1/N 0 . First we describe the gain and loss terms in the dynamics of η i . Axons born inside the basin of a fascicle contribute to the increase in η i , hence the term βλ i (we used N 0 = α/β to express axon birth rate α in terms of β). A fascicle can lose axons only by individual axon deaths, thus the loss term βη i . Any axon which is born in the inter-fascicle empty spaces ǫ i−1 and ǫ i ends up in either of the two neighboring fascicles with probability 1/2, hence the gain term (1/2)β(ǫ i−1 + ǫ i ) [53]. Next we consider the dynamics of basin size λ i . A new axon can be born in the gap ǫ i with a rate α(ǫ i −ã) and attach to the i-th fascicle with probability 1/2. If it attaches it contributes half the gap size ǫ i /2 towards the basin size λ i . Hence the gain term (α/4)ǫ i (ǫ i −ã). A similar contribution to the gain in the basin size comes from the other neighboring gap ǫ i−1 . The death of a boundary axon reduces the basin size by an amount δλ i /(η i − δ) (assuming no double occupancy at a lattice point in the y = 0 level). The contributions of this loss coming from two boundaries add up in the total loss term 2 × βδλ i /(η i − δ).
Finally, we consider the dynamics of the inter-basin gaps ǫ i . The death of boundary axons of neighboring fascicles i and i+1 that border the i-th gap ǫ i contributes to the gain in the gap size. Thus the gain terms βδλ i /(η i −δ) and βδλ i+1 /(η i+1 − δ). Birth of an axon in the gap reduces the gap size under consideration. The rate of such an axon birth is α(ǫ i −ã) and on average this event reduces the gap size by an amount ǫ i /2. Thus the loss term (α/2)ǫ i (ǫ i −ã).
We use a periodic boundary condition, such that the last fascicle is a nearest neighbor of the first fascicle. These equations obey the constraint of overall constant size i (λ i + ǫ i ) = 1. It is important to note that the cutoffã can be taken to zero meaningfully only after solving the differential equations.
For the simplest non-trivial case involving two fascicles, the steady state that follows from these equations is characterized by η 1 = η 2 = 1/2, λ 1 = λ 2 = 1/2 −ã/2 − 2δ −ãρ/16 and ǫ 1 = ǫ 2 =ã/2 + 2δ +ãρ/16. We perform a normal mode analysis for small deviations from this steady state. The constraint i=1,2 (λ i +ǫ i ) = 1 implies that there are only five independent deviations, δλ 1 , δǫ 1 , δǫ 2 , δη 1 and δη 2 . The linear stability analysis about the steady state shows that all the five possible modes are stable. Among them, four modes are shortlived. For them the deviations decay extremely fast with rates ∼ β. However, the fifth mode which in the largesize limit can be written as (δn 1 = −1, δD 1 ≈ −1/ρ, δE 1 = 0, δn 2 = 1, δD 2 ≈ 1/ρ) takes a long time to decay (see Fig. 8). It involves the loss of a boundary axon of one fascicle, which shrinks its basin size by δD 1 ≈ −1/ρ, and simultaneously a gain of a boundary axon for the other fascicle, increasing its basin size by the equal and opposite amount δD 2 ≈ 1/ρ. This operation leaves the inter-basin gap unchanged (δE 1 = 0) and can be viewed as an exchange of basin space (Fig. 8). The deviations from steady state in this mode decay over a very long time scale τ ≈n 2 /3β wheren is the steady state value of the fascicle size (n 1 = n 2 ≡n).
Using the approximate growth of mean fascicle sizen ∼ y b with b = 1/2, we find a power law growth of this time scale τ ∼ y 2b . Notice that the measured time scales τ ap and τ c obtained from MC simulations show an increase with y which approximately obeys the power law y 2b with b = 1/2 (see Fig.4). At lower densities, the simulated data agrees better with the y 2b power law. In the analytic calculation above, we assumed single occupancy of the boundary sites of a fascicle basin (removal of a boundary axon was assumed to reduce the basin size). At higher densities this assumption dose not hold, the boundary of a basin does get multiply occupied by axons and thus we see a departure from the y 2b power law.

VI. SINGLE TYPE OF AXONS, WITH DETACHMENT
In the previous section, we concentrated on the "zerotemperature", energy-minimizing dynamics, in which axons cannot detach from a fascicle once they become part of it. Now we extend our analysis to "finitetemperature" Monte-Carlo dynamics, in which the detachment of growth cones from fascicles become possible. It is important to analyze this general case as in the experimental studies of fasciculation dynamics [21,29,30], defasciculation events are clearly observed. An additional reason is that the detachment from fascicles is crucial for the formation of pure fascicles in a system containing multiple axon types -see Section VIII.
As described in Sec. III, in our Monte Carlo simulations (with effective temperature set to unity), a randomly attempted move to left (right) is accepted with probability p L = min[1, exp(−δE l )] (p R = min[1, exp(−δE r )]) where δE l and δE r are evaluated based on the axon occupancy numbers and on the additive interaction energy per axon E h (< 0). The detachment rate of a growth-cone following a fascicle of size n is exp(nE h ). Note that for weaker E h , the detachment rate is larger.

A. Impact of detachment on n∞ and τap
The approach-to-steady-state time scale τ ap decays with increasing detachment rates (decreasing inter-axon attraction |E h |). We perform the approach to steady state data analysis in the same manner as we did for the purely energy-minimizing dynamics discussed in the previous section. This analysis gives an estimate of the time-asymptotic mean fascicle size n ∞ as well as the approach-to-steady-state time scale τ ap as a function of y. We perform this analysis for a system of N 0 = 100 and L = 200 at various strengths of inter-axon interaction E h .  Strictly speaking, the presence of detachment invalidates our earlier mean field argument for the growth of fascicle size with y. While the concept of the basin of a fascicle is still meaningful, the basin can "leak", i.e., newly growing axons can escape from it through the detachment process. Despite this, for most interaction strengths our numerical results obey n ∞ = c + 2ρ eff y 1/2 exp(−βy) with the fitting parameter ρ eff taking the place of ρ [ Fig. 9(a)]. The value of ρ eff is the smallest for the weakest attractive interaction plotted in Fig. 9(a). For extremely weak attraction (e.g., E h = −0.1) almost no fasciculation can occur (ρ eff ≃ 0) and therefore n ∞ remains independent of y (data not shown). As shown in Fig. 9(b), the time scale τ ap for the ap-proach to steady state decays with reduced inter-axon attraction (i.e., with increased detachment rate). When E h ≤ −2, the growth of τ ap with y is similar to the one we observed for the "zero-temperature" dynamics. At E h = −1, however, we find that τ ap becomes comparable to the mean axonal lifetime T , even though at this value of E h the power-law growth of the mean fascicle size with y is still maintained. The reduction of τ ap with increased detachment rate may be understood as a result of enhanced interaction between neighboring fascicles. As we noted before, the slowest time scale for the case of purely energyminimizing dynamics is due to the very slow process of exchange of basin size between neighboring fascicles. In presence of finite detachment rates, axons from one fascicle can detach and connect to a neighboring fascicle, thereby opening up a new and faster mode of interaction between neighboring fascicles.

B. Impact of detachment on the fascicle size distribution
In this subsection we show that even in presence of an appreciable detachment rate (E h ≤ −1), the scaling of fascicle size distribution persists and the scaling function retains its overall functional form. However, with increasing detachment rates (decreasing |E h |), the small n portion of the fascicle size distribution gains at the cost of bigger fascicles and deviates from the scaling form. Fig. 10 shows the steady-state fascicle size distributions for a system of L = 800 and N 0 = 100. As the attraction is decreased from E h = −4 to E h = −1, we observe an overall increase in the small n portion of the distribution. Except for the lowest values of n, we obtain data collapse implying the scaling form P s (n, y) = n −r φ(n/ n ) where r = 2.1 for E h = −4, −2, and r = 2 for E h = −1.
Similarly to the case of strictly energy-minimizing dynamics, this scaling is observed only in the intermediate range of y values, 10 3 < y < 10 4 , and the scaling function is of the form φ(u) = N u exp(−νu−λu 2 ) with u = n/ n . The small-n part of the distribution does not scale. P s (n, y) is large at n = 1, and drops to lower values with increasing n, before it increases again to follow the scaling function. The functional form of the initial decay of P s (n, y) with n depends on y and also E h , as we show in detail in Fig. 10(d).
The discussion in this section shows that there are parameter regimes, e.g., the E h = −1 case discussed above, where the emergent time scales are comparable to the mean lifetime T of single axons, and at the same time the steady state statistics (the mean fascicle size and the fascicle size distribution) obey the overall features demonstrated by the energy-minimizing dynamics. This parameter regime might be utilized to attain a fasciculation pattern of this type in a relatively short time. . Note that the curve obtained at y = 10 4 is the closest to the n −5/2 line at small n. For y > 500, the initial decrease of Ps(n, y) with n is followed by a subsequent increase which merges with the scaling function. At y ≤ 20, the tail of the distribution is exponential (non-interacting limit). Inset of (d): The solid (red) lines are the fascicle size distributions for E h = −0.1 at all y-levels (10 < y ≤ 10 4 ). It shows a clear single exponential decay characterizing the (almost) non-interacting axons. This should be compared with the dotted (blue) lines that show approximate exponential decay of the distributions in the range y ≤ 20 for E h = −1.

VII. RELATION TO PARTICLE AGGREGATION AND COALESCENCE IN ONE DIMENSION
As we pointed out in Sec. IV, in a steady-state configuration at fixed time t, the axon fasciculation with increasing y may be formally viewed as the evolution of a onedimensional (1d) reaction-diffusion process, where the y coordinate takes the meaning of time. As we show in this section, this limited analogy can be used to approximately understand some steady state properties, e.g., the distributions of fascicle sizes and of inter-fascicle spatial separations. We stress, however, that the full dynamics of our system can not be mapped on to a 1d reactiondiffusion system. The process of axon turnover, which is crucial for the dynamical properties of our system, does not have any analog in the 1d models we discuss in this section.

A. In absence of detachment
In this subsection we discuss the relation of our basic model, in which axons cannot detach from fascicles, to irreversible aggregation and coalescence processes.
Interpreting the fasciculation of axons with increasing y (in a steady state configuration) as 1d irreversible aggregation of particles (mA + nA → (m + n)A) [17,46], we find a prediction u exp(−λu 2 ) (see equation 8.4.24 in [17]) for the fascicle-size distribution, which is similar to the true distribution u exp(−νu − λu 2 ) (Sec. V B) [11]. Note that the distribution obtained from the mapping to irreversible aggregation lacks the exponential part exp(−νu). As we explained in Ref. [11], having ν = 0 in our model would require τ ap = ∞. The absence of the exponential part is therefore consistent with the absence of turnover-based dynamics in the 1d analogy.
The merging of fascicles with growing y may alternatively be interpreted as an irreversible coalescence process A + A → A [16], viewing each fascicle as a particle A. The pattern formation in 1d irreversible coalescence had been quantified by the inter-particle distribution function (IPDF) [16], the distribution of distance between neighboring particles. The steady state of this process is trivial, a completely empty space. The IPDF is obtained at finite time t before this steady state arrives, and has the form (x/4Dt) exp(−x 2 /8Dt) where D denotes the particle diffusion constant [16]. The change in the inter-fascicle separation distribution with increasing y in our model may be viewed as equivalent to the time evolution of IPDF in irreversible coalescence. This leads us to a prediction of the distribution of inter-fascicle separation H(∆x, y) ≈ (∆x/y) exp(−µ ∆x 2 /y). Using the same stochastic simulation that we used to obtain Fig. 6, we calculated the histogram of spatial separations ∆x between fascicles identified at various y-levels. Note that the separation between two neighboring fascicles ∆x is measured at the same y-level at which the fascicles are identified, and is different from the gap between fascicle basins E (shown in Fig. 2). A rescaling of H(∆x) by a factor A(y) and ∆x by B(y) leads to a data collapse ( Fig. 11(a)), and approximate power law dependences A ∼ y p and B ∼ y −q with p = 0.96 and q = 0.47 ( Fig. 11(b)). This result is in reasonable agreement with the above-mentioned form of H(∆x, y) that predicts a scaling function P∆x exp(−Q∆x 2 ) (Fig. 11(a)), and scaling exponents p = 1 and q = 1/2.

B. In presence of detachment
In this subsection we discuss analogies of our model of axon fasciculation in presence of detachment to reversible aggregation and coalescence processes in 1d. Due to specific features of our model, only qualitative analogies to models from the 1d literature can be made.
Detachment events are partially captured when the axon fasciculation with increasing y in a fixed-time configuration is formally viewed as a 1d diffusion with reversible aggregation, the chipping model: mA + nA → (m + n)A and mA → (m − 1)A + A [18], denoting each axon by a particle A. (Note, however, that steady state configurations of our model show splitting of fascicles with increasing y (e.g., see Fig. 2(b)) into two fascicles containing multiple axons. The chipping model does not include the analog of such a process). The chipping model posseses a non-trivial steady state (t → ∞, corresponding to y → ∞ within our model). It shows a dynamic phase transition associated with particle density ρ [18,19]. The steady state distribution of clusters of size n is predicted to be P s (n) ∼ exp(−n/n * ) at ρ < ρ c = √ 1 + w − 1, where w denotes a constant single particle chipping rate. At the critical density ρ = ρ c , the distribution changes its shape to P s (n) ∼ n −5/2 . At density above ρ c this power-law distribution remains unaltered, and in addition to the power law distributed clusters one gets a single cluster of diverging size [18]. Fig. 10 shows fascicle size distributions obtained from our model at various interaction strengths E h . Note that, in contrast to the chipping model, in our model the rate w = 1/[1 + exp(−nE h )], with which a single axon detaches from a fascicle, depends on the fascicle size n. Only for the weakest interaction E h = −0.1 (inset of Fig. 10(d)), the detachment rate from a two-axon fascicle w = 0.45 corresponds to a critical density ρ c = 0.2 which is greater than the axon density ρ = 1/8. The corresponding fascicle size distribution shows a form con-sistent with ∼ exp(−n/n * ) (inset of Fig. 10(d)). For E h = −1 the detachment rate from a two-axon fascicle w = 0.12 corresponds to a critical density ρ c = 0.06 < ρ (= 1/8). The fascicle size distribution in the region of small n, for E h = −1, shows rough agreement with the power law ∼ n −5/2 (Fig. 10(d)). This change in shape of the fascicle size distribution from an exponential decay to a power-law decay at small n, thus, is consistent with the dynamical phase transition predicted by the chipping model. Since in our model the detachment rate of axons exponentially decays with fascicle size n, for larger fascicles the detachment rate w gets so small that the analogy with the reversible chipping model breaks down, and the behavior of the system becomes analogous to irreversible aggregation. The fascicle size distribution in the region of larger n (Fig. 10) becomes indistinguishable from axon fasciculation in absence of detachment.
We note that, for our model in presence of detachment, the change in the distribution of inter-fascicle spatial separation with increasing y can not be easily understood in terms of reversible coalescence A + A ⇋ A [16,48] (denoting each fascicle as a particle A). The main reasons are: (i) the reverse reaction A → A + A allows for splitting of a single axon into two, a mechanism not allowed in our model; (ii) in contrast to reversible coalescence, the probability of splitting of a fascicle in our model decays rapidly with increasing fascicle size. Thus the y-independent distribution of inter-fascicle separation c s exp(−c s ∆x), expected from IPDF of reversible coalescence [16], is never reached. At large y, we find a distribution of inter-fascicle separation that conforms more to P∆x exp(−Q ∆x 2 ) (data not shown), consistent with irreversible coalescence (see previous subsection).
We note again that time-dependent quantities in our model have no analog in the mapping to the 1d models we discussed above. The emergence of density-dependent long time scales in reversible coalescence [16,48,49] therefore has no relation to the long time scales in our model, which are due to a very slow reorganization of fascicle basins (Sec. V D) -a consequence of axon turnover.

VIII. MIXED POPULATION OF MULTIPLE AXON TYPES
The axons of olfactory sensory neurons expressing distinct odorant receptors [36] are believed to have shortrange interactions with interaction strengths that are correlated with the type of receptors the neurons express [41,42,50,51]. In the framework of our model, this corresponds to the introduction of multiple types of random walkers, with type-dependent probabilities for attachment to / detachment from fascicles. Neuronal systems containing multiple types of axons are known to achieve pure and stable connections. In the olfactory system, it may be expected [54] that when the growth cone is located in a relatively pure environment (i.e., when it is in contact with axons mainly of its own type), it leads to reduced axonal turnover. While we do not include such effects in our model, in this section we evaluate the mean purity of axon environment S to characterize the sorting dynamics.

A. Mean purity
Let us consider a system containing two types of axons named r and b. Assume that a fascicle contains n r of raxons and n b of b-axons. If the i-th axon in the fascicle is of type r, the purity of environment that this axon encounters within the fascicle is s i = (n r − n b )/(n r + n b ), while if the axon is of type b, the purity of environment is s i = (n b − n r )/(n r + n b ). Then the mean purity of environment obtained by averaging over all N (y) axons is S = ( N (y) i=1 s i )/N (y). The partial sum within a fascicle gives n r (n r − n b )/(n r + n b ) + n b (n b − n r )/(n r + n b ) = (n r − n b ) 2 /(n r + n b ). Thus, the degree of sorting can be quantified as the mean purity of axon environment Notice that 0 ≤ S ≤ 1; S = 1 corresponds to completely pure fascicles containing only one type of axons, whereas S = 0 describes fascicles containing an equal mixture of the two axon types.

B. Approach to steady state
We consider a system of size L = 800 having N r 0 = 50 axons of type r and N b 0 = 50 axons of type b at the y = 0 level (ρ = (N r 0 + N b 0 )/L = 1/8). We use the homotypic interaction energy (interaction between r-r or b-b) E h = −4 and the heterotypic interaction energy (between the two different types r-b) E o = −0.1. We monitor the time evolution of the mean purity of environment S and the mean fascicle sizen. A typical configuration is shown in Fig. 2(b). The mean fascicle size reaches the steady state value at each y-level with two characteristic time scales, similarly to the case of a system containing only one type of axons (Fig. 3). The shorter time scale is intrinsic, the mean lifetime of a single axon T , and the other one is the emergent approach-to-steady-state time scale τ n (equivalent to τ ap defined in the caption of Fig. 3). As in the case of a system containing only one type of axons, this time scale τ n can be orders of magnitude larger than T (Fig. 12(d)).
The intrinsic time scale T does not appear in the dynamics of S, however. S approaches its steady state value S ∞ with a single time scale τ s . The measured mean purity S fits to the form S = S ∞ + r exp(−t/τ s ) ( Fig. 12(a)). The steady-state mean purity S ∞ initially grows with y, reaching a maximum beyond which S ∞ decreases ( Fig. 12(b)). This non-monotonic behavior is seen both for low (ρ = 1/8) and high (ρ = 1/2) density systems ( Fig. 12(b)). Thus for the combination of strong homotypic and weak heterotypic attractive interaction, the fascicles achieve the highest purity at a particular distance from their starting point, beyond which the typical fascicle keeps on losing purity. From energetic considerations, one may understand this behavior as follows. The high strength of the homotypic attraction compared to the heterotypic one leads to sorting and thus the initial growth in purity at lower y-levels. Once one obtains highly pure and sufficiently large fascicles, however, the heterotypic interaction will merge the r-dominated and b-dominated fascicles to form larger mixed fascicles, and thereby lead to a loss in mean purity. The effect of the weak heterotypic interaction becomes significant only at higher y-levels, where large fascicles are formed (recall that in our model, the interactions are additive).
Similar to the case of a system containing only a single axon type in presence of detachment, we find that the steady state mean fascicle size grows as n ∞ = c + 2ρ eff y 1/2 exp(−βy). Fig. 12(c) shows that the subtracted mean fascicle size n s = (n ∞ −c) exp(βy)/2ρ eff follows the power law y b with b = 1/2. The effective density ρ eff and offset fascicle size c are treated as fitting parameters.
The time scales for approach to steady state τ n and τ s grow with y following an approximate power law y 2b with b = 1/2 [ Fig. 12(d)]. This behavior is seen to be independent of density ρ, in contrast to the single-type case, where we found reliable y 2b growth of the emergent time-scales only at lower densities. Recall that, reduced inter-axon interaction strengths lead to lower effective densities ρ eff ( Fig. 9(a)). In the present case of mixed population of axons of two types, the heterotypic interaction is very weak and may have lead to the effectively low-density (y 2b ) power-law growth. We note that τ n and τ s turn out to be approximately equal to each other at all y-levels.

C. Distribution of fascicle composition
In this section we briefly discuss the distribution of fascicle composition P (n r , n b ), measured as the number of fascicles with n r r-axons and n b b-axons at a specified y-level and time t. In Fig. 13 we show this distribution obtained by collecting data within the time interval 10T ≤ t ≤ 15T at various y-levels for a system of N r 0 = N b 0 = 50 and L = 200 (ρ = 1/2). From the approach-to-steady-state data (Fig. 12) it is clear that at t = 10T , steady state is reached only up to y ≈ 100. Thus the distributions obtained at y = 10 3 , 10 4 in Fig. 13 are far from steady state [55]. All the plots in Fig. 13 show a pronounced maximum for evenly mixed fascicles, meaning that most of the fascicles we obtain are mixed by type. However, a careful look at the plots reveals off-peak features of the distribution with reasonable weight that reflect the presence of fascicles with highly asymmetrical composition (e.g., along the n b = 1 line for the plot at y = 100 in Fig. 13, the maximum of the distribution is at n r = 5).
Finally, we comment on repulsive vs. attractive heterotypic interactions. The attractive heterotypic interaction generates larger fascicles asymptotically, but induces a reduction of purity at large y. A repulsive heterotypic interaction (combined with attractive homotypic interaction) would generate enhanced sorting, however it would be at the cost of a decreased mean fascicle size at all y levels.

IX. SUMMARY
In this paper, we provided a simple model of the dynamics of axon fasciculation and sorting. To allow us to concentrate on the collective effects that arise from axon-axon interactions in a large population of axons, we chose a particularly simple implementation of single axon growth. In our model, each growing axon is represented as a directed random walk (Sec. III A). The common preferred growth direction may arise, e.g., from the influence of a spatially distributed guidance cue emitted by a distant target. Other than this common influence, we do not include in our model the guidance of axon growth cones by graded guidance cues, and restrict our attention to axon-axon interactions. The interaction of a growth cone with other axons is modeled as a short-range attractive interaction between a random walker and the trails of other random walkers (Sec. III A). In addition, we incorporated neuronal turnover (characteristic of, e.g., the mammalian olfactory system) by assigning a finite lifetime to each growing axon (Sec. III B).
The strength of the axon-axon interaction was parametrized as an effective energy in a Monte Carlo update. In the energy-minimizing dynamics (corresponding to zero effective temperature in the Monte Carlo scheme), once a growth cone attaches to an axon fascicle, it will never detach (become a free random walker) again. For such dynamics, we extended our previous numerical and analytical results of Ref. [11] (Sec. V). In the general dynamics (corresponding to unit effective temperature), axons may detach from fascicles, with a rate that increases with decreasing axon-axon interaction strength. We systematically studied how such detachment events modify the basic dynamics with no detachments (Sec. VI). Finally, we investigated a system with two types of axons, and analyzed the sorting dynamics arising from a strong homotypic interaction (between axons of same type) combined with a weak heterotypic interaction (between axons of different types) (Sec. VIII). Our main findings were as follows.
The tendency to fasciculate is reflected in the growth of the mean fascicle sizen with the distance y in the preferred growth direction of the axons. Using a mean field argument we showed that for the energy-minimizing dynamics the mean fascicle sizen should grow asn ≃ 2y 1/2 ρ exp(−βy). This agrees with the numerical results of Sec. V. In Sec. VI and VIII, we showed that this growth law persists even in presence of detachment, with the average axon density ρ replaced by a fitting parameter ρ eff .
A more detailed characterization of the steady state is a position-dependent fascicle size distribution. Within the scaling regime L ≪ y ≪ (L/2) 2 , this distribution obeys a scaling law P s (n, y) = n(y) −r φ(n/ n(y) ) with r = 2.1 and the scaling function φ(u) = N u exp(−νu − λu 2 ) (Sec.V). At higher y-levels the distribution becomes bimodal, a new maximum arises which is characteristic of the complete fasciculation. Even in the presence of detachment, the scaling behavior of fascicle size distribution remains valid over a wide range of interaction strengths (Sec. VI and VIII).
The dynamics of reorganization of fascicles at high ylevels was found to be extremely slow. The emergent time scales, e.g., the approach-to-steady-state time τ ap or the auto-correlation time at steady state τ c can be orders of magnitude larger than the mean lifetime of an axon T . In Sec. V D, using an analytical model of effective dynamics involving two neighboring fascicles, we showed that the slowest mode of this dynamics corresponds to the exchange of basin space between the two fascicles, and grows with distance as τ ∼ y. This behavior of time scales survives even in the presence of detachment (shown in Sec. VI A).
While our model is two-dimensional, some limited analogies can be made to one dimensional (1d) models of aggregation, coalescence, and chipping (Sec. VII). We introduced a mapping in which the progressive fasciculation with increasing y (at a fixed time) in our model is mapped on to the time evolution within a 1d system of interacting particles. Using results from the literature on 1d models, we then obtained predictions for stationary quantities in our model. Thus interpreting each axon as a particle A in the irreversible aggregation model mA + nA → (m + n)A [17,46], we obtained the prediction u exp(−λu 2 ) for the distribution of fascicle sizes, which is similar to the true steady-state distribution u exp(−νu − λu 2 ) in our model with energyminimizing dynamics. Likewise, interpreting each fascicle of axons as a particle A in the irreversible coalescence model A + A → A [16], we obtained the prediction (∆x/y) exp(−µ∆x 2 /y) for the distribution of separations between fascicles, which agrees approximately with numerical results from our model (Sec. VII). A limited analogy can also be made between the 1d chipping (Ref. [18,19]) or reversible coalescence (Ref. [16,48]) models and our model in the presence of detachment. Since in our model, the rate of detachment decreases with the fascicle size, the reversible interaction models are relevant only at low axon-axon interaction strengths. In this range of parameters, we were able to relate the observed changes of distribution of fascicle sizes in our model to the phase transition that occurs in the chipping model of Ref. [18] (see Sec. VII). We stress again, however, that the mapping to these 1d models can say nothing about the time-dependent quantities in our model (time in our model has no analog in the 1d models). Therefore, for example, the slow time scales discussed in the reversible coalescence model of Refs. [16,48,49] are unrelated to the slow time scales present in our model.
In Sec. VIII we analyzed a system with two types of axons and type-specific interactions. In this system, axons sort into fascicles according to axon type. We quantified the degree of sorting by introducing the mean purity S of axon environment within fascicles. For the case of strongly attractive homotypic interaction and weakly attractive heterotypic interaction, we showed that the degree of sorting S varies with distance y in a nonmonotonic manner and has a single maximum.

X. OUTLOOK
In this paper we analyzed a model aiming to describe the formation of axon fascicles and the sorting of fascicles by neuronal types in the mammalian olfactory sys-tem. Our goal was to systematically investigate the general non-equilibrium statistical mechanics aspects of the model, leaving the task of building quantitative connections with physiology for the future.
To conclude, we discuss possible generalizations of the basic model defined in this article, and the applicability to biological data on axon fasciculation. First, we note that in our discussion of the dynamical properties of the system, it was essential that the random walkers moved in two (rather than three) spatial dimensions, and therefore cannot cross each other without interacting. In contrast, in a three dimensional system, the concept of fascicle basins would lose its validity. The resulting fascicle dynamics in three dimensions is expected to be significantly different from the two-dimensional dynamics, which we showed to be governed by the competition of fascicles for basin space. It will be necessary to examine to what extent the assumption of two dimensionality is satisfied in the olfactory system. However, this assumption is effectively satisfied in studies of growth in neuronal cell culture [21,43,52], in which the axons move on a plane surface, and interact when crossing each other. In Ref. [43], a fluorescence-based method is proposed for extracting the distribution of fascicle sizes; such experiments would permit a direct test of our model. Note that in some of these studies, our assumption of very strong adhesion of axons to the substrate is not satisfied. Events not included in our model, such as the gradual straightening of axon shafts or the local zippering / unzippering of fasciculated axons (as observed in Ref. [52]) may therefore occur. Note also that natural boundary conditions in such cell cultures are either free boundary [43] or confining channels [21,52], in contrast to the periodic boundary condition used in our simulations.
Recall that the mouse olfactory system contains about 1000 axon types. In this paper, we considered only up to 2 types of axons. It would be interesting to examine if qualitatively new features emerge in systems with many axon types and a range of heterotypic interactions. This would require, however, significantly longer simulations.
Recently, Ref. [13] examined the role of axon-axon interactions in the fasciculation and sorting of axons belonging to mouse olfactory sensory neurons. This showed that inter-axon repulsive interactions arising from Neuropilin-Semaphorin signaling play an important role in the axon sorting. In the pre-target region (before reaching the olfactory bulb), the amount of sorting grows with distance from olfactory epithelium [13]. The pretarget axon sorting is shown to affect the topographic map formation by the neurons in the olfactory bulb [13][14][15]. Using a mutant mouse, Ref. [13] further showed that heterotypic axons sort even in absence of the olfactory bulb, i.e., in complete absence of axon-target interactions. This kind of experiments forms a suitable ground for the application of our model. Notice that in our model, we have shown that a weaker heterotypic attraction combined with a stronger homotypic attraction already leads to sorting. An effective repulsion between the two types would enhance the amount of sorting, however, it would be at the expense of the size of the fascicles formed.
Our immediate future goal is to extract the parameter values of our model from controlled in vitro experiments and then use the model (and its possible extensions) to analyze in vivo data on olfactory pattern formation in mice.