Abstract

Injections of neural tracers into many mammalian neocortical areas reveal a common patchy motif of clustered axonal projections. We studied in simulation a mathematical model for neuronal development in order to investigate how this patchy connectivity could arise in layer II/III of the neocortex. In our model, individual neurons of this layer expressed the activator–inhibitor components of a Gierer–Meinhardt reaction–diffusion system. The resultant steady-state reaction–diffusion pattern across the neuronal population was approximately hexagonal. Growth cones at the tips of extending axons used the various morphogens secreted by intrapatch neurons as guidance cues to direct their growth and invoke axonal arborization, so yielding a patchy distribution of arborization across the entire layer II/III. We found that adjustment of a single parameter yields the intriguing linear relationship between average patch diameter and interpatch spacing that has been observed experimentally over many cortical areas and species. We conclude that a simple Gierer–Meinhardt system expressed by the neurons of the developing neocortex is sufficient to explain the patterns of clustered connectivity observed experimentally.

Introduction

Focal injections of neural tracers into the neocortex result in quasi-periodic patchy labeling of somata and axonal terminations that may extend over several square millimeters within a cortical area (Rockland and Lund 1982, 1983; Rockland et al. 1982; Gilbert and Wiesel 1989; Lund et al. 2003; Fig. 1A). These patches are most prominent in the superficial layers where pyramidal neurons give rise to laterally extending intralaminar axons that support the clusters of boutons (Kisvarday and Eysel 1992; Buzás et al. 2006; Binzegger et al. 2007), suggesting that the laminar patches are composed of neuronal clusters (Muir and Douglas 2011).

Figure 1.

(A) Patchy labeling of somata and terminals in layer II/III in the macaque visual cortex following bulk injection of the tracer cholear toxin B subunit (from Lund et al. 2003, with permission). (B) Scaling of interpatch distance with patch diameter observed across various cortical areas and species (from Muir et al. 2011, with permission).

Figure 1.

(A) Patchy labeling of somata and terminals in layer II/III in the macaque visual cortex following bulk injection of the tracer cholear toxin B subunit (from Lund et al. 2003, with permission). (B) Scaling of interpatch distance with patch diameter observed across various cortical areas and species (from Muir et al. 2011, with permission).

The “patch system” (Rockland and Lund 1982; Rockland et al. 1982) or “daisy architecture” (Douglas and Martin 2004) is observed across many cortical areas and species. The remarkably regular scaling of patch diameter to interpatch distance across these many areas and species (Fig. 1B) suggests that patches are a fundamental motif of cortical organization and function (Douglas and Martin 2004; Muir et al. 2011). Indeed, there is good agreement between the spatial patterns of the functional domains observed by optical imaging of the intrinsic signal associated with cortical neuronal activity and the spatial patterns of the anatomical patches (Muir et al. 2011). The interneighbor distances and angles indicate that the patches form an hexagonal lattice that is relatively periodic and isotropic in the visual cortex, but may be less so in other areas (Muir and Douglas 2011; Muir et al. 2011).

Various explanations have been offered for how this patchy organization could arise. For example, Mitchison and Crick (1982) and Buzás et al. (2006) proposed mechanisms that depend on the developing functional relationship between neurons in the visual cortex that are already able to respond to visual stimuli. However, these models are incomplete because patches are found also in areas other than in visual cortex and can be observed in a coarse form before the afferents carrying structured electrical signals arrive in the superficial layers of the cortex (Price 1986; Callaway and Katz 1990; Ruthazer and Stryker 1996). These observations imply that the patchy organization is broadly prespecified, at least on a coarse scale, and that they are refined later by the suitable patterns of afferent electrical activity (Luhmann et al. 1986; Wong 1999; Grossberg and Seitz 2003; Liets et al. 2003; Kanold 2004).

In this paper, we explore the possibility of genetic specification of the patch system. Our hypothesis is that a precursor of the patchy organization develops as early as in the cortical preplate and that the pattern is preserved during corticogenesis and lamination. Neurons preserve clonal features of their precursors in a columnar fashion due to their migration along radial glial cells in a manner consistent with Rakic's protomap hypothesis (Rakic 1988).

In our model of patch formation, neuronal precursors are genetically disposed to secrete a set of morphogens that are able to diffuse through the extracellular matrix. These morphogens are transcription factors whose interactions with the “genome” follow Gierer–Meinhardt reaction–diffusion dynamics (Turing 1952; Gierer and Meinhardt 1972). Consequently, the precursors of the preplate come to express a two-dimensional periodic profile of morphogens that provides the basis for clusters of neurons expressing similar profiles. This periodic identity is inherited from the precursors by their daughter neurons, which migrate radially to form the superficial cortical layers. When these migrating neurons come to rest, they extend lateral axons whose growth cones seek distant target neurons with similar morphogen expression profiles to their own, so generating the observed patchy organization.

We have explored this hypothesis using 2 different simulation approaches. For most of the work we used Cx3D (Zubler and Douglas 2009), in which the detailed physical mechanisms of neuronal growth and cortical development can be simulated. We used this software to demonstrate how a reaction–diffusion system in a layer of progenitor cells leads finally to the formation of patterned axonal lateral connections compatible with the observations.

Our investigations show that the superficial patch system could be specified in the very early stages of cortical development. In particular, we show how a coarse form of this connectivity pattern can be generated without any instructive electrophysiological activity. Our simulation is in agreement with and offers an explanation for several experimental findings.

Materials and Methods

Simulations were performed using the open-source Java Package Cx3D (Zubler and Douglas 2009; Zubler et al. 2011), available from http://www.ini.uzh.ch/projects/cx3d/, and with MATLAB (Mathworks Inc.). The computer used to run simulations had two 12-core processors (1.9 GHz, 64 GB of RAM).

Reaction–Diffusion Parameters

We used well established mathematical tools such as linear stability analysis to provide estimates of the regimes in the parameter space that lead to the desired pattern class. The parameters used for the various simulations of the pattern formation are listed in Table 1.

Table 1

Reaction–diffusion parameters

Simulation in Framework Parameters 
Figure 2 Cx3D Da = 8, Dh = 8000, μa = 0.08, μh = 0.16, ρa = 50, ρh = 50, ρ = 20 
Figure 3 Cx3D Da = 8, Dh = 9000, μa = 0.08, μh = 0.16, ρa = 50, ρh = 50, ρ = 600 in A, ρ = 35 in B and ρ varying between 20 and 600 in C 
Figure 4 Cx3D Same parameters and patterns as used for Figure 3 
Figure 6B,C Cx3D Da = 8, Dh = 9000, μa = 0.08, μh = 0.16, ρa = 50, ρh = 50, ρ = 850, repulsion constant r among activators = 550 
Figure 7A,B MATLAB Da = 0.6, Dh = 230, μa = 0.25, μh = 1, ρa = 0.2, ρh = 0.2, ρ = 0.9 
Figure 7C,D MATLAB Da = 0.6, Dh = 230, μa = 0.25, μh = 1, ρa = 0.2, ρh = 0.2, ρ = 0.9, repulsion constant r = 0.5 in C and 2.7 in D, respectively. Repulsion constant from inhibitor morphogen = 0 in C and 0.0075 in D, respectively 
Simulation in Framework Parameters 
Figure 2 Cx3D Da = 8, Dh = 8000, μa = 0.08, μh = 0.16, ρa = 50, ρh = 50, ρ = 20 
Figure 3 Cx3D Da = 8, Dh = 9000, μa = 0.08, μh = 0.16, ρa = 50, ρh = 50, ρ = 600 in A, ρ = 35 in B and ρ varying between 20 and 600 in C 
Figure 4 Cx3D Same parameters and patterns as used for Figure 3 
Figure 6B,C Cx3D Da = 8, Dh = 9000, μa = 0.08, μh = 0.16, ρa = 50, ρh = 50, ρ = 850, repulsion constant r among activators = 550 
Figure 7A,B MATLAB Da = 0.6, Dh = 230, μa = 0.25, μh = 1, ρa = 0.2, ρh = 0.2, ρ = 0.9 
Figure 7C,D MATLAB Da = 0.6, Dh = 230, μa = 0.25, μh = 1, ρa = 0.2, ρh = 0.2, ρ = 0.9, repulsion constant r = 0.5 in C and 2.7 in D, respectively. Repulsion constant from inhibitor morphogen = 0 in C and 0.0075 in D, respectively 

Reaction–Diffusion Mechanism

We used the Gierer–Meinhardt reaction–diffusion system (Gierer and Meinhardt 1972) for simulating the pattern formation mechanism. The system is described by the following partial differential equations:  

(1)
$$\eqalign{& \displaystyle{{\partial a} \over {\partial t}} = \rho \displaystyle{{a^2 } \over h} - \mu _a a + \rho _a + D_a \displaystyle{{\partial ^2 a} \over {\partial x^2 }} \cr & \displaystyle{{\partial h} \over {\partial t}} = \rho a^2 - \mu _h h + \rho _h + D_h \displaystyle{{\partial ^2 h} \over {\partial x^2 }}} $$
where a(x,t) and h(x,t) indicate the activator and inhibitor substance concentrations at location $$x$$ and time t, $$D_a $$ and $$D_h $$ are the morphogen diffusion coefficients, and ρ is the production coefficient. $$\mu _a $$ and $$\mu _h $$ are the substance decay coefficients, and $$\rho _a $$ and $$\rho _h $$ the basal production coefficients. Appropriate parameters allow this system to converge toward a spatially patterned equilibrium of morphogen concentrations.

The equations for the extended secretion model in the case of multiple Gierer–Meinhardt systems can be written as follows:  

(2)
$$\eqalign{& \displaystyle{{\partial a_i } \over {\partial t}} = \rho \displaystyle{{a_i^2 } \over {h_i }} - \mu _a a_i + \rho _a + D_a \displaystyle{{\partial ^2 a_i } \over {\partial x^2 }} - \sum\limits_{\mathop j\limits_{\,j \ne 1} } {ra_j } \cr & \displaystyle{{\partial h_i } \over {\partial t}} = \rho a_i^2 - \mu _h h_i + \rho _h + D_h \displaystyle{{\partial ^2 h_i } \over {\partial x^2 }}} $$
where the repulsion term $$\sum\limits_{{\scriptstyle j \atop \scriptstyle j \ne i} } {ra_j } $$ is subtracted from the activator term. The indices i,j represent different activator types.

Multiagent Simulation

Cx3D provides a sophisticated framework for simulating neural development in 3-dimensional (3D) space (Zubler and Douglas 2009). Neurons are composed of discrete physical components such as spheres (somata) and chains of cylinders (neurites), each located at a particular point in 3D space. Extracellular space is also discretized into small nonoverlapping domains, allowing for the diffusion of signaling molecules. Each cell element is the center of such a space domain. Additional space domains are added to increase the precision of diffusion. To provide the desired cellular functionality, Cx3D users can write biological modules, Java classes that are inserted into the cells or cell components of the simulation. These modules lend the cells their functional attributes.

The main difference to the classical approach of reaction–diffusion on a regular grid is that in Cx3D we decouple the mechanisms occurring in the cells (production) from the mechanisms occurring in the extracellular space (diffusion and degradation). For instance, for the implementation of the system given by equations (1), the terms $$\rho a^2 /h$$, ρa, ρa2, and ρh are coded in a “biological module” inside the cells, releasing the amounts of substances a and h according to the extracellular concentrations of these substances felt locally by each cell. Diffusion and degradation of these substances are automatically performed by the physics engine of Cx3D and occur even in the absence of cells.

Cellular somata are modeled as spheres with a diameter of 20 μm, placed in an irregular constellation, to form a sheet of 100 × 100 cells. Physical interactions ensure that they redistribute themselves after initialization, such that the average intersomatic distance converges to roughly 22 μm. We add 25 000 additional space nodes for a more accurate simulation of diffusion. For a detailed description of the numerical procedure in Cx3D, we refer to the original description in Zubler and Douglas (2009).

Regular Grid

We performed some simulations of pattern formation on regular grids, which provide an advantage in numerical simulation performance over the irregular environment of Cx3D. Consequently, some investigations requiring numerous trials or larger-scale simulation were based on MATLAB simulations on regular grids.

The differential equations were solved numerically on a discretized grid (in 2-dimensional (2D) in contrast to the 3D simulations in Cx3D) of 150 × 150 cells with the explicit Euler method for calculating the time evolution of the substance concentrations. Unlike in Cx3D, we did not make a distinction between concentration and quantity, because in this case, all domains have the same dimension. At each time step, we solved first the diffusion, that is, we updated the value of each bin ij based on the Fickian flux with its 4 neighbors (given by the concentration difference multiplied by the diffusion coefficient multiplied by the time step, dt = 0.001). Then we updated the value to account for the reaction component, by computing the change in the substance quantities based on the current activator and inhibitor values, and adding it to each bin ij (explicit Euler method).

Simulation was initialized with concentration values drawn from a uniform distribution between 0.3 and 0.7. An example script of a Gierer–Meinhardt activator–inhibitor system with parameters leading to patchy patterns is included in the Supplementary Material (Supplementary Scripts 1).

Pattern Characterization

Quantitative analysis of the outcome of reaction–diffusion simulations was performed in MATLAB. We defined a “patch” as an area of increased concentration of a morphogenetic substance, identified by an automated method based on the fitting of Gaussians, as described in Muir et al. (2011). Subsequently, a graph connecting neighboring patch locations was constructed using a Gabriel graphing technique (Gabriel and Sokal 1969; Tanigawa et al. 2005; Muir et al. 2011). This graph provided estimates of the patch locations, number of patches, mean distances between patches and angles between neighboring patches in a given pattern. These measures were used to compare simulation results against the experimental data. The estimates depend on the expected size of the patches, which is given as an input to the algorithm. This input was determined using the method described in Shen and Jung (2005). This method yields estimates of the average size and number of patches by first segmenting the concentration patterns using binary thresholding and subsequently calculating standard image processing measures. Our implementation is included in Supplementary Material (Supplementary Scripts 2).

We developed an algorithm that uses the methods described above to classify the patterns of substance concentrations as patchy or nonpatchy. We first assessed the positions of the patches of the activator and inhibitor substance patterns. In the ideal case, the patches of the activator morphogen should match in position with the inhibitor patches (see Results). If the estimated locations of the corresponding patches in those 2 patterns deviated too much from one another, that is, if more than 25% of the patches are separated by more than 6 grid cells, then the pattern was classified as nonpatchy. The parameters of the algorithm were determined heuristically to optimize its performance (Supplementary Fig. 2 for 3 example patterns that were correctly classified).

Axonal Growth

Axonal growth was simulated as a growth cone able to sense substance concentrations and gradients (in Cx3D). When an axon bifurcates, its daughter branches receive a copy of the parent growth cone (and its rules) resulting in a recursive ramification of the terminal axon.

We used a growth rule that uses the activator substance in the Gierer–Meinhardt system as a guidance cue. Growth direction was defined at each time step as the weighted vectorial sum of a random direction, the gradient of the activator substance and some measure of the previous movement direction: d= wnoisednoise + wgraddgrad + wprevdprev. In all simulations, we used wnoise = 0.34, wgrad = 0.03, and wprev = 0.63. These weights have been chosen to fit the curvature of long-range layer II/III projections as observed in the cat visual cortex (Ruesch 2011). Outgrowth speed was constant (35 μm/h, independent of activator concentration). The probability of an axon to bifurcate was computed at each time step as $$p_{{\rm bifurcate}} = u\cdot c_a + v$$, where $$c_a $$ is the concentration of the activator. The constant u determines how strongly bifurcation probability is influenced by the extracellular substance concentration, and v is added to allow bifurcation to occur also in the space between patches. Axon diameter reduces during elongation and at branch points. When a specified minimal diameter is reached, growth ceases. For the growth rule of a single-patch system (see Results for a description of the 2 growth rules), u and v were chosen to be 2.6 × 10e−4 and 2.1 × 10e−5, respectively. In the multiple patch systems, these values were 5 × 10e−5 and 2.1 × 10e−4. Diameter reduction parameters due to growth and bifurcation were −0.1% and −10% (at each step) in the first scenario, and −0.05% and −10% in the second. In both cases, the initial diameter was 3 μm and the minimal, diameter criteria to stop outgrowth was 0.5 μm.

Simulations with these parameters resulted in axonal arborizations with distinct patchy projections that matched our experimental observations. For example, the number of clusters formed by a single axon lay in the range of 1–7, and cluster diameters were between 90 and 950 μm (Binzegger et al. 2010).

Results

The Gierer–Meinhardt Model

We used the reaction–diffusion model proposed by Gierer and Meinhardt (1972) (equation 1), because it offers the necessary stability to generate patchy patterns, and because its dynamics can be interpreted as the interaction of 2 morphogens. Their system consists of a local autocatalytic activator a and a rapidly diffusing inhibitor h; a promotes its own production, as well as the production of the inhibitor h, whereas h suppresses the production of a. These 2 antagonistic substances can lead to regions (patches) of high activator and inhibitor concentration (because the activator drives production of the inhibitor), separated by regions of low morphogen concentration. With a grasp of the various regimes of this system (see Materials and Methods), we could demonstrate the required patchy patterns in continuous space.

Single Global Pattern in a Sheet of Cells

Patchy patterns were computed using standard numerical methods on a regular grid. However, our question of whether sheets of individual cells can generate a similar patterning of their behavior brings complications that need not be confronted in the more simple numerical/continuous case (Bonabeau 1997). First, the generators and consumers of morphogens are now individual cells, whose volume and separation effectively quantize the space in which the reaction–diffusion process will play out. Of course, morphogens can diffuse between cells. This diffusion through extracellular space is quantized more finely than the discretization imposed by the cellular positions alone. Secondly, cells are not restricted to a regular grid. They interact mechanically with their neighbors, and so cell locations within a layer are perturbable. And finally, these cells must maintain the patchy morphogen distributions during the developmental process, when cells are migrating radially from their birthplace in the subventricular zone to become superficial layer neurons. It is this more complex problem that we have explored by simulations in Cx3D, rather than in MATLAB.

We simulated the development of a pattern of patches on a 2D sheet of cells that secrete morphogens according to the Gierer–Meinhardt kinetics. Simulations begin with a random distribution of morphogens across random cell body locations. Figure 2 shows a sequence of images derived from a Cx3D simulation in which the process of pattern formation evolves over time. In the initial stages, symmetry-breaking drives the initial random pattern (A) toward blurred patches consisting of high levels of both the activator and inhibitor morphogens. These patches become more distinct over time and decrease in number due to the growing lateral influence of the inhibitor substance (BE). For appropriately chosen parameters, the final pattern converges to a steady state (F) that resembles strongly the observed pattern of cortical patches. The parameters of the pattern formation kinetics (equations 1) are listed in Table 1.

Figure 2.

Sequence of images from a Cx3D simulation showing the development of a Gierer–Meinhardt pattern over time. The neuronal cell bodies (somata) are colored spheres. Black indicates low morphogen concentration near the cell body. Red and green indicate the presence of high levels of activator and inhibitor substance concentrations, respectively (in the yellow patches, both morphogens are strongly concentrated). Initially, the morphogen concentrations and the locations of cell bodies are randomized over the 2-dimensional plane (A). The overabundance of patches at the earlier stages of pattern formation converges, by mutual competition, to a stable pattern of fewer and more separated morphogen concentration peaks (BF).

Figure 2.

Sequence of images from a Cx3D simulation showing the development of a Gierer–Meinhardt pattern over time. The neuronal cell bodies (somata) are colored spheres. Black indicates low morphogen concentration near the cell body. Red and green indicate the presence of high levels of activator and inhibitor substance concentrations, respectively (in the yellow patches, both morphogens are strongly concentrated). Initially, the morphogen concentrations and the locations of cell bodies are randomized over the 2-dimensional plane (A). The overabundance of patches at the earlier stages of pattern formation converges, by mutual competition, to a stable pattern of fewer and more separated morphogen concentration peaks (BF).

Supportive Evidence for the Gierer–Meinhardt Model

Linear Scaling of Spatial Characteristics

Analysis of experimental data from the cortical areas of many species has shown that there is a linear relationship between the patch diameter and mean interpatch distance (Lund et al. 1993; Binzegger et al. 2007; Muir et al. 2011), the spacing being approximately twice the average width of the patches. This as yet unexplained relationship can be understood in the context of our model.

In their mathematical analysis of the Gierer–Meinhardt model, Page et al. (2005) explored the effect of certain parameters on the resulting patterns. They show analytically that the average patch size is proportional to $${D_a } $$ and that the interpatch distance is proportional to $${D_h } $$. If we consider a common scaling factor D for Da and Dh, such that Da = D·Da′ and Dh = D·Dh′, then this result explains the linear relationship between patch diameter and interpatch distance when D and therefore both the diffusion coefficients vary. However, the assumption of different diffusion coefficients in different cortical areas and species is not biologically plausible. The magnitudes of physical variables influencing the diffusion coefficients, such as viscosity or temperature, are not thought to vary by large amounts. We therefore investigated what other parameters influence the spatial profile in a similar manner. We found both analytically and by simulation that the linear scaling is also obtained when the production coefficient (parameter ρ in equations 1) is varied. This effect can be shown as follows: In the converged state $$(\partial a/\partial t = \partial h/\partial t = 0)$$, the decay and the basal production (parameters μa, ρa, μh, and ρh in equations 1) can be neglected as they are usually relatively small:  

(3)
$$\eqalign{ - \rho \displaystyle{{a^2 } \over h} & \approx D \cdot D^{\,\prime} _a \cdot \displaystyle{{\partial ^2 a} \over {\partial x^2 }} \cr - \rho a^2 & \approx D \cdot D^{\,\prime} _h \cdot \displaystyle{{\partial ^2 h} \over {\partial x^2 }}} $$

Together with the result from Page et al. (2005), our results show that varying the diffusion coefficients or the production coefficient changes the scale of the patch size and spacing. ρ and D linearly scale the right or left-hand side term of equations (3), so they influence the spatial characteristics of the pattern by varying the proportion of the production versus the diffusion term. We find that the linear relationship for the case when the production coefficient is varied is valid also under conditions where the basal production coefficients are not negligible, as for example in the simulations in Figure 3. The slope of the linear relationship is similar for our simulations and the experimental data (approximately 2), but depends on the parameters of the model: Different parameters yield different relations. The appropriate linear scaling was also confirmed in simulations done in MATLAB (Supplementary Fig. 1). The MATLAB simulations showed a more linear scaling behavior, which can be explained through the larger simulation grid (150 × 150 cells instead of 100 × 100 in Cx3D) and the regular arrangement of substance producing cells. Figure 3 shows 2 examples of patterns (A,B) produced with 2 different production coefficients (ρ = 600 and ρ = 35) and the scaling behavior across 120 patterns with different production coefficients (B). Also note the higher density of patches with small patch diameter and spacing, a feature that is shared with biological data (Fig. 1B). The Gierer–Meinhardt parameters used in these simulations are listed in Table 1.

Figure 3.

(A) Two Gierer–Meinhardt patterns generated in Cx3D having different characteristic patch diameters and interpatch distances. The only difference between the 2 cases is the production coefficient: ρ = 600 for the left pattern and ρ = 35 for the right pattern. (B) Influence of the production coefficient ρ on the patch size and interpatch spacing. Each sample point (blue) indicates a separate simulation with different ρ. Average patch diameters and interpatch distances are linearly related. The equation of the fitted linear interpolation (red) is y = 1.72 x + 0.15 (R2 = 0.46, P-value <1 × 10−4), approximating closely the slope of 1.65 obtained from the experimental data (Muir et al. 2011). Note the increased density of sampling points toward the lower left, which has also been observed in biological systems.

Figure 3.

(A) Two Gierer–Meinhardt patterns generated in Cx3D having different characteristic patch diameters and interpatch distances. The only difference between the 2 cases is the production coefficient: ρ = 600 for the left pattern and ρ = 35 for the right pattern. (B) Influence of the production coefficient ρ on the patch size and interpatch spacing. Each sample point (blue) indicates a separate simulation with different ρ. Average patch diameters and interpatch distances are linearly related. The equation of the fitted linear interpolation (red) is y = 1.72 x + 0.15 (R2 = 0.46, P-value <1 × 10−4), approximating closely the slope of 1.65 obtained from the experimental data (Muir et al. 2011). Note the increased density of sampling points toward the lower left, which has also been observed in biological systems.

Our interpretation is that the observed relationship of patch size versus interpatch distance arises, because the production rate of morphogens secreted by cells varies between cortical areas and species.

Distribution of Interpatch Angles

We compared the distribution of angles formed between neighboring patches in our simulations with the distributions measured in biological data. Angles formed between sets of neighbors within the Gabriel graph provide a sensitive statistical measure of the spatial structure inherent in the patch layout. For example, Muir et al. (2011) have used the distribution of interpatch angles to differentiate between artificial sets of points and experimental measurements. We measured this distribution over our simulation results and compared it with data from the primary visual cortex in the cat and macaque monkey (Fig. 4). Our analysis confirms that the spatial arrangement of patches produced by our model matches experimental observations very well. The model distribution was inferred from the same patterns composed of 100 × 100 cells that were used to demonstrate the linear scaling relationship in the previous section.

Figure 4.

Comparison of the probability density functions of the measured interpatch angles from experimental data of the cat and macaque monkey (blue) and Cx3D simulations (red). Shading corresponds to the 90% confidence interval, as estimated by bootstrap analysis. The distributions are not different according to a 2-sample Kolmogorov–Smirnov test at the 5% significance level.

Figure 4.

Comparison of the probability density functions of the measured interpatch angles from experimental data of the cat and macaque monkey (blue) and Cx3D simulations (red). Shading corresponds to the 90% confidence interval, as estimated by bootstrap analysis. The distributions are not different according to a 2-sample Kolmogorov–Smirnov test at the 5% significance level.

Robustness of Patchy Patterns to Parameter Variation

An important condition for the biological plausibility of the pattern formation process is its robustness against parameter variations. We investigated this aspect by assessing the sensitivity of the patchy patterns to the various parameters (diffusion, decay, production, and basal production) of the Gierer–Meinhardt model. For each simulation, we chose a different set of parameters by modifying a template parameter set. This was done as follows: P = PO + ɛσ, where P is a parameter from the Gierer–Meinhardt model, PO is the original template parameter before modification, and ɛσ is the modification variable. For each parameter, this modification was drawn from a Gaussian distribution N(0,σ2). The variance σ2 of the modification was increased stepwise. For each level of variation, we conducted 30 simulations in MATLAB. To treat all the parameters in the same way, the standard deviation σ was chosen to be proportional to the parameter magnitude. This implies that the magnitude of the coefficient of variation cv = σ/μ (where σ is the standard deviation and μ the mean, which is equal to the original parameter magnitude) was the same for all parameters. Our results (Supplementary Fig. 3) demonstrate that the patchy pattern formation process is indeed robust. The patches form successfully for the large levels of variation, that is, 100% of all patterns were classified as patchy even when the standard deviation of the parameters was at approximately 25% of the parameter magnitude itself. From the point of view of biological stability, the robustness of the pattern formation process is also indicated by the fact that the irregular spatial arrangement of the neuronal somata in Cx3D does not prevent the formation of patches.

Development From a Single Precursor Cell

So far, we have assumed that the sheet of cells across which patches are formed is provided externally to our model. Now, we will consider the question whether these systems can develop from a single precursor cell in a self-organizing way. In this case, the whole structure has to develop from rules encoded only in the initial cell. This precursor and all subsequent daughter cells divide in the horizontal plane and give rise to a circularly arranged sheet of cells able to secrete the necessary morphogens.

The simulated developmental process is divided into 2 phases: First, a phase of cell division that produces a sheet of cells by mitotic replication of the mother cell, followed by a phase during which the cells of the sheet secrete the morphogens that control patch formation. These 2 phases are determined by the concentration of an intracellular signaling molecule, produced at a constant rate by each cell. This molecular marker is passed on at each cell division by the mother cell to its daughter cells, which continue producing it at a constant rate. This substance therefore acts as an internal clock. When the concentration of this signal exceeds a chosen threshold, cellular division stops and the soma begins releasing the activator and inhibitor substances according to the Gierer–Meinhardt kinetics as described by equations (1). Figure 5 confirms that the patchy patterns can also be generated under these self-organized conditions. A video of this self-organizing process is included as Supplementary Movie 5.

Figure 5.

Self-organized development of a patchy pattern. The initial precursor cell divides in the 2D plane, giving rise to several daughter cells (A: early division phase, B: progressed division phase). Cell division continues until an intracellular concentration reaches a predefined threshold, which stops the cell division. Based on this preplate, the Gierer–Meinhardt reaction–diffusion mechanism leads to a patchy pattern (C). A video of the self-organization of multiple patchy patterns is included in Supplementary Material (Supplementary Movie 5).

Figure 5.

Self-organized development of a patchy pattern. The initial precursor cell divides in the 2D plane, giving rise to several daughter cells (A: early division phase, B: progressed division phase). Cell division continues until an intracellular concentration reaches a predefined threshold, which stops the cell division. Based on this preplate, the Gierer–Meinhardt reaction–diffusion mechanism leads to a patchy pattern (C). A video of the self-organization of multiple patchy patterns is included in Supplementary Material (Supplementary Movie 5).

Multiple Patch Systems

It remains an open question whether the patchy architecture of the cortex is a single global pattern, or whether multiple nonoverlapping or overlapping patch systems are present (Muir and Douglas 2011). The results of our prepatterning simulations above can explain the emergence of a single, static patch system. However, those simulations use only a single activator–inhibitor pair. Multiple patch systems are possible when multiple pairs are present. Different patch systems develop as a consequence of the separate activator–inhibitor pairs, and if the substances of the activator–inhibitor pairs inhibit each other (Fig. 6A), then patches of the different systems are prevented from emerging at the same location.

Figure 6.

(A) Interactions of 2 Gierer–Meinhardt systems. In addition to the basic antagonism of the activator and inhibitor morphogens (as used in the single-patch system scenario), the activator substances of different Gierer–Meinhardt systems inhibit each other. Sufficiently strong inhibition between different activator morphogens ensures that patches belonging to distinct systems do not overlap. (B) Simulation of 4 interacting patch systems from a Cx3D simulation. Each system is based on 2 morphogens interacting according to the Gierer–Meinhardt model. Red, green, blue, and purple indicate the concentrations of distinct activator substances. Repulsion constants between activators were chosen to be identical for all of the activator pairs. (C) Multiple patch systems defined by inhibitory substances (same simulation as in Fig. 6B). Red, green, blue, and purple indicate the concentrations of distinct inhibitor types. In contrast to the scenario where patches are defined by the activator substance concentration, here, there is a smooth overlap between patches because there is no mutual inhibition between inhibitory substances. Patch systems are still guaranteed to develop at distinct locations, because the different activator substances compete with each other, such that also the patches of the inhibitor concentrations are at different locations.

Figure 6.

(A) Interactions of 2 Gierer–Meinhardt systems. In addition to the basic antagonism of the activator and inhibitor morphogens (as used in the single-patch system scenario), the activator substances of different Gierer–Meinhardt systems inhibit each other. Sufficiently strong inhibition between different activator morphogens ensures that patches belonging to distinct systems do not overlap. (B) Simulation of 4 interacting patch systems from a Cx3D simulation. Each system is based on 2 morphogens interacting according to the Gierer–Meinhardt model. Red, green, blue, and purple indicate the concentrations of distinct activator substances. Repulsion constants between activators were chosen to be identical for all of the activator pairs. (C) Multiple patch systems defined by inhibitory substances (same simulation as in Fig. 6B). Red, green, blue, and purple indicate the concentrations of distinct inhibitor types. In contrast to the scenario where patches are defined by the activator substance concentration, here, there is a smooth overlap between patches because there is no mutual inhibition between inhibitory substances. Patch systems are still guaranteed to develop at distinct locations, because the different activator substances compete with each other, such that also the patches of the inhibitor concentrations are at different locations.

The step from single to multiple patch systems is a small one, because the interactions between activator and inhibitor morphogens are unchanged. The developmental cost of increasing the number of patch systems is also low in terms of the extra information required in the genetic code: The code length increases only linearly with the number of patch systems.

We extended the secretion model to permit multiple patch systems by including a term that reduces the activator production in the presence of activator morphogens from other Gierer–Meinhardt systems at the same location (equations 2). High activator levels produce high inhibitor concentrations, and inhibitor substances can also cross-inhibit the activators of other systems. Both mechanisms prevent the patches of different activator–inhibitor systems from appearing at the same location (as for example in Fig. 6B, C).

Experimental data obtained from adjacent injections in inferior temporal (IT) cortex of the macaque monkey suggest that there is very little or no overlap in patch locations (Tanigawa et al. 2005). This situation occurs in our model when the mutual inhibition is sufficiently strong. On the other hand, the patches of several patch systems can overlap if, instead of the activator substance, the concentration of the inhibitor substance defines the patches. In this case, the appearance of patches is modified, in 2 ways: They are larger because of the larger diffusion constant (a necessary condition for stationary patterns); and they are smoother because the different inhibitory substances do not directly interact with one another. In this model, cells can exist inside multiple different patch systems. Figure 6C shows the inhibitory morphogen patches from the same pattern as in Figure 6B.

Topography of Patch Systems

Functional maps, for example maps of orientation selectivity in the visual cortex, are quasi-periodic and have a topographic neighborhood relation (Hubel et al. 1978; Bosking et al. 1997). If we assume similar characteristics for the superficial patch system, then the locations of the patches have to fulfill certain conditions. The targets of projections from nearby neurons should share common features. Neurons that are located close to each other should make target patches in locations that are only slightly offset (Lund et al. 2003; Muir and Douglas 2011; Muir et al. 2011). For example, in a scenario with 3 patch systems (U, V, and W), patches of system U could be on average closer to the patches of system V than to the patches of system W. We investigated how such a topography could be achieved with multiple Gierer–Meinhardt systems that interact with each other.

Our approach is as follows. We assume that there is Gaussian variation in the cellular secretion dynamics. Consequently, some cells have higher production and basal production coefficients than others. This variation is sufficient to predispose some regions to express the high levels of morphogen concentration, so constituting a patch. Therefore, with a sufficiently high variance in the cellular production, the patches of the different activator–inhibitor systems will fall predominantly into the same specific regions that are composed of highly productive cells (Fig. 7B).

Figure 7.

(A and B) Effect of variation in the production efficacy on patch locations. (A) Superimposition of the activator morphogen concentrations (red and blue label the 2 activator types) of 2 different noninteracting Gierer–Meinhardt systems. There is no variation, and the repulsion constant between the 2 activators is set to 0. Therefore, patches of both systems emerge at random and independent locations. For clearer visualization, overlapping patches are indicated with white circles. (B) Obtained when there was a nonzero variance on both the cellular production coefficient ρ and the basal production coefficients ρa and ρh (equations 1). The variance allows some patchy organizations to be preferred over others, because those patches will fall predominantly onto the same highly productive locations. High concentrations of both activator types are colored purple, because red and blue are added. (C and D) The magnitude of the repulsion coefficients determines the distance of the different patch systems. (C) Two superimposed concentration patterns of the activator morphogens from interacting Gierer–Meinhardt systems with small repulsion coefficients. In addition, production coefficients across the cells vary in a Gaussian manner. Due to the low repulsion, the patches of the 2 systems end up quite close to each other. (D) Two interacting Gierer–Meinhardt systems with strong repulsion. The patches of the different patch systems repel another more from the preferred locations than in (C). This was confirmed by calculating the average distances between patches belonging to distinct patch systems.

Figure 7.

(A and B) Effect of variation in the production efficacy on patch locations. (A) Superimposition of the activator morphogen concentrations (red and blue label the 2 activator types) of 2 different noninteracting Gierer–Meinhardt systems. There is no variation, and the repulsion constant between the 2 activators is set to 0. Therefore, patches of both systems emerge at random and independent locations. For clearer visualization, overlapping patches are indicated with white circles. (B) Obtained when there was a nonzero variance on both the cellular production coefficient ρ and the basal production coefficients ρa and ρh (equations 1). The variance allows some patchy organizations to be preferred over others, because those patches will fall predominantly onto the same highly productive locations. High concentrations of both activator types are colored purple, because red and blue are added. (C and D) The magnitude of the repulsion coefficients determines the distance of the different patch systems. (C) Two superimposed concentration patterns of the activator morphogens from interacting Gierer–Meinhardt systems with small repulsion coefficients. In addition, production coefficients across the cells vary in a Gaussian manner. Due to the low repulsion, the patches of the 2 systems end up quite close to each other. (D) Two interacting Gierer–Meinhardt systems with strong repulsion. The patches of the different patch systems repel another more from the preferred locations than in (C). This was confirmed by calculating the average distances between patches belonging to distinct patch systems.

The system of equations (2) implies that substances of different Gierer–Meinhardt systems mutually inhibit one another. Figure 7C,D are 2 examples in which the repulsion strength between patch systems varies (red and blue indicate the 2 different systems). The repulsion between the morphogens in Figure 7C is weaker than in D. This is reflected in the distances between the red and blue patches, because the red patches in (C) are usually closer to the blue patches than in the pattern (D). Different mutual repulsion strengths influence the average spacing that will result between patches from different systems. Our simulations indicate that the distance between patch systems is more consistent when the inhibitory morphogen also inhibits the activator production of other types, because the fast diffusion of the inhibitor substance leads to smoother patches. The pattern in Figure 7D was generated by the same model parameters, but with a stronger repulsion parameter r and an additional repulsion from the inhibitory morphogen.

Analogously, one can extend the model to include several types of activators and inhibitors. The set of repulsion coefficients imposes neighborhood relationships between the several patch systems. Theoretically, this kind of topographic relationship could also be imposed by the repulsion coefficients alone, without introducing variance in cellular production coefficients. Cellular inhomogeneity causes the different patch systems to compete for the most productive regions. In addition, cells that produce very little morphogen will not express patches and so will not participate in the superficial patch system. There are several lines of evidence that support the existence of such regions. For example, some locations in the visual cortex, namely the orientation pinwheel centers, do not collectively form long-range projections (Sharma et al. 1995; Yousef et al. 2001; Mariño et al. 2005). An early proposal for the patch system was that neurons in some regions make long-range projections, while neurons in other regions may not (Rockland and Lund 1982; Rockland et al. 1982). These points can be explained by our model by incorporating variation in morphogen production efficacy across the cells. As described later, our (simplified) model of axonal growth assumes that the long-range projections are made only from cells that lie inside a patch. Because the regions that are composed of cells with low morphogen production efficacy do not induce patches, these regions do not participate in the superficial patch system.

Laminated Structure

Although the overall laminar organization of the patch system is still unclear, the characteristic long-range projections have been observed mainly in the superficial layers. We investigated whether our hypothesis is consistent with patch expression only in specific layers, by simulating the development of a laminated structure with periodic long-range clustered projections in only the superficial layers. We achieved this by following the protomap concept of Rakic (1988). First, the cells of the cortical plate generate a patchy system as described above. If the activator concentration in a cell exceeds a given threshold, it becomes imprinted as belonging to a patch. This model can be extended to multiple patch systems on the same plate, as will be explained later.

All the cells of the initial patterned sheet divide horizontally, so that daughter cells form the next higher sheet (for simplicity we do not model here the detailed process of cortical lamination). Each division of a cell causes the decay of an intracellular lamina signal substance, so allowing the cells to track their radial position. Some of these cells have been imprinted as potential secretors as a result of the earlier Gierer–Meinhardt process, and those that reach layer II/III will secrete the appropriate activator morphogen and so contribute to the formation of a patch. The patchy expression pattern then provides the basis for establishing the appropriate axonal connectivity, as shown in Figure 8.

Figure 8.

Daisy architecture in a laminated cortical structure. Multiple patch systems (A: diagonal top-down view, B: straight top-down view) generated from an initial 2D sheet, which is the cortical plate before lamination of the cortex. Information that establishes which system a cell belongs to is preserved in the vertical direction. Red, green, blue, and purple indicate cells that belong to the 4 different patch systems. Supragranular layers are labeled gray, infragranular are black. Cells labeled as belonging to a patch system, and located in layer II/III, continue producing the activator substance that is used by the axonal growth cones to establish the long-range connectivity, as shown in Figure 9B. The simulation consists of 100 × 100 cells in a slab 18 cells thick.

Figure 8.

Daisy architecture in a laminated cortical structure. Multiple patch systems (A: diagonal top-down view, B: straight top-down view) generated from an initial 2D sheet, which is the cortical plate before lamination of the cortex. Information that establishes which system a cell belongs to is preserved in the vertical direction. Red, green, blue, and purple indicate cells that belong to the 4 different patch systems. Supragranular layers are labeled gray, infragranular are black. Cells labeled as belonging to a patch system, and located in layer II/III, continue producing the activator substance that is used by the axonal growth cones to establish the long-range connectivity, as shown in Figure 9B. The simulation consists of 100 × 100 cells in a slab 18 cells thick.

Axonal Growth Model

With the patchy morphogen concentrations in place, the neurons can go on to connect (for example) to like patches by intralaminar axonal projections under active guidance of their growth cones. Growth cones follow the concentration gradient of relevant morphogenic guidance cues, bifurcating recursively as a probabilistic function of the cue concentration (see Materials and Methods for a detailed description of the axonal branching rules). Consequently, laterally extending axons ramify in target patches.

We studied axon growth in the case of single and multiple patch systems. In the first case, the lamina of cells generates only a single global patch system. In this scenario, all (or a randomly selected percentage of) the cells that comprise the cortical sheet extend their axons. The axons follow the gradient of the activator morphogen to grow toward the patches and make their arborizations, leading to the clustered axonal long-range projections shown in Figure 9A.

Figure 9.

(A) Development of axonal long-range projections on a single-patch system. Growth cones of randomly selected cells follow the gradient of an extracellular cue secreted by cells located inside the patches. This same cue promotes more dense arborization by increasing the probability of the axon to bifurcate. (1) Initially, cells extend their axons in random directions. The yellowish somata are the strongest producers of the guidance cue. Newly extended axons and their cell bodies are colored red. (2 and 3) The region of high guidance cue concentration (yellow cell bodies) already has denser arborizations. (4) Arborization pattern at progressed stage. For the sake of clearer visualization, we set the bifurcation probability outside the patches to be quite low. (B) Sequence of axonal growth on multiple patch systems from the pattern in Figure 6B. This axonal growth rule makes use of many patchy patterns generated by the Gierer–Meinhardt model. In contrast to the simpler growth rule with only 1 global patch system, in this case, the location of the soma determines the targeted activator type (different colors indicate different types). Axons that arise from cell bodies lying in a patch will make long-range projections only to the patches of a similar identity. Cell bodies that are not located inside a region of high activator and inhibitor concentration levels make only nonspecific local connections. The specificity of the resulting connections can be controlled by the parameters of the growth model. For example, as mentioned in the Materials and Methods section, the probability to bifurcate is given by $$p_{{\rm bifurcate}} = u\cdot c_a + v$$, where $$c_a $$ is the concentration of the activator substance and u, v are variable parameters. Large v increases the probability to branch outside of the patches, leading to less specific connectivity. In the scenario with multiple patch systems, a growth cone could also target activator patches of multiple types with different scaling parameters u, such that the communication between different patch systems is increased. Additionally, the bifurcation could be guided by the inhibitor substance, which increases the overlap between different patch systems (as shown in Fig. 6B,C).

Figure 9.

(A) Development of axonal long-range projections on a single-patch system. Growth cones of randomly selected cells follow the gradient of an extracellular cue secreted by cells located inside the patches. This same cue promotes more dense arborization by increasing the probability of the axon to bifurcate. (1) Initially, cells extend their axons in random directions. The yellowish somata are the strongest producers of the guidance cue. Newly extended axons and their cell bodies are colored red. (2 and 3) The region of high guidance cue concentration (yellow cell bodies) already has denser arborizations. (4) Arborization pattern at progressed stage. For the sake of clearer visualization, we set the bifurcation probability outside the patches to be quite low. (B) Sequence of axonal growth on multiple patch systems from the pattern in Figure 6B. This axonal growth rule makes use of many patchy patterns generated by the Gierer–Meinhardt model. In contrast to the simpler growth rule with only 1 global patch system, in this case, the location of the soma determines the targeted activator type (different colors indicate different types). Axons that arise from cell bodies lying in a patch will make long-range projections only to the patches of a similar identity. Cell bodies that are not located inside a region of high activator and inhibitor concentration levels make only nonspecific local connections. The specificity of the resulting connections can be controlled by the parameters of the growth model. For example, as mentioned in the Materials and Methods section, the probability to bifurcate is given by $$p_{{\rm bifurcate}} = u\cdot c_a + v$$, where $$c_a $$ is the concentration of the activator substance and u, v are variable parameters. Large v increases the probability to branch outside of the patches, leading to less specific connectivity. In the scenario with multiple patch systems, a growth cone could also target activator patches of multiple types with different scaling parameters u, such that the communication between different patch systems is increased. Additionally, the bifurcation could be guided by the inhibitor substance, which increases the overlap between different patch systems (as shown in Fig. 6B,C).

In the case of multiple patch systems, axons follow the gradients of different morphogens for each patch system. Axonal outgrowth is triggered when a neuronal soma senses a suprathreshold activator concentration, subsequently extending an axon that follows the positive gradient of only this morphogen. This simple rule leads to clustered axonal projections over separate patch systems, as shown in Figure 9B. A video of a simulation implementing this growth rule is included as Supplementary Movie 6.

Additional modifications of the axonal growth patterns could be made in the scenario of multiple patch systems. For example, the patches could be defined by high inhibitor morphogen concentrations, as discussed before. In this case, the target regions of neurons belonging to the different patch systems can overlap. Cells lying in an overlapping location could randomly switch to one of the patch identities before the initiation of neurite outgrowth. This would increase the potential electrophysiological communication between different patch systems (although in the separated patch systems scenario there is no strict boundary, so there can always be communication between different patch systems). Another possible extension that would increase the degree of intercommunication could involve the growth rule itself, so that some additional patch systems could also be connected to one another. For example, neurons in the patches of a patch system U could project to patches belonging to systems U and V, while the patches of system V project to systems U and W. A mathematical analysis of the topology and activity dynamics in different scenarios for the superficial patch system has been performed by Voges et al. (2010).

Our goal here was not to model the detailed biological processes of the growth cone and the resulting morphology, but simply to incorporate principles that lead to patchy connectivity patterns and that match experimental observations (such as the distribution of angles between neighboring patches, addressed in the Results section).

Discussion

We have demonstrated by detailed simulation of neuronal development how the experimentally observed superficial patch system of the neocortex can arise on the basis of the quasi-periodic spatial patterns produced by a Gierer–Meinhardt type reaction–diffusion system. This entire process arises out of the genetic-like specification inserted into a single precursor cell. This progenitor cell undergoes cell division, distributing its genetic instructions to its offspring, who finally give rise to the cortical plate. The individual genetic codes instruct the cells of the plate to secrete a small number of different morphogens. These morphogens have either localized or extensive extracellular diffusion; and they either activate or suppress the production of other morphogens. Overall, these morphogenetic network interactions express Gierer–Meinhardt reaction–diffusion dynamics that converge to a spatial pattern of interlaced hexagonal lattices of similarly labeled (by their morphogen profile) neurons. These labeled neurons then spread radially to become a more superficial layer of neurons and then begin to extend long-range intralaminar axons. The direction of growth of these axons, and their propensity to bifurcation, are determined by the gradient and concentration of their target morphogen. Consequently, axons have a bias to grow toward adjacent regions of like-labeled neurons and to ramify there (Fig. 9), so giving rise to the experimentally observed patchy connectivity.

Overall, our model demonstrates how the patchy connectivity pattern of cortex can construct itself on the basis of only local cell–cell chemical interactions and each cell's (similar) genetic code, without the support of external, organizing electrical spike patterns.

Apparently Complex Structure Arising From a Simple Mechanism

The concept that the complicated spatio-temporal patterns observed in animal morphogenesis could be due to a relatively simple underlying molecular mechanism was proposed by Alan Turing in his seminal paper of 1952. He demonstrated that a suitably constrained interaction between 2 fundamental physical processes—reaction and diffusion—could give rise to many kinds of pattern formation observed in nature. This mechanism is based on the secretion, decay, and diffusion dynamics of different substances that compete for control of their sources. Small random fluctuations in the initial concentrations can lead to global symmetry-breaking and the emergence of stationary spatial “Turing” patterns. Counter-intuitively, symmetry-breaking is driven by diffusion: A necessary condition for the symmetry-breaking is that the diffusion coefficients differ between competing morphogens. Unfortunately, Turing's original formulation of the reaction–diffusion system is not molecularly plausible, because it allows negative concentrations. However, in 1972, Gierer and Meinhardt proposed the feasible formulation that we have chosen as a basis for inducing spatial organization. Other reaction–diffusion models that are also based on interacting antagonistic morphogens (e.g. the Gray–Scott [Gray and Scott 1983; McGough and Riley 2004] or Schnakenberg model [Schnakenberg 1979; Iron et al. 2004]) may yield similar results.

Grossberg (1976) explored the use of reaction–diffusion systems to model neural map formation. Since then, several reaction–diffusion systems have been used to model the development of neural tissue. For example, Cartwright (2002), Striegel and Hurdal (2009), Lefèvre and Mangin (2010), and Garzón-Alvarado et al. (2011) investigated the problem of the generation of gyri and sulci during cortical development. Other models based on reaction–diffusion dynamics have been used to explain ocular dominance and orientation maps (Swindale 1980; Grossberg and Olson 1994), as well as disparity selectivity (Siddiqui and Bhaumik 2011). All of those models were essentially analytical and so differ from the detailed physical cellular growth simulations that we report here.

Experimental Support for Reaction–Diffusion Processes

The relevance of Turing patterns as an explanation of natural morphogenesis, and specifically cortical development, is matter of debate. The mechanism is elegant and appealing, but experimental verification is difficult because it requires not only identification of the different morphogens, but also their in vivo diffusion or cellular reaction constants, data which are very difficult to obtain. Nevertheless, there is substantial circumstantial evidence for the role of reaction–diffusion processes in development. Turing patterns have been confirmed experimentally in chemical reactions (Castets et al. 1990; De Kepper et al. 1991), and there is experimental evidence that these systems play a role in animal development. For example, 2 antagonistic pattern-forming substances have been shown to be involved in the generation of the oral–aboral axis in sea urchins (Duboc et al. 2004), and in mouse hair follicle formation (Mou et al. 2006). It can predict the time evolution of the skin patterns of the angelfish Pomacanthus (Kondo and Asai 1995) and explain head regeneration in the freshwater coelenterate hydra (Gierer and Meinhardt 1972), or stripes on the shells of molluscs (Meinhardt and Klingler 1987).

Morphogens that follow dynamics of the activator–inhibitor type do occur in development. For example, transforming growth factor (TGF)-β signals can act as activator–inhibitor morphogens (Schier 2009), as demonstrated in the zebrafish blastula (Meno et al. 1999; Chen and Schier 2002; Müller et al. 2012). Müller et al. (2012) have shown that the inhibitor has a larger effective diffusion constant than the activator, which is a necessary condition for stationary Turing patterns. There is evidence that TGF-type morphogens have an autocatalytic role in vertebrate limb formation (Newman and Bhat 2007). The sonic hedgehog (SHH) signaling molecule, which is involved in various pattern formation mechanisms in the central nervous system, acts together with fibroblast growth factor as an activator–inhibitor system that is responsible for the formation of ridges in the mammalian palate (Economou et al. 2012). The WNT/bone morphogenetic protein activator–inhibitor pair has been implicated in hair and feather follicle development (Sick et al. 2006; Plikus et al. 2011). Other possible candidates of morphogens are Semaphorins (Matsuoka et al. 2011). Moreover, these morphogens have been shown to act also as guidance cues (Webber et al. 2003; Schnorrer and Dickson 2004; Charron and Tessier-Lavigne 2007; Zou and Lyuksyutova 2007; Sánchez-Camacho and Bovolenta 2009; Gordon et al. 2010).

It is difficult to estimate the parameter values of the morphogenetic substances in our simulations, because these values depend also on assumptions concerning somatic diameter, density of morphogen producing cells, computation of diffusion, etc. Consequently, we have not assigned units to the parameters. However, we are able to make qualitative estimates, such as the proportions that some coefficients have with respect to one another. For example, the inhibitor diffusion coefficient is approximately 3 orders of magnitude larger than the activator diffusion coefficient. The decay rate of the inhibitor substance in our simulations was always larger than the rate of the activator, by about a factor of two.

Agreement with Experimental Data, and Concepts of Cortical Organization and Development

The spatial patterns produced by our model are in good agreement with experimental observations, such as the remarkable observation of the linear relationship between patch diameter and spacing across cortical areas and species (Lund et al. 1993; Binzegger et al. 2007; Muir et al. 2011). It also exhibits the correct distribution of angles between neighboring patches (Muir et al. 2011). Differences in the morphogen production coefficients in the different cortical areas and species are the only assumptions required to explain this intriguing relationship. Our hypothesis could be confirmed by identifying and modifying the cellular secretion efficacies of the hypothetical morphogens.

Our reaction–diffusion model is consistent with Rakic's protomap concept (Rakic 1988), which is likely a central organizing principle of the much disputed “cortical column” (Mountcastle 1957; Hubel and Wiesel 1977; da Costa and Martin 2010). Postmitotic cells, by their migration outwards along the radial glial cells and their formation of the cortical lamination, project the patch lattice columnwise out of the 2-dimensional neural plate. Several studies provide evidence for this similar ontogenetic unity along the radial dimension (Luskin et al. 1988; Kornack and Rakic 1995; Rakic 1995), which can be seen in the laminar organization of the patch system (Rockland 1985; Galuske and Singer 1996; Angelucci and Sainsbury 2006).

Of course, the original protomap hypothesis was made in the context of the rodent brain, and there is evidence that the superficial patch system does not exist in rodents (Van Hooser et al. 2006; Muir and Douglas 2011). However, the important principle here is that organizational patterns laid down in the ventricular zone propagate radially and are elaborated in the developing cortical laminations (Kennedy and Dehay 1993). The reason why certain species do exhibit patchy connectivity may be due to the combinations of production and diffusion rates in vivo that do not fulfill the conditions for stationary hexagonal patterns.

It remains unknown whether the superficial patches reflect a single lattice, or whether there are several distinct or overlapping patch systems (Muir and Douglas 2011). Studies based on adjacent injections using multiple tracers point in the direction of multiple patch systems (Matsubara et al. 1987; Tanigawa et al. 2005). Tanigawa et al. (2005) used adjacent injections of multiple tracers to demonstrate that intra-areal long-range projections cluster at nonoverlapping locations. Our simulations show that nonoverlapping and partially overlapping patch systems can be easily realized with multiple interacting Gierer–Meinhardt systems. In this scenario, the regions targeted by an axon depend on the location of it's soma: Somata located in a given patch extend axons only toward patches belonging to the same system. We have also shown that patches in our model can overlap, if the patch identity is conferred by the concentration of the inhibitor substance.

These anatomical results offer a link to features of functional maps [periodicity and topographic mapping, described for example by Mitchison and Crick (1982) and Bosking et al. (1997)]. For example, the orientation tuning of neurons inside connected patches tends to be similar (Bosking et al. 1997; Schmidt et al. 1997). Further work in this direction could investigate how the transitions between patches of different reaction–diffusion systems could provide a scaffold for establishing the connectivity of functional maps such as orientation selectivity.

There is also evidence for a relationship between the patch system and metabolically derived patterns. For example, cytochrome oxidase–reactive regions and the superficial patch system share some features (Rockland and Lund 1983; Muir and Douglas 2011; Muir et al. 2011). We did not address ocular dominance columns or cytochrome–oxidase (CO) blobs in this work, because they stem from thalamocortical projections, and we restrict ourselves to corticocortical connections. However, it is intriguing that CO blobs in the cat visual cortex can develop in the absence of visual input (Murphy et al. 2001), raising the possibility that molecular markers such as those used in our model are involved in patterning geniculocortical afferents. Bressloff and Oster (2010) have studied just this scenario, using a patchy scaffold of markers which modulate synaptic plasticity and promote the formation of ocular dominance columns. Our work complements their model since it provides an explanation for the development of this periodic specification.

Axonal Arborization Patterns

The orderly axonal projection patterns depend on simple and biologically plausible axonal growth rules, which make use of morphogen concentration profiles for targeting arborization locations over long distances. This direct sensitivity to an organizational signal makes the arborization quite specific, rather than having an initial widespread exuberance of projections followed by pruning (Callaway and Katz 1990).

Importantly, in our model, axonal outgrowth and organization do not depend on electrical signaling that has particular sensory statistics (as for example proposed by Grabska-Barwińska and von der Malsburg (2008), suggesting that activity waves and locally periodic response patterns shape the connectivity). Electrical signaling might be required for refinement of the crude clusters established at the early developmental stage. This does not contradict the experiments of Ruthazer and Stryker (1996), where suppression of spiking activity using the sodium channel blocker tetrodotoxin prevented the formation of clustered long-range projections in the ferret. We argue that electrical activity might be permissive rather than instructive, with growth at the very early stages of development based on the intrinsic pattern, so accounting for the presence of the superficial patch system in cortical areas that are remote from the potential organizing properties of sensory afferents.

There is indeed evidence that neural activity has a large impact on axonal branching (Uesaka et al. 2006). In this vein, it has been shown in the cat visual cortex that at the very early stages of connectivity development, long-range intra-areal projections are not yet clustered (Price 1986; Callaway and Katz 1990). Only at the end of the second postnatal week is extensive clustering observed. Our interpretation of these data is that axonal outgrowth is under the control of a 2-step process. Axons initially search for their target patches without extensive ramification. A subsequent axonal cluster-forming process is invoked only when an axon finds a suitable target region and the appropriate electrical activity conditions are met. Due to retinal waves and recurrent thalamocortical loops, animals are exposed to relatively structured activity possibly very early in development, which could provide enough correlation structure for Hebbian mechanisms to play a role (Wong 1999; Weliky and Hall 2000; Butts 2002; Huberman et al. 2008). Studies in other species also show that a crude form of the patch system is present at very early stages. In the ferret patches develop before eye opening (Durack and Katz 1996); and in macaque monkey, the system can even be observed prenatally (Coogan, Van Essen, et al. 1996). Consequently, we predict that as long as the branching and axon guidance machinery is not affected, disturbance of the pattern of electrical input to the neurons will not abolish the development of the superficial patch system. This has been shown to be the case for retinal activity (Callaway and Katz 1991; Ruthazer and Stryker 1996).

Simulations of Cellular Processes in Physical 3D Space

The simulation methods used in the paper demonstrate a substantial advance in computational neuroscience. In classical mathematical models, the environment is abstracted as a regular 2-dimensional grid in which each discrete element of space participates in both diffusion and reaction. In our work, we have simulated the behavior of physical cells, behaving in a physically realistic 3D space. Each cell is an independent agent able to interact only with its local environment; there is no global organizer. Such physical simulations in Cx3D (for example) offer a novel and powerful tool for the investigation of anatomical and developmental processes that is likely to be at least as influential as the theoretical electrophysiological investigations performed in simulators such as NEURON (Hines and Carnevale 1997). Of course, dense multiagent simulations such as ours come with a high computational cost and also bring some interesting new questions. One of these, relevant to the present paper, is the question of stability.

Turing patterns have been shown to be very robust against noise (Pena and Perez-Garcia 2001; Leppänen et al. 2003). This robustness is especially important in a biological context, because there is inherent variation in the cellular secretion. We have shown that Turing patterns that emerge in our simulations are stable to disturbances, noise or the initial distribution of morphogen concentrations. However, it is not only the stability of the reaction–diffusion system that is at issue here. It is important to note that the developmental process arises by the cell division from a precursor, in which only local morphogenetic production rules and axonal growth rules are specified. Thus, overall stability of the patchy lattice depends also on the developmental strategy (Zubler et al. 2011) in which the reaction–diffusion processes are embedded.

Our goal in this paper was to describe only the principles of developmental construction of the patch system, rather than attack the more difficult problem of how these principles play out when embedded in a detailed simulation of cortical development. However, we have previously shown (Zubler and Douglas 2009) how to produce a laminated cortical structure in an unprepared environment, and how to make layer specific branching patterns in Cx3D. It remains to bring these 2 lines of work together, to demonstrate how the patchy protomap promotes columnar organization of inter- and intralaminar connections in the context of a large-scale simulation of cortical development, a task in which we are presently engaged.

Conclusion

In the modern interdisciplinary approach to neuroscience, theoretical models and their simulation, such as we present here, complement experimental neuroscience in 2 important ways: They explore and they explain. Our model explores implications of the assumption that the development of the general and rather complex patterning of the cortical patch system depends on a simple underlying molecular mechanism: The presence of a cooperative–competitive interaction of the morphogen secretion between cells. Of course, we have not yet confirmed the existence of the crucial proposed cooperative–competitive morphogens. We have, however, presented evidence from the literature for the probable existence of such substances. And our simulations have shown how, on the basis of this simple assumption, the model is able to explain a wide variety of experimental observations such as the occurrence of the patches, their spatial distribution, their linear scaling over areas and species, and their independence (initially) of electrical activity in well-formed neuronal circuits that occur only later in development. Our model also explains how this process can unfold from a single precursor cell, so demonstrating how complex organization of connections could be simply genetically encoded. Importantly, our simulations are of physical cells in a physically realistic environment, in which cells act as individual agents able to interact only with their local environment. This means that our results extend beyond the ideal of mathematical models, toward what is physically realizable in biology and future technology.

Supplementary Material

Supplementary material can be found at: http://www.cercor.oxfordjournals.org/

Funding

This work was supported by European Union Project grant FP7-216593 “SECO.” Funding to pay the Open Access publication charges for this article was provided by the Institute of Neuroinformatics, UZH/ETH Zürich.

Notes

The authors thank Michael Pfeiffer, Ueli Rutishauser, Adrian Whatley, Matthew Cook, Kevan Martin, Elisha Ruesch, Jiri Kriz, and Boris Wrobel for valuable feedback and for the very supportive critical readings of the manuscript. Conflict of Interest: None declared.

References

Angelucci
A
Sainsbury
K
Contribution of feedforward thalamic afferents and corticogeniculate feedback to the spatial summation area of macaque V1 and LGN
J Comp Neurol
 , 
2006
, vol. 
498
 (pg. 
330
-
351
)
Binzegger
T
Douglas
R
Martin
K
Stereotypical bouton clustering of individual neurons in cat primary visual cortex
J Neurosci
 , 
2007
, vol. 
27
 (pg. 
12242
-
12254
)
Binzegger
T
Douglas
R
Martin
K
Feldmeyer
D
Lübke
JHR
An axonal perspective on cortical circuits
New aspects of axonal structure and function
 , 
2010
New York (NY)
Springer
(pg. 
117
-
139
)
Bonabeau
E
From classical models of morphogenesis to agent-based models of pattern formation
Artif Life
 , 
1997
, vol. 
3
 (pg. 
191
-
211
)
Bosking
W
Zhang
Y
Schofield
B
Fitzpatrick
D
Orientation selectivity and the arrangement of horizontal connections in tree shrew striate cortex
J Neurosci
 , 
1997
, vol. 
17
 (pg. 
2112
-
2127
)
Bressloff
P
Oster
A
Theory for the alignment of cortical feature maps during development
Phys Rev E
 , 
2010
, vol. 
82
 pg. 
021920
 
Butts
D
Retinal waves: implications for synaptic learning rules during development
Neuroscientist
 , 
2002
, vol. 
8
 (pg. 
243
-
253
)
Buzás
P
Kovács
K
Ferecskó
A
Budd
J
Eysel
U
Kisvárday
Z
Model-based analysis of excitatory lateral connections in the visual cortex
J Comp Neurol
 , 
2006
, vol. 
499
 (pg. 
861
-
881
)
Callaway
E
Katz
L
Effects of binocular deprivation on the development of clustered horizontal connections in cat striate cortex
Proc Natl Acad Sci USA
 , 
1991
, vol. 
88
 (pg. 
745
-
749
)
Callaway
E
Katz
L
Emergence and refinement of clustered horizontal connections in cat striate cortex
J Neurosci
 , 
1990
, vol. 
10
 (pg. 
1134
-
1153
)
Cartwright
J
Labyrinthine Turing pattern formation in the cerebral cortex
J Theor Biol
 , 
2002
, vol. 
217
 (pg. 
97
-
103
)
Castets
V
Dulos
E
Boissonade
J
De Kepper
P
Experimental evidence of a sustained standing Turing-type nonequilibrium chemical pattern
Phys Rev Lett
 , 
1990
, vol. 
64
 (pg. 
2953
-
2956
)
Charron
F
Tessier-Lavigne
M
The Hedgehog, TGF-β/BMP and Wnt families of morphogens in axon guidance
Adv Exp Med Biol
 , 
2007
, vol. 
621
 (pg. 
116
-
133
)
Chen
Y
Schier
A
Lefty proteins are long-range inhibitors of squint-mediated nodal signaling
Curr Biol
 , 
2002
, vol. 
12
 (pg. 
2124
-
2128
)
Coogan
T
Van Essen
D
Development of connections within and between areas V1 and V2 of macaque monkeys
J Comp Neurol
 , 
1996
, vol. 
372
 (pg. 
327
-
342
)
da Costa
N
Martin
K
Whose cortical column would that be?
Front Neuroanat
 , 
2010
, vol. 
4
 pg. 
16
 
De Kepper
P
Castets
V
Dulos
E
Boissonade
J
Turing-type chemical patterns in the chlorite-iodide-malonic acid reaction
Physica D
 , 
1991
, vol. 
49
 (pg. 
161
-
169
)
Douglas
R
Martin
K
Neuronal circuits of the neocortex
Annu Rev Neurosci
 , 
2004
, vol. 
27
 (pg. 
419
-
451
)
Duboc
V
Rottinger
E
Besnardeau
L
Lepage
T
Nodal and BMP2/4 signaling organizes the oral-aboral axis of the sea urchin embryo
Dev Cell
 , 
2004
, vol. 
6
 (pg. 
397
-
410
)
Durack
J
Katz
L
Development of horizontal projections in layer 2/3 of ferret visual cortex
Cereb Cortex
 , 
1996
, vol. 
6
 (pg. 
178
-
183
)
Economou
A
Ohazama
A
Porntaveetus
T
Sharpe
P
Kondo
S
Basson
M
Gritli-Linde
A
Cobourne
M
Green
J
Periodic stripe formation by a Turing mechanism operating at growth zones in the mammalian palate
Nat Genet
 , 
2012
, vol. 
44
 (pg. 
348
-
351
)
Gabriel
K
Sokal
R
A new statistical approach to geographic variation analysis
Syst Biol
 , 
1969
, vol. 
18
 (pg. 
259
-
278
)
Galuske
R
Singer
W
The origin and topography of long-range intrinsic projections in cat visual cortex: a developmental study
Cereb Cortex
 , 
1996
, vol. 
6
 (pg. 
417
-
430
)
Garzón-Alvarado
D
Martinez
A
Segrera
D
A model of cerebral cortex formation during fetal development using reaction-diffusion-convection equations with Turing space parameters
Comput Mehods Programs Biomed
 , 
2011
, vol. 
104
 (pg. 
489
-
497
)
Gierer
A
Meinhardt
H
A theory of biological pattern formation
Biol Cybern
 , 
1972
, vol. 
12
 (pg. 
30
-
39
)
Gilbert
C
Wiesel
T
Columnar specificity of intrinsic horizontal and corticocortical connections in cat visual cortex
J Neurosci
 , 
1989
, vol. 
9
 (pg. 
2432
-
2442
)
Gordon
L
Mansh
M
Kinsman
H
Morris
A
Xenopus sonic hedgehog guides retinal axons along the optic tract
Dev Dyn
 , 
2010
, vol. 
239
 (pg. 
2921
-
2932
)
Grabska-Barwińska
A
von der Malsburg
C
Establishment of a scaffold for orientation maps in primary visual cortex of higher mammals
J Neurosci
 , 
2008
, vol. 
28
 (pg. 
249
-
257
)
Gray
P
Scott
S
Autocatalytic reactions in the isothermal, continuous stirred tank reactor: isolas and other forms of multistability
Chem Eng Sci
 , 
1983
, vol. 
38
 (pg. 
29
-
43
)
Grossberg
S
On the development of feature detectors in the visual cortex with applications to learning and reaction-diffusion systems
Biol Cybern
 , 
1976
, vol. 
21
 (pg. 
145
-
159
)
Grossberg
S
Olson
S
Rules for the cortical map of ocular dominance and orientation columns
Neural Netw
 , 
1994
, vol. 
7
 (pg. 
883
-
894
)
Grossberg
S
Seitz
A
Laminar development of receptive fields, maps and columns in visual cortex: the coordinating role of the subplate
Cereb Cortex
 , 
2003
, vol. 
13
 (pg. 
852
-
863
)
Hines
M
Carnevale
N
The NEURON simulation environment
Neural Comput
 , 
1997
, vol. 
9
 (pg. 
1179
-
1209
)
Hubel
D
Wiesel
T
Ferrier lecture: functional architecture of macaque monkey visual cortex
Proc R Soc Lond B Biol Sci
 , 
1977
, vol. 
198
 (pg. 
1
-
59
)
Hubel
D
Wiesel
T
Stryker
M
Anatomical demonstration of orientation columns in macaque monkey
J Comp Neurol
 , 
1978
, vol. 
177
 (pg. 
361
-
379
)
Huberman
A
Feller
M
Chapman
B
Mechanisms underlying development of visual maps and receptive fields
Annu Rev Neurosci
 , 
2008
, vol. 
31
 (pg. 
479
-
509
)
Iron
D
Wei
J
Winter
M
Stability analysis of Turing patterns generated by the Schnakenberg model
J Math Biol
 , 
2004
, vol. 
49
 (pg. 
358
-
390
)
Kanold
P
Transient microcircuits formed by subplate neurons and their role in functional development of thalamocortical connections
Neuroreport
 , 
2004
, vol. 
15
 (pg. 
2149
-
2153
)
Kennedy
H
Dehay
C
Cortical specification of mice and men
Cereb Cortex
 , 
1993
, vol. 
3
 (pg. 
171
-
186
)
Kisvarday
Z
Eysel
U
Cellular organization of reciprocal patchy networks in layer iii of cat visual cortex (area 17)
Neuroscience
 , 
1992
, vol. 
46
 (pg. 
275
-
286
)
Kondo
S
Asai
R
A reaction–diffusion wave on the skin of the marine angelfish Pomacanthus
Nature
 , 
1995
, vol. 
376
 (pg. 
765
-
768
)
Kornack
D
Rakic
P
Radial and horizontal deployment of clonally related cells in the primate neocortex: relationship to distinct mitotic lineages
Neuron
 , 
1995
, vol. 
15
 (pg. 
311
-
321
)
Lefèvre
J
Mangin
J
A reaction–diffusion model of human brain development
PLoS Comput Biol
 , 
2010
, vol. 
6
 pg. 
e1000749
 
Leppänen
T
Karttunen
M
Barrio
R
Kaski
K
The effect of noise on Turing patterns
Prog Theor Phys.
 , 
2003
, vol. 
150
 
Suppl.
(pg. 
367
-
370
)
Liets
L
Olshausen
B
Wang
G
Chalupa
L
Spontaneous activity of morphologically identified ganglion cells in the developing ferret retina
J Neurosci
 , 
2003
, vol. 
23
 (pg. 
7343
-
7350
)
Luhmann
H
Millan
L
Singer
W
Development of horizontal intrinsic connections in cat striate cortex
Exp Brain Res
 , 
1986
, vol. 
63
 (pg. 
443
-
448
)
Lund
J
Angelucci
A
Bressloff
P
Anatomical substrates for functional columns in macaque monkey primary visual cortex
Cereb Cortex
 , 
2003
, vol. 
13
 (pg. 
15
-
24
)
Lund
J
Yoshioka
T
Levitt
J
Comparison of intrinsic connectivity in different areas of macaque monkey cerebral cortex
Cereb Cortex
 , 
1993
, vol. 
3
 (pg. 
148
-
162
)
Luskin
M
Pearlman
A
Sanes
J
Cell lineage in the cerebral cortex of the mouse studied in vivo and in vitro with a recombinant retrovirus
Neuron
 , 
1988
, vol. 
1
 (pg. 
635
-
647
)
Mariño
J
Schummers
J
Lyon
D
Schwabe
L
Beck
O
Wiesing
P
Obermayer
K
Sur
M
Invariant computations in local cortical networks with balanced excitation and inhibition
Nat Neurosci
 , 
2005
, vol. 
8
 (pg. 
194
-
201
)
Matsubara
J
Cynader
M
Swindale
N
Anatomical properties and physiological correlates of the intrinsic connections in cat area 18
J Neurosci
 , 
1987
, vol. 
7
 (pg. 
1428
-
1446
)
Matsuoka
R
Chivatakarn
O
Badea
T
Samuels
I
Cahill
H
Katayama
K
Kumar
S
Suto
F
Chédotal
A
Peachey
N
Class 5 transmembrane semaphorins control selective mammalian retinal lamination and function
Neuron
 , 
2011
, vol. 
71
 (pg. 
460
-
473
)
McGough
J
Riley
K
Pattern formation in the gray–scott model
Nonlinear Anal
 , 
2004
, vol. 
5
 (pg. 
105
-
121
)
Meinhardt
H
Klingler
M
A model for pattern formation on the shells of molluscs
J Theor Biol
 , 
1987
, vol. 
126
 (pg. 
63
-
89
)
Meno
C
Gritsman
K
Ohishi
S
Ohfuji
Y
Heckscher
E
Mochida
K
Shimono
A
Kondoh
H
Talbot
W
Robertson
E
, et al.  . 
Mouse lefty2 and zebrafish antivin are feedback inhibitors of nodal signaling during vertebrate gastrulation
Mol Cell
 , 
1999
, vol. 
4
 (pg. 
287
-
298
)
Mitchison
G
Crick
F
Long axons within the striate cortex: their distribution, orientation, and patterns of connection
Proc Natl Acad Sci USA
 , 
1982
, vol. 
79
 (pg. 
3661
-
3665
)
Mou
C
Jackson
B
Schneider
P
Overbeek
P
Headon
D
Generation of the primary hair follicle pattern
Proc Natl Acad Sci USA
 , 
2006
, vol. 
103
 (pg. 
9075
-
9080
)
Mountcastle
V
Modality and topographic properties of single neurons of cat's somatic sensory cortex
J Neurophysiol
 , 
1957
, vol. 
20
 (pg. 
408
-
434
)
Muir
DR
Da Costa
NMa
Girardin
CC
Naaman
S
Omer
DB
Ruesch
E
Grinvald
A
Douglas
RJ
Embedding of cortical representations by the superficial patch system
Cereb Cortex
 , 
2011
, vol. 
21
 (pg. 
2244
-
2260
)
Muir
DR
Douglas
RJ
From neural arbors to daisies
Cereb Cortex
 , 
2011
, vol. 
21
 (pg. 
1118
-
1133
)
Müller
P
Rogers
K
Jordan
B
Lee
J
Robson
D
Ramanathan
S
Schier
A
Differential diffusivity of nodal and lefty underlies a reaction-diffusion patterning system
Sci Signal
 , 
2012
, vol. 
336
 (pg. 
721
-
724
)
Murphy
K
Duffy
K
Jones
D
Mitchell
D
Development of cytochrome oxidase blobs in visual cortex of normal and visually deprived cats
Cereb Cortex
 , 
2001
, vol. 
11
 (pg. 
122
-
135
)
Newman
S
Bhat
R
Activator-inhibitor dynamics of vertebrate limb pattern formation
Birth Defects Res Part C
 , 
2007
, vol. 
81
 (pg. 
305
-
319
)
Page
K
Maini
P
Monk
N
Complex pattern formation in reaction-diffusion systems with spatially varying parameters
Physica D
 , 
2005
, vol. 
202
 (pg. 
95
-
115
)
Pena
B
Perez-Garcia
C
Stability of Truing patterns in the Brusselator model
Phys Rev E
 , 
2001
, vol. 
64
 pg. 
056213
 
Plikus
M
Baker
R
Chen
C
Fare
C
de la Cruz
D
Andl
T
Maini
P
Millar
S
Widelitz
R
Chuong
C
Self-organizing and stochastic behaviors during the regeneration of hair stem cells
Sci STKE
 , 
2011
, vol. 
332
 (pg. 
586
-
589
)
Price
D
The postnatal development of clustered intrinsic connections in area 18 of the visual cortex in kittens
Dev Brain Res
 , 
1986
, vol. 
24
 (pg. 
31
-
38
)
Rakic
P
A small step for the cell, a giant leap for mankind: a hypothesis of neocortical expansion during evolution
Trends Neurosci
 , 
1995
, vol. 
18
 (pg. 
383
-
388
)
Rakic
P
Specification of cerebral cortical areas
Science
 , 
1988
, vol. 
241
 (pg. 
170
-
176
)
Rockland
K
Anatomical organization of primary visual cortex (area 17) in the ferret
J Comp Neurol
 , 
1985
, vol. 
241
 (pg. 
225
-
236
)
Rockland
K
Lund
J
Widespread periodic intrinsic connections in the tree shrew visual cortex
Science
 , 
1982
, vol. 
215
 (pg. 
1532
-
1534
)
Rockland
K
Lund
J
Humphrey
A
Anatomical banding of intrinsic connections in striate cortex of tree shrews (tupaia glis)
J Comp Neurol
 , 
1982
, vol. 
209
 (pg. 
41
-
58
)
Rockland
KS
Lund
JS
Intrinsic laminar lattice connections in primate visual cortex
J Comp Neurol
 , 
1983
, vol. 
216
 (pg. 
303
-
318
)
Ruesch
E
Functional architecture of superficial layer pyramidal neurons in the cat primary visual cortex
2011
ETH Zuerich
 
PhD dissertation
Ruthazer
E
Stryker
M
The role of activity in the development of long-range horizontal connections in area 17 of the ferret
J Neurosci
 , 
1996
, vol. 
16
 (pg. 
7253
-
7269
)
Sánchez-Camacho
C
Bovolenta
P
Emerging mechanisms in morphogen-mediated axon guidance
Bioessays
 , 
2009
, vol. 
31
 (pg. 
1013
-
1025
)
Schier
A
Nodal morphogens
Cold Spring Harb Perspect Biol
 , 
2009
, vol. 
1
 pg. 
a003459
 
Schmidt
K
Kim
D
Singer
W
Bonhoeffer
T
Löwel
S
Functional specificity of long-range intrinsic and interhemispheric connections in the visual cortex of strabismic cats
J Neurosci
 , 
1997
, vol. 
17
 (pg. 
5480
-
5492
)
Schnakenberg
J
Simple chemical reaction systems with limit cycle behaviour
J Theor Biol
 , 
1979
, vol. 
81
 (pg. 
389
-
400
)
Schnorrer
F
Dickson
B
Axon guidance: morphogens show the way
Curr Biol
 , 
2004
, vol. 
14
 (pg. 
R19
-
R21
)
Sharma
J
Angelucci
A
Rao
S
Sur
M
Relationship of intrinsic connections to orientation maps in ferret primary visual cortex: iso-orientation domains and singularities
1995
Presented at Society for Neuroscience
San Diego, CA
Shen
J
Jung
Y
Geometric and stochastic analysis of reaction-diffusion patterns
Int J Pure Appl Math
 , 
2005
, vol. 
19
 (pg. 
195
-
244
)
Sick
S
Reinker
S
Timmer
J
Schlake
T
WNT and DKK determine hair follicle spacing through a reaction-diffusion mechanism
Sci STKE
 , 
2006
, vol. 
314
 (pg. 
1447
-
1450
)
Siddiqui
M
Bhaumik
B
A reaction-diffusion model to capture disparity selectivity in primary visual cortex
PloS One
 , 
2011
, vol. 
6
 pg. 
e24997
 
Striegel
D
Hurdal
M
Chemically based mathematical model for development of cerebral cortical folding patterns
PLoS Comput Biol
 , 
2009
, vol. 
5
 pg. 
e1000524
 
Swindale
N
A model for the formation of ocular dominance stripes
Proc R Soc Lond B Biol Sci
 , 
1980
, vol. 
208
 (pg. 
243
-
264
)
Tanigawa
H
Wang
Q
Fujita
I
Organization of horizontal axons in the inferior temporal cortex and primary visual cortex of the macaque monkey
Cereb Cortex
 , 
2005
, vol. 
15
 (pg. 
1887
-
1899
)
Turing
A
The chemical basis of morphogenesis
Philos Trans R Soc Lond B Biol Sci
 , 
1952
, vol. 
237
 (pg. 
37
-
72
)
Uesaka
N
Ruthazer
E
Yamamoto
N
The role of neural activity in cortical axon branching
Neuroscientist
 , 
2006
, vol. 
12
 (pg. 
102
-
106
)
Van Hooser
S
Heimel
J
Chung
S
Nelson
S
Lack of patchy horizontal connectivity in primary visual cortex of a mammal without orientation maps
J Neurosci
 , 
2006
, vol. 
26
 (pg. 
7680
-
7692
)
Voges
N
Guijarro
C
Aertsen
A
Rotter
S
Models of cortical networks with long-range patchy projections
J Comput Neurosci
 , 
2010
, vol. 
28
 (pg. 
137
-
154
)
Webber
C
Hyakutake
M
McFarlane
S
Fibroblast growth factors redirect retinal axons in vitro and in vivo
Dev Biol
 , 
2003
, vol. 
263
 (pg. 
24
-
34
)
Weliky
M
Hall
M
Correlated neuronal activity minireview and visual cortical development
Neuron
 , 
2000
, vol. 
27
 (pg. 
427
-
430
)
Wong
R
Retinal waves and visual system development
Annu Rev Neurosci
 , 
1999
, vol. 
22
 (pg. 
29
-
47
)
Yousef
T
Tóth
É
Rausch
M
Eysel
U
Kisvárday
Z
Topography of orientation centre connections in the primary visual cortex of the cat
Neuroreport
 , 
2001
, vol. 
12
 (pg. 
1693
-
1699
)
Zou
Y
Lyuksyutova
A
Morphogens as conserved axon guidance cues
Curr Opin Neurobiol
 , 
2007
, vol. 
17
 (pg. 
22
-
28
)
Zubler
F
Douglas
R
A framework for modeling the growth and development of neurons and networks
Front Comput Neurosci
 , 
2009
, vol. 
3
 pg. 
25
 
Zubler
F
Hauri
A
Pfister
S
Whatley
A
Cook
M
Douglas
R
An instruction language for self-construction in the context of neural networks
Front Comput Neurosci
 , 
2011
, vol. 
5
 pg. 
57
 
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by-nc/3.0/), which permits non-commercial reuse, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com.