Creative Destruction: A Basic Computational Model of Cortical Layer Formation

Abstract

One of the most characteristic properties of many vertebrate neural systems is the layered organization of different cell types. This cytoarchitecture exists in the cortex, the retina, the hippocampus, and many other parts of the central nervous system. The developmental mechanisms of neural layer formation have been subject to substantial experimental efforts. Here, we provide a general computational model for cortical layer formation in 3D physical space. We show that this multiscale, agent-based model, comprising two distinct stages of apoptosis, can account for the wide range of neuronal numbers encountered in different cortical areas and species. Our results demonstrate the phenotypic richness of a basic state diagram structure. Importantly, apoptosis allows for changing the thickness of one layer without automatically affecting other layers. Therefore, apoptosis increases the flexibility for evolutionary change in layer architecture. Notably, slightly changed gene regulatory dynamics recapitulate the characteristic properties observed in neurodevelopmental diseases. Overall, we propose a novel computational model using gene-type rules, exhibiting many characteristics of normal and pathological cortical development.

Introduction

Biological tissue often exhibits complex anatomical features that are shared across many species. One well-known example is the layered distribution of different cell types seen in the mammalian neocortex, which has been a central research topic in neuroscience. A number of experimental studies have led to insights into the well-orchestrated processes that form the characteristic layers seen in the adult (Molyneaux et al. 2007; Gaspard 2011). While the pioneer of modern neuroscience Ramon y Cajal and others postulated seven layers in most neocortical regions (Kuhlenbeck 1967), Bevan Lewis (1878) mentioned the existence of six cortical layers, which was later adopted by Brodmann (1909).

Experimental work has demonstrated that the synaptic connectivity and electrical activity dynamics of cortex are strongly shaped by its layer architecture. Along those lines, the presence of connections between two cortical areas in cats and macaque monkeys can be related to the similarity of their laminar elaborations (García-Cabezas and Zikopoulos 2019; Goulas et al. 2019). Moreover, interareal connections obey generic rules dependent upon laminar identity (Markov et al. 2014). From a physiological perspective, the main prevailing model of information processing within a cortical column tightly links its functionality to the laminar architecture (Douglas et al. 1989; Douglas and Martin 2004). Although it is unknown whether layer architecture is necessary for higher order cognition, its omnipresence across different species and neural systems suggests it to be a critical factor for neural computation (Striedter 2005; Adesnik and Naka 2018). A number of experimental studies established cortical layer-specific functionalities (Bortone et al. 2014; Lur et al. 2016; Quiquempoix et al. 2018). Notably, important developments in neural imaging techniques have enabled key insights into the intricate relationship between laminar organization and electrical activity (Larkum et al. 2018).

To this day, the question of how gene regulatory dynamics give rise to the cortical neural layer architecture remains a key matter of investigation in neuroscience. Already, Wilhelm His and Ramon y Cajal studied how proliferating cells in the ventricular zone (VZ) generate neurons migrating radially toward the pia and stop before the marginal zone (MZ) (Cajal 1891; His 1904). A better understanding of how layers form during cortical development is relevant from a biological as well as medical point of view. In particular, many brain pathologies, such as in autism (Crawley 2012), schizophrenia (Nawa et al. 2000), fetal alcohol spectrum disorder (FASD) (Riley et al. 2011), or epilepsy (Bozzi et al. 2012) have been associated with developmental origins of pathological layer formation.

In the last years, in vitro models have become a viable approach to conduct research on neural development (Eiraku and Sasai 2012; Lancaster et al. 2013; Lancaster and Knoblich 2014). However, despite rapid progress along those lines, a number of serious challenges in terms of recapitulating in vivo conditions and functional behavior remain. A complementary, computational approach can benefit this research direction in multiple ways; it enables to combine various data into a coherent framework and allows to formulate and assess novel hypotheses in a quantitative manner. In particular, mechanistic computational models constitute a promising approach because they facilitate the incorporation of various genetic, molecular, and imaging data and so constitute a strong link between experiment and theory.

Only few computational/mathematical studies provide a quantitative framework for cortical layer formation. Cahalane et al. show that a mathematical model of the dynamics of neural fate determination can yield anatomical properties of superficial and deep layers of cortex across different species (Cahalane et al. 2014). However, since this framework is formulated on an abstract, analytical level, it does not take into account how factors in the physical extracellular neighborhood act on cells, such as extracellular substances and their gradients, membrane-bound proteins, or mechanical forces. Hence, it is problematic to test more direct, mechanistic hypotheses in such an approach, for example, on the involvement of certain pathways or extracellular cues. A different computational approach is taken by Caffrey et al., who employ an agent-based model-to-model cortical layer formation (Caffrey et al. 2014). This model does not comprise the developmentally crucial aspect of proliferation and relies on highly simplified intercellular interactions. Along those lines, other computational studies that model axonal and dendritic trees do not focus on the origins of a layered architecture (Koene et al. 2009; Cuntz et al. 2010). By contrast, Zubler et al. (2013) present a mechanistic computational model for cortical lamination with physical interactions. Using the “Cx3D” software (Zubler 2009), the authors studied cortical layer formation from a conceptual perspective of self-construction (Zubler et al. 2013) and put it into a gene-type coding scheme (Zubler et al. 2011; Hauri 2013). Importantly, this approach relies solely on modeling local information exchange; hence, it maintains the biologically realistic condition that cells can only experience their local 3D environment.

While the studies of Zubler et al. (2013) and Hauri (2013) propose working agent-based models for cortical layer formation, they do not demonstrate the generation of appropriate layer-specific cell numbers. Hence, it is currently unclear how a core gene regulatory network (GRN) can produce various different layer architectures. Here, we go beyond a qualitative demonstration of cortical layer formation and provide a computational model that reproduces quantitative measurements of various cortical layer architectures. In particular, we study how neocortical layer architectures observed in health and disease, across different species and cortical regions, arise. To this end, we devise a canonical GRN model, which individually plays out from a homogeneous pool of progenitor cells. In other words, we study how a core GRN, by differential alterations to its physiological interactions (Macneil and Walhout 2011), can generate a wide range of cortical layer architectures. In order to account for the evolutionary conservation of this basic structure, we reproduce experimental measurements from the human, macaque, rat, and mouse cortices. In particular, we investigate how cortical layers comprising different numbers of neurons can arise and evolve in a biologically plausible way. Importantly, this GRN model comprises the key ingredients of cortical development, that is, it incorporates the interplay between neuronal proliferation, differentiation, and apoptosis. Moreover, it respects detailed characteristics of cortical layer formation, such as an initial exponential proliferation phase that is followed by a sequential differentiation phase (Lui et al. 2011).

The overall computational model is tested in a so-called agent-based approach: initially, a homogeneous pool of neuronal precursors is represented, where each cell’s behavior is guided by its individual GRN model as predicted by the protomap hypothesis (Rakic et al. 2009). This model constitutes the genetic rules that specify the cellular dynamics of the individual cells (hence, each cell acts as an autonomous “agent”). Based on this initial configuration, simulations are conducted to demonstrate how, from the proliferating cell population, the genetically specified dynamics play out and yield laminar neuron numbers in accordance with experimental observations.

Materials and Methods

Layer-Specific Neuron Numbers

Experimentally obtained measurements of layer-and species-specific neuron numbers were obtained from (O’Kusky and Colonnier 1982; Schüz and Palm 1989; DeFelipe et al. 2002). Layer thicknesses from human cortex were obtained from von Economo (2009). These data were used to assess the impact of genetic alterations in neurodevelopmental disorders (Fig. 7, Supplementary Figs 3 and 4) and to analyze the correlations between layers (Supplementary Fig. 2).

Images

The histology sections shown in Figure 7 were obtained from Casanova (2007) and Golden and Harding (2010).

Simulation Environment

The simulations were conducted using the agent-based open-source framework “Cx3Dp,” the parallelized version of Cx3D (Zubler 2009; available at http://www.ini.uzh.ch/projects/cx3d/).

Our computational model includes a GRN that is instantiated in a few precursor cells at the beginning of the simulation. In our model, each of the cells follows internal, gene-type rules, while physically interacting with the local, extracellular environment. The GRN implementation is based upon the Cx3D LocalBiologyModule class, where each module defines specific cellular behaviors. The symmetric division phase is governed by a dedicated module, incorporating recursive cell division and growth. Upon termination of symmetric division, this module invokes the differentiation module, which can also probabilistically invoke A1-type apoptosis. After differentiation has occurred, the migration module is initiated, which finally invokes the A2-type apoptosis module. Mechanical forces also act between neighboring cells that are in physical contact. These cell–cell interaction forces are computed based on the diameter of spherical somata, their relative distance, and their adhesive properties. Moreover, cells can secrete chemicals that diffuse extracellularly in 3D space in a gradient fashion and that can be sensed by other cells. For more information on the framework, please refer to the relevant Supplementary Material Section 1 and Zubler (2009). Notably, our GRN implementation is independent from the GRNs described in Hauri (2013), Zubler et al. (2013).

Usage of Experimental Data

In order to compare layer-specific neuron numbers with the computer simulation results, the number of neurons per layer under a cortical surface area of 300 μm × 300 μm was computed. Moreover, since these numbers reflect the neuronal statistics in the adult cortex, the number of neurons dying postnatally (which was not included in our study), estimated to amount to approximately 25% of all neurons (Heck et al. 2008), was added to the adult-stage neuron numbers. Hence, the numbers listed in Table 1 are multiplied by |$0.3\cdot 0.3\!\Big/ \!0.75\Big.$|⁠. Layer thicknesses used to inform the phenomenological model were obtained from von Economo (2009). Time stamps in Figure 3 are based on Saito et al. (2011), Workman et al. (2013).

To estimate the biological time spans for the developmental phases of the model, the start and end times were determined based on the study of Workman et al. (2013). In particular, the progenitor amplification phase was estimated to start at embryonic day 30 in humans (Bystron et al. 2008). Associated times in rat, mouse, and macaque monkey were estimated based on the work of Workman et al. (2013). The selected start and end points for the differentiation phase were determined based on the “Neurogenesis cortical layer VI start” and “Neurogenesis cortical layer II/III end” events, respectively.

Simulation Setting

Initially, 175 precursor cells are arranged as a 2D layer comprising an area of 300 μm × 300 μm, where the cell body locations are given by (normally distributed) random displacements from regular grid points with a jitter of 0.3 μm standard deviation. The simulation takes place within a 3D environment with columnar boundaries of size 300 μm × 300 μm × 1400 μm, the latter being the radial dimension.

Standard Cx3Dp cell body parameters are chosen for these. Each precursor cell is initialized with an intracellular substance quantity of 100 (unitless), which decays by 1% at each time step. The end of the initial, symmetric proliferation phase is invoked when this substance quantity falls below a prespecified threshold T_SYM. 50 000 times steps were simulated to generate the final laminar arrangement. The model produced systems containing approximately 7000, 11 000, 16 000, and 24 000 cell bodies in human, rat, mouse, and macaque cortex simulations, respectively.

Model Parameters

The maximum cell body diameter during proliferation and differentiation was set to 10 μm. If this value is reached during cell body growth, the cell body either divides and subsequently increases its diameter during growth (the diameter growth per time step was 0.2 μm), or terminates its differentiation pathway. However, after a cell has terminated migration and has taken up its final location within the layer architecture, it grows its diameter to 15 μm.

Certain physical parameters change dynamically during the simulation in order to improve the segregation between the individual layers. In particular, MZ cell bodies, that have reached their prespecified external substance concentration and stopped their migratory mode, adopt the Cx3D parameters “interobject-force” coefficient of 1, “adherence” of 0.0001, and a “mass” value of 0.0001. The (standard) interobject-force coefficient of 1.0 renders the MZ cells responsive to mechanical forces from migrating cells. The small adherence and mass values enable them to be easily pushed, for if their mass were large, they would remain at their original position. In other words, these model parameters allow for the migratory cells to push the MZ cells along the radial direction. If not stated otherwise, all other cells had the standard Cx3D values for these three model parameters, that is, 1.

The GRN model of differentiation has 12 parameters. These comprise six probabilities to commit to layer-specific neuron types (C_M, C₂, C₃, C₄, C₅, and C₆). Moreover, the six parameters P_M, P₂, P₃, P₄, P₅, and P₆ determine the probability to commit to apoptosis of type A1. For simplicity, C₂ + P₂ = 1 was used, so only 11 of the 12 parameters were varied to match the experimental data. Parameters were determined in a semi-heuristic process: All parameters were initially set to 0.5. In accordance with the temporal order of layer formation, associated commitment and apoptosis parameters were changed sequentially in order to generate appropriate neuron numbers. The mathematical analysis described in the Results section was used to guide this parameter selection process. The (fixed) parameter T_SYM determines the number of symmetric cell divisions in the early proliferation phase and is set identical for all the layer architectures here due to lack of empirical data.

MZ cells that have terminated migration each create a physical bond with a neighboring cell, with the Cx3D physical bond parameters “breaking-point” set to 500%, the “dumping constant” set to 0, and the “spring constant” set to 1. This physical bonding enables MZ cells to remain connected with one another. By contrast, cell bodies of the other layers (2–6) assume a mass of 1 during migration and increase their mass to 10 after the termination of migration. Moreover, migratory cells’ interobject-force coefficients assume a value of 0.0001. Overall, these physical parameter settings facilitate the right developmental operation: Heavy cells push the lighter MZ cells before taking up their final position. Preliminary simulations indicate that limited variations of these parameters do not significantly alter the resulting cell distributions.

The stopping of cell migration is induced when certain conditions are fulfilled. In the case of layers 2, 3, 4, and 5 cells, two conditions need to be fulfilled: The number of MZ cells in the local neighborhood, that is, in the physical space where physical interactions can take place, must exceed one. Moreover, the number of neighboring cells of the other, earlier developing layers (e.g., in the case of layer 4 cells, these would be layers 5 and 6) must be zero. In the case of layer 6 cells, at least one MZ cell must be detected. Finally, in the case of MZ cells, a prespecified extracellular cue concentration needs to be reached. To this end, an extracellular substance gradient that increases in concentration toward the pia, is instantiated when simulations start. These conditions enabled cells to recognize their appropriate locations and establish the associated layer thicknesses, while pushing the MZ toward the pia.

Once cells have stopped, they initiate their apoptosis modules that trigger cell death if a given condition is fulfilled. Neurons belonging to a specific layer sense the number of local, neighboring neurons (distance <20 μm) belonging to a different layer. If at least three neurons belong to neighboring layers, or one belongs to a more distant layer (e.g., layer 6 for neurons in layer 3), this condition is fulfilled, and (A2-type) apoptosis is triggered. Since the number of layer 1 cells is very low in all species, an additional probabilistic apoptotic rate was included only for layer 1 cells, which does not depend on extracellular conditions. This apoptotic probability was set to P = 0.9 for all species, which yielded agreement with the experimentally measured neuron numbers in layer 1.

The parameters for the GRN are shown in Table 2. They were manually adjusted by taking into account that increased apoptosis at a given stage reduces the number of subsequently developing neurons and increased commitment of a proliferative state to a given neuronal layer type decreases the number of neurons left.

Analysis of Layer Properties

After simulations terminated, the cell positions and information on their cell identity (i.e., which layers they belong to) were exported as a matrix into a Matlab-readable file. Since (occasionally) cells migrated through the MZ and ended up at the topmost radial position, these were discarded in the analysis by thresholding.

The histograms shown in Figures 4 and 7E,F were computed using Matlab. In order to better account for the variation of layer-specific cell numbers, the number of bins were adapted for improved visibility and accuracy. Layer-specific neuron numbers that are compared in Figure 5 were also computed using Matlab. Agreement between experimental and simulated neuron numbers was determined if the mean simulated number was within the 95% confidence interval of the experimentally measured distribution.

Results

The Model

We studied two phases during cortical development, that is, the progenitor amplification phase, when the neural progenitor cell pool increases exponentially due to recursive cell division. In the second phase also differentiation and migration of neurons is modeled. Figure 1A visualizes these phases.

The Progenitor Amplification Phase

The simulation starts with the progenitor cells that are arranged randomly on a 2D sheet, which represents the proliferative zone. Each of these cells follows their own gene regulatory dynamics, specifying recursive cell growth and division (Fig. 1B). The initial state of the GRN is S1. In this first stage of proliferation, cells divide and grow recursively, producing an exponentially increasing number of neural progenitor cells (S2) (Lui et al. 2011). Squares colored yellow indicate cell states with proliferation capacity.

In our model, the proliferating cells possess an internal clock that is realized by means of an intracellular substance. This substance decays over time, and once the progenitor cells sense that it is below a prespecified threshold T0, they terminate this first proliferative stage of symmetric division. These progenitor cells then transition to stage S2, which is the first GRN state that leads to neuronal differentiation. Hence, this intracellular substance acts as a trigger for the neuron specification phase.

In terms of biological time, in human, rat, mouse, and macaque, cortical development the progenitor amplification process starts approximately at 30 days, 9.9 days, 8.4 days, and 34 days postconception, respectively.

The Differentiation Phase

After the exponential proliferation phase (instantiated within the cell state S1), the progenitor cells still retain the potential to proliferate, but they additionally can probabilistically commit to various other cell types (Gao et al. 2014). Once this differentiation process has occurred, neurons migrate radially to take up their final location within the cortex.

Initially, the neuronal progenitor cells (probabilistically) either differentiate into MZ cells, progress further along the differentiation path by dividing, or undergo apoptosis. States S2–S7 constitute proliferative stem cell states, where (for simplicity) only one round of proliferation can occur.

The MZ develops by MZ cells migrating along the gradient of a predefined extracellular cue that increases in concentration along the radial direction. The MZ cells continuously sense this chemical cue, and once a given threshold is reached, they stop migrating. The layer structure is not influenced by the magnitude of this threshold, which solely determines the distance between the proliferative zone and the developing cortex. It mimics the influence of the pial membrane, which in reality stops migration of the MZ cells.

For simplification, we model that layer 6 migrates separately from MZ cells, although in reality layer 6 splits the so-called preplate into MZ and subplate. In our model, layer 6 cells also commit to migration after a brief waiting period (that is specified again based on a decaying intracellular substance acting as a cellular clock). Because of the waiting period, L6 cells migration occurs after MZ cells have reached their target, and L6 cells stop migrating as soon as they sense MZ cells.

The rest of the differentiation process is a repetition of this elemental process of differentiation, proliferation, and apoptosis: Those precursor cells that did not commit to the MZ or L6 fate will continue to divide. After dividing, they again probabilistically commit to the neuronal fate (i.e., layer 5), continue with proliferation, or undergo apoptosis, and so on. This developmental process can be very efficiently encoded in the GRN because it is composed of repetitions of the same elemental genetic rule.

The parameters of the GRN are given by T1, which is a threshold that determines when the symmetric division phase ends. As soon as a cell-internal substance decays below T1, the progenitor amplification phase ends. In other words, T1 determines the number of cells before sequential differentiation is induced. In all simulations, T1 is set to 65. The probabilities C_M, C₆, C₅, C₄, C₃, and C₂ determine the probabilities for progenitor cells to commit to differentiated neuron types; P_M, P₆, P₅, P₄, P₃, and P₂ are the probabilities for progenitor cells to commit to apoptosis. To demonstrate the capabilities of the model, these parameters were adapted to produce the experimentally measured thicknesses.

In terms of biological time, in human, rat, mouse, and macaque cortical development, the differentiation process starts approximately at 48 days, 12 days, 11 days, and 45 days postconception, respectively. It ends at about 191 days, 23 days, 22 days, and 112 days postconception, respectively.

Consistency of GRN Structure

One remarkable observation across different cortical areas in humans as well as in other species is the significant variation of cell numbers in individual layers (Fig. 2; Herculano-Houzel et al. 2008; Collins 2011). It follows that an important condition for a mechanistic generic computational model of layer formation is the possibility to account for various compositions of layer-dependent cell numbers.

We first demonstrate that our model can yield a wide range of neuronal numbers across layers within different cortical areas and species. Figure 3 shows the sequential development of the layers in the human temporal cortex in the correct order, that is, the characteristic inside-out dynamics of layer formation is recapitulated (see also Supplementary Video 1). Moreover, the development of the layers in rat, mouse, and macaque cortices was simulated again by introducing changes to the model parameters. The resulting distributions of neurons across different layers are shown in Figure 4.

In Figure 5, the experimental data obtained from multiple experimental studies are shown together with the corresponding simulated layer dimensions, highlighting the explanatory power of the simple GRN structure shown in Figure 1B. Overall, these results demonstrate that inter-species variation of layer-specific cell numbers can be accounted for only by modifications to the parameters of the GRN rather than structural changes.

Figure 6 shows the cell number dynamics in a simulation of macaque cortical layer formation. First, the MZ that will subsequently become layer 1 develops and then layers 6, 5, 4, 3, and 2 form. The initial production from the stem cell pool is followed by reductions in cell numbers through the apoptotic stage A2 (Fig. 1). We find that most of the apoptosis occurs in mislocated MZ cells, predominantly in the superficial cortical layers.

The modifications of the GRN are solely changes in probabilities for transitions between certain states rather than changes to the GRN states or the connections themselves. This general applicability of the GRN structure indicates that minor evolutionary changes can yield a wide range of phenotypic differences in brain structure. In our model, these modifications can include changes to apoptotic rates in addition to changes to cellular fate commitment. However, these model parameters are underconstrained by current experimental data; hence, there are different possibilities of GRN parameter sets yielding equivalent layer-specific cell numbers. Supplementary Figure 1 shows how changes in exemplary parameters of the GRN affect the resulting layer-specific neuron numbers: Depending on where in the GRN the given model parameter takes effect, the subsequently developing layers are affected.

Differential Apoptotic Processes

In our model, changes to apoptosis at early stages affect the final cell numbers. Apoptosis later during development has less impact on the total tissue thickness because only later-developing layers will be affected. Nevertheless, all the probabilities for cells in the proliferative stages to exit and commit to apoptosis (GRN states labeled A1 in Fig. 1B) exert control over cell numbers. Overall, changes to these probabilities can yield extensive differences in the final cell numbers and layer dimensions and so the A1-type apoptotic stage is paramount in generating the right number of cells.

In a second process, apoptosis is a crucial determinant of the proper layer architecture and neuronal density (GRN states labeled A2 in Fig. 1). Our model of this second type of apoptosis is inspired by Hauri (2013): Cells can sense their local environment, and based on the neighboring cell types, commit to apoptosis. Whenever a cell senses too many neurons of a type that is different to its own, it commits to apoptosis. Hence, in this second phase of apoptosis, a different pathway is activated because it involves the binding of extracellular ligands (Guerin et al. 2002).

In our model, this A2-type apoptosis is initiated after the end of migration when the first layer 2 cells have reached their final radial position and then trigger the apoptotic process. To this end, these apoptosis-initiating layer 2 cells secrete a diffusible substance, which is then detected by all cells that have reached their final destination. Once a sufficient amount of this factor is detected in their neighborhood, these cells sense the number of cells different from their own type. If this number exceeds a prespecified threshold, the cell commits to apoptosis. This type of apoptosis that happens after migration in our model is crucial for refining the layered architecture of the final cortical structure.

These simulation results point toward apoptosis as a highly heterogeneous mechanism. In the first apoptotic stage, the fate is determined probabilistically according to a prespecified chance that is not influenced by external factors. In the second apoptotic stage, extracellular cues conveying information on the identities of neighboring cells are taken into account. In our simulations, we observe that this apoptotic behavior is particularly important for layer 1, which during development is successively pushed by the deeper layers; it occurs more often for this layer that some cells are not pushed effectively and hence remain stuck in a deeper position than the target layer.

Mathematical Analysis of Layer-specific Neuron Numbers

In our model, the number of neurons per layer is determined by overall 11 GRN parameters. These parameters specify the probabilities of either committing to a certain layer neuron type, or continuing proliferation, and then potentially committing apoptosis. These parameters are C_M, C₆, C₅, C₄, C₃, C₂, P_M, P₆, P₅, P₄, and P₃ (P₂ is given by 1 − C₂). Based on these parameters, we derived the number of neurons per layer (the GRN is inherently stochastic and so we estimate the expected number of neurons).

Again, the factor 4 (=2²) stems from the proliferative cell state (state S3 in Fig. 1B). Since only precursor cells that neither commit to layer 1 nor to apoptosis, the factors (1 − C_M) and (1 − P_M) are included. Moreover, the term is multiplied by C₆ as the probability to commit to layer 6.

To show how apoptosis increases the flexibility of the layer-specific neuron numbers, let us take N₃ as an example. N₃ is dependent upon parameters of differentiation and apoptosis, that is, the probabilities for cell fate commitment C_M, C₆, C₅, C_4, and C₃ as well as apoptosis probabilities P_M, P₆, P₅, and P₄. Hence, if there was no apoptosis, changes in any previously developing layers (e.g., layer 5) would necessarily entail changes also in the later-developing layers. In our proposed GRN model, it is possible for N₅ to change without affecting N₄, N₃, and N₂ because apoptosis and cell fate commitment rates can balance out previous changes. Hence, the inclusion of apoptosis enables layers to develop individually, allowing for a high adaptivity of the layer architecture. From an evolutionary perspective, this enables layer architectures to adapt to different requirements with regard to their inputs, outputs, and intra-areal circuitry.

In contrast to the model with apoptosis, here there is a strict dependency between individual layers. For instance, the number of layer 3 neurons decreases if the number of layer 5 neurons increases. There is no evidence for such consistent relationships between individual layers. Along those lines, Supplementary Figure 2 shows an analysis of layer thickness correlations in the human brain (von Economo 2009), without a single instance where layer thicknesses were negatively correlated. Hence, we propose that apoptosis in cortical development enables cortical layers to change individually.

Layer Architecture in Neurodevelopmental Disorders

In addition to layer architectures in healthy cortex, we also investigated whether our model could capture characteristics of certain neurodevelopmental disorders (Fig. 7A–D). Evidently, there is a large variation of the phenotypic expression across pathologies. Nevertheless, specific general properties of malformation and developmental causes are consistent within types of disorders (Barkovich et al. 2005). To investigate whether our model can generate characteristics reminiscent of certain neurodevelopmental disorders, we conducted small changes to the probability parameters of our computational model. To this end, we studied the impact on cortical thickness, which has been experimentally assessed in certain diseases (in contrast to the number of neurons). As the control, we simulated based on our GRN model a layer structure that is in agreement with the layer thicknesses of the human peristriate area 19 (Supplementary Fig. 3A). The cortical thickness of this specific area has been measured to be approximately 2.7 mm (von Economo 2009).

Polymicrogyria

Polymicrogyria is a neurodevelopmental disorder that usually entails a reduced cortical thickness and an increase in folding patterns, while a clearly layered structure is still preserved (Golden and Harding 2010; Judkins et al. 2011). We assessed the phenotypic impact of increased A1-type apoptosis, that is, the scenario where apoptotic activity at early developmental stages (A1-type) was significantly increased. As shown in Figure 7E, these simulations recapitulate a reduced cortical thickness (~1.6 mm) compared with the control thickness (~2.7 mm, see Supplementary Fig. 3A), while laminar organization is still preserved.

In addition to the scenario where apoptosis at early stages was increased, we conducted simulations where the initial phase of proliferation stopped prematurely. This decreased the neural progenitor pool significantly but did not interfere with the migration and apoptosis (Supplementary Fig. 3B). Again, a reduced cortical thickness was observed also in this scenario, while the layer thicknesses were proportionally preserved. Hence, genetic changes reducing the proliferation in the exponential proliferation stage, or increasing A1-type apoptosis, can lead to the phenotypic changes observed in polymicrogyria.

Autism Spectrum Disorder

Also autism spectrum disorders have been associated with pathological apoptosis during brain development (Gabriele et al. 2014). Along those lines, postmortem histological analysis indicates a higher cell density and less clearly structured tissue organization in terms of cell locations (Casanova 2007). Moreover, cortical thickness is increased (Khundrakpam et al. 2017). In order to compare with these observations, we simulated a model where A2-type apoptosis is not functional. As a result, a higher cell density is observed. Moreover, the layers are less clearly separated and more neurons are mislocated (Fig. 7F, Supplementary Fig. 3C). Furthermore, the cortical thickness is slightly increased (~3 vs. 2.7 mm). Overall, changes to both types of apoptosis, that is, either A1-type (in the first stage) or A2-type (in the second stage), can yield architectures reminiscent of certain neurodevelopmental disorders.

Subcortical Band Heterotopia

Finally, we conducted simulations of layer formation in subcortical band heterotopia (also known as double cortex syndrome), where neuronal migration is affected (Guerrini and Dobyns 2014). To this end, in our simulations, differentiated neurons were attributed a malfunctioning migratory behavior with a 20% chance. Once differentiated, these pathological cells did not migrate radially but remained in the proliferative zone. This resulted in the abnormal presence of neurons in the proliferative zone (Supplementary Fig. 4), which is in agreement with the observations of a deep cellular layer composed of heterotopic neurons (Manent et al. 2009).

Discussion

The presence of an abundant variety of layered structures is a common property of cortex and vertebrate brains in general. A mechanistic and in-depth understanding of how cortical layers arise is currently lacking. Computational models are a powerful tool to devise new hypotheses and to generate experimentally testable predictions. In this work, we provide a model for the development of cortical layer architectures and demonstrate that it is applicable to different cortical areas and species. To this end, our results are in accordance with the observation that the number of cortical neurons varies significantly across the cortical surface (Herculano-Houzel et al. 2008; Collins 2011). Moreover, we show that a number of key characteristics of neurodevelopmental disorders can be accounted for. Importantly, our proposed model highlights the role of apoptosis in layer formation and the adaptability of the layer architecture. As such, it provides an explanation for the high proportion of cell death during cortical development (>50%).

Conceptually, our work is different from traditional models of biological development, in that, it follows an agent-based approach. Our model does not rely on any prespecified neural architecture, and it is initialized solely from a homogeneous progenitor cell pool and a single extracellular substance gradient that enables correct cell migration along the radial direction in 3D physical space. Hence, the phenotype is generated in a self-organizing way.

It is well established that apoptosis is abundant during cortical development (Thomaidou et al. 1997; Cavallaro 2015). Already in the early days of neuroscience, the differential roles and mechanisms of apoptosis have been noted (Ernst 1926; Glücksmann 1951). Kuan et al. (2000) discuss the involvement of apoptosis not only in cell number control but also in the proper morphogenesis of the nervous system. Along those lines, Yeo and Gautier review regulatory mechanisms of neural programmed cell death (PCD) during vertebrate neural development (Yeo and Gautier 2004), which is different from neurotrophic death in differentiated neurons.

This model assumes that apoptosis during development comprises differential complexities; in an early apoptotic stage, cell death occurs independently from the local extracellular environment and helps generate roughly the final cell numbers (A1-type). After layer formation has occurred, (A2-type) apoptosis serves a different role by improving the segregation between layers, hence reducing the overlap among layers. Here, the signals from neighboring cells that convey information on the suitability of cell positions within the extracellular context are crucial. In our model, A2-type apoptosis is not necessary for the development of the layer architecture, but it improves the segregation between layers. Notably, two distinct types of apoptosis have been observed experimentally (Rakic and Zecevic 2000), which could correspond to A1-type and A2-type apoptosis. Moreover, it is well established that apoptosis can be induced during development by signaling cascades triggered by neighboring cells (Eroglu and Derry 2016) as is assumed in our proposed model of A2-type apoptosis. It should be possible to experimentally verify the existence of these two hypothetical forms of apoptosis. For instance, A2-type apoptosis could be tested by displacing individual cells during critical developmental stages. Cerebral organoids would offer a suitable testbed for such experiments (Qian et al. 2020).

Importantly, the inclusion of A1-type apoptosis in our model enables a much wider variety of cortical layer architectures. If there was no such apoptosis, changes to the differentiation of early developing layers (e.g., to parameter C₆ in Fig. 1) would necessarily impact all later-developing layers because these early differences influence the number of available progenitor cells and thereby the subsequent differentiation process. Here, we show that apoptosis during the developmental process grants independence to the later-developing layers to evolve differently from the earlier developing ones. Hence, the number of cells in all layers can increase or decrease independently from other layers thicknesses, and so, apoptosis strongly improves this flexibility of the cortical layer architecture. Moreover, our model aligns with the observations of compensatory regulation between apoptosis and proliferation during embryonic cortical development (Freret-Hodara et al. 2017).

We highlight that certain observations naturally follow from our model without tuning of model parameters. Along those lines, experimental work has demonstrated significant rates of apoptosis during proliferation (corresponding to A1-type apoptosis in our model) (Blaschke et al. 1996; Thomaidou et al. 1997). Indeed, also in our simulations, the occurrence of A1-type apoptosis is comparable to the number of generated neurons (Supplementary Fig. 5). Moreover, the rate of embryonic cell death exhibits significant variation across time. This behavior matches with the dynamics of our GRN model (Fig. 1), where apoptotic probabilities of the A1-type apoptosis depend on the GRN state. It is also in accordance with our model that apoptotic rates during development of homologous neuronal populations can vary significantly between species as well as within the same species (Williams and Herrup 1988). Finally, it is established that apoptosis plays a crucial role in determining cortical thickness (Inglis-Broadgate et al. 2005).

Our model suggests that the timing of apoptosis is a crucial factor in the generation of the respective outcome, and so, very different phenotypic properties in neurodevelopmental disorders can be generated depending on the time of apoptotic malfunction. This indicates that certain developmental brain disorders could in fact stem from causes that share many features but occur at different times. Indeed, a number of studies demonstrate the involvement of aberrant apoptosis in developmental brain disorders, such as in autism (Margolis et al. 1994; Wei et al. 2014), FASD (Creeley and Olney 2013), schizophrenia (Glantz et al. 2006), or microcephaly (Poulton et al. 2011). In addition to apoptosis, proliferative dynamics are crucial in generating the appropriate cell numbers and cell types (Lui et al. 2011). Nevertheless, distinct brain developmental disorders will affect proliferation and apoptosis to different extents and can exhibit significant heterogeneity (Donovan and Basson 2017). Hence, the genetic causes can vary depending on the subtype of the developmental disorder, the specific brain region, and the overall genetic context. Notably, specific changes to the proposed GRN model dynamics produce phenotypic outcomes that are in agreement with the observations of healthy and pathological cortical tissue: Changes to apoptosis early during development (A1-type) lead to thinner cortices (polymicrogyria), while introducing a defect to apoptosis later during development (i.e., A2-type apoptosis occurring after migration) gives rise to a more disorganized layer structure associated with autism. Migration deficits yielded pathological gray matter below deep cortical layers, which is in agreement with subcortical band heterotopia.

These findings are relevant from a clinical point of view. One common feature of developmental malformations is a difference in the cortical thickness in various regions. Our model shows the importance of apoptosis in the control of appropriate layer formation and cortical cytoarchitecture. Given this central role of apoptosis, medical intervention to redirect apoptotic behavior at suitable times during brain development is likely a fruitful direction for clinical research. Moreover, for substances that are administered during pregnancy or commonly used in early childhood, particular attention should be attributed to the apoptotic impact (Palanisamy 2012; Sinner et al. 2014). Indeed, reducing apoptosis leads to severe cortical malformation (Kuida et al. 1998). However, increased apoptosis is also a possible origin of disorder (Lainhart and Lange 2011).

Studies of polymicrogyria show increased apoptosis and a thinner cortex (Stottmann et al. 2013; Huang et al. 2016). We find that specific defects in gene regulatory dynamics controlling apoptosis yield characteristic properties associated with this disorder. Intriguingly, the gene GPR56 has been associated with polymicrogyria (Bahi-Buisson et al. 2010) as well as with apoptotic pathways (Ke et al. 2007; Yang and Xu 2012). Hence, our results are in agreement with the suggestion that polymicrogyria is not a disorder of neuronal migration (Judkins et al. 2011) and highlight the possibility that impaired apoptosis is a key player in the origins of certain neurodevelopmental disorders. However, in our simulations, we introduced a single change to apoptosis occurring before cell migration is induced. Therefore, in contrast to Judkins et al. (2011), we do not suggest that polymicrogyria should be considered a postmigration malformation of cortical development. Instead, we propose that pathological apoptosis before and during migration could be the key driving cause.

Multiple studies demonstrate a tight link between local cytoarchitecture and global brain connectivity as defined by interregional, long-range fiber tracts (Hilgetag and Amunts 2016; Wei et al. 2018). In particular, regions of similar cytoarchitecture and laminar profile tend to be more likely connected. While there are multiple possible explanations for this relationship, it suggests that alterations of the layered cytoarchitecture due to neurodevelopmental disorders will have impact on local as well as long-range connectivity. We anticipate that computational models bridging the local (Bauer et al. 2014) and global (Bauer and Kaiser 2017) spatial scales will help better understand the intricate relationship between cortical layer formation and connectivity.

Given that key neurodevelopmental genes are often involved in many functions (Manzini and Walsh 2011), a specific gene variant might operate correctly at certain developmental stages but may be wrongly expressed at others. Future treatments will need to take into account the specific timing of the genetic expression during development. Along those lines, we envisage that computational models of developmental processes will constitute an additional pillar for clinical research, allowing to anticipate fruitful experiments and formulate novel hypotheses.

Our proposed computational model gives rise to a number of experimentally verifiable and quantifiable predictions. Although the values of the model parameters are underconstrained and do not necessarily match with the true apoptotic rates, the model yields a number of experimentally verifiable predictions. In particular, apoptotic rates should be a major determinant of differences in layer-specific neuron numbers across species as well as within the same species. This is in agreement with previous studies on human and hamster cortex, suggesting that apoptosis could have a key role in shaping regional differences (Finlay and Slattery 1983; Rakic and Zecevic 2000). Along those lines, our model suggests that apoptotic rates change significantly depending on the layer that is being formed. Furthermore, our work suggests that it should be possible, by changing apoptotic events during brain development, to differentially impact cortical layers. In other words, the number of cells in individual layers should be controllable by transient changes to apoptotic rates during proliferation. Such relationship could be assessed experimentally by analyzing markers of apoptosis at different developmental stages during cortical layer formation in combination with layer-specific markers (Hevner 2007; Brunjes and Osterberg 2015). Finally, the computational model suggests that apoptosis, after cell migration has terminated, is relatively large for MZ cells. This is partially due the likelihood of not being pushed by later-developing layers, which is larger than for other layers (see Fig. 6). Additionally, in order to reproduce the small layer 1 neuron numbers, our model requires specifically high MZ cell apoptosis rates in the superficial layers. Indeed, experimental analyses are in agreement with frequent MZ cell apoptosis (Meyer and González-Gómez 2017).

From an evolutionary perspective, our work contributes to a better understanding of a longstanding conundrum on the question why there is so much apoptosis during brain development, that is, approximately one-half of all neurons (Cavallaro 2015). One possibility is that it enables cell selection, that is, inappropriate cells can be removed for the sake of a more refined cortical structure (Pompeiano et al. 2000). However, from an energetic point of view, such massive apoptosis would be surprisingly inefficient. We propose that one key function of apoptosis during brain development is to render the neural architecture more evolvable, enabling the brain to adapt faster and more efficiently to different requirements. Along those lines, it has been previously suggested that apoptosis is responsible for the rapid evolutionary reduction in neuron populations in Spanish wildcats and their domestic relatives (Williams et al. 1993). Moreover, it is in accordance with this evolutionary perspective that apoptosis seems to operate in all multicellular animals (Ameisen 2002). Finally, apoptosis is known to enable greater flexibility as certain evolutionary old structures, such as, for example, pronephric kidney tubules, which form functioning kidneys in fish and in amphibian larvae, are not active in mammals and degenerate (Meier et al. 2000). As mentioned before, it is the appropriate balance between apoptosis and proliferation that is crucial to generate healthy neural structures. Our model suggests that human brain evolution was enabled by a suitable composition of proliferation that is in line with area-specific apoptosis, which determines the properties of the ontogenetic column (Geschwind and Rakic 2013).

Any computational model is based on approximations and so does not reflect many aspects of the actual biological system. We acknowledge that neuronal survival during cortical development is impacted by electrical activity (Heck et al. 2008; Smith and Walsh 2020). Notably, postnatal electrical activity controls neuronal apoptosis in mouse cortex in an area-specific and layer-specific manner (Blanquie et al. 2017). However, it is beyond the scope of this work to cover the full developmental extent of apoptosis in cortex. Cortical layers 1–6 develop before the ingrowth of thalamic afferents, development of corticofugal projections to their targets, and formation of intracortical networks (Huttenlocher and Dabholkar 1997; Kostovic 2002).

It is known that electrical activity is also involved in cortical development during early stages (Spitzer 2006; Yu et al. 2012). Recent work indicates that electrical cellular activation can contribute to migration pausing, eventual arrest, and subsequent dendritogenesis (Medvedeva and Pierani 2020). Moreover, axonal innervation regulates the survival and proliferation of neural progenitor cells in the early developing cortex (Popovik and Haynes 2000). However, there is no evidence that these processes specify area- and layer-specific neuron numbers. We assume that electrical activity during early cortical development, that is, before axons, dendrites, and chemical synapses form, has a permissive rather than instructive role in the formation of area-specific layer architecture. In other words, our working hypothesis is that while electrical activity is necessary for normal proliferation and migration during cortical layer formation, it does not specify the layer architecture. In the future, taking into account additional developmental mechanisms will contribute toward an integrated study of the physiological implications of normal and pathological brain development. Due to their important contribution to the cortical volume, modeling glial cells, neurites, and synapses will also improve accordance with the measurements of cortical neuronal density and layer thickness.

Our model does not include radial glial fibers (RGFs), that is, the processes spanning along the full width of the developing central nervous system. RGFs play important roles during cortical development. They provide a scaffold for the migration of projection neurons that transforms into a physically discontinuous structure during layer formation (Rakic 1972; Nowakowski et al. 2016). In our model, the RGFs’ guidance function is taken over by an initial extracellular gradient that acts as a guidance cue. Moreover, the transient presence of RG subtypes is simplified by the (transient) stem cell types in our proposed core GRN (Fig. 1B). As in the case of electrical activity, we assume that RGFs are necessary but not instructive for the specification of the area-specific layer architecture. Along those lines, we did not differentiate between the excitatory and inhibitory neurons in our model. Inhibitory neurons make up between 20% and 30% of the neuronal population, with significant differences across species (Sultan and Shi 2018). However, there is no evidence that differential distribution of inhibitory neurons is the main cause for area or layer differences in neuron numbers. This is further supported by the fact that layer-specific cell numbers can differ by far more than 30% across cortical areas (Charvet et al. 2015). We did not include the development of subplate cells and preplate splitting (Marin-Padilla 1978). The role of the subplate in cortical development is complex and not yet fully understood (Molnár et al. 2019). Hence, its inclusion in the model would render the simulations more time consuming, without adding significant explanatory power relevant to the questions addressed here. A proof-of-principle for the capability to include preplate splitting in our model is shown in Supplementary Figure 6.

We are taking steps toward a comprehensive computational model of cortical development, incorporating neuronal self-organization of biologically plausible axonal arborizations (Bauer et al. 2012; Kassraian-Fard et al. 2020), complex network properties of synaptic connectivity (Bauer and Kaiser 2017), canonical microcircuit function (Bauer et al. 2014), and thalamocortical projections (Kassraian-Fard et al. 2020). We also develop modern coding techniques and software tools (Bauer et al. 2017; Gonzalez de Aledo et al. 2018; Breitwieser et al. 2020), which will enable simulations that give rise to a detailed model of neural tissue dynamics including folding patterns (Wang et al. 2016). Given that our model can account for a wide range of layered neocortical cytoarchitecture, only small amendments are likely to be necessary to explain also other layered structures in the central nervous system, such as the retina (Mellough et al. 2019). In the future, our model could be extended by integrating data on gene expression, such as obtained from public databases (e.g., Polioudakis et al. 2019; Telley et al. 2019). Such detail will enable the modeling of the relationship between genetic specification and the associated phenotypic properties of cortex. Due to the large quantities and complexities, modern machine learning techniques could be employed to produce a tight match with the experimental data. Moreover, we acknowledge the importance of reproducibility and replicability and the benefits of adhering to certain guidelines and approaches (Manninen et al. 2018; Crook et al. 2020). Hence, we provide the associated Cx3Dp code to facilitate the reproduction of our findings.

Overall, we provide a mechanistic, agent-based, and generic computational model for cortical layer formation that takes into account the interaction of intracellular dynamics with extracellular communication. In particular, our model implicates the role of apoptosis as a highly heterogeneous control mechanism, enabling the adaptability and variety of cortical layer structure. Moreover, it constitutes a platform for the generation and testing of hypotheses, such as, for example, on the malfunctioning mechanisms that give rise to developmental disorders.

Funding

Engineering and Physical Sciences Research Council (EPSRC) of the United Kingdom (EP/K026992/1) as part of the Human Brain Development Project (http://www.greenbrainproject.org/); Medical Research Council (MRC) of the UK (MR/N015037/1 to R.B.); the EPSRC (EP/S001433/1 to R.B.); Wellcome Trust (102037 to M.K.), EPSRC (NS/A000026/1, EP/N031962/1); MRC (MR/T004347/1); and the Guangci Professorship Program of Ruijin Hospital (Shanghai Jiao Tong University).

Notes

The authors wish to thank Rodney Douglas, Rob Forsyth, Evelyne Sernagor, Peter Taylor, and Yujiang Wang for stimulating discussions and useful comments. This work made use of the facilities of N8 HPC Centre of Excellence, provided and funded by the N8 consortium and EPSRC (grant no. EP/K000225/1). The Centre is co-ordinated by the Universities of Leeds and Manchester. Moreover, parts of the computational work were performed using HPC infrastructure at the Newcastle University and the University of Surrey. Conflict of Interest: None declared.

Data Availability

The simulation code can be found on ModelDB (http://modeldb.yale.edu/266895) and Zenodo (https://zenodo.org/record/4456046#.YAoAMXb7S90).

References