Random Wiring, Ganglion Cell Mosaics, and the Functional Architecture of the Visual Cortex

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.


Introduction
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][2][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 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][36][37][38][39][40][41][42] (Fig 1A, see also [23,33]). Their last common ancestor was a small shrew-like mammal [40][41][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 selforganization. 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 r ¼ 2 ffiffi ffi 3 p % 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 Common laws for the layout of iso-orientation domains in different mammalian species. A Phylogenetic relationships and macroevolution of laurasiatheria, euarchonta and glires [33-36, 52, 53]. B Key features of the thalamo-cortical pathway for cat [27,45,46,[54][55][56][57][58][59][60][61][62][63][64][65], macaque [66][67][68][69][70][71][72][73], treeshrew [19,43,44,[74][75][76] and mouse [77][78][79][80][81] at the level of retina/LGN and layer IV and II/III of V1. All species show orientation selective neurons in layer II/III, but only cat, ferret, and mouse exhibit orientation selectivity in input layer IV. Ocular dominance domain layouts differs greatly between all four species, macaque is the only species listed possessing trichromatic color vision. Only cat and macaque V1 display cytochrome oxidase blobs. Non-classical receptive fields are mediated by different circuits in cat, tree shrew and macaque. C Pinwheel density ρ in ferret (N = 82), dark-reared ferret (N = 21), cat (N = 13), tree shrew (N = 26), and galago (N = 9). Light green shading indicates one-species consistency range, dark green shading indicates common design consistency range (see text). D Illustration of the common design layout features, nearest neighbor (NN) distances, and pinwheel density in subregions of varying size. E, F Standard deviations (SD) of pinwheel densities as a function of the area A of randomly selected subregions. SD(A) is well described by a power law with variability exponent γ (F, top) and variability coefficient c (F, bottom). (G-J) Nearest neighbor distance distributions for pinwheels of arbitrary (G), opposite (H) and equal (J) topological charge in units of the column spacing. Insets indicate species means. All error bars represent 95% confidence intervals of the bootstrap distributions.
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.

A benchmark for models of orientation domains in V1
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 ( 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
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][55][56][57][58][59][60][61][62][63][64][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 Table 1. The six orientation domain layout parameters characterizing the common design. Values were calculated with the code provided in the supplemental material and intervals indicate 95% bootstrap confidence intervals. Also shown is the grand average and the associated one species and common design consistency ranges (CR). ( 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 twodimensional 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 GRF j (x) of width σ r localized at the center position x j :

Pinwheel density
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 RF y 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 w j : The synaptic weights are chosen as The parameter σ s sets the range from which a V1 neuron receives retino-thalamic inputs, x j 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 RF y (x) and the stimulus, in our case an illumination pattern L(x): Because R y 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 0 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  (52) in Methods) is shown on the right with the mosaic overlaid. Bottom: low frequency contribution of the domain layout. Black arrows indicate the lattice constant S Á r of the Moiré pattern, white hexagon indicates the unit cell of the domain layout. B Inset of the layout shown in A with RFs of three closely spaced neurons. Scale bar indicates distance on the retina. C Circular distance (see text) between the preferred angles of unfiltered and low-pass filtered domain layouts shown in A top and bottom. Bottom: Histogram of differences in preferred angles. Model parameters: σ r = 70 μm, σ s = 20 μm, lattice constants r = r 0 = 170 μm, and a relative angle Δα = 7°leading to a scaling factor of S = 8.2 (Eq (10)), as proposed in [28,29]. 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].

Orientation preference maps from crystalline RGC mosaics
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 Here, r, r 0 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 k j and fixed phase factors u j , The pattern of preferred orientations across the cortical coordinate y is given by the phase of this complex-valued field as, Fig 3A (bottom) depicts ϑ pref (y) as analytically determined. The pattern of pinwheels and isoorientation 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 2iDð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 0 ! (1 + β)r. The expression for the column spacing becomes where S is the distance between two vertices of the Moiré pattern in units of r, called the scaling factor [94-96] 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 L 2 c . Within each unit cell of area A ¼ ffi ffi one "double pinwheel" of topological charge 1, around which each orientation is represented twice, and two pinwheels of topological charge AE 1 2 . With L 2 c ¼ 3 4 ðS Á rÞ 2 and counting the pinwheel with charge 1 as two pinwheels (see below), the pinwheel density is 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 ffiffi ffi 3 p 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 .

The impact of spatially uncorrelated disorder in RGC position
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), Normalizing and radially averaging yields the so-called marginal amplitude spectrum ( Fig 5C  and 5D), The sharp peak at k c 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 ffiffi ffi 3 p 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  (6)). B As A but for a higher disorder strength η = 0.30. C Radially averaged normalized amplitude spectra of the orientation domain layouts for different disorder strengths. The fluctuation strength is color coded (legend). x-axis is given in units of k c (cf. Eq (6)). D As C but for the smoothed layouts. All other model parameters as in   . Dark green line in inset shows experimentally observed mean value ρ = 3.14. C Pinwheel density in circular regions of increasing area for η = 0.12 (blue) and η = 0.02 (red). Lines show theoretically predicted and experimentally observed values as in B inset. D The standard deviation of pinwheel density estimates for increasing subregion size. Red line shows a power law fit to the experimental data (γ = 0.5, [23]), purple line indicates a fit to the perfectly ordered hexagonal pinwheel arrangement (γ = 0.75, [23,92,93]). E Nearest neighbor (NN) distances for pinwheels irrespective of topological charge. Left: distributions for two disorder strengths (η = 0.11, blue; η = 0.02, purple) and the experimental data (red). Right: Distributions for different disorder strengths. Color encodes the (normalized) fraction of pinwheels at this distance. Blue and purple lines indicate disorder strengths shown on the left. The white line marks the theoretically predicted distance of NN pinwheels (2/3Λ c ) for vanishing disorder [92,93]. F same as E for pinwheels of equal charge, data green curve. G same as E for pinwheels of opposite charge, data blue curve. H squared deviation of NN distance distributions to the experimental estimates (shown in E-G, left) as function of disorder strength. All other model parameters as in 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.  Table 1). With increasing disorder, pinwheel density of model layouts steadily increases from  1). B The variability exponent as function of disorder strength in comparison to the common design. C As B for the variability constant. D Mean nearest neighbor distance for pinwheels independent of topological charge in comparison to common design. E As D but for pinwheels of equal charge. F As D for pinwheels of opposite charge. G Summary of ranges of disorder parameters consistent with the common design. Disorder strengths larger than 0.3 can be excluded by the lack of typical column spacing in the domain layouts (cf. Fig 5). Note that there is no disorder strengths for which all features of the common design are reproduced.  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.

Iso-orientation domain layouts from hexagonal RGC mosaics with spatially correlated disorder
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 x i was then shifted depending on its position according to x i ! x i + η y(x i ). The shift η y(x) with amplitude η for y(x) = (y 1 (x), y 2 (x)) was chosen from a Gaussian random field with vanishing mean hy 1 (x)i = hy 2 (x)i = 0, fixed standard deviation std(y 1 (x)) = std(y 2 (x)) = 1 and correlation function hyðx 1 Þyðx 2 Þi ¼ 2 exp À jx 1 Àx 2 j 2 2s 2 with correlation length σ (see Methods) where y 1 and y 2 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 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 isoorientation 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 ffiffi ffi 3 p 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 mm 2 is almost independent of either correlation length or disorder strength.

Pinwheel densities of iso-orientation domain layouts derived from PIPP mosaics
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: with a, b, β > 0. The filter function was normalized such that With this normalization, one can define the mean column spacing of the resulting layout via 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/commondesign 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.
(right). Scale bar indicates cortical distances, assuming cortical magnification factor % 1, and Λ = 0.9mm (see Eq (17) and text). Red circles indicate k c = 2π/ Λ. Black square indicates inset in A, white square indicates Λ 2 . D As C but filtered with Fermi band pass filters [23]. White square (top) indicates Λ = 0.68mm, the column spacing as measured by wavelet analysis. Red circles (bottom) indicate low pass (1.2 mm) and high pass (0.3 mm) position. Pinwheel densities are stated with standard error of the mean. E Left: Analytically predicted pinwheel density of orientation domain layouts derived from Gaussian random fields [51] as a function of filter parameter β and spatial scale (see text). Right: Pinwheel density of orientation domain layouts obtained from PIPP mosaics with the statistical connectivity model as a function of filter parameter β and spatial scale. Numbers 1-3 indicate parameter choices displayed in C. doi:10.1371/journal.pcbi.1004602.g009

Discussion
In this study, we examined whether the statistical connectivity model-a biologically plausible scheme of circuit disorder-is able to explain the common design of spatially aperiodic arrangements of orientation domains and pinwheels in the primary visual cortex. As an analytically tractable limiting case, we first considered the model with perfectly ordered hexagonal RGC mosaics (Moiré interference model). For this model we derived exact expressions for receptive fields and tuning curves as well as for unfiltered and filtered layouts. We found that unfiltered orientation domain layouts generated by Moiré-interference exhibited a fine-grained structure of subdomains with substantial and systematic variation in orientation preference on scales much smaller than the typical size of orientation domains. After smoothing, the resulting Moiré interference pattern could mathematically be expressed as the phase of a complex-valued field composed of six planar waves. The pinwheel density of this perfectly hexagonal pattern of orientation domains is r ¼ 2 ffiffi ffi 3 p % 3:46. Next, we studied the layout of numerically obtained domain layouts derived from hexagonal mosaics that are randomly distorted by spatially uncorrelated disorder. We found that pinwheel density and pinwheel nearest neighbor statistics vary substantially with the degree of randomness. Nevertheless, there was no parameter regime in which all of the common design parameters matched experimental observations. Most prominently, the pinwheel density increased monotonically with increasing disorder strength. To examine the effect of noisy RGC mosaics more broadly, we introduced a more general class of noisy hexagonal mosaics, which allows for the inclusion of spatial correlations in RGC positional disorder. We found that, while RGC dipole patterns for such mosaics are inherently aperiodic, the model still predicts domain layouts that substantially deviate from experimentally observed pinwheel layouts. Finally, we studied the model with RGC mosaics derived from Eglen's random pairwise interacting point process. The resulting layouts lacked a typical spacing between neighboring orientation domains and, after bandpass filtering, pinwheel densities were inconsistent with the values observed for any of the four species investigated.

Alternative random wiring models
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.

Spatial irregularity by disorder or optimization
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].

Hexagonal order of orientation domains
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 selforganization 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.

The impact of retinal orientation biases on visual cortical architecture
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 as benchmark for models of visual cortical development and function
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][114][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][117][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 longrange interaction model currently is the only theory known to be capable of explaining the common design of orientation domains in the primary visual cortex.

Analysis of model and experimentally obtained orientation domain layouts
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 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 and y ðx; L; φÞ ¼ ðO À1 ðφÞðy À xÞÞ: The matrix O(φ) 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 ξ 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 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.

A mathematical treatment of the Moiré-interference model
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 0 , where f = r, r 0 is the lattice constant. Describing a rotation of the lattice vectors by the rotation matrix the ON mosaic is rotated by an angle α, the OFF lattice by an angle α 0 (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 RF y 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 Evaluating the infinite Gaussian sum (Eq (3)  where α (α 0 ) and r (r 0 ) 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.

Tuning curves from receptive fields
The response R y 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 we refer to the absolute value jR y ðkÞj 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ðW; kÞ ¼ jR y ðk cos ðWÞ; k sin ðWÞÞj: ð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 where we suppressed the dependencies on α, α 0 , r, r 0 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.
Extracting preferred orientation and spatial frequency from amplitude spectra of receptive fields For simple cell receptive fields with one ON and one OFF subregion, the amplitude spectrum jRðkÞj will typically look as in Fig 10A. We follow [28,82] and define the preferred angle as ϑ pref : = arg(μ)/2, where 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 k pref = jμj [82], referred to as Center-of-Mass Method. More commonly, the circular variance CV(k) [123][124][125] CVðkÞ 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: We refer to this method as CV maximization. Finally, one can use the maximum of the amplitude spectrum 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][126][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 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.

Extracting the spatial progression of preferred orientation
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: jR y ðkÞj 2 / expðÀk 2 ðs 2 s þ s 2 r ÞÞ Á Y a;r 3 Y a;r 3 þ Y a;r 4 Y a;r 4 À Y a 0 ;r 0 3 Y a 0 ;r 0 3 À Y a 0 ;r 0 4 Y a 0 ;r 0 4 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl} with abbreviation Y a;r i Y a;r i ¼ Y i ðce ða; rÞ; tÞY i ðce r ða; rÞ; zÞ. jRðkÞj 2 is composed of a rotationally symmetric Gaussian envelope and a non-rotationally symmetric part G y (k) varying in space y. To calculate the preferred orientation ϑ pref and spatial frequency k pref , we expanded the non-rotationally symmetric part jG y ((k)j 2 to quadratic order in k: jR y ðkÞj 2 % exp Àk 2 ðs 2 s þ s 2 r Þ À Á jG 0 y j 2 þ 1 2 where H y is the Hessian matrix H y ¼ @ 2 jG y j 2 @k 2 1 @ 2 jG y j 2 @k 1 @k 2 @ 2 jG y j 2 @k 2 @k 1 and G 0 y ¼ G y ðk ¼ 0Þ. Since G y (k) = G y (−k), this Taylor expansion only contains terms of even power in k.
Using the fact that hðyÞ ¼ cos y sin y ð Þ H y cos y sin y ! ¼ a y cos 2 y þ 2b y cos y sin y þ c y sin 2 y ð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, 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] jR y ðkÞj 2 ¼ X j C j y exp À 1 2s 2 ðk À a j y Þ 2 þ C j y exp À 1 2s 2 ðk þ a j y Þ 2 where C j y and a j y 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 mðyÞ ¼ X j C j y Á ja j y j exp ð2i argða j y ÞÞ ð54Þ with C j y and a j y defined as in Eq (53). The preferred orientation then is arg(μ(y))/2. Rewriting this sum and substituting the respective expressions for C j y and a j y , we obtained mðyÞ ¼ X m;n;o;p f ðm; n; o; pÞ exp ð2iyðne þ me r À oe 0 r À pe 0 ÞÞ; ð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 ðaþa 0 Þ e ia 0 r þ e ia r 0 À Á e ia r þ e ia 0 r 0 ð56Þ is associated with an overall shift of all preferred orientations. Note that ju 0 j 2 = 1.

Correlated noise on hexagonal RGC mosaics
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 haðkÞaðk 0 Þi ¼f ðkÞ Á d k;k 0 , wheref ð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 r i = (x i , y i ) according to x i ! x i þ Z <ðzðr i ÞÞ and y i ! y i þ Z 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.

Generating PIPP RGCs mosaics
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].
Supporting Information S1 Code. Source code for pinwheel statistics analysis and simulating statistical connectivity model layouts with a variety of RGC mosaics. Numerical implementation of the statistical Table 2. Lowest frequency contributions of the orientation domain layout predicted by the Moiré interference model. The vectors k i are the wave vectors and jk i j their absolute values. u i denotes the phase factors of the complex-valued coefficient f(m, n, o, p) in Eq (55) (see also Eq (7)). The constant phase factor u 0 is given by Eq (56 wiring model. We provide all necessary code to calculate single neuron properties and orientation domain layouts along with a Mathematica program which contains the analytical solution for both, an example single neuron and the domain layout obtained from a perfect and infinite lattice. The C-program 'calculate_single_neuron.cpp' calculates the same for a single neuron numerically. The C-program is provided to illustrate the use of the rfanalyzer class. The c-program 'calculate_map.cpp' calculates the same properties as 'calculate_single_neuron.cpp' but for a whole array of cells. After finishing a run, this program generates a set of ascii files in which the output is stored. These files are read in and analyzed by the Matlab program 'plot_results.m'. It calculates the pinwheel density, pinwheel distance distributions, mean pinwheel distance and pinwheel density fluctuations as a function of subregion size. We compiled the code with gcc [g++ (Ubuntu 4.8.2-19ubuntu1) 4.8.2] and the gsl: g++ ./calculate_single_neuron.cpp -lgsl -lgslcblas -O3 -march = native g++ ./calculate_map.cpp -lgsl -lgslcblas -O3 -march = native For this article, we have calculated orientation domain layouts with aspect ratio 22x22Λ, sampled with 4096x4096 pixels. This corresponds to % 6.5 μm per cortical unit for our standard combination of parameters (r = r 0 = 170 μm and Δα = 7°). Experimental data The folder 'map_data' contains a data folder with single condition layouts and various ROIs for four ferrets cases. It also contains two Matlab files, 'run_analysis.m' and 'plot_results.m' to run the analysis and display the results. The full data set used in the present study is available on the neural data sharing platform http://www.g-node.org/. (ZIP)