Simulating Cortical Development as a Self Constructing Process: A Novel Multi-Scale Approach Combining Molecular and Physical Aspects

Current models of embryological development focus on intracellular processes such as gene expression and protein networks, rather than on the complex relationship between subcellular processes and the collective cellular organization these processes support. We have explored this collective behavior in the context of neocortical development, by modeling the expansion of a small number of progenitor cells into a laminated cortex with layer and cell type specific projections. The developmental process is steered by a formal language analogous to genomic instructions, and takes place in a physically realistic three-dimensional environment. A common genome inserted into individual cells control their individual behaviors, and thereby gives rise to collective developmental sequences in a biologically plausible manner. The simulation begins with a single progenitor cell containing the artificial genome. This progenitor then gives rise through a lineage of offspring to distinct populations of neuronal precursors that migrate to form the cortical laminae. The precursors differentiate by extending dendrites and axons, which reproduce the experimentally determined branching patterns of a number of different neuronal cell types observed in the cat visual cortex. This result is the first comprehensive demonstration of the principles of self-construction whereby the cortical architecture develops. In addition, our model makes several testable predictions concerning cell migration and branching mechanisms.


Introduction
High-throughput quantitative methods in molecular biology, such as DNA microarrays, are generating exponentially increasing information about cellular mechanisms. The need to organize these volumes of raw data, and transform them into a explanation of overall cellular function has accelerated interest in approaches to characterizing systems-level biological principles [1,2]. Generally, these methods of analysis are drawn largely from mathematical formalisms developed over decades in chemistry and biochemistry (e.g. the law of mass action, enzymatics) as well as from engineering (e.g. systems theory). They permit Systems Biologists to describe very compactly processes such as gene expression and protein interactions by using differential equations [3]. The numerical methods required to solve the resulting expressions are also well understood and widely accepted, and they are easily automated on ever more powerful computers.
While these approaches have been very successfully applied at a sub-cellular level, their application to the complex cellular interactions of tissue or organ level behavior has been more difficult and less well studied [4]. For example, the literature lacks appropriate formalisms to express the effect of specific gene expression on the mechanical properties of cells, or on their division, migration and morphological differentiation. To study the effects of genetic control at the collective cellular 'organ' level, new types of model mechanisms are required. These models should encompass the genomic and proteomic, as well as the active and passive physico-mechanical properties of cells, and offer insights into the collective synergystic behaviors of ensembles of cells at the tissue level. Large-scale agent-based simulations have been used previously to study the development of simple organisms [5,6], or specific organs (such as blood vessels [7], pancreas [8] or limb bud [9]) from a limited number of undifferentiated precursor cells.
Here we explore these questions in the context of neocortical development. The development of cortex is particularly interesting because it results in a complex yet precise architecture of connections between neurons on a wide range of spatial scales, and so provides the substrate for the meta-level of electrophysiological information processing that supports intelligent behavior. Our approach to bridging this important gap between molecular processes and cell behavior is by large-scale simulation of physical cellular mechanism.
We have previously described our simulation framework, CX3D, whereby the cellular mechanisms of brain development can be explored [10]. CX3D respects physical processes such as cell division, cell-cell interactions, movement and chemical diffusion in three-dimensional space. Each cell is an autonomous agent exerting only local actions, and using only locally available information. The behaviors of the simulated cells are determined by intracellular molecular-gene-like codes that are expressed according to intracellular or extracellular conditions. For specifying this 'genetic code' we have introduced G-code [11], a formal language based on a set of primitive neural actions, which can be combined to form networks of instructions ('G-machines'), so providing complex cellular behaviors. The same primitives can be used to specify regulating networks similar to the traditional systems biology approach, or to design models of growth cones for the elongation of long cellular processes and the formation of synaptic connections. The overall behavior of the cell at any time is the collective result of its currently instantiated G-machines (see short introduction to G-code at the beginning of Methods section).
Here we demonstrate how this instruction language can be used to control the physical development of the neocortex in a selfconstructing manner. Beginning with one progenitor cell, the code directs the mitotic sequence, differentiation, migration, and neurite formation of some 4000 neurons of different types, distributed according to the observed lamination of cortex, and with appropriate intra-and inter-laminar axonal projections. The resulting neurons approximate well the experimentally determined morphology of several cell types found in cat visual cortex (pyramidal cells of layer 2/3, spiny stellate cells of layer 4, pyramidal cells of layer 5, pyramidal cells of layer 6, and of one type of inhibitory cells, the basket cells of layer 2/3). Importantly, this developmental process is neither specified nor controlled by an external co-ordinator, but rather arises only out of local interactions between successive generations of cells. Our primary result is this first comprehensive demonstration of cortical development, based on abstract principles of biological selfconstruction. In addition, our model makes several testable predictions on cell migration and branching mechanisms, and points the way the to new types of simulations for studying the effect of biochemical parameters on the development of tissues.

Organization of the simulation
We decompose corticogenesis into three phases: (1) The formation of the preplate (composed of the marginal zone and the subplate); (2) The formation of the cortical plate (composed here of four cortical layers in addition to the marginal zone); (3) the formation of axonal and dendritic arborizations.
Corticogenesis begins from a single layer of progenitor cells at the anterior part of the neural tube called the ventricular zone (VZ). Since we do not yet model all the embryonic steps occurring before the beginning of corticogenesis, the initial phase of our simulation begins with a single progenitor cell. This cell produces a single-layered VZ through symmetric division of a precursor cell. After the formation of this single-layered structure, the mode of division becomes asymmetric, increasing the pool of precursors in the VZ cells and forming several new types of cells [12]. The next structure to emerge is the preplate (PP): asymmetric divisions of the VZ cells form first the marginal zone (MZ), that will become the future Layer 1 followed by the subplate (SP). This part of cell production occurs in an outside-in fashion, with older cells lying superficially relative to the younger deep cells [13]. In our model, the preplate is generated by a single intracellular machine that is instantiated across all cells involved in preplate formation. This machine senses its depth in the mitotic sequence though the successive dilution by cell division of two intracellular substances inherited from the original progenitor cell.
During preplate development in biology, several other types of cell appear in the VZ. We consider only two important types in this simulation: First, the radial glial cells (RGC) [14] that will be used by neuron precursors to guide their migration; and second, the intermediate progenitor cells of the sub-ventricular zone (SVZ) that provide a separate germinative compartment lying between the VZ and the PP [15,16]. Differentiating neurons arising from the VZ precursor pool will form the lower layers of the cortex, whereas the SVZ precursor pool form the neurons of the more superficial layers [17]. During this phase, the daughter cells of progenitor mitosis progressively favor differentiation over further mitosis. The resulting precursor neurons migrate radially from the VZ through the SP and stop before entering the MZ. As more precursors arrive, they form a clearly visible layer -the cortical plate (CP). This developing plate pushes the MZ further upward, splitting it away from the underlying SP.
In our simulation, the formation of these different neuronal types depends on a gene regulatory network, whose expression state triggers cell type specific behaviors. The gene regulating network produces the precursors of the different neuronal types sequentially. As each type is produced, they migrate through and past their predecessors to form the different layers of the cortex in an inside-out manner: First layer 6 neurons, then layer 5, 4, 3, and finally 2 [18]. In addition to the formation and radial migration of presumptive excitatory neurons, inhibitory interneurons are formed in the lateral and medial eminences of the future basal ganglia, and they must migrate tangentially to take up position in the cortex. We simulated this tangential migration for only one type of inhibitory cells, as proof of concept. The production of these inhbitory neurons is not yet governed by our specific gene regulatory network.
Once at their final location, or during migration, the neurons differentiate further by extending neurites (axons or dendrites) and producing cell type specific branching patterns. In real neurons, neurite extension is achieved by a sensitive motile structure, the growth cone, located in the tip of the neurite [19,20]. This growth cone decides on the basis of local internal conditions and external physical and chemical cues (e.g. diffusible or membrane-bound signaling molecules) whether to extend, retract, or bifurcate the neurite. We have designed several growth cone machines in Gcode, reproducing typical features of the branching patterns observed in various cortical cells [21,22]. We considered four types of excitatory cells: pyramidal cells of layer 2/3 (P23), spiny stellate cells of layer 4 (SS4), pyramidal cells of layer 5 (P5), pyramidal cells of layer 6 (P6), and of one type of inhibitory cell, the basket cells of layer 2/3 (B23).

Author Summary
The proper operation of the brain depends on the correct developmental wiring of billions of neurons. Understanding this process of living self-construction is crucial not only for biological explanation and medical therapy, but could also provide an entirely new approach to industrial fabrication. We are approaching this problem through detailed simulation of cortical development. We have previously presented a software package that allows for simulation of cellular growth in a 3D space that respects physical forces and diffusion of substances, as well as an instruction language for specifying biologically plausible 'genetic codes'. Here we apply this novel formalism to understanding the principles of cortical development in the context of multiple, spatially distributed agents that communicate only by local metabolic messages.

Preplate formation
An important aspect of biological development is the manner in which progenitor cells replicate to produce nearly exact copies of themselves, but are also able to form cells of different types. This process can be viewed as a lineage-tree that has the progenitor cell as its root, and recognizable phenotypic types as leaves [23]. The symmetry breaking needed for the formation of different cell types can rely on two different mechanisms: extrinsic (asymmetric exposure of daughter cells to extracellular factors) or intrinsic (asymmetric partition of intracellular markers between the two daughter cells) [24]. We make use of the latter principle for controlling the formation of the PP (Figures 1, 2), using the direction of division to regulate the asymmetric partitioning of intracellular factors [25]. In our model, the initial progenitor contains two intracellular substances, 'X' and 'Y' in precise concentrations. X is distributed symmetrically in every kind of division (whatever the axis of division is), while Y is distributed symmetrically when the axis of division is perpendicular to the internal axis of the cell, and asymmetrically if the division axis is aligned with the internal axis, the daughter on the 'south' side inheriting all of the substance. By initially choosing a division axis perpendicular to the internal cell axis, we form a homogeneous plane (the VZ). The number of divisions is controlled by the dilution of the intracellular substance X occurring at each division. When the concentration of X falls below a given threshold T 1 , the cells align the division axis with the internal axis. This results in an asymmetric division, since the substance Y is now transmitted only to one of the daughter cells (the 'southern' one). Cells which do not contain any Y are not allowed to divide further, so that only one cell continues to divide. Since the molecule X is still distributed equally across both daughter cells, its concentration further decreases, defining the total number of divisions allowed, as well as the type of cell the 'northern' daughter becomes (type L1 if the concentration of X is greater than a second threshold T 2 , type SP otherwise). The initial concentrations of X and Y and the different threshold values specify the number of cells of each type that are produced. The division rate is constant, and independent of the intracellular substances. Division cycles stop in cells with no Y. When the X concentration falls below a concentration T 3 , the Gcode machine for the preplate formation is removed, and the next machine -responsible for the cortical plate formation -is instantiated.
As already mentioned, the formation of the PP from an initial precursor is not biologically accurate, but had to be implemented because we do not simulate the stages of embryogenesis taking place before the beginning of corticogenesis. Even if there is no experimental evidence supporting the control of cell fates in the PP through dilution over several divisions, we used this model because of its simplicity and its theoretical importance [26],

Cortical plate formation
After preplate formation, the VZ progenitors begin to generate the neurons of the layered cortical plate. This developmental phase has two aspects: (1) the generation of the approriate cell lineage, i.e. the sequential generation of the correct numbers and types of neurons; and (2) their migration to the correct position so as to form the appropriately layered cortical organization.
Cell lineage. Labeling cell state and time by simple dilution and partitioning of substances may be useful for guiding the first few steps of development such as the construction of the preplate. But biological cells are likely to have practical limitations on precision of measurement, the accurate partitioning of marker substances, and the implementation of thresholds. Consequently, the dilution approach quickly becomes unreliable as the depth of the lineage tree increases. Therefore, for producing the different cortical cell types, we use a simplified gene regulatory network [27][28][29] (GRN; Figure 1). The activation profile of GRN genes defines the cell type, and can trigger the expression of specific machines (Table 1). Throughout this phase, the same cell cycle machine is used as during the PP formation.
Our GRN is composed of sequentially activated bi-stable switches, that is pairs of self enhancing and mutually inhibiting genes a i and b i coding for transcription factors that have differential distributions during division, one daughter inheriting more a i , the other receiving more b i . The competition between these paired genes accounts for a branching in the cell lineage tree, whereby the state space of gene expression is partitioned in two basins of attractions. For instance the first branch point, at which cells either become P6 neurons or stay in a proliferative state as VZ cells, is under the control of the genes a 6 and b 6 . Initially there is the same quantity of protein a 6 and protein b 6 in the cells (the expression of which is controlled by two additional genes s 1 and s 2 active at the beginning of this phase). After the first division, due to the asymmetric distribution of internal proteins, one daughter cell receives more a 6 , the other one more b 6 . This asymmetry is reinforced over time by the network dynamics and by further division cycles, driving the expression state toward one of the two attractor states (defining cell types). Cells in which the a 6 protein exceeds a given threshold will exit the cell cycle, migrate and acquire the morphological characteristics of P6 cells. In the other cells, the next pair of mutually inhibiting genes (a 5 and b 5 ) is expressed, determining which cells will become L5 neurons (if a 5 ''wins''), and which cells will divide further. The neurons of the superficial layers are known to derive from a second pool of precursor cells, in the supraventricular zone (SVZ) [30]. Similarly, in our model, cells in which b 5 exceeds a given threshold become SVZ cells. Due to the asymmetrical partition of the gene product, the threshold is often reached by division. Reaching the b 5 threshold triggers the expression of the next bistable switch a 4 /b 4 . The intensity of selfexcitation and mutual inhibition, the degree of asymmetry in distribution at cell division and the value of the threshold used regulates the number of cells produced. Using the parameters described in the Methods section the following cell numbers were produced: P6: 560 cells (18%), P5: 805 cells (25.8%), SS4: 980 cells (31.5%), P2/3: 770 cells (24.7%); see also Table 2 for comparison with biological data.
In many non-primate mammals, a large proportion of the cortical inhibitory cells are produced in the lateral and medial eminences, and migrate tangentially into the cortex [31][32][33]. Such cells are not yet produced by our GRN. As a proof of concept, we nevertheless incorporate the migration of one sort of inhibitory cell (the basket cell of layer 2/3). We insert ten cell bodies at about 300 mm from the borders of the developing cortex. These cells have a simple chemoattraction mechanism that enables them to follow the gradient of a signaling molecule produced by the P2/3 cells, and so migrate tangentially and upward toward L23.
Obviously, the organization and operation of our GRN is only an approximation to the more complex multistable transcriptional networks expressed during neural development (Sonic Hedgehog, Notch-Delta, Hes, etc.). However, there are no models of the detailed biological systems that have been shown to be capable of orchestrating development as the one reported here is. It is also worth noting that there is experimental evidence for the existence of pairs of self-enhancing and mutually inhibiting genes in hematopoetic cell lineages [34]. Pairs of self enhancing and mutually inhibiting genes form the most simple gene regulatory network motif that incorporates cooperation and competition, thus guaranteeing a bistable behavior (for a theoretical analysis see [35]).
Migration. The neuron precursors formed in the proliferative zones (ventricular zone, subventricular zone) use several modes of radial migration to reach their final position and form the cortical plate [36][37][38]. In the present work we consider only locomotion along radial glial processes (RGPs). RGPs are extensions of some of the progenitor cells attached to the pial surface, and are thus being elongated when the distance between the pial and the ventricular basal membrane increases. However, due to technical difficulties with the stability of the RGPs in the current version of our simulation framework CX3D, we had do make two simplifications before the start of the cortical plate formation: First, the RGPs do not attach to the pia but may extend beyond it. Second, we remove the subplate cells ( Figure 3A), whereas these cells normally disappear at later stages of development. More realistic mechanisms will be incorporated in future versions of our simulations.
As in biology, the newly formed neuron precursors begin to migrate randomly as soon as they stop dividing. During this The initial progenitor in the ventricular zone (VZ) contains precise amounts of the intracellular substances X and Y. Through symmetric and later asymmetric divisions with differential partition of the substances in the daughter cells, the VZ progenitor pool increases, and the future layer 1 cells (L1) as well as the the subplate (SP) are formed from cells that have left the proliferative state. Cortical plate formation: After the preplate formation, the VZ cells which are still in a proliferative stage express a gene regulatory network (GRN) for the sequential formation of neuron precursors. The GRN consists of pairs of self enhancing, mutually inhibiting genes, sequentially activated at each branch point in the lineage tree, and deciding between further division (in the VZ and later in the subventricular zone -SVZ) or exit from the cell cycle and differentiation into a cortical neuron. (X, Y: intracellular substances; T1, T2: thresholds on intracellular substances; a6, b6, a5, b5, a4, b4, a23, b23: genes. See text and Table 1 for details). doi:10.1371/journal.pcbi.1003173.g001 random motion they encounter a RGP and attach themselves to it. Then the precursors migrate along their RGPs, through the preplate, and stop before entering the marginal zone (MZ, i.e the future layer 1 - Figure 3B). Their accumulation pushes the MZ further up. The successive generations of neuron precursors migrate past their predecessors, leading to an inside-out development of the cortex ( Figure 3C-F): first layer 6, then layer 5, then layer 4, and finally layer 2/3 (in our model as in many mammalian cortices we consider layers 2 and 3 as a single layer). In the present model a stopping molecule is produced by the cells in L1. (One such molecule is the protein Reelin, produced mainly by specific cells in the MZ, and necessary for the migrating neuroblast to detach from the radial glial process. The absence of this protein leads to severe disorganization of the cortical architecture [39][40][41]). We have shown in previous simulations (see Figure 7 in ref. [10]) that a stopping mechanism due to contact with L1 is sufficient if the neuron precursors are really produced in distinct separate waves (first all the L6 cells, then all the L5 cells etc.). However our current GRN, as seems to be the case in biology [42], generates the different cell types with a certain degree of overlap. If the hypothesis of overlapping production phases holds, a stopping signal expressed only at the border of L1 is not sufficient to provide correct lamination, since later born cells belonging to deep layers could pass by earlier born cells belonging to more superficial layers. We thus postulate the existence of additional stopping mechanisms which prevent a cell from going through already settled cells of the next layer, or alternatively the existence of correction mechanisms such as a backward migrational movement toward white matter. As far as we know such mechanisms have not yet been described. We have implemented a version of the first proposed mechanism (additional stopping signals): neurons in our model express a type-specific membrane cue; in addition, cells which stop migrating express a further membrane cue that is non-specific for cell type. The additional stopping signal is triggered if the migrating neuron precursor is in contact with a sufficient number of cells of the same type that have already stopped their migration. When the MZ or the additional postulated stoping signal is detected, the neuron precursor detaches itself from the RGF, and differentiation into a cortical neuron can begin.
As in biology [43], cells in our simulation can still end up in the wrong layer. To be able to quantify this phenomenon we have used the following criterion: neurons that are in contact with less than three other cells of the same type are considered as being in a wrong layer. With this criterion, 18.4% of the neuron precursors were in the wrong layer. To explain this relatively high rate of migrational defect one can of course incriminate the oversimplification of our biological migration model. However we found that not only biological, but also mechanical aspects of the simulation have a important influence on the quality of the lamination. For instance the parameters chosen for the cell density, for the cell-cell interactions and/or the friction coefficient of cells have a crucial  Table 1. Conditional instantiation of cell type-specific machines based on concentration of intracellular substances and gene activity. influence on the passive displacement of cells which have already stopped migrating. The border conditions at the edges of the simulated volume also have a disruptive effect in allowing cells to avoid the stopping signal or because of a non-equilibrium of forces. These effects are particularly strong in our case because of the relatively small dimensions of the simulated cortical column.
To preserve the correct steering of the axonal and dendritic growth cones, the misplaced cells in our simulated cortex do not secrete any guidance cue. The biological justification is the observation that certain cells that have migrated to the 'wrong' layer do not acquire a normal phenotype [43].

Neuron differentiation
Axonal projections of cortical neurons form specific patterns of axonal and dendritic branching, often in a layer-dependent manner [44,45] (Figure 4A). For instance P23 cells have their main axonal shaft going down toward the white matter, forming tangential collaterals in L2/3 and in L5, but not in L4 or L6;  whereas P6 cells emit branches in L6 which move up to L4 where they ramify (neither the main axonal shaft nor the collaterals branch in L5) [46][47][48]. These differences are due to the layerspecific expression of signaling molecules, such as EphrinA5 which is present in L4 and in L6 and acts as a repulsive signal for P2/3 axons and a branch promoting signal for P6 [49]. In addition, guidance cues such as semaphorin3A [50,51] span the whole thickness of the cortex, and guide the movement of the axon and of the apical dendrite toward the white matter or toward the pial surface respectively). Such cell-type specific axonal projection patterns can be generated using G-code. We obtain appropriate code through the following procedure: First, we decompose the cell structure into distinct motifs ( Figure 4B); for each of these motifs we design the mechanisms (G-code machine(s)) which can generate them in a realistic manner; and finally we link the machines in a sequential order (which machine instantiates or kills which other machine). A crucial question is which are the essential features of neural arbors? Currently, the most successful methods for generating simulated branching patterns use recursive algorithms using empirically determined parameters (data driven simulation, e.g. [52]). Even though the end results often look very realistic, the principles at stake in these construction methods are so far from real development that they cannot easily be used to implement and test biologically accurate growth mechanisms. Therefore in this study we do not refine our algorithms using classical quantitative morphological methods such as Sholl analysis [53], fractal dimension [54], or asymmetry index [55]. Instead we aim at characterizing neurons qualitatively. We consider the following criteria: (1) the growth direction: does a neurite extend 'vertically' toward the pia, or toward the white matter, or does it extend more 'horizontally' by staying in a plane parallel to the layer borders? (2) the layer specificity for branch points and ramifications (crucial for the information flow in the cortex); and (3) the distinction between the 'backbone' of the tree (the major shafts of the axon and of the apical dendrite, defining the shape of a neuron at larger scale, having only a few selected branch points, with asymmetric and perpendicular branching modes, and making very sparse synaptic connections) versus terminal 'patches' (ramifications of the neurite, characterized by a symmetrical branching mode and a high density of synapses). These principles are illustrated in Figure 5.
Neurons extend their dendrites in specific directions ( Figure 5AB). Generally, the axon leaves the soma at the inferior pole and goes 'down' toward the white matter; in P23, P5 and P6 an apical dendrite extends from the upper pole of the soma, and moves 'up' to L1; and finally, basal dendrites extend more horizontally. We designed a mechanism for neurite sprouting that depends on two orientation signals. Firstly, we assume that cells have an inherent internal axis. This assumption is based on the observation that the leading process during migration becomes the apical dendrite, and the trailing process the axon [56]. We also assume a general long-range orientation gradient through cortex due to Semaphorin, which in our simulation is secreted by L1 cells, as observed in biology. The newly extended neurites are provided with cell-type specific growth cone (GC) models. In our model, neurite outgrowth and branching are independently regulated [57], with the neurite containing a mechanism for elongation (up or down the gradient of semaphorin for the apical dendrite and the axon respectively; in a random direction for the basal dendrites) as well as regulation mechanisms for producing branches.
We consider essentially two modes of branching: side-branching and tip bifurcation. The former (side branching [58]) refers to the case where a perpendicular minor branch leaves a major shaft. This defines the shape of a neuron on a larger scale, for instance to which layer it projects to, in a cell type and layer-dependent manner ( Figure 5CD). Patchy axonal arborizations( Figure 5EF), or the apical tufts of dendrites, arise by a different growth mode [59] that involves recursive bifurcation of the growth cone. The growth cone can bifurcate at each time step with a small probability; after bifurcation, the growth cone of each daughter branch acts independently, according to the same rule as the parent branch did; during extension and at bifurcation a diameter reduction occurs; when the branch diameter falls under a given threshold, the elongation stops. This mechanism has the useful property that markedly different tree structures are generated by changing only a few parameters (see also Figure 2 in ref. [11]). If diameter reduction during elongation is the major process, then the arborization will be characterized by branch tips having all the same path length, regardless of the number of bifurcation points. This pattern is used to generate basal dendrites, or the apical tuft. If the reduction occurs only at branch points, then the number of branches is balanced independently of the path length (so emulating a possible competition between growth cones). This mechanism is used for generating the axonal patches.
The simulated neural branches ( Figure 6, Video S1) respect the layer specificity described above, especially the ramifications in L2/3 and L5 for P2/3 cells and in L4 for P6 cells ( Figure 7AB). On a smaller scale the neurons are less realistic: Firstly, because we have not yet used statistical methods to obtain branching parameters that match observed geometrical properties of distal branches. Secondly because of border effects due to the limited size of the simulation: When branches leave the simulated column, the concentration gradients they are following point in the opposite direction (toward the column), resulting in a turning of the branches. Interestingly, many geometrical features of the neurons that were not specified in the genetic code, for instance the tortuosity of the branches, emerge because of the physical properties of the developing tissue. To further investigate the role of physics we performed a growth simulation with test cells containing the same G-code instructions in an environment without any mechanical obstacles (i.e. no other cells) and with noiseless distributions of guidance cues, analytically defined as gaussians for the layer specific signaling substances and as a linearly increasing concentration for the Semaphorin produced in layer 1. As expected we observed a clear decrease in branch tortuosity ( Figure 7C): The axons grown in a 'cellular' environment have an average tortuosity of 1.27, whereas the axons grown in the 'empty' environment have an average tortuosity of 1.03. According to ref. [59], tortuosity in biological axons as measured in Layer 1 of mice at P13 was on average 1.2 in Cajal-Retzius cells and 1.48 in thalamo-cortical projections.
The layer (or more precisely the concentration profile of layer specific guidance cues) as well as the plane (perpendicular to the neurite shaft) in which side branching occurs is encoded in our genome. But the exact outgrowth direction is not. It can thus happen by chance that all side branches are on the same side, as it is the case for the P5 cells in Figure 6 (see also Discussion).

Discussion
Much research has been devoted to describing and understanding the developmental process by which relatively few precursor cell types unfold into huge numbers of variously differentiated neurons whose interactions create intelligent behavior. Thus far, experimental exploration of this extraordinary process has focused on crucial attributes such as the rate of cell proliferation, migrational patterns, axonal guidance and branching behavior based on intra-and extra-cellular cues. Although aspects of this process have been studied in simulation (for a review see [60]), most models have focused on single aspects of development that are studied in great detail, but independently of one another. This bottom-up approach [35] allows quantitative local predictions to be made, but do not consider how grown patterns constrain more To gain a better understanding of its structure, we visualize selectively subparts of each cell and describe them qualitatively. (A) Orientation of the neurites of a P23 cell; the axon (black) leaves the soma toward the white matter, the apical dendrite (blue) toward the pia, the basal dendrites (red, pink, yellow, green) in a perpendicular direction. The lengths of all basal dendrites (path lengths from the tip of the branches to the soma) seem to be constant. The apical dendrite ramifies before entering layer 1. (B) Similar neurite orientation for a P6 cell; the basal dendrites make less bifurcation than in the P23 cell, but here too the path length seems to be constant. Here the apical dendrite does not ramify (no apical 'tuft'). (C) Large scale axonal structure of a P23 cell. The main trunk of the axon goes down (initially toward the white matter), and makes side branches in layer 2/3 (green, blue) and in layer 5 (pink, red). These side branches usually extend before forming a patch in the same layer (several side branches have been removed for clarity; the complete axon is displayed in E). (D) Large scale axonal structure of a P6 cell (for clarity several side branches have been removed). (E) One quantitative measure which helps to differentiate between the 'backbone' (forming the large scale shape of a neuron) and the 'patches' is the bouton density, illustrated here with a P23 cell: the segments with more than 0.2 boutons per mm are in red, the segments with fewer boutons are in blue. (F) The differentiation between backbone and patches based on the branching angles is less obvious. Segments displayed in black make an angle w80 degrees to the parent branch (corresponding in principle to side branches); for segments in red the angle with the parent branch is 20-80 degrees (corresponding in principle to bifurcation); the segments in blue are continuations of the mother branch (angle v20 degrees; continuation of a major shaft global circuit formation. At the other end of the spectrum, top-down models focus primarily on large scale properties of neural systems, at the cost of using a coarser granularity that does not allow a detailed representation of individual cell behaviors.
By contrast, our group has over the last few years been developing methods for simulation of detailed multi-scale models, in which development is an interaction between numerous independent agents, each of which may use only local information and acts only in its own local coordinate frame. This approach provides a synthesis of a top-down and bottom-up description, since it allows us to study how phenomena at higher levels of organization emerge from and impact on the different developmental phases at lower spatial and temporal scale. Using the G-code instruction language [11], we have shown in this paper how simple genomiclike codes that provide for the expression of independent cellular functions lead to the elaborate collective spatiotemporal behavior of corticogenesis. In our simulation individual cells express local physical processes such as replication, differentiation, migration, and morphogenesis. These apparently complex actions arise out of the independent behaviors of many intracellularly spatially localized G-machines, which are conditionally instantiated forms of a genomic-like G-code. This means that an individual cell with an elaborate morphology such as a neuron behaves as a collection of spatially distributed machines, constrained by their common cell membrane. For example, somata and neuritic growth cones act independently, only exchanging information by means of extracellular signaling molecules, or with diffusion of intracellular substances/proteins. This collective developmental process takes place in an unprepared environment. This means that individual machines and/or cells cannot simply assume their coordinate frame. Instead, the developmental space must be systematically labelled by the guidance cues which are produced by the successive generations of cells formed during the simulation.
In this paper we have focused on the broader problem of modeling the developmental process at all. It is by no means selfevident how this complex process should unfold from the genome inserted into a single precursor. Now that this fundamental result is possible, there are many refinements of method, and experimentally relevant manipulations to be pursued. The simulation framework CX3D offers a unique opportunity to explore these many questions; in particular, to examine the effect of genetically determined parameters on the emerging physical properties of developing neural tissues [61][62][63][64].
At the level of physical interactions, the model raises questions concerning which physical/mechanical properties and constraints are relevant for a simulated environment. For instance friction and border effects have a strong impact on the quality of the lamination, which suggests that 2-dimensional models of cell migration [65,66] might be oversimplistic. A further example is the tortuosity of neurite outgrowth. This property is usually attributed to intrinsic growth programs [59,67,68], however we find that it is dramatically influenced by mechanical constraints due to other cells. In this respect, our approach contrasts with other simulation frameworks such as L-Neuron [69] or NETMORPH [70], in which the physical properties of the neurons and of the extracellular matrix are taken less into account, and in which the versatility of the growth models are less extensive (for instance these programs do not incorporate cell division). On the other hand, those simulators are faster and less memory intensive, allowing for larger simulations on a single processor (for a comparison of NETMORPH with CX3D in the context of dissociated cultures, see [71]).
In addition to drawing attention to the role of mechanical interactions in neural development [72], our model raises several questions that could be tested experimentally. For example, if there is indeed an overlap in the production of different cell types in the ventricular zone, then additional stopping mechanisms are required during migration, in order to prevent cells of layer n from migrating past earlier born cells of layer n{1. We suggest that stopping could be triggered by contact with cells of the same type that have already stopped.
We further postulate the existence of mechanisms for alternating the direction for side branch formation. Reaction-diffusion based models for selecting the circumferential branching domain have been studied early on by theoreticians [73], and have been experimentally investigated in the context of lung development [74]. In the context of axonal growth, however, alternation of branching direction is not well studied, and mentioned only infrequently in reviews on axonal growth.
Even in areas which are actively investigated, the modeler is often obliged to make assumptions. One such example is the termination of neurite outgrowth, for which there is currently no clear consensus in the literature. Based on morphological observations on dendrites [75], several models based on diameter tapering have been proposed [52]. Other models of stopping condition have been published, for instance based on limited metabolic resource [76], a model especially used in the case of axons [77]. We could also consider other mechanisms, such as a limited time window during which axons are allowed to elongate; or the interplay between an intracellular growth signal and an extracellular retraction cue reaching an equilibrium at the desired axonal length. We decided here to use a condition on the diameter as the only branching mechanism for both axons and dendrites because it is the simplest model to implement, and because it has interesting theoretical implications [68]. This decision does not constitute a prediction based on the confirmation or refutation of which the importance of our work should be judged a posteriori. Instead, the relevance of our work here lies in offering the possibility of testing in simulation mechanisms yet to be discovered.
Obviously, the simulations presented here have many limitations. For example they contain only five different cell types, and the branching patterns are overly simplistic. More detailed implementations will be needed before we can consider that some satisfactory degree of biological plausibility has been reached. However, this verisimilitude is not restricted by the G-code language. The modularity of the G-code allows theoreticians or experimentalists to sequentially refine crucial parts of the model, without compromising the stability in the rest of the simulation. For example, we are currently exploring how the model genome can be inferred directly from cell differential data. Another limitation of our work is the absence of electrical activity dependent growth mechanisms. Although the establishment of synaptic connections is possible in the CX3D framework, the simulation of activity still requires the export into an electrophysiology simulator. Finally, because of the cost of computational resources it is possible to simulate the developmental lamination of the order of 10 4 cells using a single, current desktop computer. Only a few cells are allowed to generate their dendritic and axonal morphologies (our current cortical cells are composed of approximatively 10 3 elements each; more realistic cell architecture would require in the order of 10 5 elements). These simulations take in the order of one hour to model regional corticogenesis (Video S1), that in the mouse plays out over 5 days. Our goal is to describe the development from progenitor cells of the observed neocortical connectivity within a region [21], as well as the long range projections between some cortical areas [78]. Obviously, this is a very large computational task and cannot be performed on a single procesor. For this reason we are currently developing a new version of CX3D that can be distributed across multiple cores of multiple networked computers.

Methods
The G-code instruction language A complete specification of the G-code instruction language can be found in [11]. In this section we summarize the core features of the G-code instruction language needed to follow the implementation of the model presented in this paper. The G-code instruction language is based on a set of nine primitives that define basic neural actions: move, secrete, detect, morph, fork, attach, replicate, synapse, die, defining respectively (1) displacement of a growth cone or of a cell body, (2) substance production (intracellular, extracellular or membrane bound), (3) signaling molecule detection (concentration or gradient), (4) modification of the physical properties (e.g. change in volume or elasticity), (5) formation of a new growth cone (neurite extension from the soma or branching of an already existing neurite), (6) attachment onto the extracellular matrix or onto another cell (e.g. for fasciculation), (7) cell division, (8) synapse formation and (9) apoptosis. Two additional primitives, instantiate and kill are used for regulating the expression and removing networks of primitives (such networks are called 'G-machines', or simply 'machines') in cellular space. Most of the primitives take parameters, specifying the context of the action that has to be taken.
When two primitives are linked in a machine, the output of the first one becomes the input to the second one. For instance a machine for chemoattraction in a growth cone or a migrating soma can be derived by connecting two action primitives in the following way: detect(A,e)Rmove. The parameters (A, e) indicate that we consider the extracellular substance 'A'. The arrow indicates the direction of the signal; in this case the gradient of concentration of the substance A is fed into the move primitive, indicating the movement direction. The G-code language also contains a set of filters that modify the intensity or duration of signals exchanged between primitives. For instance, a filter can be intercalated between the two primitives of the previous machine to invert the direction, leading to chemorepulsion e.g.: detect(A,e)RiRmove, where i denotes the inversion filter. [More precisely: primitives and elements use input and output ports for the formation of links. The correct description of the chemoattraction machine described above is: detect(A,e):gradient?direction:move. The output port gradient of the primitive detect is linked to the input port direction of the primitive move. For reasons of brevity we will omit the mention of input and output ports for the rest of this paper and adopt the simplified representation detect(A,e)Rmove.] Although each G-machine typically executes only a simple behavior, they can be concatenated to generate quite elaborate functionality by calling G-code 'subroutines' for the instantiation of successors or helpers. For example after completing its assigned actions a machine can call for the instantiation of the next machine (for which it has the address on the genetic code), and then remove or deactivate itself. At the beginning of the simulation, the only information that cells have is a list of machines in an inactive form (the 'genome'), as well as the indication of the first machine that is instantiated to launch the developmental process.

Implementation of large machines in the simulation framework CX3D
All simulations were performed with the Java-based open-source framework CX3D [10], running on a DALCO personal computer (CPU: Intel Core 2 Quad Q6600 at 2.40 GHz, RAM: 3.8 GB). The technical implementation of the G-code instruction language in CX3D has been previously described in [11]. In the original G-code setting, each instance of a primitive, a filter or a link corresponds to an independent Java object (sequentially run by the CX3D scheduler). While clean, this approach comes with a prohibitive cost in term of computational time and memory in case of larger simulations (for instance a growth cone model of moderate complexity can contain up to a few tens of Java objects, and this in each terminal neurite segment). For these computational reasons, the G-code machines used in this work were rewritten as single Java classes, with the exact same functionalities as their equivalent multiobject implementation. The XML-based G-code 'genome' framework has been extended with a look-up table specifying which Java class corresponds to the requested G-code machine for instantiation and incorporation into the neurite elements at runtime. Our corticogenesis model relies on 36 main G-code machines (Table 3), often containing internal anonymous sub-machines. In the following, we describe the implementation in 'classical' G-code of the most important machines. As an illustration, Figure 8 offers a graphical representation of many of the machines used for the growth of layer 2/3 pyramidal cells (P23).

Cell lineage
The simulation starts with a single cell in an unprepared environment. The first G-code machine to be instantiated in the initial cell is PreplateFormation which, as its name suggests, is used for the formation of the preplate. This machine contains as embedded machines a cell cycle for cell replication as well as a few other internal machines detecting specific configurations of concentrations of the substances X and Y. We use the cell cycle described in Figure 3A of [11], which works as follows: the size of the soma is asserted with the primitive morph; as long as the diameter is below a given threshold, the cell increases its volume; when the cell exceeds After it has exited the cell cycle, the P23 precursor cell expresses 'RadialMigration', a machine which consists of two sequentially activated internal machines. The first internal machine ('MoveRandom') performs a random walk (with an instance of the G-code primitive 'move' without directional input) until the cell touches a radial glial fiber (RGF); this contact (recognized with the primitive 'detect') triggers the removal (primitive 'kill') of the internal machine, and the activation (primitive 'instantiate') of the second internal machine ('ClimbFiber') for migration along the RGF, until the cell enters layer 1 or touches other P23 cells. When this happens, the migration stops and the differentiation machine is launched. (B) The machine 'DifferentiationP23' is expressed when a P23 cell stops its migration. It extends (primitive 'fork') an axon, an apical dendrite and several basal dendrite and instantiates into them specific growth cone machines. The direction of sprouting is determined by the gradient of the extracellular substance S (Semaphorin). (C) The axon growth cone is composed of three sub-machines: the first one moves the branch tip down the gradient of Semaphorin; the two others produce side branches when crossing a region with high concentration of the P23 or P5 signaling molecule. (D) The growth cone machine used in the terminal branches forms a recurrent branching process, inserting copies of itself in each daughter branch after bifurcation. During elongation the direction is chosen along the gradient of the layer specific signaling cue. The diameter decreases during elongation and at bifurcation, and is constantly assessed (with the primitive 'morph'): as long as the diameter is large enough branching occurs further; when the diameter falls under a threshold, the growth cone removes itself and the elongation stops. doi:10.1371/journal.pcbi.1003173.g008 the threshold, the division occurs (replicate primitive). Since the daughter cells are smaller, they increase their volume in turn and the cycle repeats. The primitive replicate can take as an input a direction specifying the axis of division. The internal axis of each soma is used as the division axis for the asymmetric division in the preplate formation (internal axis?replicate). For symmetric divisions used at later stages the G-code filter p is used to provide a random direction perpendicular to the input internal axis?p?replicate). The concentration of the intracellular substances X and Y is assessed with the primitive detect. X and Y are present in the first cell at the beginning of the simulation. (An alternative would be to have the initial machine produce them until the desired concentration is reached.) Based on the concentration of these intracellular substances, internal machines are removed, and others are instantiated inside PreplateFormation, defining the cell type or the type of divisions. When the concentration of X falls below a threshold T 3 the PreplateFormation machine removes itself, and instantiates the next machine CorticalPlateGRN.
CorticalPlateGRN contains a cell cycle, a gene regulatory network (GRN) and several 'read-out' genes (machines for recognizing specific conditions on gene activity and selectively instantiating further machines). Details on the implementation in G-code and the parameter setting of the GRN and the read-out genes can be found in [11]. In essence, the GRN consists of a set of 10 differential equations describing how the activity of each gene evolves over time as a function of the current expression pattern. The dynamics of the system form four sequentially activated bistable switches (one for each cell-type) deciding between the exit from the cell cycle to become a neuron precursor or further division (genes a i , b i see also Figure 1). In addition two genes s 1 and s 2 are expressed during the transition preplate to cortical plate (not shown in Figure 1); the purpose of these two additional genes is to permit an expansion of the progenitor pool. The detailed dynamics are: Since the genes of the GRN are viewed as intracellular substances (and not as machines), their expression is controlled by the primitive secrete. The read-out genes (not part of the GRN, but activated by genes of the GRN) contain one or more instances of the primitive detect (for probing concentration), and an instance of the primitive instantiate for the activation of new machines linked with some filters or logical elements. The conditional activation of the machines is summarized in Table 2. When the conditions for becoming a specific neuron precursor are met, the cell cycle, the GRN and the read-out genes are removed, and the next machine for the layer specific migration is instantiated.

Migration
The RadialMigration machine used by excitatory cells contains two sequentially instantiated internal machines ( Figure 8A). The first one, (MoveRandom), contains an instance of the move primitive which does not receive any directional input, and so performs the initial random walk of the neural precursors. This machine also contains a mechanism for detecting the contact with a radial glial fiber (RGF) and triggering the expression of the second sub-machine (ClimbFiber), which is then responsible for the adherence to and migration along the RGF. This step is coded with the primitive attach, which takes as its argument the name of the substance expressed at the surface of the RGF. Contact with cells of the same cell type or with layer 1 cells results in the removal of the migration machine, and in the instantiation of the next machine. Once they stop their migration, neural precursors also express a cell type specific membrane marker, and secrete a celltype specific diffusible signaling molecule.

Growing the axonal and dendritic arbors
The differentiation is launched by a cell-type specific machine ( Figure 8B) which extends neurites -containing specific growth cone machines -in specific directions. Sprouting direction is defined according to the internal soma axis or with the help of the semaphorin gradient (a diffusible guidance cue produced by the L1 cells, and spanning the whole cortical plate), and transmitted though the input port 'direction' of the primitive fork; the growth cone machine which has to be instantiated in each new neurite is specified as the argument to this primitive. The growth cones contain (1) a movement mechanism with chemoattraction which moves up (in case of the apical dendrite) or down (for the axon) the gradient of the extracellular substance Semaphorin, or in a random direction (basal dendrite) and (2) machine(s) extending side branches based on specific extracellular conditions. For instance in P23 axons ( Figure 8C), the machine PyramidalAx-onMain extends the axon down the gradient of Semaphorin, P23Side23Outgrowth produces side branches while the axon is in a region where the extracellular concentration of P23 specific diffusible cue is higher than a given threshold, and P23Side5Outgrowth extends a side branch if the concentration of the extracellular substance produced by Layer 5 cells is higher than the concentration of the extracellular substance produced by Layer 4. After the formation of a side branch, the machines responsible for its formation have a certain probability to remove themselves. The side-branches ramify in their respective layers (Figure 8 D): the path finding is as usual done with detect and move, while the bifurcation is controlled with fork and the instantiation of a similar growth cone in both daughter branches. To control this recursive process we use some conditions on the diameter of the branches (assessed with the primitive morph).
The initial P6 machine is very similar to its P23 counterpart: it initiates the sprouting of an apical dendrite, an axon and several basal dendrites based on the extracellular gradient of the substance Semaphorin. The only differences are the types of machine instantiated into the newly produced axon and apical dendrite (the machines for the basal dendrite are identical in all pyramidal cells). The initial axonal growth cone is also very similar to the one in P23: it moves down the Semaphorin gradient, and produces collaterals (but only in layer 6, and not in layer 2/3 or layer 5). The machine instantiated into the axonal side branches in layer 6 are also layer specific. The growth cones in the P6 collaterals move up the Semaphorin gradient toward layer 5. Once in that layer, they make side branches of another type, which move tangentially within layer 5, and finally ramify. We use the same type of branching process as for the P23 side branches. We apply the same principles to the three other cell types: The basket cell of layer 2/3 (B23) has a down-going axon which can bifurcate and produce two types of collaterals. The proximal collaterals are formed close to the soma and surround it with branches. The distal collaterals are much shorter and branch much less often. These different branching patterns are achieved by changing the parameters in our bifurcating growth cone. The layer 4 spiny stellate cell (SS4) axon has first-order side-branches which travel horizontally; these fibers have then second-order collaterals which move up and ramify when they contact layer 2/3. Finally, the main trunk of the axon of the pyramidal cells of layer 5 (P5) also has two types of collaterals (similar to the B23 cells): one type creates a patch next to the soma, whereas the other type moves and ramifies in layer 2/ 3. Each growth step is performed by a specific G-code machine, which only executes a simple task, after which it removes itself, and instantiates the next machine (for which it has the address in the genetic code). Note that we have used two different modes of branching: for the general structure of the cell we use side branching from the main axonal shaft; each branch is made by a different growth cone machine; for ramification we use bifurcation; the growth cone in the two child branches are instances of the same machine type that the was acting in the parent branch. Details on the implementation in G-code of neurite outgrowth and of the branching processes can be found in [11].
Because of the memory and computational limitations on the single threaded version of the simulator CX3D we only let one cell of each type fully differentiate (see also Discussion).

Tortuosity measures
In the following we define an axonal branch as a segment delimited proximally by a soma or a branch point, and distally by a branch point or a terminal point. The tortuosity of a single branch b is defined as t b~lb =d b , where l b is the path length of b, and d b is the euclidian distance between the two ends of b. The overall tortuosity of an axon a is then defined as t a~P b[a,bw20 mm l b =d b ; branches shorter than 20 mm are excluded from the sum because their tortuosity is invariably 1 due to the discretisation scheme used in CX3D. Compared to a mere average (St b T b[a ), our definition of axonal tortuosity has the advantage of not overemphasizing the contribution per linear distance of the smaller branches. Similarly, Portera-Cailliau et al. [59] only considered branches within a certain range of length.

Robustness to parameter variation
Most G-machines operating directly on cell motility or on morphological changes have a linear dependence on parameters. For instance, doubling the rate of diameter decrease per distance in a growth cone will reduce the distance traveled before reaching the cutting threshold by a factor two. And so we can in principle control many aspects of the cell morphology with arbitrary precision. However, in practice we often decide to introduce some randomness by using probabilities (e.g. to kill a growth cone) instead of hard thresholds. This stability is also reflected at the population level: during the pre-plate formation the number of divisions can be precisely defined by choosing the appropriate thresholds (for instance for n symmetric divisions, set T 1~Xstart =2 n ). Obviously, when G-code is used to encode more biologically realistic dynamical systems (such as the set of differential equations of the GRN), the machines inherit the stability conditions of the system, independently from the implementation in G-code.

Supporting Information
Video S1 Movie of the complete simulation. For clarity only two cells differentiate: a P23 cells (axon shown in black, dendrites in grey) and a P6 cell (axon shown in orange, dendrites in green). (MP4)