Studying the role of axon fasciculation during development in a computational model of the Xenopus tadpole spinal cord

During nervous system development growing axons can interact with each other, for example by adhering together in order to produce bundles (fasciculation). How does such axon-axon interaction affect the resulting axonal trajectories, and what are the possible benefits of this process in terms of network function? In this paper we study these questions by adapting an existing computational model of the development of neurons in the Xenopus tadpole spinal cord to include interactions between axons. We demonstrate that even relatively weak attraction causes bundles to appear, while if axons weakly repulse each other their trajectories diverge such that they fill the available space. We show how fasciculation can help to ensure axons grow in the correct location for proper network formation when normal growth barriers contain gaps, and use a functional spiking model to show that fasciculation allows the network to generate reliable swimming behaviour even when overall synapse counts are artificially lowered. Although we study fasciculation in one particular organism, our approach to modelling axon growth is general and can be widely applied to study other nervous systems.

Fasciculation. Axonal fasciculation is the process of a growing axon adhering to another, potentially forming groups of axons known as bundles ("fascicles"), which follow similar growth trajectories. Gradient guidance is a common neurodevelopmental mechanism used to navigate growth cones during axonogenesis and this mechanism can include fasciculation. The fasciculation process depends on molecular interactions between proteins in axonal membranes (axolemmas), which can promote fasciculation, defasciculation (i.e. an axon leaving a bundle of fasciculated axons) or growth cone repulsion.
Fasciculation appears to be a common mechanism that has a wide range of different roles. For example, fasciculation helps glomerular formation by sensory neurons in the Drosophila melanogaster olfactory lobe, and guides retinotectal axon growth in vertebrates 5,6 . The importance of fasciculation is further illustrated by experiments in which the proteins that mediate it are mutated. For example, in Xenopus laevis, loss of the Neural Cell Adhesion Molecule (NCAM) that is expressed by sensory Rohon-Beard (RB) neurons causes a reduction in the RB population and its absence from the dorsolateral tract 7 . In Drosophila melanogaster, loss of fasciclin II, a vertebrate NCAM homologue, leads to loss of the vMP2 and MP1 fascicles, which are the first longitudinal tracts to develop in the Drosophila embryo 8 . Evidence from both C. elegans and human connectomes suggests that the formation of fascicles also helps to limit dispersion of axonal trajectories, reducing their average path length and producing more targeted connections between brain regions 9 .
Does loss of fasciculation have functional implications? The answer is unclear. The loss of fasciculation will most likely be a result of mutation of adhesive proteins, which likely have other biological roles; fasciclin II, for example, acts as a chemical signal 10 . This creates a difficulty in analysing the mechanistic contribution non-fasciculation makes towards a disease, because the disease may result from perturbation of the other functions of the mutated adhesive proteins. Nevertheless, some authors predict disease results from non-fasciculation. For example, the loss of vMP2 and MP1 fascicles in Drosophila melanogaster (due to mutated fasciclin II) predicted deleterious effects on correct synaptogenesis in the affected neuronal populations 8 . Furthermore, NCAM deficient mice exhibit mossy fibre fasciculation defects in the hippocampus, which is associated with reduced spatial learning and exploratory behaviour 11 . Aberrant synaptogenesis therefore represents a plausible mechanistic explanation for how loss of fasciculation could contribute to neurological disease. This is relevant to human diseases such as autism, which are associated with aberrant synaptogenesis 12 , and may help to explain why neurodevelopmental disorders arise in patients with NCAM mutations 13 . Given the potential importance of fasciculation to neurodevelopment, we attempt to computationally investigate the importance of its role in maintaining correct axon guidance.
Several computational models of axon growth attempt to represent the dynamic nature of the growth cone and its fasciculation behaviour. In 14 fasciculation is modelled as a culmination of short and long range attractive cues, while in 15 attractive short-range interactions between axons of different types (referring to axons which express different levels of adhesion proteins in their membranes) are considered. These models were used to investigate the likely mechanisms that influence fasciculation and the physical dynamics of the process. These models successfully reproduce the process of fasciculation and incorporate important biological features, such as different attraction strengths between different axon types or axon turnover, which matches a more realistic mechanism in development. However, neither studies the effects of fasciculation within a complex chemical growth environment involving long-range growth cues that can attract/repulse axons in multiple directions in order to generate complex axonal trajectories, nor are they based on and measured against biological data from real animals. A different approach was taken in a recent paper 16 , where the mechanism by which axons can fasciculate and defasciculate by a process termed "zippering" was studied; a biophysical model was used to determine certain properties of this process, such as the force of adhesion between axons.
When fasciculation is included in our growth model, the resulting axons have a structure that consists of multiple bundles. Varying the model parameters changes the number of bundles and their average size. This corresponds with experimental results from real animals, which show that commissural axons at earlier stages of spinal cord development (stage 28/29) are arranged into multiple bundles, which are located at certain positions along the spinal cord 17 . Our functional simulations show that when the number of synapses is reduced, networks with fasciculation are able to generate more reliable swimming activity than the networks without fasciculation.
The paper is organized in the following way: Section 2 describes the models, including our new parallel axon growth model and the spiking physiological model. Section 3 shows the results of our experiments with the models, including patterns of axonal projection and spiking activity. Finally, section 4 gives the conclusions of our modelling and compares our findings with the results of other models and available biological data.

Model description
Our model of nervous system development is an improved version of the model that we have previously described 3,4 . We will briefly summarize the original model here, but for full details see the earlier papers.
We consider the spinal cord as a 2D rectangular area upon which somata, dendrites and axons are allocated (Fig. 1). This rectangle can be imagined as the result of opening the spinal cord like a book along the dorsal midline (dashed line in Fig. 1A). We disregard the third dimension (the thickness of the spinal cord "tube") as this is very thin, with all somata, dendrites and axons located in a 10 µm thick region around the sides of the tube. The horizontal co-ordinate (x) corresponds to the distance from the tadpole's midbrain in the rostro-caudal (RC) direction, and we only consider the region μ < × < μ 500 m 2000 m. The vertical co-ordinate (y) corresponds to the distance from the ventral mid-line in the dorso-ventral (DV) axis, with positive values representing positions on the left side of the body and negative values the right. In this direction we only consider the region − μ < < μ y 145 m 145 m. Although in this paper we only consider this short 1500 µm section of the spinal cord, when a longer section is considered the behavior is similar 4 .
There are seven different types of neuron included in the model: sensory Rohon-Beard (RB) neurons that respond directly to skin stimulation; dorso-lateral ascending and commissural neurons (dlas and dlcs) that make up the touch sensory pathway; ascending, descending and commissural interneurons (aINs, dINs, cINs) that make up the locomotor central pattern generator, and motoneurons (mns). In each generated connectome, the number of neurons of each type is fixed, with the same number of neurons of each type on both sides of the body. See Table 1 for neuron counts and the colour coding used throughout this paper.
Basic axonal growth procedure and model optimization. Axon growth in the model is based on two chemical gradients: one in the rostro-caudal direction and another in the dorso-ventral direction. The growth cone is characterized by different sensitivities to these gradients. Axons grow in discrete time with elongation m 1 Δ = μ per time step. To describe the dynamics of the growing axon we use the growth angle θ A which is characterised by "stiffness", the tendency of the growing axon to grow straight, keeping the same growth angle which was used on the previous growth step. The influence of environmental cues (according to their gradients) deviates the growth cone from a straight path. The addition of a random variable at each step of growth provides an additional degree of freedom. Axon growth is guided by two gradients, resulting in a change of the growth angle value: One gradient, G RC , influences the axon growth in the rostro-caudal direction, and the other, G DV , influences axon growth in the dorso-ventral direction. Each gradient is projected to the direction perpendicular to the current growth direction to describe the change of the growth angle. The model is described by the following difference equations: x cos Here x y ( , , ) n n n θ are the current coordinates and growth angle of the axon tip at time step n and Δ is the axon elongation at each step (usually 1 µm). The term θ n A describes how the angle of growth changes according to rostro-caudal and dorso-ventral growth cues, as specified by the functions G x y ( , ) RC and G x y ( , ) DV respectively, as well as a random variable ξ n that is sampled at each step from a uniform distribution in the range θ + are used to model interactions between growing axons, and were not present in our previous model of axon growth. Parameter s, where − < < s 1 1 , represents the sensitivity of axons to other nearby axons, while the term θ p in equation 1d is an angle calculated from the growth angle of the closest nearby axon, as described below. Note that if = s 0 then θ + n B 1 has no effect on the angle of growth and the model behaves identically to the previously published version.
To start axon growth, we assume that the initial axon position coincides with the soma position. The initial value of the growth angle is randomly selected using the 2D generalization procedure described in 3 for experimental measurements of pairs θ y ( , ), where y is the DV coordinate of the soma and θ is the initial angle of axon growth near the soma. The axon length is selected using the same procedure, using experimental measurements of total axon length.
The functions that control axon growth in response to gradients, G x y ( , ) RC and G x y ( , )

DV
, contain another set of parameters that determine the chemical gradient environment, and how axons grow in response to this. The values of the parameters which describe the gradient environment are independent of cell type and were selected according to a general biological knowledge on the distribution and properties of chemical gradients. For each cell type four parameters controlling growth are required: three giving sensitivities to gradients and one giving the range of the random deviation variable (α). To find values for these parameters we use a pattern-search optimization procedure to minimize a cost function, which takes into account how much the generated axons differ from experimentally measured axons in terms of their dorso-ventral distribution and tortuosity ("wiggliness"). For details of the cost function and optimization procedure see 3 . The optimization procedure was repeated in order to find optimal parameter values for each neuron type.
It is important to note that the experimental measurements of individual axons that were used for the optimization procedure were collected from different animals and therefore represent trajectories of single axons growing independently from each other. Consequently, these data cannot be expected to reflect any possible effects of fasciculation within individual animals. In the case of the axon growth model with fasciculation (or repulsion), multiple axons from a single animal should be used for fitting, however such experimental data from tadpole spinal cord axons is not currently available. Therefore, to generate axons in the case of fasciculation/repulsion we used the same parameter values that were found by the optimization procedure in the model without fasciculation. Re-running the fitting procedure to find the optimal chemical gradient parameters with fasciculation produced worse fits to the experimental data than without axonal interactions (see section "Model optimization with fasciculation"), but this is to be expected given the limitations of the biological data. However, multiple simulations (N = 12) confirmed that using parameter values optimized for fasciculation produced axon growth that Branching points and commissural neuron growth. The primary axon of most neuron types grows from the soma in an ascending direction (from tail to head); the exceptions being the dINs and mns which have descending (head to tail) primary axons. Neurons of most types also have secondary axons, which grow from a particular branching point in the opposite direction to the primary axon. The neuron types that do not have secondary axons are the dlas and mns, as well as some dINs. The coordinates of branching points are selected in a random way using the generalization procedure described above, using available experimental data on branching point locations for each neuron type.
The axons of commissural neurons (dlcs and cINs) start to grow on one side of the body and rapidly navigate in the ventral direction due to the strong influence of DV gradients. Upon reaching the floor plate they cross the body and change their sensitivity to the DV gradients. After crossing the ventral midline, instead of moving with the gradient they continue to grow on the contralateral side of the body against the DV gradient, with deviation to the ascending direction. Their secondary axons are positioned on the contralateral side and grow in the descending direction 17,18 .
Axon growth is limited by the rectangular modelling area, and also by two longitudinal barriers. The model contains two dorsal barriers: one formed by RB somata located at DV co-ordinate = ± y 137 and spanning most of the rostro-caudal extent > x 500, and another formed by dlc/dla somata located at y 127 = ± and spanning the range > x 700. Together these barriers form a dorsal tract which exclusively contains axons of RB neurons. Ventrally, another barrier that spans the entire length of the environment with DV co-ordinate y 25 = prevents axons (except those of commissural neurons) from entering the floor plate and crossing to the other side of the body.
There is limited biological evidence on the process of branching and secondary axon growth after bifurcation. For example, it was found that in the Xenopus spinal cord commissural neurons axon branching occurs after crossing the floor plate and in response to external cues found there, including Slit-Robo signalling 19,20 . It is not known exactly when after bifurcation secondary commissural axonogenesis begins, although some data suggests it is soon after appearance of the branching point 21 . Mouse spinal cord sensory neurons possess secondary axons that similarly bifurcate in response to Slit-Robo signalling, and which grow soon after bifurcation 20 . In the model all primary axons grow first and then all secondary axons. We chose to do this because of the lack of concrete biological evidence of secondary axon growth timing, and will instead investigate the effects of varied secondary axon growth timings in future versions of the computational model.
Dendrites and synapse formation. Each neuron has a dendritic field represented as a single line oriented in the dorso-ventral direction, centred on the soma position. The dendritic extents (i.e. length of the line) for each neuron are chosen using the same generalization process as axon initial angle and length, based on experimental measurements of dendritic fields of each neuron type. If a growing axon crosses a dendrite then synaptic connection will be generated with some probability. The probabilities of synaptic connections between neurons of different types have been experimentally defined in various pairwise recording experiments (for details see 1 ). The "geography" of neurons has a strong impact on synaptogenesis. RB axons, for example, have y-coordinates in the range y 127 137 < < because they are allocated in the dorsal tract. Therefore, the axons of RB neurons have no chance to cross the motor neuron dendrites, as these are located in very ventral positions.   θ we consider the nearest existing axon (shown in orange) within a certain radius given by a parameter (r). The growth angle at the nearest point on the existing axon is denoted θ p , and the angle perpendicular to this is denoted θ p . In case of fasciculation s ( 0) > the current axon will tend to grow with the same angle as the existing axon, whereas in the case of repulsion it will try to grow perpendicularly away from it. While our previous model grew axons one at a time sequentially, the new version of the model allows axons from all neurons to grow at the same time, since this is required for axons to interact with each other. At each time step the growth cone positions of all axons are updated and each axon is elongated by the space step Δ. The growth process runs until all axons have been generated up to their full length. Different axons can begin growing at different times. For simplicity, we consider two groups of axons: "pioneers", which grow first, and "followers", which start to grow after all pioneers have finished.

Assumptions/restrictions of anatomical model implementation.
This paper presents a model based on anatomical measurements of developing Xenopus spinal cord neurons. We have attempted to incorporate as many biological details as possible into the model. However, there is no available anatomical information concerning fasciculation in animals at this stage of development, and there are still many open questions awaiting answers from experimental data. As a result, we have to formulate some assumptions for the model, some of which are motivated by logical considerations and some for the purposes of model simplification.
Assumption 1. Primary axons first, secondary after. We consider growth of primary axons and secondary axons as two consecutive processes. All primary axons grow first, and after that secondary axon growth starts. Naturally we prescribe a starting time for growth of each secondary axon. A future version of the model may include the possibility of starting secondary axon growth before the primary axon has finished growing. Assumption 2. Fasciculation/synapses of commissural neurons on the opposite side only. We assume that the commissural neurons (dlcs and cINs) can fasciculate/repulse and create synaptic contacts on the contralateral side only. At the beginning, commissural neurons grow in the ventral direction according to the axon growth equations (1) with specially adjusted parameter values. After crossing the boundary of the ventral plate on the opposite side = ± y ( 25) the axon growth is described by the same equations but with another regular set of parameter values 3 . During the initial pre-crossing stage the axons grow without fasciculation/repulsion and do not produce synapses. This assumption allows us to model axon growth on each side of the body independently of each other. Assumption 3. Interaction with axons of the same cell type only. We assume that axons can fasciculate/repulse on to/away from the other axons (either primary or secondary) of the same type only (e.g. axons of cIN neurons can only fasciculate onto axons of other cINs). Although the biological reality is likely to be more complicated than this, we hypothesise that in many cases it is important for axons of a particular type to follow a specific patterns of growth, in order to achieve the macroscopic connectivity between populations required for proper network function. We would expect that if all neuron types were able to fasciculate onto each other, network connectivity would become less type-specific, which we would expect to negatively affect the network's function. This hypothesis should be tested in future work. Assumption 4. Universal values of fasciculation sensitivities. The fasciculation/repulsion sensitivity parameter has the same value for all cell types. However, we allow primary and secondary axons to have different sensitivities, denoted s pr and s se respectively. Assumption 5. Fasciculation to the nearest point of a pioneer axon. A growing axon's growth is only affected by a single point on an existing axon -specifically the point that is closest to its current growth cone position and within distance r. This point could be part of either a primary or secondary axon. At the next step of growth a new closest point is selected independently of the previous step.
Functional model. The connectome produced by the growth model provides detailed information on connectivity in the spinal cord -specifically a list of synaptic connections between neurons. We use this information to build a functional model that simulates the spiking activity of spinal cord neurons. Briefly, this model is of single compartment Hodgkin-Huxley type, with channel properties chosen to match known electrophysiology. The functional model includes axonal and synaptic delays, as well as electrical coupling between dINs. For full details, see 4 . In contrast to the previous paper where custom software was used, in this paper the functional model simulations were performed using NEURON 22 with a step size of 0.01ms. All other details of the model were identical to those previously published.

Results
In this section we present the results of simulations of the growth and functional models. For the growth model we first show typical patterns of axon distribution for different neuron types in the case without fasciculation or repulsion. We then compare these results with the patterns that appear in the cases of fasciculation and repulsion. Similarly, we analyse how these different patterns of axonal projection affect synaptic connectivity. We then use Here, and in the following section, we focus on presenting qualitative results of our simulations. For each set of parameter values used, at least 12 connectomes were generated and the results compared to ensure they were visually similar. We look at more quantitative measures of the effects of fasciculation in the later sections. There are several parameters that control the effects of axonal interactions in the model (see Methods section). The sensitivity parameter s controls the strength of attraction/repulsion between a growing axon and an existing pioneer axon. When = s 0 axons do not interact, whereas if s 0 1 < ≤ axons will fasciculate, and if − ≤ < s 1 0 they will repulse each other. The model includes two sensitivity parameters, s pr and s se , for primary and second axons respectively. The range parameter r determines the maximum distance (in µm) from a growing axon tip over which an existing axon can exert an influence. Finally, for each neuron the time at which both its primary and secondary axons begin growing are specified. These parameters allow us to divide axons into two groups: pioneers that start to grow first, and followers which start to grow after the pioneers. Taken together, these parameters have a large effect on the resulting pattern of axon growth; for example, increasing the sensitivity parameter value results in axon bundles appearing.  (Fig. 3A), = .
s 0 1 pr (Fig. 3B) and s 0 5 pr = . (Fig. 3C); in all cases r 1 m = μ . In these simulations there are four pioneer axons (black lines in Fig. 3B,C) that start to grow simultaneously. All other axons (followers) grow after all the pioneer axons have finished growing. The time at which follower axons begin growing varies rostro-caudally, with the most rostral axon starting first and the most ventral axon last, with 200 time units between the start of each axon's growth. One time unit corresponds to the real time which is needed for an axon to elongate by 1 µm (we assume that all axons grow uniformly with the same speed). Thus, when the length of the first follower axon is 200 µm the second follower starts to grow. Figure 3B clearly demonstrates that fasciculation, even with a small sensitivity value, leads to the axons grouping. The bundles are clearly visible in Fig. 3C where the value of sensitivity is even higher. Now we vary two parameters of the model: sensitivity and range. To illustrate how these parameters influence the pattern of axon growth we consider primary RB axons only. Comparing Fig. 4(B),(D), with different sensitives and the same range, shows again that higher sensitivities give stronger groupings of axons. The same result is clear from comparing Fig. 4(C),(E). Although we expected that increasing the range at which axons can interact would increase the degree of fasciculation, this was not true for our model, as can be seen by comparing Fig. 4(B,D) (r 1 m = μ ) with (C, E) ( = μ r 3 m). The reason for this is clear: in our model, a positive value for the sensitivity means that neurons that are within range of each other will tend to grow in the same direction, not towards each other. This means that for larger values of the range (e.g. 3 μm) axons tend to grow in parallel with gaps between them, rather than bundling closely together. For the rest of this paper we avoid using such large values of the range, and instead fix = μ r 1 m.

Patterns of primary and secondary axons.
In this section we consider patterns of both primary and secondary axons for different values of the sensitivity parameters, with the range parameter fixed as 1 µm. Starting times for primary axons were defined in the same way as described in the previous section. When all primary axons have reached their full length, the secondary axons start to grow. The starting times of growth for secondary axons are chosen in a similar way to the primary axons. There are nine pioneer secondary axons which start to grow simultaneously, and after that all other secondary axons start to grow in sequence from rostral to caudal positions with a 200 time unit step. In the figures in this section, secondary axons are shown in magenta, to distinguish them from primary axons (yellow). Black lines show the pioneer axons.  . It is clear in Fig. 5(A) that without fasciculation axons spread inside the dorsal tract and uniformly fill the space between two barriers: primary axons are located preferentially in the rostral part of the body and secondary ones in the caudal part. In the case of weak fasciculation (B) the pattern of axons is less dense and a tendency for axons to organise into groups is visible. Also, the secondary axons are more concentrated in the upper part of the simulated space near the higher barrier. Increasing the sensitivity value (C) makes this feature even clearer: several bundles of primary axons are visible and the secondary axons are located in the narrow area near the upper barrier.
In Fig. 6 we show all axons for the full connectome with fasciculation sensitivities = = .  Simulations of the anatomical model show that variation of the sensitivity parameter value results in very different patterns of axon growth with fasciculation/repulsion. To illustrate this we present, for two example  pr se (Fig. 7). It is clear from these figures that in case of fasciculation (middle column) the pattern of axons is much more concentrated than without fasciculation (left column), and multiple axon bundles are visible, particularly in the case of cINs. The cIN axons (especially secondary ones) are more concentrated in the ventral region when fasciculation is present, since axons that would otherwise grow to more dorsal positions are influenced to turn to grow longitudinally by those that have already grown, such as pioneers. In the case of repulsion (right column), even a small value of the sensitivity coefficient causes growing axons to grow in extremely tortuous trajectories that occupy all the available space. We would therefore expect repulsion to be associated with a decrease in the specificity of connections between types resulting in a loss of network function.
Fasciculation reduces the number of synapses. Our simulations show that the total number of synapses decreases when fasciculation is present. For example, when fasciculation sensitivity is 0.2 and range is 1 µm, the mean total number of connections is 72, 554 ± 753 (n = 12), compared with ± 81, 877 744 (n = 300) without fasciculation. Despite the overall reduction in synapse counts, the actual impact of fasciculation varies based on the type of the pre-and post-synaptic neuron, as Table 2 shows. Each entry gives the average (rounded) number of synapses between different cell types for connectomes with fasciculation = .
= s r ( 02, 1) and without fasciculation (first and second number in each cell respectively).
We can use Table 2 to study how each part of the network is affected by fasciculation. Synapse counts in the sensory pathway (RB→dli/dlc→CPG) are relatively unchanged, although it is interesting to note that fasciculation causes an increase in the strength of the ipsilateral reflex pathway from dlas directly onto motoneurons. Within the central pattern generator network, the biggest decreases in synapse counts are seen in connections to and from ascending interneurons (aINs), but since these neurons are almost completely silent during swimming this change cannot have an effect on the swimming pattern. Within the core "active" CPG network (cINs and dINs; bold bordered area in Table 2) the reduction in synapse counts is consistent but modest. However, we know that the characteristics of the swimming pattern generated by the core CPG network is very dependent on relative synaptic strengths; in the next section we will investigate the effect of fasciculation on the network's spiking activity. Finally, fasciculation causes a marked increase in the strength of rhythmic excitation to motoneurons, thanks to an increase in dIN→mn excitation and a decrease in cIN→mn inhibition.
Fasciculation compensates for imperfect growth barriers. The axons of sensory RB neurons are confined to a narrow dorsal tract by two longitudinal barriers: one more ventrally at µ = ± y m 127 formed by the somata of dlc and dla neurons, and one more dorsally formed by the somata of RB neurons at y m 137 = ± μ . The majority of dlc and dla dendrites are located in this dorsal tract, so the barriers ensure that RB axons make many synapses onto these dendrites, resulting in a strong sensory pathway. In our model these barriers are implemented as hard barriers that axons can never cross, but in reality gaps in the barriers may be possible, for example due to spaces To answer this question we modified the growth model so that both barriers to RB growth contained 25 µm gaps at 25 µm intervals. We then grew RB axons as normal, and calculated what proportion of axon points lay outside of the dorsal tract. Figure 8A shows the results of this process in one case without fasciculation, demonstrating that RB axons are able to escape the dorsal tract via the gaps in the barrier. Figure 8B shows a simulation with fasciculation enabled s r ( 02, 1) = .
= , showing that in this case axons are generally restricted to the dorsal tract despite the gaps in the barrier. Without fasciculation, 24% of RB axon points (std. 2.6, N = 12 simulations, 63 RBs in each) were outside the dorsal tract, whereas with fasciculation only 15% of points (std. 3.3, N = 12 simulations, 63 RBs in each) were (p < 0.001; Student's unpaired t-test, two-tailed). The reason for this reduction can be seen intuitively, by considering that RB axon growth is primarily longitudinal. Since attraction causes growing axons to turn in the same direction as existing axons, pioneers and other already-grown axons all act as a form of soft barrier that prevents growing axons from growing in dorsal or ventral directions.
In most cases all of the pioneer axons remained in the dorsal tract, since the probability of any given axon escaping is relatively low. This laid down a scaffold that the follower axons followed. Fasciculation can therefore reduce the spread of axonal trajectories and help to restrict axon growth to a particular path or region.
Model optimization with fasciculation. As previously discussed, we performed the optimization process used to determine the axons' sensitivities to the chemical growth environment without the presence of fasciculation. For a range of non-zero sensitivity values (from 0.1 to 0.7, µ r m 1 = ) the resulting pattern of axonal growth did not match the experimental data as well as when fasciculation was not present = s ( 0). Figure 9 shows the effect of fasciculation on the dorso-ventral distribution of axon trajectories, which is the main component of the optimization procedure's cost function. The data shown in this figure is for aINs, but we observed similar effects for other neuron types. From this it is clear that fasciculation has two main effects: a slight ventral shift in the modal axon position and a pronounced sharpening of the peak of the distribution. The ventral shift can be explained by the fact that axons initially grow ventrally before, at some point, turning to grow longitudinally. With fasciculation present, axons that might otherwise turn longitudinally at more dorsal positions are prevented from doing so by the influence of other axons. We explain the narrower peak in the distribution by noting that the axon data used for optimization was based on traced single axons from different individuals, and hypothesising that the wider peak in the anatomical data results from natural variation in the growth field / growth cue sensitivities between animals. If data from multiple axons in the same animal were available, these may show the effects of fasciculation as a sharper peak in the distribution.
We also repeated the optimization process using the same set of fixed sensitivity values in order to obtain new optimal parameters for the chemical gradients. Unsurprisingly this did not produce axons that more closely matched the anatomical data, as there was still a relatively tall and narrow peak in the distribution of dorso-ventral positions. This peak is a consequence of fasciculation causing axons to bundle together, which is not something that the optimization procedure can affect by changing growth cue sensitivities. Since the best fit to the limited anatomical data available is that produced when there was no fasciculation present during optimization, we use those optimized parameters throughout this paper.
Effect of Fasciculation on Swimming Activity. As previously reported 4 , without fasciculation the growth model produces connectomes that very reliably generate swimming activity in response to sensory stimulation. In this section, we use a functional spiking model to investigate how fasciculation affects the activity generated by the network. = . Each data point shows the average of N = 12 connectome simulations.

Effects of Fasciculation.
Introducing fasciculation reduces the number of synapses in a connectome. If fasciculation is too strong, the number of synapses is reduced so much that the generated connectomes are unable to swim. However, connectomes generated with only modest fasciculation (sensitivity 0.2, radius 1) do reliably swim in response to sensory input (6/6 connectomes), with a similar frequency (18.6 ± 0.8 Hz) to unfasciculated connectomes. This result raised the possibility that fasciculation would allow connectomes with fewer synapses than normal to still generate swimming activity. To investigate this we artificially lowered the probability of synapse formation in order to generate connectomes with reduced numbers of synapses. Swimming, defined as a stable pattern of alternating left and right motoneuron spikes at about 18 Hz, could be seen in the majority of connectomes even when the average synapse count was reduced to around 43,000. This threshold for swimming was roughly the same regardless of whether connectomes were generated with fasciculation or not, although fasciculated connectomes were able to swim with slightly fewer synapses (Fig. 10). However, in connectomes without fasciculation but with fewer synapses there was a much greater tendency for dINs to fire mid-cycle spikes, i.e. at roughly the same time as contralateral neurons (Fig. 11A). This happens because dINs are capable of firing NMDA-driven pacemaker spikes after their normal post-inhibitory swimming spikes. Normally these spikes do not occur, since contra-lateral inhibition arrives early enough to stop them, but decreasing excitatory synaptic drive to a dIN decreases the delay to the first NMDA-driven spike. If this delay is short enough (i.e. when the number of synapses have been reduced dramatically) then spikes can happen before commissural inhibition arrives. This aberrant activity is not seen in real swimming, and we therefore consider swimming in which there is a lot of mid-cycle dIN activity to be lower quality. By this measure, fasciculation allows the same quality of swimming to be achieved with fewer synapses (Fig. 11B).
Fasciculation's effect of reducing the number of mid-cycle dINs is somewhat surprising, since fasciculation also reduces the number of synapses. This apparent contradiction is explained by the fact that dINs that do not receive any excitatory input from other dINs rarely fire mid-cycle spikes (Fig. 11C). Fasciculated connectomes tend to have many more dINs that do not receive any synapses from other dINs (non-fasciculated: 6 ± 3 dINs, fasciculated: 31 ± 6 dINs, 24 connectomes, p < 0.001, Welch's t-test). The spatially localised axon bundles that result from fasciculation mean that a dIN tends to either receive a lot of excitatory input or none at all, avoiding the intermediate levels of excitation that lead to mid-cycle spiking.  Table 1. (B) Reducing the number of synapses increases the number of dINs that fire mid-cycle spikes. For a given number of synapses, connectomes generated with fasciculation = .
= s r ( 02, 1) have fewer mid-cycle dINs. (C) A dIN's spiking pattern depends on how many synapses it receives from other dINs. Neurons that receive little excitation (but not none) have a roughly 50% chance of firing mid-cycle spikes. Inactive dINs tend to receive a lot of excitatory input, but are unable to spike due to depolarisation blockade. Data from four connectomes with fasciculation and probability reduced to give approximately 50,000 synapses.

Discussion
In this paper we present a new computational model of neuronal development in the Xenopus laevis spinal cord that includes the effects of axon-axon interactions. Our results suggest that attraction between axons (fasciculation) produces networks that have a simpler structure than networks without fasciculation, and that can generate reliable swimming behaviour with fewer synapses. Fasciculation also helps to produce appropriate patterns of axon growth when other growth cues (here the marginal zone barriers) are damaged. As discussed in section "Model optimization with fasciculation", the model was optimized to match a set of measured axons that were taken from different individuals, which potentially means that axon positions in this data are more widely spread out than they would be within a single animal. When fasciculation is included in the model the spread of axon positions becomes narrower, and we have shown that the more selective connectivity that results increases the reliability of swimming. This highlights a potential important limitation of combining measurements of single axons from different animals together into a single dataset.
We are aware of two other computational models of axon fasciculation in the literature. The model described in 15 simulates axon growth as a random walk on a lattice. This model includes the ability for axon type dependent interactions (e.g. axons may be attracted to the same type but repulsed by a different type), and while this feature was not included in our model for reasons of computational efficiency it is something we plan to investigate in future work. Their paper also includes a process whereby axons can detach from a fascicle (i.e. when they have reached a target region); our model does not explicitly include a mechanism for this but does allow axons to detach from bundles naturally, for example if the concentration of growth cues becomes strong enough to overcome the attractive force. The advantage of the model in 15 is that its reasonably simple mathematical formulation allows theoretical results (e.g. about fascicle size) to be calculated, whereas we have followed a more experimental approach. A relative strength of our model is that it is built upon on an existing model of axon growth in the tadpole spinal cord that already includes other biologically realistic features of the growth environment, such as growth cue gradients and barriers. In 14 a model is presented that is closer to ours in terms of its implementation, as it features axons that grow in continuous space in response to external growth cues (but without barriers). The model simulates the effects of cell adhesion molecules that allow touching axons to join together, and includes chemoattractants that diffuse from axons to attract other nearby axons. In this model they found that the strength of environmental growth cues was not generally strong enough to cause axons to unbundle from a fascicle, so they allowed axons to switch to repulsing each other once the external growth environment reaches a certain trigger threshold; a mechanism that is biologically motivated. We have done some preliminary work in this direction, whereby commissural neurons switch from attraction to repulsion once they have crossed to the other side of the body, but we will present these results more fully in a later paper. An advantage that our model has over both of the previous models described here is that it produces complete tadpole spinal connectomes that can be simulated using a spiking network. The results of these spiking simulations can be compared with known biological results, in order to study the functional effect of fasciculation.
As with the previous models, ours can be used to investigate physical properties of fasciculation, including sorting dynamics and the mechanics of fascicle formation. Despite the simplicity of the attraction/repulsion model, we are still able to visually observe the tendency for axons to grow together, even in the case of relatively weak attraction. With some further refinement it would be possible to use the model to study the effect of varying the growth cone size/axonal interaction distance on the formation of fascicles. This question is relevant to understanding the pathobiology of diseases of abnormal axon pathfinding and subsequent synaptogenesis, such as autism 12 , as protein mutations that alter growth cone size could alter fascicle formation patterns and subsequently impact on correct axon pathfinding. These results would also be relevant to the design of regenerative therapies that treat peripheral nerve injuries, predicting, for example, that growing neurons with large growth cones along pre-laid axon tracts would lead to the formation of fewer, larger fascicles. This would be desirable if the goal were to grow a fascicle of axons towards a single target area. To investigate this using our model, a modification could be made such that nearby axons initially grow towards each other, before growing in along parallel trajectories when very close.
The results presented in this paper show that there are at least two plausible roles for fasciculation: reducing the number of synapses required for reliable swimming and improving RB axon guidance when growth barriers are damaged. This latter result is at least partially in agreement with the finding that the formation of a longitudinal tract by RB axons is disrupted when expression of cell adhesion molecules related to fasciculation is disrupted 7 . Conversely, however, it has been suggested that fasciculation may make aberrant axonal growth worse in some situations, since in Drosophila embryos with disrupted growth barriers follower axons may grow far away from their target areas as a result of fasciculating onto pioneers 23 . It is difficult to assess whether the patterns of axon growth that our new model generates match biological spinal cord axons more closely than when fasciculation is not included. The data available to us about spinal axons for developmental stage modelling (stage 37/38) are generally from studies where only single neurons were stained, which makes it impossible to discern bundles. Evidence from tadpoles at an earlier stage of development does suggest that fasciculation is a feature of axon growth for both RB neurons 18 and commissural neurons 21 ; however it is possible that as the animal grows fasciculated axons are pulled apart, which would again lead to bundles not being present at later stages. Our model does not include any simulation of axon shaft dynamics (i.e. how axons move after their initial growth), meaning it is unable to show the effects of further body growth on axon position, or reproduce the "zippering" phenomenon studied in 16 . Furthermore, there are hypothesised advantages to fasciculation which our model has not explored, such as it leads to more efficient and targeted network formation 9 . Nevertheless, our model suggests that fasciculation could play other important roles in axon guidance in the tadpole spinal cord. It should also be noted that our approach is general and can be applied to simulating the development of other nervous systems including humans. A greater understanding of how nervous systems develop and recover from damage will hopefully one day lead to new treatments for developmental disorders and traumatic injuries. Data Availability. The code used to perform the experiments described in this paper is available from the authors on request.