Modelling the emergence of whisker barrels
Abstract
Brain development relies on an interplay between genetic specification and selforganization. Striking examples of this relationship can be found in the somatosensory brainstem, thalamus, and cortex of rats and mice, where the arrangement of the facial whiskers is preserved in the arrangement of cell aggregates to form precise somatotopic maps. We show in simulation how realistic whisker maps can selforganize, by assuming that information is exchanged between adjacent cells only, under the guidance of gene expression gradients. The resulting model provides a simple account of how patterns of gene expression can constrain spontaneous pattern formation to faithfully reproduce functional maps in subsequent brain structures.
eLife digest
How does the brain wire itself up? One possibility is that a precise genetic blueprint tells every brain cell explicitly how it should be connected to other cells. Another option is that complex patterns emerge from relatively simple interactions between growing cells, which are more loosely controlled by genetic instruction.
The barrel cortex in the brains of rats and mice features one of the most distinctive wiring patterns. There, cylindrical clusters of cells – or barrels – are arranged in a pattern that closely matches the arrangement of the whiskers on the face. Neurons in a barrel become active when the corresponding whisker is stimulated. This precise mapping between individual whiskers and their brain counterparts makes the whiskerbarrel system ideal for studying brain wiring.
Guidance fields are a way the brain can create cell networks with wiring patterns like the barrels. In this case, genetic instructions help to create gradients of proteins across the brain. These help the axons that connect neurons together to grow in the right direction, by navigating towards regions of higher or lower concentrations. A large number of guidance fields could map out a set of centrepoint locations for axons to grow towards, ensuring the correct barrel arrangement. However, there are too few known guidance fields to explain how the barrel cortex could form by this kind of genetic instruction alone.
Here, James et al. tried to find a mechanism that could create the structure of the barrel cortex, relying only on two simple guidance fields. Indeed, two guidance fields should be enough to form a coordinate system on the surface of the cortex. In particular, it was examined whether the cortical barrel map could reliably selforganize without a full genetic blueprint prespecifying the barrel centrepoints in the cortex.
To do so, James et al. leveraged a mathematical model to create computer simulations; these showed that only two guidance fields are required to reproduce the map. However, this was only the case if axons related to different whiskers competed strongly for space while making connections, causing them to concentrate into whiskerspecific clusters. The simulations also revealed that the target tissue does not need to specify centrepoints if, instead, the origin tissue directs how strongly the axons should respond to the guidance fields. So this model describes a simple way that specific structures can be copied across the central nervous system.
Understanding the way the barrel cortex is set up could help to grasp how healthy brains develop, how brain development differs in certain neurodevelopmental disorders, and how brain wiring reorganizes itself in different contexts, for example after a stroke. Computational models also have the potential to reduce the amount of animal experimentation required to understand how brains are wired, and to cast light on how brain wiring is shaped by evolution.
Introduction
Spatial patterns in neural connectivity provide clues about the constraints under which brains evolve and develop (Purves et al., 1992). Perhaps the most distinctive pattern can be found in the barrel cortex of many rodent species (Woolsey and Van der Loos, 1970). The barrels are identifiable soon after birth in layer 4 of primary somatosensory cortex as dense clusters of thalamocortical axons, which are enclosed by borders a few neurons thick from postnatal day 3 (Erzurumlu and Gaspar, 2012). In the plane tangential to the cortical surface the barrels constitute a somatotopic map of the whiskers, with cells within adjacent barrels responding most strongly and quickly to deflection of adjacent whiskers (ArmstrongJames et al., 1992). Barrel patterning reflects subcortical whisker maps comprising cell aggregates called barrelettes in the brainstem and barreloids in the thalamus (Ma, 1991; Van Der Loos, 1976).
Barrel formation requires afferent input from whisker stimulation and thalamic calcium waves (AntónBolaños et al., 2019), and depends on a complex network of axon guidance molecules such as ephrinA5 and A7 and adhesion molecules such as cadherin6 and 8 (Vanderhaeghen et al., 2000; Miller et al., 2006). This network is orchestrated by interactions between morphogens Fgf8 and Fgf17 and transcription factors Emx2, Pax6, Sp8, and Couptf1 (Shimogori and Grove, 2005; Bishop et al., 2000), which are expressed in gradients spanning the cortical sheet that mark orthogonal axes and can be manipulated to stretch, shrink, shift, and even duplicate barrels (Assimacopoulos et al., 2012).
The barrel boundaries form a Voronoi tessellation (Senft and Woolsey, 1991; Figure 1A), suggesting that barreloid topology is preserved in the projection of thalamocortical axons into the cortex, and that a barrel forms by lateral axon branching from an initial centrepoint that ceases upon contact with axons branching from adjacent centres. However, the assumption of prearranged centrepoints is difficult to resolve with the observation that axons arrive in the cortical plate as an undifferentiated bundle, prior to barreloid formation (Agmon et al., 1993). In mice, axons from the trigeminal ganglion arrive in the principal division of the trigeminal nucleus (PrV) at E12, then axons from the PrV arrive in the ventroposteromedial nucleus of the thalamus (VPM) at E17, then axons from the VPM arrive in the cortical plate at E18/P0. Distinct whiskerrelated clusters then become apparent in the PrV at P0P1, in the VPM at P2P3, and in the cortex at P3P5 (Erzurumlu and Gaspar, 2012; Sehara and Kawasaki, 2011).
Alternatively, reactiondiffusion dynamics could generate a Voronoi tessellation without prearranged centres, by amplifying characteristic modes in a noisy initial distribution of axon branches, as a net effect of shortrange cooperative and longrange competitive interactions. Accordingly, the barrel pattern would be determined by the relative strength of these interactions and by the shape of the cortical field boundary. However, intrinsic cortical dynamics alone cannot account for the topographic correspondence between thalamic and cortical domains, the irregular sizes and specific arrangement of the barrels in rows and arcs, or the influence of gene expression gradients.
The centrepoint and reactiondiffusion models are not mutually exclusive. Preorganized centres could bias reactiondiffusion processes to generate specific arrangements more reliably, and mechanisms of lateral axon branching may constitute the tension between cooperation and competition required for selforganization. However, proof that barrel patterning can emerge from an undifferentiated bundle of axons, based only on local interactions, would show that a separate stage and/or extrinsic mechanism for preorganizing thalamocortical connections need not be assumed. To this end, we ask whether barrel maps can emerge in a system with reactiondiffusion dynamics, under the guidance of signalling gradients, and in the absence of predefined centres.
Models
Karbowski and Ermentrout, 2004 developed a reactiondiffusion style model of how extrinsic signalling gradients can constrain the emergence of distinct fields from intrinsic cortical dynamics. Their model defines how the fraction of occupied synapses ${c}_{i}(x,t)$ and the density of axon branches ${a}_{i}(x,t)$ interact at time t, along a 1D anteriorposterior axis x, for N thalamocortical projections indexed by i. The model was derived from the assumption that the rates at which a_{i} and c_{i} grow are reciprocally coupled. Extending the original 1D model to simulate arealization on a 2D cortical sheet, we use ${a}_{i}(\mathbf{\mathbf{x}},t)$ and ${c}_{i}(\mathbf{\mathbf{x}},t)$, and model synaptogenesis as
Accordingly, where the total fraction of synaptic connections sums to one, connections decay at rate α. Otherwise, ${c}_{i}(\mathbf{\mathbf{x}},t)$ increases nonlinearly ($k>1$) with the density of axon branching. Axon branching is modelled as
The first term on the right describes the divergence (indicated by $\mathrm{\nabla}\cdot$) of the quantity in parentheses, which is referred to as the ‘flux’ of axonal branching. The flux represents diffusion across the cortical sheet, at rate D, and the influence of M molecular signalling fields, $\rho (\mathbf{\mathbf{x}})$. The influence of a given field (indexed by j) on a given thalamic projection (indexed by i), is determined by ${\gamma}_{i,j}$, which may be positive or negative in order that axons may branch in the direction of either higher or lower concentrations. Note that computing the divergence in simulation requires cells on the cortical sheet to communicate with immediately adjacent cells only (see Methods). Here ${\chi}_{i}=0$ is a placeholder. The second term on the right represents the coupling between axon branching and synaptogenesis, and an assumption that the spatial distribution of synaptic density across the cortical sheet is broadly homogeneous. As such, the quantity c_{i} can be thought of as the connection density.
Results
First, we verified that all results established by Karbowski and Ermentrout, 2004 for a 1D axis could be reproduced using our extension to a 2D cortical sheet. Using an elliptical domain, S, with $M=3$ offset guidance gradients aligned to the longer axis, $N=5$ thalamocortical projections gave rise to five distinct cortical fields at locations that preserved the topographic ordering defined by the original γ values. However, we found that specifying N ordered areas required $M\approx (N+1)/2$ signalling fields. This is because localization of axon densities occurs only when projections are influenced by interactions with two or more signalling gradients that encourage migration in opposing directions. As the number of guidance fields is unlikely to approach the number of individual barrels, modifications to the model were required.
We reasoned that an arbitrary number of distinct field locations may be determined by a minimum of two guidance gradients, if the concentration of the projection densities is influenced by competition between projections, and if a projection that interacts more strongly with a given guidance gradient migrates further in the direction of that gradient. Accordingly, projections that interact most strongly with a given guidance gradient would come to occupy cortical locations at which that field has extreme values, leaving adjacent locations available to be occupied by projections with the next strongest interactions, and so forth. This would in principle allow the relative locations of the fields to be specified by the relative values of the interaction parameters, γ, and hence for a topological map in the cortex to be specified by a spatial ordering of the γ values at the level of the thalamus.
Such dynamics are quite unlike those described by classic chemospecificity models (Sperry, 1963), which essentially assume centrepoints by specifying conditions in the target tissue that instruct preidentified afferents to stop growing. Consider, for example, that when simulated in isolation from oneanother, all projections in the model described would simply migrate to the extrema of the cortical guidance fields.
Testing this reasoning required increasing the strength of the competition between simulated thalamocortical projections for cortical territory, by increasing the tendency for each projection to compete for cortical space in which to branch and make connections. The major modification required was thus to introduce into the model an additional source of competition between thalamic projections. The term in parentheses in Equation 1 represents competition between thalamocortical projections for a limited availability of cortical connections. To introduce competition also in terms of axon branching, whilst ensuring that a_{i} is conserved over time, we redefined
This term contributes to the flux of axonal branching as an additional source of diffusion, scaled by $\u03f5$, which reduces the branching density for a given projection where the branches of other projections are dense. Note that this operation is local to individual afferent projections.
In addition, the model we have outlined requires that molecular guidance gradients in the cortex are complemented by graded values of the interaction strengths, γ, at the level of the thalamus. While the precise mechanisms by which thalamic and cortical gradients interact during development have not been fully characterised, the presence of complementary thalamic and cortical guidance gradients has been well established experimentally. In particular, the EphA4 receptor and its ligand ephrinA5 are distributed in complementary gradients in the somatosensory thalamus and cortex (Vanderhaeghen et al., 2000; Miller et al., 2006). Cells originating in VPM express high levels of EphA receptors and project to the lateral part of S1, which expresses low levels of ephrinA5, and cells originating in the VPL express low levels of EphA receptors and project to the medial part of S1, which expresses high levels of ephrinA5 (see Gao et al., 1998; Dufour et al., 2003; Vanderhaeghen and Polleux, 2004; Speer and Chapman, 2005; Torii et al., 2013). We assume that such patterning arises because the relative strengths of interaction with guidance molecules (e.g., ephrinA5) in the cortex are correlated with the relative concentrations of complementary molecules (e.g., EphA4) in the thalamus, and thus with thalamic position along the axis to which their gradients are aligned.
For simplicity, the two simulated thalamic interaction gradients, as well as the two cortical guidance gradients, were initially chosen to be linear and orthogonal. Hence a given pair of γ values corresponds to the coordinate of a barreloid centre in the VPM. Coordinates, in a reference plane defined by the anteriorposterior and mediallateral axes, were estimated from Figure 5d of Haidarliu and Ahissar, 2001, and scaled such that $\gamma \in \pm 2$. Note that this scaling is arbitrary because according to the model the coordinates provide relative position information only.
A cortical boundary enclosing barrels for 41 macrovibrissae was traced from a cytochrome oxidase stain from Zheng et al., 2001 (using original data kindly supplied by the authors), and Equations 1–3 were solved for $N=41$ projections on the resulting domain, S, using $M=2$ linear signalling gradients aligned with the anteriorposterior and mediallateral axes. These gradients are shown with the barrel field boundary in Figure 1A for clarity, though like ephrinA5 they may be thought of as extending across the cortical hemisphere (Miller et al., 2006). Simulations were stepped through 30000 iterations of Equations 1–3 ($\delta t=0.0001$).
Across a wide range of parameter values, random initial conditions (a uniform random distribution for $a(\mathbf{\mathbf{x}},0)\in (0.2,0.4)$, $c(\mathbf{\mathbf{x}},0)=0$) eventually yielded a clear Voronoilike tessellation of topographically organized thalamocortical projections, confirming that barrel maps can selforganize in the absence of prespecified centre points. The organization is apparent in a plot of the identity of the projection for which the connection density is maximal at each simulated cortical location, as shown in Figure 1C. Parameter values for this example simulation (see also Figure 1—video 1) were obtained by conducting a full parameter sweep and choosing a combination ($\alpha =3.6$, $\beta =16.67$, $k=3$, $D=0.5$, $\u03f5=1.2$) that scored well against the following three measures.
First, we used an algorithm introduced by Honda to measure the discrepancy of each barrel shape from a Dirichlet domain shape (Honda, 1983). Low overall values of this Hondaδ metric obtained from simulated barrels indicate a close correspondence of the simulated barrel field with a Voronoi tessellation, and thus with a biological barrel field (for mice $\delta \approx 0.054$, Senft and Woolsey, 1991, and our analysis of data from Zheng et al., 2001 indicates that the value for rats is similar). For the tessellation that is overlaid on the real barrel field in Figure 1A, $\delta =0.025$, and a reduction in δ in the example simulation over time confirmed that an equivalent ‘good’ Voronoi pattern can emerge within ≈ 20000 iterations (Figure 1D, red circles). Second, we devised a pattern difference measure that is sensitive to deviations in the component shapes and overall topographic registration between two tessellations, η, and we used this measure to compare the simulated barrel fields to the real barrel field from which the boundary shape applied to the simulation was obtained (see Methods for details). A similar reduction in η in the development of the example simulation confirmed that the shapes and arrangement of emergent connection fields came to match those of the real barrel field by around 20000 iterations (Figure 1D, black squares). Third, we measured the connection selectivity, ω, at each location on the cortical sheet, as the connection density of the most dense projection divided by the sum over all projection densities. The overall connection selectivity increased as the barrel map selforganized in the example simulation (Figure 1D, grey hexagons), and the selectivity became concentrated in regions overlapping with the emergent barrel centres (Figure 1E).
Against these three metrics we are also able to characterise the robustness of selforganization to the model parameters, and to investigate the sensitivity of the model to variation in its inputs. Figure 2A shows values of δ, η and ω obtained after 30000 iterations, from 216 independent simulations, each representing a unique combination of the model parameters D, $\u03f5$, and the ratio $\alpha /\beta $. First, we observe that selforganization is highly robust to the ratio $\alpha /\beta $, across five orders of magnitude, with respect to all three metrics. Second, the most strongly Voronoiconforming patterns (low Hondaδ) were generated by simulations in which the diffusion constant D and the strength of competition $\u03f5$ were high. Third, strongest overall connection selectivities, ω, were obtained for lower values of D. Fourth, variation in the pattern difference metric, η, indicated that the alignment between real and simulated patterns was greatest for intermediate rates of diffusion, $D\approx 0.5$. Together these results indicate that when competition is strong, the rate of diffusion determines a tradeoff such that fields emerge to be barrelshaped when diffusion is fast and they emerge to be more selectively innervated when diffusion is slow.
The parameters of the example simulation are indicated in Figure 2A using an asterisk. In Figure 2B, we also present examples of alternative patterns that emerge for different choices of D. Decreasing the rate of diffusion may be considered equivalent to increasing the overall size of the domain, S. Hence, insights into barrel development in species with a larger representation of the vibrissae, which do not have barrel fields, may be gained by studying pattern formation when D is small. In this context, it is interesting to note that for small D, the organization is predicted to be topological but highly irregular, with a general expansion in the territory occupied by the central versus peripheral domains that would presumably manifest as an absence of identifiable barrel fields (Figure 2Biii).
Next we conducted a sensitivity analysis to determine the extent to which the quality of the pattern (after $t=30000$ iterations) is affected by perturbations to (i) the magnitude and offset of the noise applied to a_{i} at $t=0$; (ii) noise applied to the interaction parameters, ${\gamma}_{i,j}$; (iii) noise (at various length scales) applied to the guidance fields; and (iv) the magnitude and orientation of one cortical guidance field relative to the other (Figure 1A).
Using the parameters of the example simulation (Figure 1C) we established baseline mean and standard deviations from ten independent simulations with initial uniform random values for $a(\mathbf{\mathbf{x}},0)\in (0.2,0.4)$, to be $\delta =0.089\pm 0.004$, $\eta =0.2108\pm 0.002$, and $\omega =0.2165\pm 0.0001$. Repeating with the variation in the initial noise doubled ($a(\mathbf{\mathbf{x}},0)\in (0.1,0.5)$), or removed altogether ($a(\mathbf{\mathbf{x}},0)=0.3$), generated distributions of δ, η, and ω that were not statistically different, as established using paired twosample ttests. Adding noise to the interaction parameters (γ) affected neither the Hondaδ or the connection selectivity measures substantially (see Figure 3A), and an increase in the pattern difference reflected an increase in the occurrence of topological defects only when perturbations became so large as to cause the ordering of γ values from neighbouring thalamic sites to be switched (see example map Figure 3Bi). Adding noise to the cortical guidance field values, ${\rho}_{1}(\mathbf{\mathbf{x}})$ and ${\rho}_{2}(\mathbf{\mathbf{x}})$, disrupted pattern formation only for high levels of noise applied at short length scales, which manifested as nonstraight edges at the domain boundaries (Figure 3Bii). Varying the slope of one linear gradient ${\rho}_{1}(\mathbf{\mathbf{x}})$ while keeping that of the other constant caused elongation of the emergent domains along the corresponding axis (Figure 3Biii), while pattern formation was not strongly influenced by relaxing the assumption that the gradients of the two cortical guidance fields are orthogonal (Figure 3Biv). Overall, the sensitivity analysis revealed that selforganization of barrellike fields in the model is highly robust to a wide range of sources of perturbation.
To further investigate the interplay of genes intrinsic to the developing neocortex and extrinsic factors such as thalamocortical input, we simulated two well known experimental manipulations of barrel development. First, we simulated a seminal barrel duplication paradigm (Shimogori and Grove, 2005; Assimacopoulos et al., 2012) in which the growth factor Fgf8, which is normally expressed at the anterior end of the cortical subplate from around E9.5 (Crossley and Martin, 1995), is ectopically expressed (by electroporation) also at the posterior pole. We assume that this results in a mirror of the primary barrel cortex boundary along the rostrocaudal axis (Assimacopoulos et al., 2012) and a mirroring of the anteriorposterior guidance gradient ${\rho}_{1}$ at the border between them (Figure 4A). The result after 30000 iterations, and otherwise using the parameters of the example simulation, was two mirrorsymmetrical barrel fields comprising $2N$ barrels (Figure 4B), consistent with the outcome of the original experiments.
Finally, to investigate the response of the model to environmental manipulation, we simulated a whisker deprivation experiment. In a critical period comprising the first postnatal days, removal of the whiskers by electrocauterization, plucking, or trimming leads to observable changes in brain structures, including the barrel field (Jeanmonod et al., 1981). Amongst other changes, deprivation of individual whiskers leads to smaller barrels (Kossut, 1992). We simulated trimming of the individual whisker C3 during the critical period by reducing the competitiveness of the C3 thalamocortical projection, $\u03f5$. As a result, the corresponding field size was smaller (Figure 3C), and the size of the fields representing the neighbours of C3 increased in size. A reduction in area to 65%, comparable to that induced by Kossut, 1992, was obtained in simulation when $\u03f5$ was reduced to 86% of its default value (Figure 3D), and the C3 barrel disappeared altogether when $\u03f5$ was less than half of its original value.
Although individual whisker trimming reduces barrel size, if an entire row is trimmed, barrel sizes for the trimmed row are not obviously changed (Land and Simons, 1985). We investigated with a simulation in which we varied $\u03f5$ for all of the row C projections. Although the barrels which formed for row C did show some reduction in area (Figure 4D), this reduction was small when compared with the individually trimmed C3 simulation. If the effect of trimming any whisker is to reduce its $\u03f5$ (the competitiveness of its projection) to 86% of its original value, then the model predicts that the cortical barrels for the trimmed C row will retain 91% of their area on average.
Together with the results of simulated misexpression, the consistency of the simulated whisker trimming results with those of the original studies demonstrates how the model can be used to investigate the contribution of intrinsic and extrinsic factors to the development of cortical fields.
Discussion
The present results suggest that the key requirements for the emergence of realistic barrel patterning are (i) at each cortical location thalamocortical projections compete for a limited number of available synaptic connections (Equations 1–2), (ii) at each location the branching rate of a given projection is reduced by the density of other projections (Equation 3), and (iii) the branch density of each projection is conserved over time.
The emergence of barrels in simulation required competition between thalamic projections in terms of synaptic connectivity and also competition in terms of cortical space, as represented by χ, with an implicit requirement for a self/other identifier amongst projections. This latter form of competition may account for the absence of barrels in rodents with larger brains, such as capybara, for which competition for space is presumably weaker (Woolsey et al., 1975). Hence, irrespective of whether barrels are necessary for adaptive whisker function, the emergence of somatotopically ordered modular structures may be an inevitable consequence of local competition for cortical territory driven by input from an array of discrete sensory organs (Purves et al., 1992).
In reality, the Voronoi tessellation is extended by scores of smaller barrels alongside the Erow barrels, which represent the microvibrissae, and presumably form via the same competitive processes. Enforcing here the same boundary condition as used to represent the true edges of the barrel field was necessary to ensure the stability of the simulation, though we acknowledge that this region of the boundary was enforced primarily to keep the number of simulated projections, and hence the overall computational complexity of the simulation, manageable (simulating an extra projection introduces 13030 new dynamical variables).
It is important to emphasize that the formulation of the model is entirely local, insofar as simulation requires no information to be communicated from a given cortical grid cell to any but those immediately adjacent (via diffusion). Hence the simulations demonstrate how a selforganizing system, constrained by genetically specified guidance cues and by the shape of the cortical field boundary, can faithfully reproduce an arrangement of cell aggregates in one neural structure as a topographic map in another.
Moreover, the present results confirm that somatotopic map formation does not require the prespecification of centrepoints by as yet undetermined additional developmental mechanisms.
Methods
We concentrated on the representation of the fortyone macrovibrissae that constitute a given barrel field, because their thalamic and cortical correlates are easily identifiable and consistently located, excluding the five rhinal whiskers as their cortical representation is isolated from the main barrel field. We excluded the representation of the microvibrissae to limit the overall complexity of the simulations.
The cortical sheet was modelled as a two dimensional hexagonal lattice, which simplifies the computation of the 2D Laplacian. Within a boundary traced around the edge of a rat barrel field (Figure 1A) we set the hextohex distance d to 0.03 mm, which resulted in a lattice containing 6515 hexes for the simulations shown in Figure 1A,C and D and 12739 hexes for the Fgf8 misexpression study shown in Figure 3. Each hex contained 82 timedependent variables: 41 branching densities (a_{i}) and 41 connection densities (c_{i}). The rate of change of each of the timedependent variables (Equations 1 and 2) was computed using a fourthorder RungeKutta method.
The most involved part of this computation is to find the divergence of the flux of axonal branching, ${\mathbf{\mathbf{J}}}_{i}(\mathbf{\mathbf{x}},t)$, the term in parentheses in Equation 2:
where ${\widehat{a}}_{i}\equiv {\sum}_{j\ne i}^{N}{a}_{j}$. Note that the sum of the guidance gradients is timeindependent and define ${\mathbf{\mathbf{g}}}_{i}(\mathbf{\mathbf{x}})\equiv {\sum}_{j=1}^{M}{\gamma}_{i,j}\nabla {\rho}_{j}(\mathbf{\mathbf{x}})$. Because the divergence operator is distributive, Equation 4 can be expanded using vector calculus identities (dropping references to $\mathbf{\mathbf{x}}$ and t for clarity):
Applying the vector calculus product rule identity yields
which has five elements to compute: (i) $D\mathrm{\nabla}\cdot \mathrm{\nabla}{a}_{i}$ (the Laplacian of a_{i}); (ii) a timeindependent modulator of a_{i} (because $\mathrm{\nabla}\cdot {\mathbf{g}}_{i}$ is a timeindependent static field); (iii) the scalar product of the static vector field ${\mathbf{\mathbf{g}}}_{i}$ and the gradient of a_{i}; (iv) the Laplacian of ${\widehat{a}}_{i}$; and (v) a term involving the gradients of a_{i} and ${\widehat{a}}_{i}$. Each of the divergences can be simplified by means of Gauss’s Theorem following Lee et al., 2014.
The computation of the mean value of the Laplacian across one hexagon of area $\mathrm{\Omega}=\frac{\sqrt{3}}{2}{d}^{2}$, located at position ${\mathbf{\mathbf{p}}}_{0}$, with neighbours at positions ${\mathbf{\mathbf{p}}}_{1}$–${\mathbf{\mathbf{p}}}_{6}$ is
(7) $\begin{array}{ll}{\displaystyle \u27e8D\mathrm{\nabla}\cdot \mathrm{\nabla}{a}_{i}({\mathbf{p}}_{0},t)\u27e9=\frac{D}{\mathrm{\Omega}}{\u222f}_{\mathrm{\Omega}}\mathrm{\nabla}\cdot \mathrm{\nabla}{a}_{i}(\mathbf{x},t)\phantom{\rule{negativethinmathspace}{0ex}}\mathrm{d}\mathrm{\Omega}}& {\displaystyle =\frac{D}{\mathrm{\Omega}}\oint \frac{\mathrm{\partial}{a}_{i}}{\mathrm{\partial}\hat{\mathbf{n}}}\phantom{\rule{negativethinmathspace}{0ex}}\mathrm{d}\gamma}\\ & {\displaystyle \approx \frac{D}{\mathrm{\Omega}}\sum _{j=1}^{6}\frac{\mathrm{\partial}{a}_{i}({\mathbf{p}}_{j})}{\mathrm{\partial}\hat{\mathbf{n}}}{}_{\mathrm{m}\mathrm{i}\mathrm{d}}v}\\ & {\displaystyle =\frac{2D}{\sqrt{3}{d}^{2}}\sum _{j=1}^{6}\frac{{a}_{i}({\mathbf{p}}_{j}){a}_{i}({\mathbf{p}}_{0})}{d}\frac{d}{\sqrt{3}}}\\ & {\displaystyle =\frac{2D}{3{d}^{2}}\sum _{j=1}^{6}({a}_{i}({\mathbf{p}}_{j}){a}_{i}({\mathbf{p}}_{0})),}\end{array}$where $v=d/\sqrt{3}$ is the length of each edge of the hexagon and $\mathrm{d}\gamma $ is an infinitesimally small distance along its perimeter.
The computation of the second term in Equation 6, $\u27e8{a}_{i}({\mathbf{p}}_{0},t)\mathrm{\nabla}\cdot {\mathbf{g}}_{i}({\mathbf{p}}_{0})\u27e9$, can be written out similarly:
(8) $\begin{array}{ll}{\displaystyle \frac{1}{\mathrm{\Omega}}{\u222f}_{\mathrm{\Omega}}{a}_{i}\mathrm{\nabla}\cdot {\mathbf{g}}_{i}\phantom{\rule{thickmathspace}{0ex}}\phantom{\rule{negativethinmathspace}{0ex}}\mathrm{d}\mathrm{\Omega}}& {\displaystyle =\frac{{a}_{i}({\mathbf{p}}_{0},t)}{\mathrm{\Omega}}\oint {\mathbf{g}}_{i}\cdot \phantom{\rule{negativethinmathspace}{0ex}}\mathrm{d}\hat{\mathbf{n}}}\\ & {\displaystyle \approx \frac{{a}_{i}({\mathbf{p}}_{0},t)}{\mathrm{\Omega}}\sum _{j=1}^{6}\frac{{\mathbf{g}}_{i}({\mathbf{p}}_{j})+{\mathbf{g}}_{i}({\mathbf{p}}_{0})}{2}\cdot \hat{\mathbf{n}}\phantom{\rule{thickmathspace}{0ex}}v}\\ & {\displaystyle =\frac{2{a}_{i}({\mathbf{p}}_{0},t)v}{\sqrt{3}{d}^{2}}\sum _{j=1}^{6}[\frac{{g}_{i}^{x}({\mathbf{p}}_{j})+{g}_{i}^{x}({\mathbf{p}}_{0})}{2}\cdot \hat{\mathbf{n}}+\frac{{g}_{i}^{y}({\mathbf{p}}_{j})+{g}_{i}^{y}({\mathbf{p}}_{0})}{2}\cdot \hat{\mathbf{n}}]}\\ {\displaystyle \Rightarrow \u27e8{a}_{i}({\mathbf{p}}_{0},t)\mathrm{\nabla}\cdot \mathbf{g}({\mathbf{p}}_{0})\u27e9}& {\displaystyle \approx \frac{{a}_{i}({\mathbf{p}}_{0},t)}{3d}\sum _{j=1}^{6}[({g}_{i}^{x}({\mathbf{p}}_{j})+{g}_{i}^{x}({\mathbf{p}}_{0}))\mathrm{cos}(\frac{\pi}{3}(j1))+({g}_{i}^{y}({\mathbf{p}}_{j})+{g}_{i}^{y}({\mathbf{p}}_{0}))\mathrm{sin}(\frac{\pi}{3}(j1))],}\end{array}$where ${g}_{i}^{x}$ and ${g}_{i}^{y}$ are the Cartesian components of ${\mathbf{\mathbf{g}}}_{i}$. Both this last expression, and the final expression of Equation 7 can be computed locally, by summing over values of the nearest neighbours.
The middle term in Equation 6 is the scalar product of two vector fields which is straightforward to compute from their Cartesian components.
The same method used to compute $\mathrm{\nabla}\cdot \mathrm{\nabla}{a}_{i}$ in term (i) is used to compute $\mathrm{\nabla}\cdot \mathrm{\nabla}{\hat{a}}_{i}$.
The final term is the scalar product of the two vector fields $\nabla {a}_{i}$ and $\nabla {\widehat{a}}_{i}$.
By separating the computation of Equation 4 into parts (i)–(v), the noflux boundary condition,
can be fulfilled. On the boundary, the contribution to $\mathbf{\mathbf{J}}$ resulting from the first term of Equation 6 can be fixed to 0 by the ‘ghost cell method’ in which, during the evaluation of (i), a hex outside the boundary containing the same value as the hex inside the boundary is imagined to exist such that the flux of $\mathbf{\mathbf{J}}$ across the boundary is 0. Then, ${\mathbf{\mathbf{g}}}_{i}(\mathbf{\mathbf{x}})$ can be tailored so that it, and its normal derivative, approach 0 at the boundary, ensuring that the second and third terms of Equation 6 also contribute nothing to $\mathbf{\mathbf{J}}$. This is achieved by multiplying ${\mathbf{\mathbf{g}}}_{i}(\mathbf{\mathbf{x}})$ by a sharp logistic function of the distance, d_{b}, from $\mathbf{\mathbf{x}}$ to the boundary, of the form $1/[1+\mathrm{exp}(100({d}_{f}{d}_{b}))]$, where ${d}_{f}=0.1$ mm $\approx 3d$ is the boundary falloff distance.
The pattern difference metric, η, incorporates information about the differences in the areas of simulated and experimentally determined (real) barrels (${\mathcal{A}}_{i}^{\mathrm{sim}}$ and ${\mathcal{A}}_{i}^{\mathrm{exp}}$), as well as information contained in the ‘adjacency vector’ for each barrel, ${\mathcal{V}}_{i}$, the jth element of which is the length of the border between barrel i and j. For a wellformed barrel pattern, ${\mathcal{V}}_{i}$ is a sparse vector. A dimensionless quantity can be obtained from the scalar product of the simulated and experimental adjacency vectors: $\frac{1}{N}\sum _{i}\frac{{\mathcal{\mathcal{V}}}_{i}^{\mathrm{s}\mathrm{i}\mathrm{m}}}{{b}_{i}^{\mathrm{s}\mathrm{i}\mathrm{m}}}\text{}\cdot \text{}\frac{{\mathcal{\mathcal{V}}}_{i}^{\mathrm{e}\mathrm{x}\mathrm{p}}}{{b}_{i}^{\mathrm{e}\mathrm{x}\mathrm{p}}}$, where for example ${b}_{i}^{\mathrm{exp}}$, is the total length of the border around real barrel i. This quantity is small when the fields that form in simulation have dissimilar neighbour relations to those of the real barrels (e.g., Figure 1C, $t=1000$), and maximal for a precise topological map ($t=10000$). A second comparison considers the mean magnitude of the difference between the simulated and experimental vectors: $\frac{1}{N}{\sum}_{i}\parallel {\mathcal{V}}_{i}^{\mathrm{sim}}{\mathcal{V}}_{i}^{\mathrm{exp}}\parallel $. This tends to 0 for a perfect match and can separate patterns with straight boundaries from those with ‘noisy edges’ (as in Figure 3Bii). We combined these terms into a single metric:
which has units of mm^{3}.
The code required to reproduce these results is available at https://github.com/ABRGModels/BarrelEmerge/tree/eLife (James, 2020, copy archived at https://github.com/elifesciencespublications/BarrelEmerge). The computations described in (i)–(v) may be found in the class method RD_James_comp2::compute_divJ() which calculates term1, term2, term3, term1_1 and term1_2, respectively. BarrelEmerge depends on the software library morphologica (RRID:SCR_018813), which must also be compiled.
Data availability
All data and code required to reproduce the results in this work are available in the public repository https://github.com/ABRGModels/BarrelEmerge/tree/eLife (copy archived at https://github.com/elifesciencespublications/BarrelEmerge).
References

Flow of excitation within rat barrel cortex on striking a single vibrissaJournal of Neurophysiology 68:1345–1358.https://doi.org/10.1152/jn.1992.68.4.1345

Fibroblast growth factor 8 organizes the neocortical area map and regulates sensory map topographyJournal of Neuroscience 32:7191–7201.https://doi.org/10.1523/JNEUROSCI.007112.2012

The mouse Fgf8 gene encodes a family of polypeptides and is expressed in regions that direct outgrowth and patterning in the developing embryoDevelopment 121:439–451.

Development and critical period plasticity of the barrel cortex: barrel cortex plasticityEuropean Journal of Neuroscience 35:1540–1553.https://doi.org/10.1111/j.14609568.2012.08075.x

Size gradients of barreloids in the rat thalamusThe Journal of Comparative Neurology 429:372–387.https://doi.org/10.1002/10969861(20010115)429:3<372::AIDCNE2>3.0.CO;23

Geometrical models for cells in tissuesInternational Review of Cytology 81:191–248.https://doi.org/10.1016/s00747696(08)623396

Model of the early development of thalamocortical connections and area patterning via signaling moleculesJournal of Computational Neuroscience 17:347–363.https://doi.org/10.1023/B:JCNS.0000044876.28268.18

Effects of sensory deprivation upon a single cortical vibrissal column: a 2dg studyExperimental Brain Research 90:639–642.https://doi.org/10.1007/BF00230950

Hexagonal grid methods with applications to partial differential equationsInternational Journal of Computer Mathematics 91:1986–2009.https://doi.org/10.1080/00207160.2013.864392

The barrelettesarchitectonic vibrissal representations in the brainstem trigeminal complex of the mouse. I. normal structural organizationThe Journal of Comparative Neurology 309:161–199.https://doi.org/10.1002/cne.903090202

EphA7ephrinA5 signaling in mouse somatosensory cortex: developmental restriction of molecular domains and postnatal maintenance of functional compartmentsThe Journal of Comparative Neurology 496:627–642.https://doi.org/10.1002/cne.20926

Iterated patterns of brain circuitry (or how the cortex gets its spots)Trends in Neurosciences 15:362–368.https://doi.org/10.1016/01662236(92)90180G

Neuronal circuits with whiskerrelated patternsMolecular Neurobiology 43:155–162.https://doi.org/10.1007/s1203501181708

Mouse barrel cortex viewed as Dirichlet domainsCerebral Cortex 1:348–363.https://doi.org/10.1093/cercor/1.4.348

Fibroblast growth factor 8 regulates neocortical guidance of areaspecific thalamic innervationJournal of Neuroscience 25:6550–6560.https://doi.org/10.1523/JNEUROSCI.045305.2005

Grading the thalamus: how can an 'Eph' be excellent?Thalamus and Related Systems 3:235.https://doi.org/10.1017/S1472928807000234

Role of EphA/ephrina signaling in the development of topographic maps in mouse corticothalamic projectionsJournal of Comparative Neurology 521:626–637.https://doi.org/10.1002/cne.23195

Barreloids in mouse somatosensory thalamusNeuroscience Letters 2:1–6.https://doi.org/10.1016/03043940(76)900367

A mapping label required for normal scale of body representation in the cortexNature Neuroscience 3:358–365.https://doi.org/10.1038/73929

Comparative anatomical studies of the SmL face cortex with special reference to the occurrence of "barrels" in layer IVThe Journal of Comparative Neurology 164:79–94.https://doi.org/10.1002/cne.901640107
Decision letter

Upinder Singh BhallaReviewing Editor; Tata Institute of Fundamental Research, India

Timothy E BehrensSenior Editor; University of Oxford, United Kingdom

Bard ErmentroutReviewer; University of Pittsburgh, United States
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Acceptance summary:
The study proposes a compact selforganization model for formation of whisker barrels in the absence of predefined centers for the barrels. Further, the authors have come up with a biologically plausible rule for competition and conservation of axonal density.
Decision letter after peer review:
Thank you for submitting your article "Modelling the emergence of whisker barrels" for consideration by eLife. Your article has been reviewed by three peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Timothy Behrens as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Bard Ermentrout (Reviewer #2).
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
As the editors have judged that your manuscript is of interest, but as described below that additional analyses are required before it is published, we would like to draw your attention to changes in our revision policy that we have made in response to COVID19 (https://elifesciences.org/articles/57162). First, because many researchers have temporarily lost access to the labs, we will give authors as much time as they need to submit revised manuscripts. We are also offering, if you choose, to post the manuscript to bioRxiv (if it is not already there) along with this decision letter and a formal designation that the manuscript is 'in revision at eLife'. Please let us know if you would like to pursue this option. (If your work is more suitable for medRxiv, you will need to post the preprint yourself, as the mechanisms for us to do so are still in development.)
Summary:
This compact paper proposes a selforganization model for formation of whisker barrels. The key idea is that reactiondiffusion dynamics can lead to the observed topology, in the absence of predefined centers for the barrels.
Essential revisions:
1) How do the authors obtain 41 pairs of γ values (Results, third paragraph)? Are these parameters or were they inferred from experiments? This must be better motivated.
2) The competition term χ_{i} requires renormalization, which seems biologically implausible. The authors may wish to try a form such as
which does not need renormalization. Several other points about this competition term are unclear as mentioned in the reviewer comments.
3) There should be more exploration of the model: some parameter exploration and sensitivity analysis, and some more predictions.
Reviewer #1:
This compact paper proposes a selforganization model for formation of whisker barrels. The key idea is that reactiondiffusion dynamics can lead to the observed topology, in the absence of predefined centers for the barrels.
The model is well presented and the motivation of mathematical choices is mostly clear. It may be worth expanding on the motivation for competition for axonal branching (Equations 3 and 4).
It is a little unclear how the misexpression experiment (Shimogori and Grove, 2005) in Figure 1E was done. The simulation approach and outcome for this section is described very tersely.
The authors also mention another easily modeled experiment, in which capybara brains lack barrels because they are big. It should be a simple matter to do this run.
Overall I feel this study presents an attractive and compact model for the formation of whisker barrels, which has good biological motivation, and does a good job of reducing assumptions and molecular guidance cues.
Reviewer #2:
This is an interesting paper that with a few assumptions shows that an old model for areal formation in cortex is sufficient to quantitatively reproduce the patterns of barrels observed in mouse S1. It would appear from the model that the key is the parameters γ_{ij} that are presumably (hypothesized) to be assigned at the level of the thalamus. I have a few questions about the paper.
1) Does the same model work with respect to projections from the brain stem (barrelettes) to the thalamus (barreloids)? This would be a good way to check the ideas. Related to this, is it true that the barrelettes (barreloids) precede the development of the barreloids (barrels)? It would seem to be necessary? Or perhaps, starting with a double gradient in the thalamus and cortex and a prepattern in the barelettes, would the correct patterns emerge simultaneously?
2) There seems to be a strong prediction in this concerning the development of the patterns over time. Figure 1C indicates that early on there are large distortions in the shape of the barrels particularly in Figure 1D, E rows. is this known to occur?
3) It seems to me that without the chi, then possible connections plus axons are conserved which is reasonable. But with the necessary competition, there seems to be a flaw in the model if they have to renormalize at each point. If axons make connections should they not be lost from the pool forever (this is the dci/dt the model). For example, since the gradient has no flux in the original K&E model, there is conservation of the total number of connections and axons of a given type (int ai+ci dx = constant). This principle seems to make sense to me. However, the competition term χ_{i} seems to disrupt this. Is there a way to introduce the axonal competition in a way the prevents the unrealistic (or biologically implausible, at least) renormalization at each step? I'd be more comfortable with the model if there were a more physiological way to renormalize. For example, I don’t know if the authors considered something like an additional flux of the form:
This makes the axons of type i move away from type j while at the same time enforcing conservation without recourse to some sort of post normalization.
Reviewer #3:
The manuscript studies a theoretical model within the framework of reactiondiffusion equations coupled to signalling gradients to possibly explain the emergence of whisker barrels in the cortex.
1) The model considered by the authors is identical to the one studied by Karbowski and Ermentrout, 2004. The only new features are the extension of the original 1D model to 2D and the addition of an extra term to represent competition in axonal branching.
2) The authors consider 2 guiding fields. What are their explicit spatial profiles? Notice that since these fields essentially guide the emergent pattern and hence their profiles, in relation to the geometry of the 2D domain, are crucial. A different profile would certainly lead to a different pattern. I feel that it is not enough to say '…linear signalling gradients aligned with the anteriorposterior and mediallateral axes….' since the domain is 2D and of nonrectangular shape.
3) The justification for the introduction of the extra term for competition amongst axons (Equation 3) is missing. Why that form? What is the reasoning for introducing axonal competition? What essential features of the resultant patterns are missed out if this term is absent? Or has a different form? In the Discussion section, the authors mention, without any justification, that the conservation of branch density in each projection is a key requirement for the emergence of barrel patterns. This is totally unclear.
4) Related to the above point, the authors mention that the axonal branch density is bounded by their dynamics. I presume that the integrations on the RHS of Equation 4 are spatial integrals over the domain. Then how come a spatial index survives in the LHS of this equation? How did the authors arrive at this equation? Is there a continuoustime version of this equation (like a conservation law), i.e., one that does not make a reference to the discrete timestepping dynamics?
5) A typical mathematical modelling study should explore the space of relevant parameters to demonstrate the possible range of behaviours that the model can exhibit. This is usually presented as a phasediagram. The authors do not explore the parameter space (or the possible spatial profiles of the guiding fields) in their study.
6) Throughout, the authors emphasize the spatiallocality of their mathematical model and conclude 'Hence the simulations demonstrate how a selforganizing system…'. A mathematical model with spatiallocality alone does not imply selforganized dynamics. With a sufficiently large number of spatiotemporal fields (N=42), and the concomitant parameters, and nonautonomous guiding fields, it is possible to reproduce any desired pattern. As such, it is crucial in the mathematical modelling of living systems to delineate the essential requirements from the incidental.
https://doi.org/10.7554/eLife.55588.sa1Author response
Essential revisions:
1) How do the authors obtain 41 pairs of γ values (Results, third paragraph)? Are these parameters or were they inferred from experiments? This must be better motivated.
We have addressed this comment through edits at several points in the manuscript. At the beginning of the revised Materials and methods section we have clarified why 41 pairs were chosen:
“We concentrated on the representation of the fortyone macrovibrissae that constitute a given barrel field, because their thalamic and cortical correlates are easily identifiable and consistently located, excluding the five rhinal whiskers as their cortical representation is isolated from the main barrel field. We excluded the representation of the microvibrissae to limit the overall complexity of the simulations.”
To the revised Introduction we have added the following motivation:
“We reasoned that an arbitrary number of distinct field locations may be determined by a minimum of two guidance gradients, if the concentration of the projection densities is influenced by competition between projections, and if a projection that interacts more strongly with a given guidance gradient migrates further in the direction of that gradient. […] Consider, for example, that when simulated in isolation from oneanother, all projections in the model described would simply migrate to the extrema of the cortical guidance fields.”
Shortly thereafter we have added the following, which identifies the experimental observations on which our inferences are based:
“In addition, the model we have outlined requires that molecular guidance gradients in the cortex are complemented by graded values of the interaction strengths, 𝛾, at the level of the thalamus. […] We assume that such patterning arises because the relative strengths of interaction with guidance molecules (e.g., ephrinA5) in the cortex are correlated with the relative concentrations of complementary molecules (e.g., EphA4) in the thalamus, and thus with thalamic position along the axis to which their gradients are aligned.”
Finally, in the revised Discussion, we have acknowledged the following:
“In reality, scores of smaller Dirichletform barrels representing the microvibrissae form alongside the Erow barrels, presumably via the same competitive processes. Enforcing here the same boundary condition as used to represent the true edges of the barrel field was necessary to ensure the stability of the simulation, though we acknowledge that this region of the boundary was enforced primarily to keep the number of simulated projections, and hence the overall computational complexity of the simulation, manageable (simulating an extra projection introduces 13030 new dynamical variables).”
Together, we feel that these additions better explain and motivate our assumptions about the interaction parameters in a transparent way.
2) The competition term χ_{i} requires renormalization, which seems biologically implausible. The authors may wish to try a form such as
$X{a}_{i}\mathrm{\nabla}\frac{1}{N1}\sum _{j}\ne i{a}_{j}$which does not need renormalization. Several other points about this competition term are unclear as mentioned in the reviewer comments.
This is a very insightful comment. Addressing it led us to make important modifications to the model, and substantial improvements to the paper. We had previously attempted to use the form for the competition term that is suggested here, but without success. However, to address point 3 below we have since developed new software tools to facilitate a largescale parameter sweep (and sensitivity analysis), and when we carried out that full parameter sweep using the form for the competition term suggested, we discovered that during our original investigations we had been exploring only a small region of the parameter space in which the dynamics happened to be unstable. Now, in the region of the parameter space in which simulations are numerically stable, we have found pattern formation via the reformulated model to be qualitatively similar and highly robust. We too prefer the suggested reformulation, which removes the need for an explicit renormalization step (original Equation 4), and in the revised paper all of the results are now obtained using the more elegant form of the model (revised Equation 3). Note that the resulting expression is identical to that suggested, but that it must be included as part of the flux of axonal branching term, hence the placeholder variable 𝜒 is included within the parentheses in Equation 2 of the revision (previously it appeared outside). Note also that concomitant edits have been made to the Materials and methods section to reflect minor changes to the method of numerical integration that was required to realize the revised form of the competition.
3) There should be more exploration of the model: some parameter exploration and sensitivity analysis, and some more predictions.
The revision now includes a report of the results of fully exploring the parameter space of the model, as well as a sensitivity analysis. In addition to the Hondaẟ measure that we used in the original submission to measure the extent to which emergent fields were Dirichletform, we have in the Revision devised additional metrics of the pattern quality that allow us to better quantify the correspondence between real and simulated barrels as the parameters are varied (in terms of barrel area and map topology; see final section of Materials and methods). Using these metrics, we analyse the results of 216 independent simulations, presented in the new Figure 2, as well as 200 new simulations used in the new Figure 3 to present an analysis of the sensitivity of the model to perturbations to the following inputs: (i) the magnitude and offset of the noise applied to a_{i} at t = 0, (ii) random noise applied to the interaction parameters, γ, (iii) random noise applied to the guidance fields, and (iv) the angle and magnitude of one of the guidance fields.
Based on these new analyses we arrive at the following prediction about how the tension of competition between thalamocortical axons and their lateral branching (diffusion) manifests in barrel organization:
“First we observe that selforganization is highly robust to the ratio ɑ/𝛽, across five orders of magnitude, with respect to all three metrics. […] Together these results indicate that when competition is strong, the rate of diffusion determines a tradeoff such that fields emerge to be barrelshaped when diffusion is fast and they emerge to be more selectively innervated when diffusion is slow.”
We now include several examples of patterns that are predicted to emerge in different regions of the parameter space. Specifically we confirm in simulation an intuition about the influence of cortical sheet size that we articulated only informally in the original submission (see new Figure 2B; specifically plot iii):
“The parameters of the example simulation are indicated in Figure 2A using an asterisk. […] In this context, it is interesting to note that for small D, the organization is predicted to be topological but highly irregular, with a general expansion in the territory occupied by the central versus peripheral domains that would presumably manifest as an absence of identifiable barrel fields (Figure 2B iii).”
We also demonstrate in the revision how the model can be used to study the effects of specific experimental manipulations that have been conducted with the whiskerbarrel system, by replicating two whisker deprivation experiments:
“Finally, to investigate the response of the model to environmental manipulation, we simulated a whisker deprivation experiment. […] Together with the results of simulated misexpression, the consistency of the simulated whisker trimming results with those of the original studies demonstrates how the model can be used to investigate the contribution of intrinsic and extrinsic factors to the development of cortical fields.”
Reviewer #1:
This compact paper proposes a selforganization model for formation of whisker barrels. The key idea is that reactiondiffusion dynamics can lead to the observed topology, in the absence of predefined centers for the barrels.
The model is well presented and the motivation of mathematical choices is mostly clear. It may be worth expanding on the motivation for competition for axonal branching (Equations 3 and 4).
Thank you. We have expanded on the motivation for competition for axonal branching as outlined in our responses to point 1 above.
It is a little unclear how the misexpression experiment (Shimogori and Grove, 2005) in Figure 1E was done. The simulation approach and outcome for this section is described very tersely.
We have expanded our description of this experiment in the main text, as follows:
“To further investigate the interplay of genes intrinsic to the developing neocortex and extrinsic factors such as thalamocortical input, we simulated two wellknown experimental manipulations of barrel development. […] The result after 30000 iterations, and otherwise using the parameters of the example simulation, was two mirrorsymmetrical barrel fields comprising 2𝑁 barrels (Figure 4B), consistent with the outcome of the original experiments.”
The authors also mention another easily modeled experiment, in which capybara brains lack barrels because they are big. It should be a simple matter to do this run.
Thank you for the suggestion. As explained in response to point 3 above, we have run this simulation as part of the parameter sweep, and the results (as anticipated in the original submission) are shown in the new Figure 2Biii.
Overall I feel this study presents an attractive and compact model for the formation of whisker barrels, which has good biological motivation, and does a good job of reducing assumptions and molecular guidance cues.
Reviewer #2:
This is an interesting paper that with a few assumptions shows that an old model for areal formation in cortex is sufficient to quantitatively reproduce the patterns of barrels observed in mouse S1. It would appear from the model that the key is the parameters γ_{ij} that are presumably (hypothesized) to be assigned at the level of the thalamus. I have a few questions about the paper.
1) Does the same model work with respect to projections from the brain stem (barrelettes) to the thalamus (barreloids)? This would be a good way to check the ideas. Related to this, is it true that the barrelettes (barreloids) precede the development of the barreloids (barrels)? It would seem to be necessary? Or perhaps, starting with a double gradient in the thalamus and cortex and a prepattern in the barelettes, would the correct patterns emerge simultaneously?
To answer the reviewers questions directly – yes, our thinking is that this model could be applicable for thalamic pattern formation, though we note that there is a more intricate 3D (peppershaped) structure to the thalamic barreloids that would also need to be considered in an extended version of the model, and there is less data available that describes these structures insofar as would be necessary to validate such a model. Yes, the formation of discrete whiskerrelated structures (or at least their visibility using current experimental techniques) precedes from brainstem barrelettes to thalamic barreloid to cortical barrels, though the (undifferentiated) projections from brainstem to thalamus to cortex are in place just before the barelettes become apparent. We are keen in future work to extend this modelling framework to investigate pattern formation along the neuraxis.
2) There seems to be a strong prediction in this concerning the development of the patterns over time. Figure 1C indicates that early on there are large distortions in the shape of the barrels particularly in Figure 1D, E rows. is this known to occur?
Care should be taken when relating the timeevolution of the simulation to the timeevolution of the development of barrels. The model considers axon branching density and connection density in the cortex without explicit consideration of which layers these quantities may be evolving within. It is known (e.g., from Agmon, 2005) that the axons initially form an undifferentiated mass in the subplate, and it appears to be here that they find their location before growing up, fairly straight, into layer IV where they then form the barrels. Thus, distortions in the pattern seen early on in the model development may be observable in the subplate, without necessarily being observable in the barrel field, and we are not aware of an appropriate dataset for subplate pattern formation that would be applicable.
3) It seems to me that without the chi, then possible connections plus axons are conserved which is reasonable. But with the necessary competition, there seems to be a flaw in the model if they have to renormalize at each point. If axons make connections should they not be lost from the pool forever (this is the dci/dt the model). For example, since the gradient has no flux in the original K&E model, there is conservation of the total number of connections and axons of a given type. (int ai+ci dx = constant). This principle seems to make sense to me. However, the competition term χ_{i} seems to disrupt this. Is there a way to introduce the axonal competition in a way the prevents the unrealistic (or biologically implausible, at least) renormalization at each step? I'd be more comfortable with the model if there were a more physiological way to renormalize. For example, I don’t know if the authors considered something like an additional flux of the form:
$X{a}_{i}\mathrm{\nabla}\frac{1}{N1}\sum _{j}\ne i{a}_{j}$This makes the axons of type i move away from type j while at the same time enforcing conservation without recourse to some sort of post normalization.
We thank reviewer 2 for this very insightful comment. We have addressed this comment in our response to point 2 above. We now present the model with the competition expressed in the form suggested, and we certainly agree that the model is much more elegant in this form.
Reviewer #3:
The manuscript studies a theoretical model within the framework of reactiondiffusion equations coupled to signalling gradients to possibly explain the emergence of whisker barrels in the cortex.
1) The model considered by the authors is identical to the one studied by Karbowski and Ermentrout, 2004. The only new features are the extension of the original 1D model to 2D and the addition of an extra term to represent competition in axonal branching.
This is correct. Though note that in the original simulations of Karbowksi and Ermentrout, 2004, a given field becomes localised around a specific point that is specified by a unique combination of two (or more) guidance fields, whereas in our simulations we show that the locations of many fields can be reliably influenced by only two fields.
2) The authors consider 2 guiding fields. What are their explicit spatial profiles? Notice that since these fields essentially guide the emergent pattern and hence their profiles, in relation to the geometry of the 2D domain, are crucial. A different profile would certainly lead to a different pattern. I feel that it is not enough to say '…linear signalling gradients aligned with the anteriorposterior and mediallateral axes….' since the domain is 2D and of nonrectangular shape.
Minor edits to the Introduction text underscore that gradients in the guidance molecules are not restricted to the barrel field domain but extend across the entire cortical sheet (e.g., see Miller et al., 2006, Figure 11, for a clear illustration), and in the Results we have clarified:
“These gradients are shown with the barrel field boundary in Figure 1A for clarity, though like ephrin5 they may be thought of as extending across the cortical hemisphere (Miller et al., 2006)”.
Note that the new Figure 3, which accompanies the sensitivity analysis, includes an analysis of the sensitivity of the model to various perturbations to the profile of the chemical cues, and that the specific profile is clarified in the corresponding caption:
“The effect of setting the rotational angle of the linearly varying guidance field, 𝜌1 (x), to 𝜙𝜌1, and modifying its overall gain to 𝐺𝜌1, whilst keeping the parameters of 𝜌2(x) unchanged from those used in the example simulation, for which 𝜙𝜌2 = 84°and 𝐺𝜌2 = 1.”
Note also that in the revision we have included an analysis of how the relationship between the two guidance fields affects pattern formation (new Results section on sensitivity analysis and new Figure 3A).
3) The justification for the introduction of the extra term for competition amongst axons (Equation 3) is missing. Why that form? What is the reasoning for introducing axonal competition? What essential features of the resultant patterns are missed out if this term is absent? Or has a different form? In the Discussion section, the authors mention, without any justification, that the conservation of branch density in each projection is a key requirement for the emergence of barrel patterns. This is totally unclear.
This comment is addressed through our response to point 2 above.
4) Related to the above point, the authors mention that the axonal branch density is bounded by their dynamics. I presume that the integrations on the RHS of Equation 4 are spatial integrals over the domain. Then how come a spatial index survives in the LHS of this equation? How did the authors arrive at this equation? Is there a continuoustime version of this equation (like a conservation law), i.e., one that does not make a reference to the discrete timestepping dynamics?
We understand and appreciate the points raised here, as they continue on from the previous point. As explained in full in our response to point 2 above, the post normalization mechanism (original Equation 4) has indeed been replaced by a continuoustime version (revised Equation 3).
5) A typical mathematical modelling study should explore the space of relevant parameters to demonstrate the possible range of behaviours that the model can exhibit. This is usually presented as a phasediagram. The authors do not explore the parameter space (or the possible spatial profiles of the guiding fields) in their study.
Thank you for encouraging us to fully explore the parameter space. This comment is addressed in response to point 3 above. Briefly, we have included results from a full parameter sweep (new Figure 2), and also include as part of a sensitivity analysis (new Figure 3) an investigation of how the profiles of the guiding fields affect pattern formation.
6) Throughout, the authors emphasize the spatiallocality of their mathematical model and conclude 'Hence the simulations demonstrate how a selforganizing system…'. A mathematical model with spatiallocality alone does not imply selforganized dynamics. With a sufficiently large number of spatiotemporal fields (N=42), and the concomitant parameters, and nonautonomous guiding fields, it is possible to reproduce any desired pattern. As such, it is crucial in the mathematical modelling of living systems to delineate the essential requirements from the incidental.
We agree that a mathematical model with spatial locality alone does not imply selforganized dynamics. However, nothing in the formulation of the model locally at the resolution of cellcell interactions was designed explicitly (in terms of any global supervisory mechanism present in the model) to produce the key aspects of the emergent patterns as they may be described at the global scale (i.e., a Voronoi tesselation of fields that each comprise hundreds of cells). The guidance gradients (and the overall boundary shape) are indeed defined globally, but these gradients vary linearly across the domain, and they do not directly specify the Dirichletform shapes that emerge from cellcell interactions. We trust that this point is clearer in the light of the results of the parameter sweep and sensitivity analysis, which demonstrate that the form of the patterns that develop in simulation is highly robust to variation in any of the component parameters.
https://doi.org/10.7554/eLife.55588.sa2Article and author information
Author details
Funding
James S. McDonnell Foundation (220020516)
 Sebastian S James
 Leah A Krubitzer
 Stuart P Wilson
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
The authors thank Jason Berwick at the University of Sheffield for advice and for access to the rat barrel stains used to construct Figure 1A. This work was supported by a Collaborative Activity Award, Cortical Plasticity Within and Across Lifetimes, from the James S McDonnell Foundation (grant 220020516).
Senior Editor
 Timothy E Behrens, University of Oxford, United Kingdom
Reviewing Editor
 Upinder Singh Bhalla, Tata Institute of Fundamental Research, India
Reviewer
 Bard Ermentrout, University of Pittsburgh, United States
Publication history
 Received: January 29, 2020
 Accepted: August 26, 2020
 Version of Record published: September 29, 2020 (version 1)
Copyright
© 2020, James et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,253
 Page views

 112
 Downloads

 3
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Computational and Systems Biology
 Neuroscience
Biological motor control is versatile, efficient, and depends on proprioceptive feedback. Muscles are flexible and undergo continuous changes, requiring distributed adaptive control mechanisms that continuously account for the body's state. The canonical role of proprioception is representing the body state. We hypothesize that the proprioceptive system could also be critical for highlevel tasks such as action recognition. To test this theory, we pursued a taskdriven modeling approach, which allowed us to isolate the study of proprioception. We generated a large synthetic dataset of human arm trajectories tracing characters of the Latin alphabet in 3D space, together with muscle activities obtained from a musculoskeletal model and modelbased muscle spindle activity. Next, we compared two classes of tasks: trajectory decoding and action recognition, which allowed us to train hierarchical models to decode either the position and velocity of the endeffector of one's posture or the character (action) identity from the spindle firing patterns. We found that artificial neural networks could robustly solve both tasks, and the networks'units show tuning properties similar to neurons in the primate somatosensory cortex and the brainstem. Remarkably, we found uniformly distributed directional selective units only with the actionrecognitiontrained models and not the trajectorydecodingtrained models. This suggests that proprioceptive encoding is additionally associated with higherlevel functions such as action recognition and therefore provides new, experimentally testable hypotheses of how proprioception aids in adaptive motor control.

 Computational and Systems Biology
Correlation between objects is prone to occur coincidentally, and exploring correlation or association in most situations does not answer scientific questions rich in causality. Causal discovery (also called causal inference) infers causal interactions between objects from observational data. Reported causal discovery methods and singlecell datasets make applying causal discovery to single cells a promising direction. However, evaluating and choosing causal discovery methods and developing and performing proper workflow remain challenges. We report the workflow and platform CausalCell (http://www.gaemons.net/causalcell/causalDiscovery/) for performing singlecell causal discovery. The workflow/platform is developed upon benchmarking four kinds of causal discovery methods and is examined by analyzing multiple singlecell RNAsequencing (scRNAseq) datasets. Our results suggest that different situations need different methods and the constraintbased PC algorithm with kernelbased conditional independence tests work best in most situations. Related issues are discussed and tips for best practices are given. Inferred causal interactions in single cells provide valuable clues for investigating molecular interactions and gene regulations, identifying critical diagnostic and therapeutic targets, and designing experimental and clinical interventions.