Prairie voles (Microtus ochrogaster) are one of the few mammalian species that are monogamous and engage in the biparental rearing of their offspring. Biparental care impacts the quantity and quality of care the offspring receives. The increased attention by the father may translate to heightened tactile contact the offspring receives through licking and grooming. In the current study, we used electrophysiological multiunit recording techniques to define the organization of the perioral representation in the primary somatosensory area (S1) of prairie voles. Functional representations were related to myeloarchitectonic boundaries. Our results show that most of S1 is occupied by the representation of the contralateral mystacial whiskers and the lower and upper lips. The mystacial vibrissae representation encompassed a large portion of the caudolateral S1, while the representation of the lower and upper lips occupied a large portion of the rostrolateral aspect of S1. We found that neuronal populations representing the perioral structures tended to have small receptive fields relative to other body part representations on the head. The representation of the mystacial whiskers and perioral structures was coextensive with cytoarchitectonically defined barrel fields that extend from the caudolateral to a rostrolateral aspect of S1. We discuss our findings in the context of the magnification of behaviorally relevant sensory surfaces in other rodents, the ubiquity of the barrel systems in rodents, and behaviors associated with specialized sensory surfaces.

The organization of the extant mammalian brain is influenced by development, evolutionary history and the environment. Ecological adaptations specifically have had a major role in shaping the structures and associated functions of the mammalian brain. Although general organization of the brain is relatively conserved in modern mammals, throughout millions of years of evolution mammals have acquired diverse sensory and nervous system adaptations as they invaded new ecological niches. Here, we synthesize palaeontological and neurobiological evidence on mammalian brain structure evolution, the mechanisms behind the observed variation in the size and organization of brain structures, and the effect of behavioural ecology on the evolution of brain functions and associated structures. Neuroecology has advanced greatly over the past 40 years and is now unravelling the complex relationship between specific behaviours and brain organization and function. Relying on different types of data, comparative neurobiologists and palaeontologists strive to answer similar questions about brain evolution, benefiting from a synergistic approach. We conclude this Review by outlining outstanding questions regarding the relationships between structure, function, behaviour and evolution that deserve future research attention, and propose methodologies and approaches to help to resolve these problems. Palaeontologists and comparative neurobiologists share a common interest in the evolution of the mammalian brain, but often fail to realize the benefits of this shared interest. This Review draws these fields together, demonstrating the utility of a cross-disciplinary, synergistic approach.

INTRODUCTION:The gray short-tailed opossum, Monodelhis domestica (M. domestica), is a widely used marsupial model species that presents unique advantages for neurodevelopmental studies. Notably their extremely altricial birth allows manipulation of postnatal pups at timepoints equivalent to embryonic stages of placental mammals. A robust literature exists on the development of short-tailed opossums, but many researchers working in the more conventional model species of mice and rats may find it daunting to identify the appropriate age at which to conduct experiments. METHODS:Here, we present detailed staging diagrams taken from photographic observations of 40 individual pups, in 6 litters, over 25 timepoints across postnatal development. We also present a comparative neurodevelopmental timeline of short-tailed opossums (M. domestica), the house mouse (Mus musculus), and the laboratory rat (Rattus norvegicus) during embryonic as well as postnatal development, using timepoints taken from this study and a review of existing literature, and use this dataset to present statistical models comparing the opossum to the rat and mouse. RESULTS:One aim of this research was to aid in testing the generalizability of results found in rodents to other mammalian brains, such as the more distantly related metatherians. However, this broad dataset also allows the identification of potential heterochronies in opossum development compared to rats and mice. In contrast to previous work, we found broad similarity between the pace of opossum neural development with that of rats and mice. We also found that development of some systems was accelerated in the opossum, such as the forelimb motor plant, oral motor control, and some aspects of the olfactory system, while the development of the cortex, some aspects of the retina, and other aspects of the olfactory system are delayed compared to the rat and mouse. DISCUSSION:The pace of opossum development is broadly similar to that of mice and rats, which underscores the usefulness of this species as a compliment to the more commonly used rodents. Many features that differ the most between opossums and rats and mice were either clustered around the day of birth and were features that have functional importance for the pup immediately after or during birth, or were features that have reduced functional importance for the pup until later in postnatal development, given that it is initially attached to the mother.

Motor cortex and anterior and posterior parietal cortex form a sensorimotor integration network. We tested the extent to which parietal areas could initiate movements independent of M1. Our findings support the contention that, although areas 2, 5L, PF, and PFG are highly dependent on M1 to produce movement, area 1 may constitute a parallel corticospinal pathway that can function somewhat independently of M1. A similar functional architecture may underlie dexterous tool use in humans.

A central question in comparative neurobiology concerns how evolution has produced brains with expanded neocortices, composed of more areas with unique connectivity and functional properties. Some mammalian lineages, such as primates, exhibit exceptionally large cortices relative to the amount of sensory inputs from the dorsal thalamus, and this expansion is associated with a larger number of distinct cortical areas, composing a larger proportion of the cortical sheet. We propose a link between the organization of the neocortex and its expansion relative to the size of the dorsal thalamus, based on a combination of work in comparative neuroanatomy and experimental research.

In the following review, we describe the types of phenotypic changes to the neocortex that occur over the longer time scale of evolution, and over the shorter time scale of an individual lifetime. To understand how phenotypic variability emerges in the neocortex, it is important to consider the cortex as part of an integrated system of the brain, the body, the environment in which the brain and body develops and evolves, and the affordances available within a particular environmental context; changes in any part of this brain/body/environment network impact the neocortex. We provide data from comparative studies on a wide variety of mammals that demonstrate that body morphology, the sensory epithelium, and the use of a particular morphological structure have a profound impact on neocortical organization and connections. We then discuss the genetic and epigenetic factors that contribute to the development of the neocortex, as well as the role of spontaneous and sensory driven activity in constructing a nervous system. Although the evolution of the neocortex cannot be studied directly, studies in which developmental processes are experimentally manipulated provide important insights into how phenotypic transformations could occur over the course of evolution and demonstrate that relatively small alterations to the body and/or the environment in which an individual develops can manifest as large changes to the neocortex. Finally, we discuss how these phenotypic alterations to the neocortex impact an important target of selection – behavior.

Advances in sequencing techniques have made comparative studies of gene expression a current focus for understanding evolutionary and developmental processes. However, insights into the spatial expression of genes have been limited by a lack of robust methodology. To overcome this obstacle, we developed methods and software tools for quantifying and comparing tissue-wide spatial patterns of gene expression within and between species. Here, we compare cortex-wide expression of RZRβ and Id2 mRNA across early postnatal development in mice and voles. We show that patterns of RZRβ expression in neocortical layer 4 are highly conserved between species but develop rapidly in voles and much more gradually in mice, who show a marked expansion in the relative size of the putative primary visual area across the first postnatal week. Patterns of Id2 expression, by contrast, emerge in a dynamic and layer-specific sequence that is consistent between the two species. We suggest that these differences in the development of neocortical patterning reflect the independent evolution of brains, bodies, and sensory systems in the 35 million years since their last common ancestor.

In the present investigation, we examined the role of different cortical fields in the fronto-parietal reaching and grasping network in awake, behaving macaque monkeys. This network is greatly expanded in primates compared to other mammals and coevolved with glabrous hands with opposable thumbs and the extraordinary dexterous behaviors employed by a number of primates, including humans. To examine this, we reversibly deactivated the primary motor area (M1), anterior parietal area 2, and posterior parietal areas 5L and 7b individually while monkeys were performing two types of reaching and grasping tasks. Reversible deactivation was accomplished with small microfluidic thermal regulators abutting specifically targeted cortical areas. Placement of these devices in the different cortical fields was confirmed post hoc in histologically processed tissue. Our results indicate that the different areas examined form a complex network of motor control that is overlapping. However, several consistent themes emerged that suggest the independent roles that motor cortex, area 2, area 7b, and area 5L play in the motor planning and execution of reaching and grasping movements. Area 5L is involved in the early stages and area 7b the later stages of a reaching and grasping movement, motor cortex is involved in all aspects of the execution of the movement, and area 2 provides proprioceptive feedback throughout the movement. We discuss our results in the context of previous studies that explored the fronto-parietal network, the overlapping (but also independent) functions of different nodes of this network, and the rapid compensatory plasticity of this network. NEW & NOTEWORTHY This is the first study to directly compare the results of cooling different portions of the fronto-parietal reaching and grasping network (motor cortex, anterior and posterior parietal cortex) in the same animals and the first to employ a complex, bimanual reaching and grasping task that is ethologically relevant. Whereas cooling area 7b or area 5L evoked deficits at distinct task phases, cooling M1 evoked a general set of deficits and cooling area 2 evoked proprioceptive deficits.

Bats have evolved behavioral specializations that are unique among mammals, including self-propelled flight and echolocation. However, areas of motor cortex that are critical in the generation and fine control of these unique behaviors have never been fully characterized in any bat species, despite the fact that bats compose similar to 25% of extant mammalian species. Using intracortical microstimulation, we examined the organization of motor cortex in Egyptian fruit bats (Rousettus aegyptiacus), a species that has evolved a novel form of tongue-based echolocation.(1,2) We found that movement representations include an enlarged tongue region containing discrete subregions devoted to generating distinct tongue movement types, consistent with their behavioral specialization generating active sonar using tongue clicks. This magnification of the tongue in motor cortex is comparable to the enlargement of somatosensory representations in species with sensory specializations.(3-5) We also found a novel degree of coactivation between the forelimbs and hindlimbs, both of which are involved in altering the shape and tension of wing membranes during flight. Together, these findings suggest that the organization of motor cortex has coevolved with peripheral morphology in bats to support the unique motor demands of flight and echolocation.

Our understanding of the neural basis of somatosensation is based largely on studies of the whisker system of mice and rats and the hands of macaque monkeys. Results across these animal models are often interpreted as providing direct insight into human somatosensation. Work on these systems has proceeded in parallel, capitalizing on the strengths of each model, but has rarely been considered as a whole. This lack of integration promotes a piecemeal understanding of somatosensation. Here, we examine the functions and morphologies of whiskers of mice and rats, the hands of macaque monkeys, and the somatosensory neuraxes of these three species. We then discuss how somatosensory information is encoded in their respective nervous systems, highlighting similarities and differences. We reflect on the limitations of these models of human somatosensation and consider key gaps in our understanding of the neural basis of somatosensation.

This protocol presents a workflow for detecting differences in kinematics between experimental conditions. It is tailored for short-tailed opossums but can be applied to any species capable of completing the ladder rung task. There are four phases of this protocol: (1) data collection, (2) pose tracking, (3) analysis of single trials, and (4) cross-condition comparisons. This pipeline implements aspects of machine learning and signal processing, allowing for rapid data analysis that provides insight into how animals perform this task.For complete details on the use and execution of this protocol, please refer to Englund et al. (2020).

Advances in sequencing techniques have made comparative studies of gene expression a current focus for understanding evolutionary and developmental processes. However, insights into the spatial expression of genes have been limited by a lack of robust methodology. We therefore developed a set of algorithms for quantifying and comparing tissue-wide spatial patterns of gene expression within and across species. Here we apply these algorithms to compare cortex-wide expression of Id2 and RZRβ mRNA in early postnatal mice and voles. We show that neocortical patterns of Id2 expression are moderately conserved between species, but that the degree of conservation varies by cortical layer and area. By comparison, patterns of RZRβ expression are highly conserved in somatosensory areas, and more variable between species in visual and auditory areas. We consider if these differences reflect independent evolution in the 35 million years since the last common ancestor.

ABSTRACT Behavioral strategies that depend on sensory information are not immutable; rather they can be shaped by the specific sensory context in which animals develop. This behavioral plasticity depends on the remarkable capacity of the brain to reorganize in response to alterations in the sensory environment, particularly when changes in sensory input occur at an early age. To study this phenomenon, we utilize the short-tailed opossum, a marsupial that has been a valuable animal model to study developmental plasticity due to the extremely immature state of its nervous system at birth. Previous studies in opossums have demonstrated that removal of retinal inputs early in development results in profound alterations to cortical connectivity and functional organization of visual and somatosensory cortex; however, behavioral consequences of this plasticity are not well understood. We trained early blind and sighted control opossums to perform a two-alternative forced choice texture discrimination task. Whisker trimming caused an acute deficit in discrimination accuracy for both groups, indicating the use of a primarily whisker-based strategy to guide choices based on tactile cues. Mystacial whiskers were important for performance in both groups; however, genal whiskers only contributed to behavioral performance in early blind animals. Early blind opossums significantly outperformed their sighted counterparts in discrimination accuracy, with discrimination thresholds that were lower by ∼75 μm. Our results support behavioral compensation following early blindness using tactile inputs, especially the whisker system.

ASSOCIATE EDITORS DEANNA L. BENSON Icahn School of Medicine at Mount Sinai EDWARD M. CALLAWAY Salk Ins tute THOMAS E. FINGER University of Colorado School of Medicine JEFFREY H. KORDOWER Rush University Medical Center ANDREW D. HUBERMAN Stanford University School of Medicine, Stanford, CA, USA KATHLEEN S. ROCKLAND Boston University School of Medicine JOHN L.R. RUBENSTEIN University of California, San Francisco UWE HOMBERG University of Marburg

Brain development relies on an interplay between genetic specification and self-organization. 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 self-organize, 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.

ABSTRACT The early loss of vision results in a reorganized visual cortex that processes tactile and auditory inputs. Recent studies in the short-tailed opossum ( Monodelphis domestica) found that the connections and response properties of neurons in somatosensory cortex of early blind animals are also altered. While research in humans and other mammals shows that early vision loss leads to heightened abilities on discrimination tasks involving the spared senses, if and how this superior discrimination leads to adaptive sensorimotor behavior has yet to be determined. Moreover, little is known about the extent to which blind animals rely on the spared senses. Here, we tested early blind opossums on a sensorimotor task involving somatosensation and found that they had increased limb placement accuracy. However, increased reliance on tactile inputs in early blind animals resulted in greater deficits in limb placement and behavioral flexibility when the whiskers were trimmed.

The early loss of vision results in a reorganized neocortex, affecting areas of the brain that process both the spared and lost senses, and leads to heightened abilities on discrimination tasks involving the spared senses. Here, we used performance measures and machine learning algorithms that quantify behavioral strategy to determine if and how early vision loss alters adaptive sensorimotor behavior. We tested opossums on a motor task involving somatosensation and found that early blind animals had increased limb placement accuracy compared with sighted controls, while showing similarities in crossing strategy. However, increased reliance on tactile inputs in early blind animals resulted in greater deficits in limb placement and behavioral flexibility when the whiskers were trimmed. These data show that compensatory cross-modal plasticity extends beyond sensory discrimination tasks to motor tasks involving the spared senses and highlights the importance of whiskers in guiding forelimb control.

One of the hallmarks of human evolution is the extraordinary degree to which we can manipulate the physical world with our hands or with tools that extend or amplify our limbs. This manual dexterity coevolved with an expansion of posterior parietal cortex (PPC), which contains areas involved in programming voluntary movements, coding reach targets in multiple reference frames, and decision-making. To enable our body to interact with our physical surroundings, these fields must also construct an internal model of the physical self: our body's configuration, the boundary between our body and external physical objects, and the temporary expansion of that self as we wield a tool that extends our reach and manual capabilities. Such comprehension of where and what the self is and even the ability to manipulate objects and use them as tools did not evolve de novo in humans, but rather emerged from simple networks present in early mammals. In this chapter, we consider not only the structure and function of PPC in primates, but also include data from nonprimate mammals in an effort to understand the basic processing networks that were present in our early ancestors, and how these networks evolved and expanded in primates.

Article Figures and data Abstract eLife digest Introduction Results Discussion Methods Data availability References Decision letter Author response Article and author information Metrics Abstract Brain development relies on an interplay between genetic specification and self-organization. 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 self-organize, 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 whisker-barrel 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 centre-point 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 self-organize without a full genetic blueprint pre-specifying the barrel centre-points 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 whisker-specific clusters. The simulations also revealed that the target tissue does not need to specify centre-points 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 (Armstrong-James 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ón-Bolaños et al., 2019), and depends on a complex network of axon guidance molecules such as ephrin-A5 and A7 and adhesion molecules such as cadherin-6 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 Coup-tf1 (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 centre-point that ceases upon contact with axons branching from adjacent centres. However, the assumption of pre-arranged centre-points 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 whisker-related clusters then become apparent in the PrV at P0-P1, in the VPM at P2-P3, and in the cortex at P3-P5 (Erzurumlu and Gaspar, 2012; Sehara and Kawasaki, 2011). Figure 1 with 1 supplement see all Download asset Open asset The emergence of whisker barrels. (A) Left shows a cytochrome oxidase (CO) stain obtained from rat S1 by Zheng et al., 2001, with black lines to delineate barrels and to measure departure (Honda-δ; see Senft and Woolsey, 1991) from a perfect Voronoi tessellation. Right shows the initial distribution of axon branching density (a) for one thalamocortical projection, and two molecular guidance fields (ρ), where the domain S has been traced from the CO stain. (B) The strengths of interaction γ with fields ρ1 and ρ2 are indicated for each of 41 projections by the lengths of green and blue arrows respectively, assuming that similar fields aligned to the posterior-anterior and medial-lateral axes in the ventroposterior medial nucleus of the thalamus are sampled at the locations of putative barreloid centres (reconstructed from Haidarliu and Ahissar, 2001, their Figure 5b). (C) Results for the example simulation, with parameters N=41, α=3.6, β=16.67, k=3, D=0.5, γ∈±2, ϵ=1.2 and δ⁢t=0.0001. Colours indicate the thalamic projection for which the connection density is maximal, barrel labels are located at the centroid of each region and black lines delineate boundaries (see Figure 1—video 1). (D) Red dots show the Honda–δ⁢(t) metric obtained from the simulation approaching that obtained from the real barrels in A (dotted line); black squares show the pattern difference metric η⁢(t), and reveal the emergence of a correspondence between the real and simulated barrel shapes (units mm3); grey hexagons show how selectively each cortical site is innervated; ω(t)=∯Sμ(x,t)dS, where μ⁢(𝐱)≡maxi⁢(ci⁢(𝐱,t))/∑j=1Ncj⁢(𝐱,t). (E) Plotted across the cortical sheet, the selectivity develops to reveal an alignment with the emergent barrel boundary shapes. Greyscale colour indicates values of μ⁢(𝐱). All scale bars 1 mm. Alternatively, reaction-diffusion dynamics could generate a Voronoi tessellation without pre-arranged centres, by amplifying characteristic modes in a noisy initial distribution of axon branches, as a net effect of short-range cooperative and long-range 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 centre-point and reaction-diffusion models are not mutually exclusive. Pre-organized centres could bias reaction-diffusion processes to generate specific arrangements more reliably, and mechanisms of lateral axon branching may constitute the tension between cooperation and competition required for self-organization. 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 pre-organizing thalamocortical connections need not be assumed. To this end, we ask whether barrel maps can emerge in a system with reaction-diffusion dynamics, under the guidance of signalling gradients, and in the absence of pre-defined centres. Models Karbowski and Ermentrout, 2004 developed a reaction-diffusion 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 ci⁢(x,t) and the density of axon branches ai⁢(x,t) interact at time t, along a 1D anterior-posterior axis x, for N thalamocortical projections indexed by i. The model was derived from the assumption that the rates at which ai and ci grow are reciprocally coupled. Extending the original 1D model to simulate arealization on a 2D cortical sheet, we use ai⁢(𝐱,t) and ci⁢(𝐱,t), and model synaptogenesis as (1) ∂⁡ci∂⁡t=-α⁢ci+β⁢(1-∑j=1Ncj)⁢[ai]k. Accordingly, where the total fraction of synaptic connections sums to one, connections decay at rate α. Otherwise, ci⁢(𝐱,t) increases non-linearly (k>1) with the density of axon branching. Axon branching is modelled as (2) ∂ai∂t= ∇ ⋅ (D∇ai−ai∑j=1Mγi,j∇ρj(x)+χi)−∂ci∂t. The first term on the right describes the divergence (indicated by ∇⋅) 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, ρ⁢(𝐱). The influence of a given field (indexed by j) on a given thalamic projection (indexed by i), is determined by γ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 χ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 ci 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≈(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 centre-points by specifying conditions in the target tissue that instruct pre-identified afferents to stop growing. Consider, for example, that when simulated in isolation from one-another, 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 ai is conserved over time, we redefined (3) χi(x,t)=ϵaiN−1∇∑j≠iNaj. This term contributes to the flux of axonal branching as an additional source of diffusion, scaled by ϵ, 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 ephrin-A5 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 ephrin-A5, 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 ephrin-A5 (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., ephrin-A5) 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 anterior-posterior and medial-lateral axes, were estimated from Figure 5d of Haidarliu and Ahissar, 2001, and scaled such that γ∈±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 anterior-posterior and medial-lateral axes. These gradients are shown with the barrel field boundary in Figure 1A for clarity, though like ephrin-A5 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 (δ⁢t=0.0001). Across a wide range of parameter values, random initial conditions (a uniform random distribution for a⁢(𝐱,0)∈(0.2,0.4), c⁢(𝐱,0)=0) eventually yielded a clear Voronoi-like tessellation of topographically organized thalamocortical projections, confirming that barrel maps can self-organize in the absence of pre-specified 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 (α=3.6, β=16.67, k=3, D=0.5, ϵ=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 δ≈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, δ=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 self-organized 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 self-organization 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, ϵ, and the ratio α/β. First, we observe that self-organization is highly robust to the ratio α/β, across five orders of magnitude, with respect to all three metrics. Second, the most strongly Voronoi-conforming patterns (low Honda-δ) were generated by simulations in which the diffusion constant D and the strength of competition ϵ 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≈0.5. Together these results indicate that when competition is strong, the rate of diffusion determines a trade-off such that fields emerge to be barrel-shaped when diffusion is fast and they emerge to be more selectively innervated when diffusion is slow. Figure 2 Download asset Open asset Exploring the parameter space. (A) Colour indicates the quality of the pattern at t=30000 steps, against three measures: Low values of Honda–δ (top row) suggest a Voronoi-like pattern of fields. The pattern difference, η (middle row) measures the difference in the area and arrangement of barrels in real and simulated fields. The connection selectivity, ω (bottom row), measures the specificity with which the cortical sheet is innervated. Colour maps are chosen so that lighter (orange and yellow) colours indicate higher quality patterns against each measure. White squares indicate combinations of parameters for which simulations were numerically unstable. The parameter space explored is three dimensional with the ratio α/β varying between plots, and the competition parameter ϵ and the diffusion constant D varying within plots. An asterisk (*) marks the parameters used in Figure 1. Boxes (i), (ii) and (iii) mark parameter sets for which corresponding patterns are shown in B (for t=30000 steps). (B) Varying the diffusion constant D generates qualitatively different patterns. Higher values cause expansions of the peripheral barrels and a corresponding compression of the inner barrels (i). Lower values instead cause an expansion of the central barrels and compression of the peripheral barrels (ii). Further reducing the rate of diffusion (iii), which is equivalent to increasing the size of the domain and hence simulating development in an animal with a larger cortex, causes a large area to be occupied by projections with intermediate interaction parameters; those with strong interaction parameters are compressed around the edge of the domain, and consequently a barrel pattern fails to form. 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 ai at t=0; (ii) noise applied to the interaction parameters, γ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⁢(𝐱,0)∈(0.2,0.4), to be δ=0.089±0.004, η=0.2108±0.002, and ω=0.2165±0.0001. Repeating with the variation in the initial noise doubled (a⁢(𝐱,0)∈(0.1,0.5)), or removed altogether (a⁢(𝐱,0)=0.3), generated distributions of δ, η, and ω that were not statistically different, as established using paired two-sample t-tests. 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, ρ1⁢(𝐱) and ρ2⁢(𝐱), disrupted pattern formation only for high levels of noise applied at short length scales, which manifested as non-straight edges at the domain boundaries (Figure 3Bii). Varying the slope of one linear gradient ρ1⁢(𝐱) 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 self-organization of barrel-like fields in the model is highly robust to a wide range of sources of perturbation. Figure 3 Download asset Open asset Sensitivity analysis. (A) Metrics of map quality δ⁢(t) (top row), ω⁢(t) (middle row) and η⁢(t) (bottom row), were evaluated at t=30000 steps. The y-axes and colour scales have identical ranges to the colour scales in Figure 2, for easy comparison. Left column: The effect of adding noise drawn from a uniform distribution, (γm⁢a⁢x-γm⁢i⁢n)⁢U⁢(0,νγ), to the values of the interaction parameters, γ used in Figure 1. Middle column: The effect of changing the magnitude and the length scale of noise applied to the guidance fields. Uniform random noise (ρm⁢a⁢x-ρm⁢i⁢n)⁢U⁢(0,νρ) was added to each (hexagonal) element of ρ1⁢(𝐱) and ρ2⁢(𝐱) and the result was smoothed by convolution with a symmetric 2D Gaussian kernel of width σρ. Right column: The effect of setting the rotational angle of the linearly varying guidance field, ρ1⁢(𝐱), to ϕρ1, and modifying its overall gain to Gρ1, whilst keeping the parameters of ρ2⁢(𝐱) unchanged from those used in the example simulation, for which ϕρ2=84° and Gρ2=1. (B) Four ways in which the perturbations in A affect the patterns. (i) The effect of significant interaction parameter noise is to introduce topological defects. (ii) High magnitude, short length scale noise in ρ⁢(𝐱) leads to non-straight edges between adjacent barrels. (iii) Reducing the slope of gradient ρ1 (by a factor of 10) causes barrel rows B, C and D to become ‘crushed’ down the centre line, and edge barrels to dominate. (iv) Rotating ρ1 by 20° causes a slight distortion of the pattern, resulting in an overall anticlockwise rotation of the field locations. 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 anterior-posterior guidance gradient ρ1 at the border between them (Figure 4A). The result after 30000 iterations, and otherwise using the parameters of the example simulation, was two mirror-symmetrical barrel fields comprising 2⁢N barrels (Figure 4B), consistent with the outcome of the original experiments. Figure 4 Download asset Open asset Simulating altered barrel development. Guidance fields (A) and emergent barrel pattern (B) in a Fgf8 misexpression experiment (Assimacopoulos et al., 2012), simulated by reflecting ρ1 from Figure 1A at the join of the original boundary with its mirror. All other model parameters match those in Figure 1C. C Simulating whisker trimming by reducing the competitiveness, ϵ, of one projection. For C3 only, ϵ was multiplied by m∈(0,1). The pattern shown is that formed after 30000 steps with m=0.86, which reduces the size of the C3 field to 65% of its original size, matching the average barrel area reduction observed by Kossut, 1992. D The area of the C3 barrel (black squares) reduces as m is reduced, whereas the mean area of neighbouring barrels (B3, C2, C4, D2 and D3, grey circles) increases. The dotted grey line indicates m=0.86 for comparison with panel C. If the ϵ value is instead reduced for all row C projections (including the interstitial γ projection), the mean row C barrel area is only slightly reduced (light grey triangles). In this case, the mean area of a simulated row C barrel at m=0.86 is 91% of that for m=1. Finally, to investigate the response of the model to environmental man