Models of Purkinje cell dendritic tree selection during early cerebellar development

We investigate the relationship between primary dendrite selection of Purkinje cells and migration of their presynaptic partner granule cells during early cerebellar development. During postnatal development, each Purkinje cell grows more than three dendritic trees, from which a primary tree is selected for development, whereas the others completely retract. Experimental studies suggest that this selection process is coordinated by physical and synaptic interactions with granule cells, which undergo a massive migration at the same time. However, technical limitations hinder continuous experimental observation of multiple cell populations. To explore possible mechanisms underlying this selection process, we constructed a computational model using a new computational framework, NeuroDevSim. The study presents the first computational model that simultaneously simulates Purkinje cell growth and the dynamics of granule cell migrations during the first two postnatal weeks, allowing exploration of the role of physical and synaptic interactions upon dendritic selection. The model suggests that interaction with parallel fibers is important to establish the distinct planar morphology of Purkinje cell dendrites. Specific rules to select which dendritic trees to keep or retract result in larger winner trees with more synaptic contacts than using random selection. A rule based on afferent synaptic activity was less effective than rules based on dendritic size or numbers of synapses.


Introduction
The cerebellum is involved in coordinating motor functions as well as in cognition and emotion [1][2][3][4][5]. The cerebellar volume comprises only 10% of the whole brain but holds more than 70% of all neurons. Most of these neurons are cerebellar granule cells whose population seems too numerous to be counted accurately (estimated densities range from 500,000 [6] to 6,562,500 [7] cells per mm 3 in mice, compared in [8]). The adult mammalian cerebellum consists of a cortex and nuclei, with the cortex containing three layers. A middle layer with cell bodies of its output neuron, the Purkinje cell, an outer molecular layer with Purkinje cell dendrites and granule cell axons and an inner granule cell layer with granule cell somata ( Fig  1A2). But at birth cerebellar cortex is largely undeveloped, with a different layering structure. Granule cells proliferate in an external granular layer (Fig 1A1), most of them postnatally, and then migrate through a rapidly expanding molecular layer to arrive in the granular layer. This migration is guided by radial fibers of a specialized glial cell, the Bergmann glia [9][10][11][12][13]. The granule cell axons, parallel fibers and descending axons, are formed during the migration through the molecular layer and establish synaptic contacts with the rapidly growing Purkinje cell dendrites. Adult Purkinje cells have a uniquely flat dendrite, which is optimal to maximize connectivity with perpendicular parallel fibers [14][15][16][17], with usually a single root segment. But during early postnatal development, each Purkinje cell grows multiple dendritic trees (defined as the arbor connected to a dendritic root) ( Fig 1B1) and then selects a primary tree among them by retracting most of the newly grown dendrites (Fig 1 B2) [18][19][20][21][22]. This is followed by additional growth and development of a more monoplanar tree in the third postnatal week [23].
In isolated in vitro environments, Purkinje cells do not fully retract excessive dendritic trees, resulting in persistent multiple primary dendrites [21,24]. Granule cells and their parallel fibers are strong environmental factors to regulate the dendritic arborizations and retractions in terms of physical and synaptic interactions with Purkinje cells [25][26][27][28]. However, the mechanisms by which granule cells select the final primary dendritic tree remains unclear.
Because it is challenging to investigate the dendritic selection stage of Purkinje cells in vivo using time-lapse imaging [29], resulting in very small data sets, we decided on a computational approach. Our primary goal is to investigate the interaction between the massive granule cell migration and growth and retraction of Purkinje cell dendritic trees, and to compare different hypothetical selection rules for selection of the winner dendritic tree.
We use the NeuroDevSim software [30], a parallelized agent-based simulation environment to model neural development. NeuroDevSim simulates development as consecutive cycles, in this study we simulate from P1 to P14 and a cycle corresponds roughly to 2.2 hours of growth. Neural growth always starts with a soma, that can be 'born' at any cycle and is either stationary (Purkinje cells, Bergmann glia) or migrates (granule cells). Growth of dendrites and axons is represented by cylindrical 'fronts' that act as independent agents, fronts are usually active during one cycle only. As is standard for neural growth algorithms [14,[31][32][33], new fronts can at the next cycle form one (extension) or two (branching) new fronts, that will become their children in the tree-like structure, or they can terminate growth. The choice is determined by neuron specific growth rules that compute probabilities for extension, branching or termination and by the local environment (see Materials and methods). An important effect of the local environment is crowding, fronts cannot overlap with other fronts of the same or different neurons. In addition, attraction and repulsion is also used for specific structures.
Using this computational approach, we first build a model of granule cell migration along Bergmann glia processes and of parallel fiber development that will provide an environment in which Purkinje dendrites can grow. Next, we introduce growth of multiple Purkinje cell dendrites starting at P4. Before we compare possible selection mechanisms, we investigate control scenarios with random selection, no selection, or no environment.

Granule cell migration model
In order to provide a 3D environment for Purkinje dendrites to grow, a granule cell migration model was built first. The model simulates a small volume of space in which we estimate about 4,500 granule cells and 24 Purkinje cells [34] will be present at adult stage. As we only simulate P1 to P14, about 3,000 granule cells migrate in this volume and form parallel fibers (Fig 2A) and Purkinje cell dendrites do not grow. The model does not simulate tissue expansion. To approximate this, granule cells are born in 12 consecutive phases from bottom to top (Fig 2E), simulating the expansion of the molecular layer pushing the external granular layer upward. Parallel fibers are much longer than the width of this volume, therefore, we simulate an additional 10,000 parallel fibers growing in from neighboring volumes that do not simulate Purkinje cell or granule cell growth. Fig 2B summarizes the sequences simulated during granule cell migration. First, a granule cell soma (small green dot) migrates horizontally to one of the closest Bergmann glia processes (light blue, more distant glia processes of same cell are dark blue) by extending a cylindrical leading process (gray). The granule cell will migrate following its leading process till it gets close enough to the target Bergmann glia process. There, it changes direction and starts a radial migration by extending the leading process downwards, staying close to the process. It also extends a short axon (pink) from its tail which will further bifurcate into parallel fibers. The soma keeps migrating down and parallel fibers further extend along the yaxis in both directions. A line of axonal fronts is laid down along the path of the radial migration of the soma.
In a representative model case, 80.59% of the granule cells succeeded in the radial migration, and 67.73% of parallel fibers extended through the complete volume of y-axis ( Fig 2C). Because of crowding, only 29.31% of incoming parallel fibers completed extensions through the main volume. Total parallel fibers in the central volume was 3,487 towards front and 3,430 to back. Distributions of success rates from 10 simulations are shown in Fig 2D.

Purkinje cell growth model
Purkinje cell growth in the model corresponds to development in a mouse from P4 to P14. It starts by growing five new dendrites on the upper sphere of each soma, of which four will be retracted in two phases, followed by a few days of growth of the winner tree ( Fig 1B). The initiation of dendritic growth waits till early granule cell migrations are completed, and starts at initiation of phase 6 granule cell migrations (Fig 2E, cycle 65). To reduce overall run-times, 10-20 simulations of the early granule cell migration phases were saved and used as the We assumed that a major goal for Purkinje cell dendrite growth is to establish synaptic contacts with parallel fibers [16,35,36]. Therefore, dendrites grow toward nearby parallel fibers. Once they reach a parallel fiber, a synaptic contact is made and growth continues in a direction perpendicular to the parallel fiber [37,38]. In addition, dendrites are repelled by other nearby dendrites and have a low probability of branching. The specific growth rules used are described in Materials and Methods. In NeuroDevSim, synapses do not have a physical structure but, because each front can only have one synapse, they can never be very close to each other. In most of the models, synapses are passive and used only to compare the effectiveness of dendritic growth on connectivity.
Two phases of retraction are implemented in all models [29]. In the first phase, three trees retract and the remaining two trees continue to grow. In the final phase a single winner tree is selected and the other one is retracted. Before exploring possible mechanisms that control the selection of winning dendrites, we performed three control scenarios to investigate the influence of dendritic selection processes and physical presence of granule cells on dendritic tree morphology.

Purkinje cell retraction controls
The control scenarios investigate how the growth rules affect Purkinje cell dendritic growth in the absence of specific retraction mechanisms. We simulated random selection of a winner, no retractions at all, or growth and retraction of dendrites in the absence of granule cells. Each of these scenarios is described in more detail below, and examples of dendritic growth and general comparisons are shown in Figs 3 and 4. Here only C2 shows 2 phases of retraction. In C1 and C3, all dendritic trees at the first screening phase were larger than the retraction threshold and kept growing until the second screening phase. (D) Average number of synapses on each winner dendritic tree in the three scenarios with error bars representing standard deviations. (E) Average maximum length of each winner tree in the three scenarios with standard deviations. (F) Average number of branch points in each winner tree in the three scenarios with standard deviations. P values: **** indicates p < 0.00005, Welch's t-test. Actual p values, t statistics, and degrees of freedom for comparing each data set are summarized in S1 Table. https://doi.org/10.1371/journal.pcbi.1011320.g003 The RandomRetract scenario randomly chooses a single winner dendritic tree out of five candidates at 15 hours after start of growth (cycle 72) and retracts all others (visible at the beginning in each plot in Fig 3A1-3A3). The NoRetract scenario never initiates retractions and all the dendritic trees in each cell keep growing until they reach z = 100 μm (Fig 3B1-3B3). The NoGranCells scenario does not use the granule cell migration model. Instead, Purkinje cell somata are initiated in an empty volume of the same dimensions and dendritic growth is initiated at cycle 5. After 2.3 days of growth (cycle 30), cells that achieve a growth threshold retract small dendrites. 15 hours later (cycle 37), the tree with the largest number of fronts is selected as winner (Fig 3C1-3C3).
We next compare properties of the winner dendritic trees in the three control scenarios. The average number of synapses with parallel fibers in winner trees shows differences ( Fig  3D). Most clearly, dendritic trees in NoGranCells did not have any synapses because there are no granule cells. Each winner dendritic tree in RandomRetract got a much larger number of synapses than in NoRetract, probably because growing all candidate trees in NoRetract led to more intense competition over limited numbers of synapse locations with parallel fibers. Moreover, in NoRetract, crowding by more dendritic branches in the simulation cube led to a decrease in the number of fully extended parallel fibers, which made the competition even more severe. Although experimental data to compare number of synapses on dendritic trees at the age of interest is not available, it seems the model provided enough opportunity for dendrites to form synapses with parallel fibers. Because retraction opens up space for growth of winner dendrites, they have a larger max path length and number of branch points (40.89 ± 45.2 versus 20.64 ± 25.2) in RandomRetract than in NoRetract ( Fig 3E). This may also have contributed to the increase in number of synapses. Winner trees in NoGranCells had an artificially large increased number of total fronts (Fig 3F), due to how the Purkinje cell growth rule deals with an absence of nearby parallel fibers.
Examples of the morphology of the winner trees and a comparison of their dimensions are shown in Fig 4. Flat dendritic trees, measured by the maximum y-distance of winner trees, depend on the presence of granule cells (Fig 4A1-4C1 and 4D). Only RandomRetract and NoRetract show well separated planar morphologies, conversely NoGranCells dendritic trees intermingle each other. This indicates that repulsion by nearby dendrites [39], present in all three scenarios, was not sufficient to separate dendritic trees. Instead, growth along directions perpendicular to close by parallel fibers was more effective [37,38] (see dendritic tree growth algorithm in Materials and Methods). Comparing the maximum y-distance of RandomRetract and NoRetract, RandomRetract has slightly thinner trees. Therefore, selection of a single winner tree seems to contribute to planar dendrites [20] and further growth of the winner tree did not result into expansion along the Y axis.
Although differences in maximum x-distance were less obvious, they were statistically significant ( Fig 4A2-4C2 and 4E). The RandomRetract and NoRetract scenarios resulted in much greater variation between samples than NoGranCells. Physical hindrance from granule cells possibly contributed to enlarge the divergence in the tree morphologies.
In adult control mice, y-distance of dendritic trees was 13.9 ± 1.5 μm and x-distance was 82.9 ± 2.5 μm [40]. Y-distance from NoRetract and RandomRetract are slightly larger than this observation, though Purkinje cell dendritic trees attains single planar morphology only after P22 in mice [23] which is later than the simulated growth. For growth in x-axis direction, the dendritic tree reaches its full x-width roughly by P13 in mice [20]. The x-distance of the three scenarios was marginally wider than experimental observation [18]. Confirming the importance of lack of granule cells in NoGranCells, larger y-distances (28.7 ± 4.8 μm) were observed in reeler mice where the structure of cerebellar cortex is disrupted due to a failure of granule cell migration [40].
The simulation results from these scenarios suggest that excess trees growing from the same cell lead to a reduced dendritic tree morphology by competition for space and synapses (NoRetract scenario), implying an important role for the retraction stage. Also, the physical presence of granule cells contributes to the distinct planar structure of Purkinje cell dendritic trees.

Possible mechanisms for Purkinje dendrite retraction
The mechanisms controlling the timing and selection of dendrites to retract are at present unknown, but neuronal and synaptic activity are required and retraction occurs in two stages [29]. In this study we simulated four simple retraction mechanisms and compare their outcomes, assuming that the number of parallel fiber synapses made by the winner trees is an appropriate measure of success [36]. In the FixedRetract scenario all Purkinje cells retract dendritic trees at the same times and winners are selected based on number of synapses on each tree. In the other scenarios, individual cells need to achieve a growth threshold before retractions are triggered, resulting in different retraction times, and winner trees are selected based on different criteria. In SizeRetract, the size of the entire dendritic tree, measured by its number of fronts, is used to control the timing of retractions and winners are selected based on the size of each tree. In SynapseRetract the number of synapses controls both timings of retraction and selection of winners. In the first three scenarios, the role of granule cells in promoting retraction is only structural by providing axons that can have synaptic contacts with Purkinje cells. In the final scenario, InputRetract, we assume that granule cell network activity also contributes. Specifically, granule cells have random firing rates, with larger firing rates for cells that completed migration to the granular layer and can receive synaptic input by mossy fibers (not simulated). The amplitude of synaptic input to Purkinje cells is granule cell firing rate dependent and integrated over time as a synaptic signal (see Materials and methods). The summed synaptic signals are used in InputRetract to control both retraction timings and selection of winner dendritic trees.
Each of the four retraction scenarios has two to three control parameters. In first instance, we performed parameter space evaluations to select the control parameter values that resulted in winner trees with the largest number of parallel fiber synapses (Fig 5 and Table 1).
There was a high variance in the number of synapses generated for different parameter combinations for all scenarios. While for scenarios FixedRetract and InputRetract clear trends were visible in the color maps, this was less the case for SizeRetract and SynapseRetract, resulting in noisier maps. Values for the best parameter sets (red stars in Fig 5) and average number of synapses on winners are listed in Table 1. Data from these best parameter sets were further analyzed in the following Figures.
Maximum x-distance and y-distance of winner trees were compared. Although especially y-distances showed great variability between cells, this was consistent for all parameter combinations (data for SynapseRetract shown in Fig 5E and 5F), and there were no differences with  Table 1. (A) FixedRetract scenario: the control parameters are cycles at which each retraction phase occurs. X-axis shows the first retraction cycle and y-axis shows number of additional cycles to trigger second retraction. Each point represents 10 or 20 samples. (B) SizeRetract scenario: has 3 control parameters. X-axis shows number of fronts required to trigger first retraction and y-axis to trigger second retraction. In addition, during first retraction, a parameter (named '1st_f_th' in Table 1) sets minimum number of fronts required to survive, 20 for this figure (range 20-160 evaluated). Each point represents 20 samples, surviving trees without synapses excluded. (C) SynapseRetract scenario: has 3 control parameters. X-axis shows total number of synapses required to trigger second retraction, and y-axis is minimum number of synapses required to survive first retraction. Total number of synapses required to trigger first retraction (named '1st_s_th' in Table 1)   the measurements for other scenarios. Compared to experimental data [40], the same conclusions apply as for the measurements in the control scenarios (Fig 4E and 4F). An important difference between the retraction scenarios is that FixedRetract forces two phased retractions at fixed cycles, while other scenarios can have retractions occur at different cycles for each cell (Fig 6).
Measurements that represent sizes of the winner trees, their number of fronts, branch points, and terminals, were strongly correlated with each other and therefore we show only the number of branch points (Fig 7A). InputRetract had significantly smaller trees than the other scenarios, which correlates with its smaller number of synapses (Table 1). Compared with control scenarios RandomRetract and NoRetract, winner trees from retraction scenarios have significantly more branch points (S3 Table), which suggests that proper selection of which trees to retract promotes growth of the winner dendritic trees. We found that the size of winner dendritic trees strongly correlates with the total volume of the parallel fibers in their direct neighborhood (Fig 7B-7E). This may be expected from the dendritic growth rule that promotes growth towards parallel fibers [38], but it also implies that the positive influence of parallel fibers is much larger than their negative crowding effect on growth.
Finally, we examined the variability of the number of synapses for the different retraction scenarios. In Fig 8 we show examples of how this variability develops over simulation time for different cells in each scenario. There were obvious differences between the scenarios in the range of number of synapses in the winner trees. The SizeRetract scenario clearly has a much lower variability (Fig 8B1), while that of InputRetract seems largest. Randomly selected Purkinje cells from these simulations are shown in A2-3, B2-3, C2-3, and D2-3. It is difficult to detect morphological differences between scenarios from the figures, but variation in morphology of Purkinje cell dendrites from the same simulations can be observed. Movies for the dendritic growths and retractions of selected cells in each scenario can be found in supplementary materials.  Table 1, except D1 and D2 are from a parameter set with Whole_signal_1st = 10,000 and Whole_signal_2nd = 13,000. https://doi.org/10.1371/journal.pcbi.1011320.g006 In Fig 9, the distributions of winner trees with a given number of synapses at the end of growth in all simulations are compared (See S4 Table for comparative statistics). All dendritic selection scenarios result in winner dendrites with significant higher number of synapses than RandomRetract, conversely NoRetract tends to have more synapses because it has 5 trees.
The variability in number of synapses of the InputRetract scenario ( Fig 9D, Table 1) is larger than that of other retraction scenario (Fig 9A-9C) and this is mostly due to a larger proportion of winner trees with small numbers of synapses. This is caused by the input firing rate  https://doi.org/10.1371/journal.pcbi.1011320.g008 dependent selection rule of InputRetract, which favors winner trees that mainly synapse onto parallel fibers extending in the bottom of the molecular layer. The somata of the granule cells of these parallel fibers reach the granular layer first, strongly increasing their average firing rates.

Summary and discussion of the main results
While it was challenging to achieve reasonable success rates for granule cell migration and parallel fiber growth in the very crowded environment of the molecular layer in the model (Fig 2), crowding was less of an issue for Purkinje cell dendritic growth. In fact, the presence of numerous parallel fibers was essential for planar growth of the Purkinje cell dendrite (Fig 4) and promoted the growth of larger dendritic trees (Fig 7). Conversely, dendrite-dendrite repulsion was not sufficient to generate flat dendrites. It should be noted, however, that these observations are emergent properties of the growth algorithms used, which were not varied systematically.
Comparisons of the control scenarios with the retraction ones allows to suggest possible roles of the dendritic selection stage. Winner trees based on selection were larger and had more synapses than those generated by RandomRetract (Figs 4, 7 and 9). They were also larger than those obtained in NoRetract (Figs 4 and 7), but no retraction did result in a larger number of synapses on the dendrite (Fig 9F).
The model then explored detailed mechanisms of dendritic selection comparing four retraction scenarios. In the absence of experimental data about the exact mechanisms controlling the selection of a winner dendrite, we assumed that an optimal developmental sequence generates the largest number of parallel fiber synapses onto a Purkinje cell dendritic tree. This is based on the uniquely high number of synaptic inputs on these neurons [36,41] and on the effects of up-or down-regulation of parallel fiber synapses on dendritic morphology [23,28]. We then devised scenarios that assume that the retraction process is primarily controlled by a genetic clock identical to all cells (FixedRetract), by size or synapse dependent Purkinje cell metabolic mechanisms resulting in variable retraction times (SizeRetract and SynapseRetract) or by network activity in the granule cell layer (InputRetract). Based on studies of the development of other brain regions [42,43], the afferent network activity InputRetract scenario seemed the more plausible one. This would also match the  Table. https://doi.org/10.1371/journal.pcbi.1011320.g009 late, postnatal development of cerebellar structures, which has been hypothesized to reflect the requirement of structured synaptic input to establish correct synaptic connectivity in cerebellum [44]. However, in our simulations, InputRetract was less successful than the other three scenarios, resulting in smaller winner trees with less synapses. This was due to a specific property of cerebellar development: the input granular layer is still forming while output Purkinje neurons grow their dendrites. A large fraction of early parallel fiber synapses will, therefore, be made onto axons of granule cells that have not yet reached their final position or had time to grow dendrites and establish mossy fiber contacts. In other words, around P6-P8 when dendritic tree retractions occur, the input layer is structurally very heterogenous, and this is probably reflected also in its firing rates. As a result, in the model, the Inpu-tRetract scenario selected winner trees that primarily made synapses in the lower region of the molecular layer. Conversely, the other retraction scenarios used number of synapses (Fix-edRetract and SynapseRetract) or dendritic size (SizeRetract) to select the winner tree and this resulted in less variable outcomes with more synapses (Fig 9). Therefore, we postulate that, in cerebellum, the dendritic tree selection cannot be based on afferent synaptic activity alone. Postsynaptic metabolic mechanisms that are sensitive to tree size and/or number of synapses are expected to play an important role, though we cannot include an additional effect of afferent input.

Limitations of the early cerebellar development model
A major limitation for model development was the lack of experimental data. There are only a few reconstructions of normal Purkinje cell morphologies at P10-P14 available [45], insufficient for a systematic statistical comparison with the model morphologies [33,46,47]. Neither is data about the development and properties of parallel fiber synapses in normal mice at this critical age available.
The largest difference between the model and biological systems is that cerebellar expansion during development was not implemented, because volume expansion cannot be simulated in NeuroDevSim or any other neural development simulator [31,48]. The simulator also does not support growth of front diameters. The cerebellum undergoes a massive expansion, especially along the anterior-posterior plane (x-axis of the simulation cube) [49] during the first two postnatal weeks in mice [50]. While the anterior-posterior plane in an embryonic days 17.5 mouse has a length of half the width of its molecular layer axis (z-axis of the cube), it becomes in postnatal mice 7.8 times longer than the molecular layer which also expands significantly [49]. The lack of volume expansion strongly increases crowding of cells, making it more difficult to simulate migration and growth of large cell populations. Growth may also stretch and remodel existing neural structures [31,51], an effect not included in the simulations. It has been suggested that retraction of Purkinje cell branches may contribute to making the dendrite more planar, but this occurs at later ages (P18 to P25) than simulated here [23].
There is a difference in the process of extending parallel fibers during granule cell migration between the model and real neurons. In murine cerebellum, parallel fibers already extend during horizontal migration. However, the model extended parallel fiber axons after horizontal migration (Fig 2B). The purpose of horizontal migration in biological system is unclear, but it might be related to granule cell allocation upon expansion of the cortex [52], which is not implemented in the model. The model ignores the formation of synapses onto the descending part of the granule cell axon [53] and did not include development of inhibitory interneurons [54].
For Purkinje cell models, apoptosis and reorganizations of soma positions in Purkinje cell layer were not simulated. Also, the distinction between spiny dendrites and smooth dendrites is not made in the model. Climbing fibers and mossy fibers, two afferents to cerebellar cortex forming synapses with developing Purkinje cells [55,56] are not included. Climbing fibers are involved in organizing the final shapes of Purkinje cell dendrites [57]. A reduced branching pattern in dendritic arborizations of Purkinje cells is observed in rats when climbing fibers were removed by thermocoagulation in vivo [58]. Similarly, transient contacts of mossy fibers with developing Purkinje cells are likely involved in assembly of zonal circuit maps of the cerebellum rather than dendritic selections [59,60]. However, there is no strong evidence that they are involved in early dendritic development, and inclusions of these axons may be considered in future models. Finally, development was simulated only until P14, excluding later growth of the winning dendritic trees from analysis.

Potential improvements and future directions of the models
In the model, all granule cells are homogeneous in that there is no zonal patterning by molecular phenotype or birth order [52]. Such inhomogeneity in granule cells can be represented by, for example, 2 different types of molecular diffusions in the model to facilitate the zoning. Also, we can set different types of granule cell objects and assign different preferences in making synapses with distal or proximal regions of Purkinje cells' dendritic trees as observed [53,61]. It is technically also possible to include granule cell proliferation, but in the absence of tissue expansion in the model this did not seem useful.
For the dendritic structure of Purkinje cells, an additional growth algorithm to facilitate elaboration of the proximal region of dendritic trees should be introduced to match observed morphologies. Also, while dendrite to dendrite repulsion was included in the growth algorithm [39], it was not studied systematically nor did we distinguish auto-repulsion from repulsion by other cells. Neither was the possible attraction of dendrites by Bergmann glia implemented [62], though one cannot exclude that the anatomical proximity reported was due to nearby descending granule cell axons instead of glia itself. Recent work on local versus global effects on dendritic morphology [63] may also provide a useful framework to analyze the typical branching patterns of Purkinje cell dendrites. In that regard, our Purkinje cell morphologies are quite similar to those of Luczak [64], which suggests that the parallel fibers in our model act as 'neurotrophic particles'.
The present model is the first attempt to reconstruct early development of the main neurons in the cerebellar cortex with detailed interactions between several types of cells. The real cerebellum has evolved to be very complex, and experimental approaches to reveal how it develops are obviously important but limited by technical constraints. Modeling helps to reorganize and synthetize the experimental data and facilitates to understand how development works at different time scales.

Materials and methods
The NeuroDevSim software version 1.0 [30] was used to construct and simulate the models (https://github.com/CNS-OIST/NeuroDevSim). Simulations were executed on the high-performance computing cluster Deigo at OIST, using AMD Epyc 7702 CPUs with 128 cores and 512GB memory. 64 cores were used for each simulation, and it took about 3 hours to finish a single simulation. Python scripts for each model can be found in the Model Database (https:// modeldb.science/267591). Python notebooks and scripts and NeuroDevSim database files are available at https://zenodo.org/record/7771504.
The model embodies a Purkinje cell layer and a molecular layer in a central volume of (x1, y1, z1), (x2, y2, z2) = (-20, -20, -20), (180, 160, 140) μm 3 , representing a part of a murine cerebellar cortex (lobule VI) during postnatal days 1 to 14. The x-axis of the cube represents the sagittal plane, y-axis represents the transversal plane (long axis of folium), and z-axis is the depth of the cerebellar cortex (Fig 2A and 2E). The volume of the cube is fixed, therefore cerebellar tissue expansion is not implemented in the model.
The Purkinje cell layer has somas of Purkinje cells and Bergmann glia. The molecular layer has simplified Bergmann processes for guiding granule cells which gradually fill up the layer with their descending axons and parallel fibers ( Fig 1A) and growing Purkinje cell dendrites (Fig 1B). Additionally, we simulated the in-growth of parallel fibers from neighboring regions, resulting in a total simulation volume of (-20, -160, -20), (180, 300, 140) μm 3 .
All growth follows a standard NeuroDevSim procedure: the end position of a new front is calculated and this is created as a new child, provided it does not collide with existing structures (Algorithm 1 in S1 Fig). In case of collision different solutions are used, dependent on which structure grows.

Bergmann glia
On first cycle each Bergmann soma generates seven 4 μm long root processes that are spread along the x-axis. These then grow upward in 4 μm steps till they reach z = 131 μm (width of molecular layer at P11 [65] at cycle 35). Because Bergmann glia are first to grow, collisions are not an issue.

Granule cell migration model
About 3,000 granule cells are gradually initiated in 12 consecutive phases in the central volume ( Fig 2E) along with incoming parallel fibers of 10,000 granule cells located outside this central volume (orange structures in Fig 2E). Each phase initiates around 250 granule cells in the central volume and about 850 in-growing parallel fibers at the neighboring volumes. Upon initiation, granule cells migrate horizontally to a nearby Bergmann glia process ( Fig 2B) and start downward radial migrations after they get close enough to the process (Fig 2B) (Algorithms 2 in S1 Fig and 3 in S2 Fig). During the radial migration, each granule cell soma extends an axon which further bifurcates as parallel fibers (Fig 2B). The cells keep migrating until they reach the internal granule layer (z = -15μm), while parallel fibers extend up to the ends of y-axis ( Fig  2B). Granule cells start with a low random firing rate (0-0.2 / cycle) that strongly increases when they arrive in the internal granule layer (0-10 / cycle). The incoming parallel fiber structures simulate only parallel fibers themselves as a non-migrating soma with fibers extending along the positive or negative y-axis.
Granule cells in the model often experience collisions with surrounding structures during their migrations and extensions of axonal fibers. Built-in methods in NeuroDevSim are used to manage collisions; providing a bypass (solve_collision()) for migration, and finding an alternative space around a colliding structure (alternate_location()) for the fiber extensions. Also, the granule cell soma radius in the model is smaller than the actual value since physical shoving between cells cannot be simulated in the model, the smaller diameter representing the 'squeezed' soma structure of actual cells. The basic growth algorithm (Algorithm 6 in S4 Fig) simulates growth of a tip of dendrite towards a nearby synapse free parallel fiber segment. Once it gets close to the target, it makes a synapse and continues to grow perpendicular to the parallel fiber as observed in biological systems [37,38]. If no free parallel fiber segment is nearby, the dendrites continue growth along the direction of their current heading with a small upward force assuming phenomenological involvement of interneurons in the molecular layers [66]. When dendrites experience collisions, the code use the built-in solve_collision() method to find a detour to reach a destined position. Repulsion force by the closest dendritic tree of other cells or of the same cell also affect the growth direction, with strength of repulsion depending on distance along the y-axis to the closest dendrite segment (Algorithm 6 in S4 Fig). During the elongation process, random branching events occur with a small probability (Algorithm 7 in S5 Fig). A dendrite tip also initiates a branch towards its target parallel fiber if the direction to it strongly deviates from its current heading direction (Algorithm 8 in S5 Fig).

Control scenarios
Three control scenarios were simulated: • RandomRetract: Purkinje cells skip the dendritic selection stage by randomly selecting a single primary tree out of 5 candidate trees for each Purkinje cell 7 cycles after initiating growth.
• NoRetract: never triggers retractions of the candidate trees resulting in a physically more crowded environment and increased competition for synaptic locations.
• NoGranCells: Purkinje cells grow in an empty cube without granule cell migrations. The retraction mechanism for this scenario used that of the SizeRetract scenario so that they can initiate retractions without granule cells. At cycle 30 any cell with more than 250 fronts will retract trees that have less than 60 fronts. At cycle 37 cells having at least 400 fronts will select their largest tree and retract all others.

Retraction Scenarios
To select the surviving dendrite, four different retraction scenarios are implemented in the model: • FixedRetract triggers dendritic retractions at fixed cycles. It uses two fixed simulation cycles to trigger retractions. At first cycle, three out of 5 trees are retracted and later the two remaining compete to become winner(s). Winner trees are determined by relative numbers of synapses for each dendritic tree in a cell. Parameters: the two cycles at which to retract, same for all cells.
• SizeRetract uses an internal parameter of growth, the size of each branch measured as its number of fronts. Granule cells merely acts as physical obstacles to dendritic growth. Cycles at which to retract are determined by cell thresholds for the summed sizes of all trees and vary between cells. At first retraction any tree smaller than a tree threshold is retracted, as a consequence a variable number of trees remains. Parameters: two size cell thresholds to start retraction and tree threshold.
• SynapseRetract uses the presence of synaptic connectivity. Cycles at which to retract are determined by cell thresholds for the number of synapses with parallel fibers of all trees and vary between cells. At first retraction any tree with less parallel fiber synapses than a tree threshold is retracted, as a consequence a variable number of trees remains. Parameters: two synapse cell thresholds to start retraction and tree threshold.
• InputRetact uses synaptic input determined by granule network activity levels. Synaptic input is integrated as an input signal proportional to afferent firing rate (equation in Algorithm 4 in S3 Fig) and relative total signal per tree is used to determine winner trees. Cycles at which to retract are determined by cell thresholds for total signal summed over all trees and vary between cells. Parameters: two signal cell thresholds. 24 Purkinje cells enclosed by outer Purkinje cells were sampled for analysis (Fig 2F). For all scenarios, morphology of resulted dendritic trees was analyzed and compared. Pairwise statistics were made using the Welch's t-test.

Data analysis
The following properties were computed: • number of winner trees • size of winner trees: computed as number of fronts because dendritic fronts tend to have uniform length of~5 μm • number of synapses • maximum x-and maximum y-distance: computed for each tree as difference between minimal and maximal x or y coordinate value for the end point of all fronts • number of branch and terminal points of winner trees.
The correlation analysis in Fig 7B-7E used the following algorithm: • compute bounding box around Purkinje cell as minimal and maximal x or y coordinate value for the origin and end points of all dendritic fronts, with z in 17-100 μm • sum volume of all cylindrical Purkinje dendrite fronts • detect all parallel fiber fronts that have origin and/or end point in bounding box, compute their volume for the part that is contained in the box and sum • correlate the two summed volumes for all cells.
All movies show only a few Purkinje cells: cells 21 in red, 24 in blue, 29 in green, and 33 in yellow (see Purkinje cell location map on Fig 4D). Granule cells, parallel fibers and Bergmann glia are not drawn. For each Purkinje cell, the dendritic trees in color are final primary trees, while ones that will be retracted are drawn in black. The movies start from simulation cycle 61 and ends at cycle 140. S1 Movie. Dendritic growths and retractions movie for control scenario RandomRetract (grow one tree). Seed 2 of scenario RandomRetract simulation was used to generate the movie.