The architecture of iso-orientation domains in the primary visual cortex (V1) of placental carnivores and primates apparently follows species invariant quantitative laws. Dynamical optimization models assuming that neurons coordinate their stimulus preferences throughout cortical circuits linking millions of cells specifically predict these invariants. This might indicate that V1’s intrinsic connectome and its functional architecture adhere to a single optimization principle with high precision and robustness. To validate this hypothesis, it is critical to closely examine the quantitative predictions of alternative candidate theories. Random feedforward wiring within the retino-cortical pathway represents a conceptually appealing alternative to dynamical circuit optimization because random dimension-expanding projections are believed to generically exhibit computationally favorable properties for stimulus representations. Here, we ask whether the quantitative invariants of V1 architecture can be explained as a generic emergent property of random wiring. We generalize and examine the stochastic wiring model proposed by Ringach and coworkers, in which iso-orientation domains in the visual cortex arise through random feedforward connections between semi-regular mosaics of retinal ganglion cells (RGCs) and visual cortical neurons. We derive closed-form expressions for cortical receptive fields and domain layouts predicted by the model for perfectly hexagonal RGC mosaics. Including spatial disorder in the RGC positions considerably changes the domain layout properties as a function of disorder parameters such as position scatter and its correlations across the retina. However, independent of parameter choice, we find that the model predictions substantially deviate from the layout laws of iso-orientation domains observed experimentally. Considering random wiring with the currently most realistic model of RGC mosaic layouts, a pairwise interacting point process, the predicted layouts remain distinct from experimental observations and resemble Gaussian random fields. We conclude that V1 layout invariants are specific quantitative signatures of visual cortical optimization, which cannot be explained by generic random feedforward-wiring models.
Processing high-dimensional external stimuli and efficiently communicating their essential features to higher brain areas is a fundamental function of any sensory system. For many sensory modalities, this task is implemented via convergent and divergent neural pathways in which information from a large number of sensors is compressed into a smaller layer of neurons, transmitted, and then re-expanded into a larger neuronal layer. When sensory inputs are sparse, compression of the inputs through random convergent feedforward projections has been shown to retain much of the information present in the stimuli [1–3]. On the other hand, random expanding projections can lead to computationally powerful high-dimensional representations of such compressed signals, which combine separability of the inputs with high signal-to-noise ratio to facilitate downstream readouts [4]. Given these computational benefits, one might expect randomness to be a fundamental wiring principle employed by different sensory systems. The most striking example of a random expansion so far has been observed in the olfactory system of Drosophila melanogaster. Kenyon cells in the fly brain’s mushroom body were shown to integrate input from various olfactory glomeruli in combinations that are consistent with purely random choices from the overall distribution of glomerular projections to the mushroom body [5]. What is the role of random projections between neural layers in mammalian sensory systems? Sompolinsky and others have argued that the human visual system, for instance, implements a compression-transmission-expansion strategy [3,4]. In fact, visual stimulus information is transmitted from about 5 million cone photoreceptors [6,7] to 1 million retinal ganglion cells (RGCs) [7] and then via the optic nerve to about 1 million lateral geniculate relay cells [8] to on the order of 100 million neurons in the primary visual cortex (V1) [9,10]. We note that, while the overall connectivity indeed suggests compression for peripheral retinal regions [11], close to the fovea RGC density is higher than the density of photoreceptors [12,13]. How much does randomness contribute to shaping the functional architecture of early visual cortical areas? Projections between individual layers of the early mammalian visual pathway are clearly not entirely random. Visual information is mapped visuotopically from the retina to V1 such that neighboring groups of V1 neurons process information from neighboring regions in visual space. Yet, it has long been realized that many features of the spatial progression of receptive fields across V1 layer IV naturally emerge if random feedforward connections from groups of RGC cells to layer IV neurons (via the lateral geniculate nucleus (LGN) ) are assumed (see [14] for an early example). The most important of such features is orientation selectivity, i. e. the selective response to edge-like stimuli of a particular orientation. In carnivores, primates and their close relatives, orientation selectivity is arranged in patterns of iso-orientation domains. Iso-orientation domains (orientation domains for short) in V1 exhibit a continuous, roughly repetitive arrangement. A distance in the millimeter range, called the column spacing, separates close-by domains preferring the same orientation. The continuous progression of preferred orientations is interrupted by a system of topological defects, called pinwheel centers, at which neurons selective to the whole complement of stimulus orientations are located in close vicinity [15–20]. These topological defects exhibit two distinct topological charges, indicating that preferred orientations change clockwise or counterclockwise around the defect center [15,18,21–23]. More than 25 years ago, Soodak [24,25] (see also [26]) proposed random wiring between irregularly positioned retinal ganglion cells (RGCs) and layer IV neurons in V1 via the thalamus as a candidate mechanism defining the pattern of iso-orientation domains. According to this statistical wiring hypothesis, a V1 neuron randomly samples feedforward inputs from geniculate projections in the immediate vicinity of its receptive field center (see e. g. [27]). The neuron then is likely to receive the strongest inputs from a central pair of ON/OFF RGCs, forming a so-called RGC dipole [28–30]. In this scheme, one ON and one OFF subregion dominate the receptive field (RF) of the V1 neuron and its response is tuned to the orientation perpendicular to the dipole axis. Thus, the preferred orientation of the neuron in this case is determined by the orientation of the RGC dipole. Consequently, the key prediction of the statistical wiring hypothesis is that the spatial arrangement of ON/OFF RGC cells in the retina essentially determines the spatial layout of orientation preference domains in V1. Recently, Paik & Ringach showed that the statistical wiring hypothesis—when constructed with a hexagonal grid of RGCs—predicts a periodic orientation domain layout with a hexagonal autocorrelation function [28]. Moreover, it predicts that orientation preference is differently linked to the visuotopic map around pinwheels of positive or negative topological charge [29]. Qualitative signatures of both predictions were reported to be present in experimentally measured patterns [28,29]. Thus the statistical wiring model has conceptual appeal and is a mechanistically particularly transparent candidate explanation for V1 functional architecture (see however [31,32]). Does the predictive power of the random wiring hypothesis for the early visual pathway reach beyond this qualitative agreement? The recent discovery of species-invariant quantitative layout laws for the arrangement of pinwheel centers in tree shrews, galagos and ferrets [23] provides a unique opportunity to address this question. Kaschube et al. demonstrated that in these species, the statistics of pinwheel defect layouts is quantitatively invariant, with potential deviations in geometrical layout parameters of at most a few percent [23]. Specifically, the overall pinwheel density, defined as the average number of defects within the area of one square column spacing Λ2 was found to be virtually identical. Subsequently, orientation domain layouts from cat V1 were shown to exhibit pinwheel densities very close to those of the three species previously studied [33]. Additionally, Kaschube et al. found an entire set of local and non-local quantitative pinwheel layout features to be species-invariant (see below). Following [23], we refer to this overall layout of orientation domains as the common design. During mammalian evolution, the common design most likely arose independently in carnivores and euarchontans and potentially even in scandentia [23,34]. This is suggested by two lines of evidence: (i) The four species in which the common design has been observed so far are widely separated in terms of evolutionary descent, belonging to distinct supra-ordinal clades that split already during basal radiation of placentals [35–42] (Fig 1A, see also [23,33]). Their last common ancestor was a small shrew-like mammal [40–42] that is unlikely to have possessed a columnar V1 architecture [23,34]. (ii) Distinct neuronal circuits underlie the generation of orientation selectivity in galago, ferret, tree shrew, and cat (Fig 1B). Tree shrews, for instance, lack orientation selectivity in the input layer IV of V1 [43,44] and use intracortical circuits to compute contour orientation. In contrast, cats exhibit both, orientation selectivity and organization of selectivity into orientation domains already in layer IV and thus first generate orientation selectivity by thalamo-cortical circuits [45,46] (see Fig 1B for further differences). Kaschube et al. used a dynamical self-organization model with long-range suppressive interactions, the long-range interaction model, to explain all features of the common design [23]. The hypothesis that randomness of feedforward connections between the retina/LGN and V1 could explain the common design is conceptually diametrically opposed to large-scale self-organization. In the long-range interaction model, the orientation preference of a neuron is chosen from an, in principle, unlimited afferent repertoire of potential receptive fields. Single neurons dynamically select a particular preferred orientation as a result of large-scale circuit interactions involving millions of other cortical neurons. In the statistical connectivity model, to the contrary, the preferred orientation of a cortical neuron is essentially imposed by the alignment of only one pair of neighboring ON-OFF RGCs, a local process involving in principle not more than 5 cells. Can the invariant layout laws of iso-orientation domains and pinwheels be explained as the generic outcome of a locally stochastic feedforward wiring of the early visual pathway? More generally, do iso-orientation domains and pinwheels in different species adhere to identical layout laws because any mechanism that generates a retinotopic random feedforward circuit will automatically set up a layout that adheres to the common design? Here, we systemically investigate the arrangements of iso-orientation domains generated by the statistical connectivity model and assess their consistency with the experimentally observed common design invariants. First, we consider the statistical wiring model with perfectly hexagonal mosaics of RGCs, its most tractable form. We derive closed-form expressions for cortical neuron receptive fields and orientation domain layouts resulting from the Moiré interference effect of hexagonal ON and OFF ganglion cell mosaics [28,29]. The pinwheel density of these pinwheel layouts is ρ = 2 3 ≈ 3. 46, substantially larger than experimentally observed. We then characterize the orientation domain layouts resulting from spatially disordered hexagonal mosaics. We find that parameters of RGC position disorder can not be tuned such that the statistical wiring model’s layouts match the quantitative invariants of the common design. Next, we examine a generalized class of noisy hexagonal mosaics that allows for spatially correlated disorder of RGC positions. This correlated retinal disorder induces local variations in column spacing, mimicking column spacing heterogeneity in the visual cortex [47,48]. With these mosaics, Moiré interference persists to larger disorder strength. Pinwheel densities, however, are unaffected by low and intermediate levels of disorder and increase from a lower bound of 3. 5 for stronger disorder. Finally, we characterize the statistical connectivity model with RGC mosaics generated by Eglen’s pairwise interaction point process (PIPP), the most realistic model for RGC mosaics currently available [31,32,49]. The resulting arrangements of iso-orientation domains and pinwheels are identical to those predicted by Gaussian random field models [22,50,51]. Their pinwheel densities can be tuned by applying band pass filters of different bandwidths. However, for all plausible filter shapes, pinwheel densities are substantially larger than experimentally observed. Our findings demonstrate that the mechanism for seeding patterns of iso-orientation domains described by the stochastic wiring model predicts column arrangements substantially different from the long-range interaction model and distinct from the experimentally observed invariant common design.
Our overall goal was to assess whether the layout of orientation domains predicted by the statistical wiring model are consistent with the observed common design invariants. To achieve this, we first sought to establish a benchmark for models of orientation domain layouts in general, to which predicted layouts can then be compared. To this end, we re-analyzed the data set used in [23] using the fully automated method described in the same study. The data set contains optical imaging of intrinsic signal experiments from tree shrew (N = 26), ferrets (N = 82), dark-reared ferrets (N = 21) and galagos (N = 9). Because many previous studies used the statistical wiring model with parameters optimized to mimic the early visual pathway of the cat, e. g. [82], we additionally analyzed data from 13 cat V1 hemispheres. Following [23], we first computed the average pinwheel densities (Fig 1C). Pinwheel densities of all four species, including cat were statistically indistinguishable from each other and statistically indistinguishable from π (dark-reared ferrets excluded) —the value predicted for the average pinwheel density by the long-range interaction model [23]. As a measure of pinwheel position variability, spanning all scales from single hypercolumn to the entire imaged region, we calculated the standard deviation, SD, of pinwheel density estimates in circular subregions of area A (see Fig 1D for an illustration). For all species, the function SD (A) was well described by SD (A) = c ρ A γ (1) (Fig 1E) with ρ denoting the average pinwheel density. The variability exponents γ and variability coefficients c were similar in all four species (Fig 1F). As a measure of relative pinwheel positioning on the hypercolumn scale, we computed the nearest neighbor (NN) distance statistics for pinwheels of same or opposite topological charge as well as independent of their topological charge (see Fig 1D for an illustration). Distance distributions were unimodal and very similar (Fig 1G–1J). Importantly, the distributions obtained from cat V1 were indistinguishable from the other three species. Mean NN distances, when measured in units of hypercolumns, were statistically indistinguishable (Fig 1G–1J, insets). These findings confirm the results of [23,33] and show that cat primary visual cortex follows the same quantitative layout laws as in tree shrew, galago and ferret. From the above results, we extracted two types of consistency ranges that can be used as a benchmark for models of orientation domains in V1. To be consistent with an observed layout of orientation domains, a model’s predictions should not be significantly different from experimental observations in at least one species. We thus defined one species consistency ranges spanned by the minimal lower and maximal upper margin of the single species confidence intervals for each parameter. If a model’s predicted layout parameters are located outside one or more of the one species consistency ranges, data from every species rejects this model at 5% significance level. This criterion is thus conservative in nature and does not assume that there is in fact one species invariant common design. If such a truly universal common design for orientation domains in fact exists, it would be appropriate to pool data from different species and to consider the more precisely defined confidence intervals of the grand average statistics as the relevant benchmark. To perform this more demanding test of model viability, we also defined common design consistency ranges as the 95% bootstrap confidence intervals obtained from the whole data set. If the layout parameters predicted by a model are within all of the common design consistency ranges, the model offers a quantitative account of the bona fide universal common design. If one or more layout parameters have predicted values outside the common design consistency ranges the model is inconsistent with the common design. With the current data set, if a model is common design consistent, it is also one species consistent. One species (common design) consistency ranges are shaded in light (dark) green in Fig 1C, 1E–1J and summarized in Table 1. The tests of model viability defined above are most simply performed if the parameter values predicted by a model are determined exactly or with a numerical error that is much smaller than the empirical uncertainties. For models that can be solved accurately numerically, this can in principle always be achieved by a sufficiently large sample size of simulations. In the following, through analytical and numerical calculations, we will perform a comprehensive search through the statistical wiring model’s parameter space to identify regimes in which the model is one species consistent or common design consistent. The statistical wiring model formalizes the hypothesis that the spatial progression of orientation preference domains arises from the spatial distribution of RGC receptive fields on the retina via feedforward wiring. Fig 2A shows a simplified schematics of the early visual pathway in the cat [27,45,46,54–65], from the retina to layer IV of V1. A stimulus is focussed onto the retina through the cornea and lens, is sampled by RGC RFs and transmitted to the LGN. LGN neurons project to stellate cells in layer IV of V1, whose responses are orientation tuned. Orientation tuning varies smoothly across the cortical surface. In the model, RGCs are assumed to be mono-synaptically connected one-to-one to relay cells in the LGN. Thus, the receptive fields of LGN neurons are similar to those of RGCs and the spatial arrangement of ON/OFF receptive fields of relay cells in the LGN mirrors the RGC receptive field mosaic. Neurons in the model visual cortex linearly sum inputs of LGN neurons (Fig 2B). Their spatial receptive fields and orientation preferences are assumed to solely depend on the spatial arrangement of their afferent inputs. V1 neurons are assumed to receive dominant inputs from a small number of geniculo-cortical axons. Most of them sample from a single pair of ON/OFF RGCs, a so-called RGC dipole (Fig 2B). The neuron’s receptive field then consists of one ON and one OFF subregion and its response to edge-like stimuli is tuned to an edge orientation orthogonal to the RGC-dipole vector (Fig 2B). Within a mosaic of ON and OFF center RGCs, many such dipoles are present and the spatial arrangement of dipoles on the retina determines how tuning properties, e. g. the preferred orientation, change along a two-dimensional sheet parallel to the layers of the visual cortex. If ON and OFF RGCs are positioned on hexagonal lattices, the model predicts that a hexagonal pattern of orientation preference can arise through Moiré interference (MI) between the two lattices (Fig 2C). Following [28,83], we model RGC receptive fields using a Gaussian function GRFj (x) of width σr localized at the center position xj: GRF j (x) = ± exp - (x j - x) 2 2 σ r 2, (2) where x indicates position in retinal space. All subsequent results remain qualitatively unchanged if a biologically more realistic difference-of-Gaussians (see [84]) is used. A plus or minus sign in Eq (2) indicates an ON or OFF center cell, respectively. The receptive field RFy of a visual cortical neuron at position y in the two-dimensional cortical sheet is obtained by summing several ganglion cell receptive fields with positive synaptic weights wj: RF y (x) = ∑ j w j (y) GRF j (x). (3) The synaptic weights are chosen as w j (y) = exp - (x j - y) 2 2 σ s 2. (4) The parameter σs sets the range from which a V1 neuron receives retino-thalamic inputs, xj denotes the center of an RGC receptive field. According to Eq (3) the spatial distribution of RGC locations determines how response properties change across cortex. For σs smaller than the lattice spacing, each cortical cell receives substantial input only from a very small number of ganglion cells. Inputs received by most cortical cells are dominated by one ON and one OFF center RGCs (see inset in Fig 2A), forming an RGC dipole. The small σs regime is thus generally referred to as the dipole approximation of the model. While the dipole approximation leads to the robust emergence of simple-cell receptive fields with one (ON, OFF) or two (ON-OFF) subfields in the model V1 layer, it is worth mentioning that simple cells in cat and macaque monkey sometimes have more than two aligned, regularly spaced subfields (e. g. ON-OFF-ON or OFF-ON-OFF) (see [85,86]). In the dipole approximation of the statistical wiring model, such simple-cell RFs almost never occur. While the model as defined above implements a deterministic wiring scheme, it represents a simplification of a more detailed formulation of the statistical connectivity model proposed in [83]. In the more detailed formulation, the synaptic weights between the cortical units and the retina/LGN are chosen at random from a Gaussian distribution with the shape given in Eq (4). Ringach established in [83] that the spatial structure of the resulting domain layouts for the detailed and simplified model are nearly identical. We therefore refer to the model as statistical connectivity model. We used the linear response assumption [87,88] to determine cortical stimulus responses. A response R of a cortical neuron is modeled by the inner product between its receptive field RFy (x) and the stimulus, in our case an illumination pattern L (x): R y = ∫ d 2 x RF y (x) L (x). (5) Because Ry can become negative, a firing rate f of the cortical neuron is then defined through a static nonlinearity, e. g. half-wave rectification [87]. For the purposes of the present study, this nonlinearity can be neglected assuming that it does not alter core properties of the receptive field such as orientation preference and spatial frequency preference [82,83]. We derived close-form expressions for the pattern of cortical receptive fields across V1 that arises through Moiré interference in the case that ON and OFF center cells are localized on perfectly hexagonal lattices with different lattice constants r and r′ and relative angle α between the lattices (see Fig 2C). Detailed derivations are provided in Methods, along with closed-form expressions for receptive fields, the frequency response of orientation selective neurons, and their spatial organization. Fig 3A depicts the analytically calculated orientation preference pattern generated through Moiré interference between two hexagonal ON/OFF RGC mosaics. Iso-orientation domains are organized in fine-grained parcellations on small scales and repeat in a hexagonal pattern on a larger scale Λ. The larger scale is the predicted column spacing of the orientation domain layout (see Methods), and model parameters are chosen such that Λ ≈ 1mm as experimentally measured (see Methods and [28,29] for details). The scale of the small parcels therefore is < 200μm. A magnified view of a small region of the domain layout is provided in Fig 3B along with three analytically determined cortical receptive fields at closely spaced locations roughly 100μm apart from each other. These receptive fields highlight that individual parcels contain highly tuned units with vastly different preferred orientations. This means that orientation preference changes abruptly on scales < 200μm in the predicted patterns. Clearly, these features distinguish the obtained pattern of orientation preferences from the experimentally observed domain layouts. While orientation selectivity in V1 exhibits some small scale scatter within orientation domains [89,90], two-photon imaging suggests that orientation preferences progresses rather smoothly across the cortical surface [65,78]. Paik & Ringach implicitly assumed that random feedforward wiring from the retina/LGN to V1 effectively results in a smoothed version of the dipole layout (see Fig 2C). To extract this smooth pattern of orientation preferences from the statistical connectivity model, they adopted a two-step procedure to suppress the small-scale variation in the Moiré interference pattern: First, locations with orientation selectivity index (OSI) larger than a threshold value are determined [28,29,82,83]. Second, the orientation selectivity of all other location is set to zero. The resulting layout is then filtered with a Gaussian lowpass filter resulting in continuous and smooth array of iso-orientation domains [28,29,82]. We find that the thresholding/smoothing procedure effectively extracts the dominant lowest spatial frequency Fourier components of the Moiré interference pattern. As derived in Methods, the lowest spatial frequency contribution to the Moiré interference pattern consists of six Fourier modes with identical amplitude and wave number k c = 4 π 3 r r ′ r 2 + r ′ 2 - 2 r r ′ cos (Δ α), (6) Here, r, r′ denote to the lattice constants and Δα the angle between the hexagonal ON/OFF lattices (Fig 2C). The smooth orientation domain layout resulting from Moiré interference can therefore be summarized in a complex-valued field z (y) composed of six planar waves with wave numbers kj and fixed phase factors uj, z (y) = ∑ j = 1 6 exp (i k j y) · u j. (7) The pattern of preferred orientations across the cortical coordinate y is given by the phase of this complex-valued field as, ϑ pref (y) = 1 2 arg z (y). (8) Fig 3A (bottom) depicts ϑpref (y) as analytically determined. The pattern of pinwheels and iso-orientation domains is organized into a smooth hexagonal crystalline array. Interestingly, an identical layout of iso-orientation domains was constructed by Braitenberg et al. [91] based on an the idea that orientation preference is generated by discrete centers of inhibition in V1. It was also found by Reich et al. to solve a symmetry defined class of models for the self-organization of iso-orientation domains [92,93]. Fig 3C shows the differences in preferred angle between the unfiltered domain layout of the Moiré interference pattern and its low frequency contribution, together with a histogram of the differences. With Δ (x) = ϑ1 (x) − ϑpref (x), the difference d (x) between the two preferred angles is defined as d (x) = 1 2 abs (arg (e 2 i Δ (x) ) ). The bimodal shape of the histogram indicates that the orientation preference of a large fraction of cortical locations differs substantially between unfiltered and smoothed layout. Roughly one fifth of all locations exhibit differences of orientation preferences of more than 45°. To compare our mathematical expression for the column spacing of the orientation domain layout to previous results, Eq (6) can be rewritten by introducing a parameter β representing the detuning between the two lattice constants in units of the lattice constant r′ → (1 + β) r. The expression for the column spacing becomes Λ c ≡ 2 π k c = 3 2 · S · r, (9) where S is the distance between two vertices of the Moiré pattern in units of r, called the scaling factor [94–96] S = 1 + β β 2 + 2 (1 - cos (Δ α) ) (1 + β). (10) The difference between S ⋅ r and Λc is displayed in Fig 3A (bottom). Eqs (9) and (10) are identical to previous results for the spacing of hexagonal Moiré patterns derived via geometrical considerations [95,96]. Using these explicit expressions for the iso-orientation domain layout and its column spacing Λc, we first evaluated the central quantity of the common design—the pinwheel density, i. e. the number of pinwheels per unit area Λ c 2. Within each unit cell of area A = 3 2 (S · r) 2, there is one “double pinwheel” of topological charge 1, around which each orientation is represented twice, and two pinwheels of topological charge ± 1 2. With Λ c 2 = 3 4 (S · r) 2 and counting the pinwheel with charge 1 as two pinwheels (see below), the pinwheel density is ρ = 2 + 2 · 1 · Λ c 2 A = 2 3 ≈ 3. 46. (11) Notably, this value is outside of both, the common design consistency range and the single species consistency range for the experimentally measured pinwheel densities (cf. Table 1). Since the statistical connectivity model for perfectly hexagonal RGC mosaics results in a periodic array of pinwheels, all three nearest neighbor distance distributions of pinwheels are sharply peaked (see also Supplementary Material of [23]) and, thus, in disagreement with the distributions experimentally observed (cf. Fig 1). We compared these analytical results to numerically evaluated Moiré interference patterns (Fig 4). The fine-grained layouts of numerically and analytically obtained unfiltered layouts are almost indistinguishable (cross-correlation coeff. 0. 9, Fig 4A top). This confirms the analytical treatment and indicates accuracy of the numerical implementation. A hierarchy of discrete spatial frequency contributions is apparent in amplitude spectra of both domain layouts (Fig 4A (bottom) ). The peaks at larger spatial frequencies in Fig 4A and 4B are localized at 3 k c as analytically predicted (see Methods). To numerically generate the smoothed array of orientation domains, the layout in Fig 4A (top) was thresholded (OSI > 0. 25, see Methods) and subsequently smoothed with a Gaussian lowpass filter (Fig 4B top) [28,82]. In general, strongly tuned locations are those exactly between ON-OFF RGC pairs (Fig 4B, inset). Fig 4C depicts the crystalline pinwheel arrangement of the analytically calculated smoothed layout (Eq (8) ) as well as the six dominant low frequency Moiré modes in the amplitude spectrum. While the numerically obtained layouts and its analytical approximation are similar (cross-correlation coeff. 0. 6), one major difference can be observed: the pinwheel of topological charge 1 is replaced by two pinwheels of topological charge 1 2 in the numerically obtained layouts, along with subtle deformations of adjacent orientation domains (compare Fig 4B and 4C). To see why this is the case, we note that the pinwheel with charge 1 in the analytically calculated pattern (Eq (8) ) arises from a zero of the field z (x) (Eq (7) ) with multiplicity two. A phase singularity of a complex-valued field arising from a zero with multiplicity N > 1 is structurally unstable and unfolds upon generic infinitesimal perturbations into N closely spaced singularities of multiplicity one [97]. The numerical procedure of discretizing V1 unit positions on a numerical grid, OSI thresholding, and smoothing realizes such a perturbation and this explains why in the numerical solutions the pinwheel of charge 1 unfolds into two adjacent pinwheels of charge 1 2. So far, we have studied the idealized situation of iso-orientation domains induced by perfectly ordered hexagonal RGC mosaics. RGC mosaics in the eye, however, are not perfectly hexagonal but exhibit substantial spatial irregularity [26,98]. Therefore, we next turned to numerically investigate the statistical connectivity model with hexagonal RGC mosaics subject to Gaussian disorder in RGC position as previously described [28,29]. The effect of ganglion cells displaced by Gaussian distributed offsets with standard deviation σ = η ⋅ r is illustrated in Fig 5. The parameter r is the lattice constant and η is the disorder strength. Fig 5 shows the unfiltered orientation domain layout (far left), the layout thresholded for cells with an OSI > 0. 25 (left), the smoothed thresholded layout (right) as well as its amplitude spectrum (far right), numerically obtained for η = 0. 12. As in the perfectly ordered case, the unfiltered layout of the noisy Moiré interference model exhibits a substantial scatter of orientation preferences across small scales. For a disorder strength of η = 0. 12 (Fig 5A), the domain layout is still dominated by the six lowest spatial frequency Moiré modes also present in the perfectly ordered system. For a disorder strength of η = 0. 3 (Fig 5B), the amplitude spectrum (Fig 5B, far right) lacks any indication of theses Moiré modes indicating that Moiré interference no longer takes place. As a consequence the resulting layouts of iso-orientation domains lack a typical column spacing. To characterize the model orientation domain arrangements, we first calculate amplitude spectra for both, unfiltered and smoothed layouts (Fig 5A and 5B), | R (k) | = ∫ d 2 x z (x) e i k x where z (x) = e 2 i ϑ pref (x). (12) Normalizing and radially averaging yields the so-called marginal amplitude spectrum (Fig 5C and 5D), f (k) = ∫ 0 2 π d ϑ | R (k cos (ϑ), k sin (ϑ) ) | max k ∫ 0 2 π d ϑ | R (k cos (ϑ), k sin (ϑ) ) |. (13) The sharp peak at kc corresponds to the dominant Moiré mode indicating that orientation domain layouts exhibit a typical column spacing. For increasing disorder, the relative levels of peak height to background decreases while the peak width remains small. As expected, marginal amplitude spectra of unfiltered and the smoothed layout mainly differ in the strength of background components. The flat amplitude spectrum of the unfiltered iso-orientation domain layouts for large disorder strength is transformed into a Gaussian amplitude spectrum by the lowpass filtering. Based on this assessment, the disorder strength η has to be smaller than 0. 3 to ensure that layouts exhibit a typical spacing between adjacent iso-orientation domains. We next systematically evaluated the core layout parameters of the common design—pinwheel density and pinwheel nearest neighbor distance distributions for the statistical wiring model with disordered hexagonal RGC mosaics. To compare the model predictions with experiments, we estimated the column spacing of the model orientation domain layouts as well as pinwheel layout parameters using the exact same methods that we applied to the experimental data (see Methods). For weak disorder, column spacing estimates closely match the theoretical prediction Λc (Fig 6A), confirming the accuracy of the wavelet method. For disorder strengths larger then 0. 12, Moiré modes are no longer the dominant spatial frequency contribution in the model layouts and the estimated column spacing increases with disorder strength. Having estimated the column spacing, we analyzed model orientation domain layouts with respect to the common design parameters (Fig 6B–6D). As expected, pinwheel densities approach the analytical predicted value of 2 3 for weak disorder (Fig 6B) and increase with increasing disorder strength. This increase is largely caused by the increase in the estimated column spacing (Fig 6A) and does not involve a massive generation of additional pinwheels for larger disorder strength. We next calculated the standard deviation of pinwheel densities as a function of the area A of randomly selected subregions of the iso-orientation domain layouts. Generally, the standard deviation’s decay with subregion size followed a power law with increasing area size, with larger exponents for weak disorder (Fig 6D). Fig 6E and 6F show a complete characterization of pinwheel nearest neighbor (NN) distance distributions of the noisy Moiré interference model. Histograms for NN distances for arbitrary charge are bimodal for weak disorder (Fig 6E). The peak at smaller NN distances results from the unfolding of pinwheels with topological charge 1 into two adjacent pinwheels of topological charge 1/2 for finite disorder strength (see above and Figs 3 and 4). For the same reason, the NN distance histogram for pinwheels of identical topological charge is also bimodal (Fig 6F). With increasing disorder strength, both distributions become unimodal (Fig 6E and 6F left). The NN distance distribution for pinwheels of opposite sign is unimodal for all parameter values, indicating that only very few additional piwheel pairs are added to the pinwheels of the Moiré layout for small and intermediate disorder strengths (Fig 6G). The overall decay in mean NN distance for strong disorder in all three histograms mostly reflects the increase in measured column spacing (see Fig 6A). Based on these results, we attempted to identify a disorder strength for which all NN pinwheel distance distributions resembled the experimental data. To this end, we calculated the squared error between the calculated histograms and the experimental data as a function of disorder strength (Fig 6H). Smallest deviations from experimental data were obtained around η ≈ 0. 11 for all three NN distance distributions. Fig 7 summarizes all common design features determined for disordered Moiré interference model layouts as a function of disorder parameter and compares them to the experimentally observed values in tree shrew, galago, ferret, dark-reared ferrets, and cats. Light (dark) green shaded areas indicate the single species (common design) consistency ranges (see Fig 1, cf. Table 1). With increasing disorder, pinwheel density of model layouts steadily increases from ρ = 2 3 ≈ 3. 46 (Fig 7A) and always lies above the single species consistency range. Thus, the pinwheel density of model orientation domain layouts is inconsistent with pinwheel densities observed in all species. Next, we fitted the empirically observed power law, Eq (1), to the standard deviation of the pinwheel density estimate in increasing subregions of area A (see Fig 6D) [23]. The variability exponent γ is consistent with experiments for disorder levels exceeding η = 0. 15. The variability constant c is monotonically increasing up to η ≈ 0. 15 at which point the model domain layouts lose their typical column spacing (cf. Fig 5) and the increasing pinwheel density ρ causes a drop (Fig 7C). Fig 7D–7F displays the mean pinwheel NN distances as function of the disorder strength, all of which substantially decrease with increasing disorder strength. This can be attributed to the increasing mean column spacing of the domain layouts under increasing disorder (see Fig 6A). Mean NN distances for weak disorder strength are close to the experimental data, but NN distance distributions for pinwheels of different topological charge and independent of topological charge are bimodal for weak disorder (Fig 6E and 6F). The latter is clearly distinct from the experimental data (cf. Fig 1, [23]). Fig 7G shows an overview of the consistency of model orientation domain layout parameters with the data for various disorder strengths. As can be seen, no strength of disorder results in layouts that are consistent with the common design for all layout parameters. Perhaps, even more surprising, pinwheel density and NN distance for pinwheels of the same sign are inconsistent with the individual values obtained for each species, no matter how the strength of disorder is chosen. The above results show that the statistical connectivity model with hexagonal mosaics is unable to reproduce all features of the common design, even if spatially uncorrelated position disorder is imposed on the RGC positions. Whatever the source of disorder that causes the irregularity in the RGCs’ positions, it is plausible to assume that it is correlated on scales spanning several RGCs. Such spatial correlations would preserve the Moiré effect locally, yet generate spatial irregularity in orientation domain layouts. To test whether correlated positional disorder can produce model arrangements of orientation domains that match experimental observations, we generalized the noisy hexagonal mosaics proposed in [28,83] to include spatial correlations. To obtain noisy hexagonal RGC mosaics with spatial correlations, we started with a hexagonal array of RGC positions. The position of each lattice point xi was then shifted depending on its position according to xi → xi + η y (xi). The shift η y (x) with amplitude η for y (x) = (y1 (x), y2 (x) ) was chosen from a Gaussian random field with vanishing mean 〈y1 (x) 〉 = 〈y2 (x) 〉 = 0, fixed standard deviation std (y1 (x) ) = std (y2 (x) ) = 1 and correlation function 〈 y (x 1) y (x 2) 〉 = 2 exp (- | x 1 - x 2 | 2 2 σ 2) with correlation length σ (see Methods) where y1 and y2 are statistically independent. The two parameters, correlation length σ and amplitude η were expressed in units of the lattice constant r. Fig 8A and 8B illustrates this procedure. RGCs are shifted in a coordinated manner across the plane, correlated in both direction and magnitude of the shift. The determinant of the Jacobian det J (x) = det ∂ y 1 (x) ∂ x 1 ∂ y 1 (x) ∂ x 2 ∂ y 2 (x) ∂ x 1 ∂ y 2 (x) ∂ x 2, (14) measures the local change of RGC lattice constant. In regions of negative det J, RGCs are closer together than average, in regions of positive det J, RGCs are further apart (contour lines in Fig 8B). In primates, regions of higher density are predicted to have higher cortical magnification and vice versa [12]. The RGC mosaics with correlated positional noise therefore imply local fluctuations in the cortical magnification factor on the scale of the noise correlation length. Fig 8C-8E display unfiltered and smoothed model layouts obtained with spatially correlated noisy hexagonal mosaics as well as their amplitude spectra. As expected, orientation domain layouts exhibit a typical column spacing up to higher disorder strengths, when the position disorder was correlated (compare Fig 8F with Fig 5D). Locally, Moiré interference leads to a roughly hexagonal layout of columns that is distorted on larger scales. Both, the orientation of the hexagons as well as column spacings change continuously across the layout. For weak disorder, the amplitude spectrum still exhibits six peaks, indicating a globally hexagonal layout (Fig 8C, right). For intermediate disorder local column spacing and direction of the hexagons varies to the extend that peaks can hardly be identified in the amplitude spectrum of the resulting domain layout. In particular, the spatially varying local column spacing leads to a broader peak in the radially averaged amplitude spectrum with increasing disorder strength (Fig 8F). This is in contrast to the case of uncorrelated disorder (cf. Fig 5D). Note that experimental iso-orientation domain layouts exhibit a similarly broad peak in their marginal amplitude spectra [99]. We quantified the pinwheel density of orientation domain layouts obtained with correlated noisy hexagonal RGCs (Fig 8G) as a function of disorder correlation length and disorder strengths. Independent of the disorder correlation length, pinwheel densities plateau around 2 3 for weak disorder and monotonically increase above a critical disorder strength. This critical disorder strength is higher, the larger the correlation length. Thus, model pinwheel densities are inconsistent with the individual values obtained for each species, no matter what the strength of disorder or correlation length is. Fig 8H illustrates that the pinwheel density increases with increasing disorder strength largely because the overall measured column spacing increases, not because additional pinwheels appear in the layouts. In fact, the number of pinwheels per mm2 is almost independent of either correlation length or disorder strength. Finally, we examined whether the statistical connectivity model could reproduce the common design invariants with RGCs distributed in space according to a pairwise interacting point process (PIPP). The PIPP developed by Eglen et al. [49] is currently the experimentally best supported model for RGCs mosaics and was shown by several studies to generate RGC positions which accurately reproduce a variety of spatial statistics of RGC mosaics [31,32,49,82]. The PIPP model generates samples from a statistical ensemble of RGC mosaics by iteratively updating RGC positions to maximize a target joint probability density, specified by pairwise interactions between neighboring RGCs (for details see Materials & Methods). Each PIPP mosaic represents a random realization of a regularly-spaced RGC mosaic with radially isotropic autocorrelograms [31] and lacks long-range positional order. Fig 9A depicts a realization of a PIPP with parameters choosen to reproduce cat RGC mosaics (for details see Materials & Methods). We generated thresholded and smoothed iso-orientation domain layouts from PIPP RGC mosaics as from the ordered mosaics (Fig 9A and 9B left). The lack of long-range positional order in ON and OFF mosaics prevents any Moiré interference between them. Thus, no typical spacing between adjacent columns preferring the same orientation is set in the model layouts (Fig 9B, right, see also [31,32]). Spectral power in these layouts is broadly distributed and monotonically decays with increasing spatial frequency. As a consequence, one needs to apply bandpass filtering to obtain orientation domain layouts that are at least qualitatively resembling the experimental data. We used a flexible band-pass filtering function f (k) of the following form to the amplitude spectrum: f (k) = a | k | β exp (- | k | 2 b), (15) with a, b, β > 0. The filter function was normalized such that ∫ d 2 k f (k) = 2 π. (16) With this normalization, one can define the mean column spacing of the resulting layout via Λ = 2 π / k ¯ with k ¯ = ∫ 0 ∞ d k 2 π k f (k) = 1. (17) By increasing the parameter β, the shape of the filter can be changed from Gaussian lowpass (β = 0, Fig 9C, left) to wide bandpass (β = 2, Fig 9C, middle) and to narrow bandpass (β = 10, Fig 9C, right). There is an additional degree of freedom in this filter definition, namely how the filter is scaled relative to the absolute physical units mm−1 of the amplitude spectrum. To scan a wide range of filter shapes and column spacings, we varied β between 0 and 10 and choose the scaling such that Λ varied between 0. 6 mm and 1. 2 mm, i. e. covering the entire range of experimentally observed mean column spacings in tree shrew, galago, cat, and ferret [23,33]. We then measured the pinwheel densities of the resulting statistical connectivity model layouts (Fig 9E, right), where pinwheel density was defined as the number of pinwheels within an area Λ2. Pinwheel densities were independent of the scale Λ and increased monotonically with increasing spectral width (decreasing β). They are in general substantially larger than the experimentally observed value of 3. 14 (see also Fig 9C) and outside of the single-species/common-design consistency range. Qualitatively, iso-orientation domain layouts generated with the PIPP RGC mosaics resembled those generated from Gaussian random field (GRF) [22,50,51] (Fig 9C). In fact, we find that this resemblance is quantitative. Fig 9E depicts the analytical prediction [51] for the pinwheel density of orientation domains obtained with GRFs with a marginal amplitude spectrum corresponding to the filter function in Eq (15) (Fig 9E, left). Pinwheel densities for GRF layouts and layouts obtained from PIPP mosaics with the statistical connectivity model are indistinguishable. For the pinwheel density to be consistent with at least the single-species consistency range (ρ < 3. 42), amplitude spectra had to be much more peaked (β ≥ 17) than experimentally observed [99]. Finally, we filtered statistical connectivity model layouts with the Fermi-Filter function as used in [23,33] with cut-off wavelengths of 0. 3 mm and 1. 2 mm (Fig 9D). Again, pinwheel densities were much larger than those observed in the experimental data and outside of the single-species and common-design consistency range. In summary, iso-orientation domain layouts generated by the statistical connectivity model using PIPP RGC mosaics quantitatively resemble layouts derived from Gaussian random fields. Their statistics is distinct from the statistics of experimentally measured layouts.
The statistical wiring model analyzed in the present study is only one representative of possible random wiring schemes. One could argue that alternative, perhaps more realistic, schemes might do a better job at reproducing the experimentally observed pinwheel layouts. There is good evidence that the spatial statistics of RGC mosaics is well approximated by Eglen’s PIPP [31,32,49], and, hence, there is little freedom of choice at the retinal level. In contrast, at the next network layer, two main modifications or extensions of the statistical wiring model could be considered: (i) adding an additional layer to the feedforward network implementing the transformation of the retinal input structure by the lateral geniculate nucleus (LGN) (ii) choosing different probabilistic connectivity rules between the retinal/LGN layer and the primary visual cortex. We argue that both modifications of the random wiring approach are unlikely to improve the consistency of the model with experimental data. Regarding (i), Martinez et al. [100] have recently tried to infer the mapping between RGC inputs and LGN relay cells using a statistical connectivity approach. In their model, ON and OFF cell types were homogeneously distributed and their polarity (ON or OFF) was inherited from the nearest retinal input. Connection probability between RGCs and LGN neurons was modeled as an isotropic Gaussian function of the relative distance between the RF centers of the presynaptic and postsynaptic partners. With this simple wiring scheme, together with similar connectivity rules for the population of inhibitory interneurons, several spatiotemporal properties of LGN RFs robustly agreed with the experimental data. In the architecture between the retina and the LGN proposed by these authors, the dipoles of ON- and OFF-center cells that characterize the retinal mosaic are transformed into small clusters of same-sign relay cells. The LGN ON-OFF dipoles occur at the boundaries of these clusters with LGN dipole orientations strongly correlating but not necessarily matching dipole orientations in the retina. Notably, dipole density and dipole angle correlation length in the LGN is not increased compared to the retina. The data and modeling by Martinez et al. suggest that the LGN mosaics do not systematically alter the spatial structure of RGC inputs beyond providing an additional source of dipole disorder. Additional disorder imposed by the LGN mosaics would likely add a uniform level of disorder to all spatial frequency components in the unfiltered orientation domain layouts predicted by RGC dipole structure. When hexagonal mosaics are considered, the disorder strength that has to be assumed to match the spatial distribution of RGC cells found in experiment is already rather high [28,29,31,32,83]. Additional noise is likely to obstruct any remaining Moiré interference. We therefore speculate that when considering an additional LGN layer, after smoothing, domain layouts for both, the noisy Moiré interference model and the model with PIPP mosaics would be similar to those obtained with PIPP mosaics [31,32,82]. As we have shown in the present study, the spatial statistics of these layouts resembles those derived from Gaussian random fields and is inconsistent with the data obtained for any of the four species analyzed (cf. Fig 9, see also [23,50]). A similar argument can be made for alternative probabilistic connectivity rules between the retinal/LGN layer and V1. As our analysis shows, the dipole structure emerging from “realistic” RGC mosaics (be it very noisy hexagonal mosaics or PIPP mosaics) is spatially fine-grained because dipole angles vary over short distances in cortical space relative to the typical size of an iso-orientation domain. For this reason, the statistical connectivity model requires an additional smoothing step (cf. Figs 4,5 and 9) to yield smooth orientation domain layouts as observed in experiment [65,78]. Unless the connectivity rule is assumed to specifically select dipoles with a similar angle from a larger spatial region of the retina, or neurons within an iso-orientation domain are assumed to choose one particular dipole to receive the input from and ignore all other dipoles in the vicinity, such spatial averaging within the cortical layer will always be required no matter what the actual probabilistic connectivity rule is. Domain layouts resulting from such spatial averaging of weakly correlated dipole angles (see also [32]) are likely to follow layout statistics that resemble those of Gaussian random fields, independently of the connectivity rule assumed. If neurons are assumed to select specific dipoles out of the repertoire of “available” ones, then the overall spatial layout of RGC dipoles is not informative about the resulting domain layout, which contradicts the main hypothesis of the statistical wiring model. Ultimately, the key experiment to provide support for the random wiring approach consists of determining both, the orientation domain layout and the retinotopic map in a single animal and, in a second step, correlate these with the spatial arrangement of RGCs in the same animal. This challenging experimental task is still awaiting its completion. So far, the only model class able to robustly reproduce all common design parameters, describes the formation of orientation domain layouts as a deterministic optimization process converging to quasi-periodic pinwheel-rich orientation domain patterns associated with and stabilized by a matching system of intrinsic horizontal connections [23,101,102]. Irregular layouts of orientation domains dynamically emerge as a consequence of large-scale circuit optimization of domain patterns and intrinsic circuits. Is this agreement between model and data good evidence for global circuit optimization or are there simple alternative explanations such as the random feedforward wiring hypothesis that can explain the invariant statistical properties of orientation domains? Qualitatively, it is in fact tempting to attribute the spatial irregularity and apparent randomness of pinwheel layouts in V1 to some general kind of “biological noise”. In this view, the quantitative laws of pinwheel organization that Kaschube et al. found [23] would then be conceived as outcome of a largely random process underlying the emergence of orientation domains. By now, however, all proposals based on the assumption of disorder as the determinant of spatial irregularity have failed to reproduce the common design parameters and laws that have now been observed in four divergent species. Orientation domain layouts obtained from statistical ensembles of Gaussian random fields [22,50,51] as well as phase randomized layouts derived from experimental data [23], exhibit pinwheel densities that are substantially higher than experimentally observed. Importantly, most dynamical models for the development of orientation domains produce such Gaussian random domain layouts during the initial emergence of orientation selectivity [22,50]. Therefore approaches based on “frozen” early states of such models are also ruled out by the existing data (see [103]). The present study shows that a mechanistic and biologically plausible feedforward model of the early visual pathway based on (i) noisy hexagonal placement of RGCs or (ii) a more realistic semi-regular positioning of RGCs generated by the PIPP also generates layouts distinct from experimental observations. These findings illustrate that orientation domains and pinwheels positions, although spatially non-periodic and irregular, follow a rather distinct set of layout laws. These laws cannot easily be accounted for by a spatial irregularity or randomness in the structure of afferent projections to visual cortical neurons. A further conceivable and potentially critical source of stochasticity that is often overlooked is randomness within intracortical circuits. The field approach employed in various models for orientation domain layouts, such as the long-range interaction model, represents an idealization of a complex network, in which every neuron is characterized by its own set of inputs and outputs. These inputs and outputs may, at least to some extent, be stochastic. How and to what extent randomness in intracortical connections can affect and shape orientation domain layouts is not well understood. In that respect, it is interesting to note that model networks for largely stochastic intracortical circuits are able to generate and robustly maintain orientation selective responses to afferent inputs and can lead to highly coherent orientation domains [104,105]. Paik & Ringach reported indications of hexagonal order in visual cortical orientation domain patterns of tree shrew, ferret, cat, and macaque monkey [28,29]. Two recent studies have casted doubt on the hypothesis that this hexagonal order echoes hexagonal or quasi-hexagonal arrangements of ganglion cell mosaics in the retina [31,32]. Hore et al. showed that noisy hexagonal lattices do not capture the spatial statistics of RGC mosaic. Moreover, the positional correlations in measured mosaics extends to only 200–350μm, far less than required for generating Moiré interference [31]. More generally, Schottdorf et al. studied the spatial arrangement of RGC dipole angles in cat beta cell and primate parasol RF mosaics [32]. According to the statistical wiring hypothesis, dipole angle correlations should follow the spatial correlations of preferred orientations in the primary visual cortex, i. e. be positively correlated on short scales (0–300 μm) and negatively correlated on larger scales (300–600 μm) in the retina. By introducing a positive control point process that (i) reproduces both, the nearest neighbor spatial statistics and the spatial autocorrelation structure of parasol cell mosaics and (ii) exhibits a tunable degree of spatial correlations of dipole angles, they were able to show that, given the size of available data sets, the presence of even weak angular correlations in the data is very unlikely. If not from the structure of RGC mosaics, where does the apparent hexagonal organization in orientation domains come from? A variety of self-organization models on all levels of biological detail have been shown to generate orientation domains with hexagonal arrangements, e. g. [14,92,93,103,104,106,107] (notably including the earliest theory for the self-organization of orientation preference by von der Malsburg in 1973). Thus, hexagonal order, even if present, would not provide specific evidence in favor of Moiré interference between RGCs. Future work will have to elucidate whether the long-range interaction model for orientation domains [23,101,102] can not only explain the common design statistics, but at the same time account for the observed hexagonal order in the visual cortex. Compared to the dense sampling of stimulus space by cortical neurons, the repertoire of detectors on the retina that input into a given cortical area is limited. For the cat visual pathway, Alonso et al. estimated the number of LGN X-relay cells converging onto a single simple cell in V1 to be ∼ 20–40 [27], based on measuring the probability of finding a connection between individual geniculate and cortical neurons with overlapping receptive fields. This estimate was later confirmed by directly measuring population receptive fields of ON and OFF thalamic inputs to a single orientation column [108]. With an expansion of around 1. 5–2. 0 from X-cells in the retina to X-relay cells in the LGN [109,110], each simple cell in V1 receives on average input from only ∼ 10–25 RGCs. This not only implies that random afferent inputs to cortical neurons might seed groups of V1 neurons with similar orientation preferences but also that they might in fact impose substantial biases on the preferred orientation that can be adopted by the cortical neurons. The postnatal development of orientation columns could then be imagined as a dynamical activity-dependent process which refines and remodels an initial set of small biases provided by the RGC mosaic model through Hebbian learning rules and other mechanisms of synaptic plasticity. The up to now most striking experimental evidence that retinal organization can impose local biases in V1 function architecture was revealed by the finding that the pattern of retinal blood vessels can specifically determine the layout of ocular dominance columns in squirrel monkey [111] (for a modeling study see [112]). Dynamical models of orientation column formation generally assume no a priori constraints or biases as to which preferred orientation a given position in the cortical surface can acquire. Usually random initial conditions determine which instance from the large intrinsic repertoire of stable potential domain layouts is adopted. Including seeds and biases derived from RGC mosaics in such models for the dynamical formation of V1 orientation domain layouts may elucidate the potentially complex interplay between a sparse set of subcortical feedforward constraints and self-organization in a dense almost continuum-like intracortical network. The present study provides a detailed description of a candidate set of such subcortical biases and, therefore, can serve as a foundation for such future investigations. The common design invariants comprise four distinct functions in addition to the apparently invariant pinwheel density. As such, they represent a rather specific quantitative characterization of orientation domain layouts. It is, thus, not surprising that entire model classes have been rejected based on whether their predictions match these invariants. Since the discovery of visual cortical functional architecture more than fifty years ago, a large number of models based on a variety of circuit mechanisms has been proposed to account for their postnatal formation (see [113–115] for reviews). Many of these models are explicitly or implicitly based on optimization principles and attribute a functional advantage to the intriguing spatial arrangement of orientation domains. Because many, even mutually exclusive, models could qualitatively account for main features of orientation domain layouts such as the presence of pinwheels or the roughly periodic arrangement of columns, theory could not provide decisive evidence on whether to favor one hypothesis over another. It is only in recent years that the large available data sets have started to allow for a rigorous quantitative analysis of visual cortical architecture and its design principles in distinct evolutionary lineages. To date, abstract optimization models have been analyzed most comprehensively, providing in many case the complete phase diagram. For instance, Reichl et. al. systematically evaluated energy-minimization-based models for the coordinated optimization of orientation preference and ocular dominance layouts [92,93]. By quantitatively comparing model solutions to the common design, they were able to rule out a whole variety of otherwise intuitive and plausible principles for their emergence. It is desirable to obtain a similarly quantitative understanding of more detailed models for the formation of orientation domains. In this regard, the analysis of abstract models is informative because there is a many-to-one relationship between detailed models of the visual cortical pathway and those abstract formulations. Abstract models often can be shown to be representative of an entire universality class and, once comprehensively characterized, the questions becomes whether more complex modeling schemes are simply complicated instantiations of such a class. For models of an intermediate degree of realism, semi-analytical perturbation methods can be employed to explicitly provide this mapping. Using this approach, Keil & Wolf studied orientation domain layouts predicted by a widely used representative of a general optimization framework [103]. According to this framework, the primary visual cortex is optimized for achieving an optimal tradeoff between the representation of all combinations of local edge-like stimuli, i. e. all positions in the visual field and all orientations, and the overall continuity of this representation across the cortical surface [116]. While this framework has successfully explained a variety of qualitative aspects of orientation domain design, e. g. [116–118], the authors found quantitative disagreement with the common design in all physiologically realistic parameter regimes of the representative model [103]. Their analysis enabled an unbiased comprehensive search of the model’s parameter space for a match to the experimental data and indicated alternative more promising optimization hypotheses to explain the experimentally observed V1 functional architecture. Although the statistical wiring model is still rather simplistic, it is hard to make analytical or semi-analytical progress as soon as RGC mosaics with the necessary degree of realism are considered. In this case, the question of whether models account for the cortical architecture can only be answered with the approach we have pursued here, i. e. by systematic comparison between their solutions, experimental data, as well as predictions from minimal approaches. The results presented here show that this approach can indeed be successfully applied to rule out candidate mechanisms as sufficient explanations for the emergence of V1 functional architectures. We expect a re-examination of the quantitative predictions of other modeling approaches to be highly informative about candidate mechanisms for the formation of V1 functional architecture. This present study provides the first systematic assessment as to whether the common design of orientation domains could result from an inherently random process, as realized through the local feedforward structure of the early visual pathway rather than an optimization process coordinated on large scales. Given the disagreement between the layouts predicted by the statistical wiring model and the data, global circuit optimization as proposed by the long-range interaction model currently is the only theory known to be capable of explaining the common design of orientation domains in the primary visual cortex.
Column spacing and pinwheel statistics of both data and simulation were analyzed using the wavelet method introduced in [47,119]. This method specifically takes into account that experimentally measured domain layouts often exhibit local variations in column spacing (as opposed to most model layouts) and is thus well-suited to unbiasedly compare the pinwheel layouts of model layouts and experimental data. Matlab source code for preprocessing of experimentally obtained layouts, column spacing analysis, and the analysis of pinwheel layouts can be found in the Supplemental Material, along with four example cases from ferret V1 to test the code. The full data set used in the present study is available on the neural data sharing platform http: //www. g-node. org/. For comparison between model orientation domain layouts and experimentally obtained layouts, both were analyzed with the exact same wavelet parameters settings. Raw difference images obtained in the experiments were Fermi bandpass filtered as described in [23]. Filter parameters were adapted to the column spacings of the different species such that structures on the relevant scales were only weakly attenuated (see [23]). To determine the local column spacing of the layouts, we first calculated wavelet coefficients of an image I (x), averaged over all orientations Ψ (y, Λ) = ∫ d φ π ∫ d 2 x I (x) · ϕ y (x, Λ, φ) (18) where y is the position, φ the orientation and Λ the scale of the wavelet ϕy (x, Λ, φ). We used complex-valued Morlet wavelets composed of a Gaussian envelope and a plane wave ϕ (x) = 1 σ exp - x 2 2 σ 2 · exp (i k ϕ x) (19) and ϕ y (x, Λ, φ) = ϕ (Ω - 1 (φ) (y - x) ). (20) The matrix Ω (φ) is the two dimensional rotation matrix (Eq (25) ). To compute the wavelet orientation average in Eq (18), 16 equally spaced wavelet orientations were used. For a given Λ, the parameters of the Morlet wavelet were chosen as k ϕ = 2 π Λ 1 0 (21) σ = ξ Λ 2 π. (22) ξ determines the size of the wavelet and was chosen to be ξ = 7, as in [23]. This captures column spacing variations on scales larger than 4 hypercolumns while at the same time enabling robust column spacing estimation. To obtain the map of local column spacing Λlocal (y), we calculated the scale Λ with the largest wavelet coefficient Λ local (y) = argmax Λ Ψ (y, Λ) (23) for every position y. To estimate the pinwheel density and other pinwheel layout parameters, we used a fully automated procedure proposed in [23]. We refer to the Supplemental Material accompanying [23] for further details. Here, we examine the analytically most tractable variant of the statistical wiring model, in which ON and OFF center cells are localized on perfectly hexagonal lattices L, (Fig 2B), that may exhibit different lattice constants r and r′, L = 1 0 k + 1 2 1 3 l f ∀ k, l ∈ Z, (24) where f = r, r′ is the lattice constant. Describing a rotation of the lattice vectors by the rotation matrix Ω (α) = cos (α) - sin (α) sin (α) cos (α), (25) the ON mosaic is rotated by an angle α, the OFF lattice by an angle α′ (Fig 2B). Paik & Ringach found in numerical simulations that in this case, Moiré interference between two hexagonal RGC mosaics results in a hexagonal layout of orientation domains [28,29,83]. We now first derive an explicit expression for cortical receptive fields RFy spatially varying with y predicted by the model. Calculating the preferred orientation of these receptive fields then provides an explicit expression for the hexagonal domain layouts. The sum in Eq (3) can be evaluated analytically for rectangular lattices using Jacobi Theta functions [120]. To solve the Moiré interference model, we used the fact that every hexagonal lattice can be written as sum L = L 1 + L 2 of two rectangular lattices with orthogonal base vectors by separating even and odd numbers in l and shifting the l-sum so that the x-component is equal to zero. The two rectangular lattices are L 1 = 1 0 k + 0 3 l f ∀ k, l ∈ Z, L 2 = 1 0 k + 0 3 l + 1 2 1 3 f ∀ k, l ∈ Z. (26) Evaluating the infinite Gaussian sum (Eq (3) ) yields the result for a single sub lattice (either ON or OFF) RF α, r, y ON/OFF (x) = T Θ 3 b e ϕ, τ Θ 3 b e r, κ + Θ 4 b e ϕ, τ Θ 4 b e r, κ (27) where b = x σ s 2 + y σ r 2 σ s 2 + σ r 2 (28) e r = - π r cos (α) sin (α) (29) e ϕ = - π 3 r - sin (α) cos (α) (30) τ = e - 2 π 2 σ r 2 σ s 2 3 r 2 (σ r 2 + σ s 2) (31) κ = e - 2 π 2 σ r 2 σ s 2 r 2 (σ r 2 + σ s 2) (32) T = 2 π σ r 2 σ s 2 3 r 2 (σ r 2 + σ s 2) exp - (x - y) 2 2 (σ s 2 + σ r 2) (33) and Θ3 and Θ4 are the third and fourth Jacobi theta functions [120]. The cortical receptive fields RFy (x) are obtained by summing the ON and OFF sublattices RF y (x) = RF α, r, y ON (x) - RF α ′, r ′, y OFF (x), (34) where α (α′) and r (r′) are the angle and the lattice spacing of the ON (OFF) lattice. The inset of a receptive field in Fig 3B shows a plot of Eq (34). These receptive fields resemble simple cell receptive fields in V1 with a size of about 1° for our choice of parameters [121,122]. An implementation/visualization of the equations for receptive fields can be found in the Supplemental Material. The response Ry of a neuron with receptive field RF (x) to a sine wave grating can be calculated using L (x) = exp (−ik x) as a stimulus in Eq (5). Evaluating the integral then corresponds to Fourier transforming the receptive field RF (x). Denoting the Fourier transform of the receptive field as R y (k) = 1 2 π ∫ d 2 x RF y (x) e - i k x, (35) we refer to the absolute value | R y (k) | as the amplitude spectrum of the receptive field. Given the above definition, the amplitude spectrum represents the response to a sine wave grating with wave vector k = (k cos (ϑ), k cos (ϑ) ), where ϑ is the grating orientation and k its spatial frequency. A tuning curve for spatial frequencies k and orientations ϑ is given by TC (ϑ, k) = | R y (k cos (ϑ), k sin (ϑ) ) |. (36) We calculated the Fourier transform of Eq (34) by transforming Eq (2) and subsequently summing over the two rectangular lattices L 1 and L 2 in Eq (26). Interchanging summation and integration is valid because all infinite sums are uniformly convergent. The result is R (k) α, r, y ON/OFF = U Θ 3 (c e ϕ, ν) Θ 3 (c e r, ζ) + Θ 4 (c e ϕ, ν) Θ 4 (c e r, ζ), (37) where U = 2 π σ r 2 σ s 2 3 r 2 exp - i k y - 1 2 k 2 σ s 2 + σ r 2 (38) c = (y - i σ s 2 k) (39) ν = exp - 2 π 2 σ s 2 3 r 2 (40) ζ = exp - 2 π 2 σ s 2 r 2. (41) The Fourier transform of cortical receptive fields is given by the sum of the ON and OFF sublattice Fourier transforms R y (k) = R (k) α, r, y ON - R (k) α ′, r ′, y OFF, (42) where we suppressed the dependencies on α, α′, r, r′ on the left hand side. Receptive fields depend on the two scales σr and σs. With increasing σs, more RGCs are pooled to form the cortical receptive field. If the cortical receptive field is dominated by more than two RGCs, it can exhibit multiple ON and OFF subregions. The spatial arrangement of these ON and OFF subregion mirrors the hexagonal lattices of the ON and OFF center RGCs. The parameter σs must be of a minimal size since for very small values of many cortical cells are connected with only a single, dominant RGC input and exhibit no orientation selectivity. Varying σr does not qualitatively change the shape of cortical receptive fields. For simple cell receptive fields with one ON and one OFF subregion, the amplitude spectrum | R (k) | will typically look as in Fig 10A. We follow [28,82] and define the preferred angle as ϑpref: = arg (μ) /2, where μ = ∫ d 2 k | R (k) | · e 2 i arg (k) k ∫ d 2 k | R (k) |. (43) While there is consensus about the definition of the preferred orientation, methods for extracting the preferred spatial frequency differ within the literature. Ringach proposed to use kpref = ∣μ∣ [82], referred to as Center-of-Mass Method. More commonly, the circular variance CV (k) [123–125] CV (k) = ∫ 0 2 π d ϑ TC (ϑ, k) e 2 i ϑ ∫ 0 2 π d ϑ TC (ϑ, k), (44) is first computed as a measure of orientation selectivity at a given spatial frequency k. Maximizing the circular variance across all spatial frequencies is then performed to obtain an estimate of preferred spatial frequency: k pref = argmax { k } CV (k). (45) We refer to this method as CV maximization. Finally, one can use the maximum of the amplitude spectrum k pref = argmax { k } | R (k) |, (46) as an estimate of the preferred spatial frequency (Maximum method). We argue that Maximum method and CV maximization in most cases yield similar results. They extract preferred spatial frequencies that one would obtain when searching for the “strongest response” by presenting a set of gratings of varying orientation and spatial frequency to a subject [125–127]. In contrast, estimates made with the center-of-mass method are usually substantially smaller then this intuitive measure (Fig 10C and 10D). As a consequence, the orientation selectivity index defined as OSI = ∫ 0 2 π d ϑ TC (ϑ, k pref) e 2 i ϑ ∫ 0 2 π d ϑ TC (ϑ, k pref) (47) for all three methods, will usually be substantially smaller, when estimated with the Center-of-Mass method compared to the other two methods. Among all three methods, the Maximum method has the advantage that its estimates are unaffected by monotonic nonlinearities applied to R commonly used to convert it to a firing rate of a neuron. For this reason, the Maximum method is our method of choice for extracting the preferred spatial frequency from amplitude spectra of receptive fields. Using the equations for the receptive fields of cortical neurons in the Moiré interference model, we extracted the spatial progression preferred orientation and spatial frequency from their squared amplitude spectrum: | ℛ y (k) | 2 ∝ exp (− k 2 (σ s 2 + σ r 2) ) · | Θ 3 α, r Θ 3 α, r + Θ 4 α, r Θ 4 α, r − Θ 3 α ′, r ′ Θ 3 α ′, r ′ − Θ 4 α ′, r ′ Θ 4 α ′, r ′ ︸ G y (k) | 2, (48) with abbreviation Θ i α, r Θ i α, r = Θ i (c e ϕ (α, r), τ) Θ i (c e r (α, r), ζ). | R (k) | 2 is composed of a rotationally symmetric Gaussian envelope and a non-rotationally symmetric part Gy (k) varying in space y. To calculate the preferred orientation ϑpref and spatial frequency kpref, we expanded the non-rotationally symmetric part ∣Gy ( (k) ∣2 to quadratic order in k: | R y (k) | 2 ≈ exp - k 2 (σ s 2 + σ r 2) | G y 0 | 2 + 1 2 k 1 k 2 H y k 1 k 2, (49) where H y is the Hessian matrix H y = ∂ 2 | G y | 2 ∂ k 1 2 ∂ 2 | G y | 2 ∂ k 1 ∂ k 2 ∂ 2 | G y | 2 ∂ k 2 ∂ k 1 ∂ 2 | G y | 2 ∂ k 2 2 k = 0 ≡ a y b y b y c y (50) and G y 0 = G y (k = 0). Since Gy (k) = Gy (−k), this Taylor expansion only contains terms of even power in k. Using the fact that h (θ) = cos θ sin θ H y cos θ sin θ = a y cos 2 θ + 2 b y cos θ sin θ + c y sin 2 θ (51) yields the second directional derivative in the direction of (cos θ, sin θ), ϑpref can be found as the maximum of h (θ), i. e. the direction of largest increase in amplitude spectrum, ϑ pref (y) = atan (a y - c y) 2 + 4 b y 2 - a y + c y 2 b y. (52) This formula for ϑpref (y) represents an expression for the orientation domain layouts produced by the Moiré interference model for hexagonal RGC lattices. To calculate the smoothed domain layout of the Moiré interference model analytically, we identified the low frequency components of our analytical solution. To this end, we expanded the Jacobi theta functions in Eq (42) [120] | R y (k) | 2 = ∑ j C y j exp - 1 2 σ 2 (k - a y j) 2 + C y j exp - 1 2 σ 2 (k + a y j) 2 (53) where C y j and a y j are determined by Eq (42). According to this equation, the power spectra of receptive fields in the Moiré interference model is represented by an infinite sum of Gaussians, each mirrored at the origin (0,0) of Fourier space. The preferred orientation of a receptive field represented by such an infinite sum is set by the direction in which the “center-of-mass” of the Gaussians is located. Due to the symmetry of the power under spatial inversion, there are two peaks located at ϑ (y) and π + ϑ (y). The direction towards the center-of-mass of the peak is obtained through the complex number μ (y) = ∑ j C y j · | a y j | exp (2 i arg (a y j) ) (54) with C y j and a y j defined as in Eq (53). The preferred orientation then is arg (μ (y) ) /2. Rewriting this sum and substituting the respective expressions for C y j and a y j, we obtained μ (y) = ∑ m, n, o, p f (m, n, o, p) exp (2 i y (n e ϕ + m e r - o e r ′ - p e ϕ ′) ), (55) with coefficients f (m, n, o, p). This is a decomposition of the orientation domain layout of the Moiré interference model into Fourier modes, indexed by four numbers m, n, o, p = 0, ±1, ±2, …. Table 2 lists the first terms of this series in ascending spatial frequency order. By rewriting μ → z and selecting only the lowest contributing spatial frequencies, we obtain Eq (7). The phase factor u 0 = e i (α + α ′) e i α ′ r + e i α r ′ e i α r + e i α ′ r ′ (56) is associated with an overall shift of all preferred orientations. Note that ∣u0∣2 = 1. Realization of a random distortion field y (x) were generated by finding a complex-valued field z (x) of which real and imaginary part correspond to dislocations in x and y direction, respectively. We constructed such a field using established methods (e. g. [51]) in the Fourier domain. In short, we drew complex-valued amplitudes a (k) from a Gaussian distribution satisfying 〈 a (k) a (k ′) 〉 = f ˜ (k) · δ k, k ′, where f ˜ (k) was a chosen power spectrum, in our case a Gaussian with width 1/σ, σ being the desired correlation length. The corresponding amplitude spectrum was then inversely Fourier transformed to obtain a complex-valued field z (x). Real and imaginary part of this field constitute two independent real-valued Gaussian random fields, both with the desired spatial statistics. We then transformed the coordinates of the hexagonal ON/OFF lattice points ri = (xi, yi) according to x i → x i + η ℜ (z (r i) ) and y i → y i + η ℑ (z (r i) ). For the displacements of ON and OFF lattices, we used two independent complex-valued Gaussian random field realizations. Source code to generate hexagonal RGC mosaics with correlated spatial noise along with Matlab code for visualization is part of the supplementary material to this manuscript. We generated RGC mosaics with a pairwise interacting point process using the code published by Schottdorf et al. [32] derived from the method developed in [49,82].
In the primary visual cortex of primates and carnivores, local visual stimulus features such as edge orientation are processed by neurons arranged in arrays of iso-orientation domains. Large-scale comparative studies have uncovered that the spatial layout of these domains and their topological defects follows species-invariant quantitative laws, predicted by models of large-scale circuit self-organization. Here, we ask whether the experimentally observed layout invariants might alternatively emerge as a consequence of random connectivity rules for feedforward projections from a small number of retinal cells to a much larger number of cortical target neurons. In this random wiring framework, the semi-regular and spatially granular arrangement of retinal ganglion cells determines the spatial layout of visual cortical iso-orientation domains—a hypothesis diametrically opposed to cortical large-scale circuit self-organization. Generalizing a prominent model of the early visual pathway, we find that the random wiring framework does not reproduce the experimentally determined layout invariants. Our results demonstrate how comparison between theory and quantitative phenomenological laws obtained from large-scale experimental data can successfully discriminate between competing hypotheses about the design principles of cortical circuits.
Random Wiring, Ganglion Cell Mosaics, and the Functional Architecture of the Visual Cortex