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 (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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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).

https://doi.org/10.1371/journal.pcbi.1004602.t001

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–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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Early visual pathway, RGC dipoles, and Moiré interference of RGC mosaics.

A Schematic illustration of the early visual pathway following the organization in the cat (see text for details). B Orientation selective receptive fields can arise through summation of two adjacent rotationally symmetric retinal/LGN receptive fields (RGC dipole). Shown are ON and OFF center mosaics from cat retina [26]. Colors indicate ON (red) and OFF (blue) regions of a receptive field in the LGN (left and middle) and V1 (right). For illustration, the RGC mosaic is overlaid and the two RGCs whose RFs are summed are shown as black and white dots. C Left: Moiré interference between a hexagonal ON (white dots) and OFF (black dots) RGC lattice with relative orientation Δα and lattice constants r and r′ (black bars) creates a Moiré pattern with lattice constant S ⋅ r. Middle: sampling from this RGC mosaic as described in B (and text) yields a periodic orientation preference pattern through Moiré interference. Right: Model layout predicted by the statistical wiring model, obtained by thresholding and smoothing the pattern in the middle (see text).

https://doi.org/10.1371/journal.pcbi.1004602.g002

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 GRF_j(x) of width σ_r localized at the center position x_j: (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 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: (3) The synaptic weights are chosen as (4) 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): (5) 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′ 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].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Receptive fields and iso-orientation domains in the Moiré interference model.

A Top: Moiré interference between two RGC mosaic (left) with ON and OFF center RGCs illustrated as white and black dots. The corresponding orientation domain layout (Eq (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′ = 170μm, and a relative angle Δα = 7° leading to a scaling factor of S = 8.2 (Eq (10)), as proposed in [28, 29].

https://doi.org/10.1371/journal.pcbi.1004602.g003

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 (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 k_j and fixed phase factors u_j, (7) The pattern of preferred orientations across the cortical coordinate y is given by the phase of this complex-valued field as, (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) = ϑ₁(x) − ϑ_pref(x), the difference d(x) between the two preferred angles is defined as . 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 (9) where S is the distance between two vertices of the Moiré pattern in units of r, called the scaling factor [94–96] (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 . Within each unit cell of area , there is one “double pinwheel” of topological charge 1, around which each orientation is represented twice, and two pinwheels of topological charge . With and counting the pinwheel with charge 1 as two pinwheels (see below), the pinwheel density is (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 as analytically predicted (see Methods).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Comparison of analytically and numerically obtained solutions of the Moiré interference model.

A Top: unfiltered Moiré interference patterns. Black line separates analytical (left, see Eq (52)) from numerical result (right). Bottom: amplitude spectrum of numerically obtained Moiré interference patterns. Red circle marks k_c (cf. Eq (6)), blue circle indicates . Note the high frequency contributions, indicated by the small yellow circles. B Top: Numerically obtained Moiré interference pattern after thresholding for cells with OSI > 0.25 (left, see Methods for the OSI definition used) and subsequent smoothing (right, see text). Inset shows a magnified region of the OSI-filtered domain layout together with the RGC mosaic from which the neurons sample. Bottom: amplitude spectrum of the numerically obtained thresholded and smoothed layout. Red circle indicates k_c (cf. Eq (6)). C Orientation domain layout (top) and amplitude spectrum (bottom) obtained by calculating the lowest spatial frequency contributions of the layout in A (Eq (8)). All model parameters as in Fig 3.

https://doi.org/10.1371/journal.pcbi.1004602.g004

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 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 .

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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Spatially uncorrelated position disorder in hexagonal RGC mosaics induce broadband noise in iso-orientation domain layouts.

A Numerically calculated orientation domain layouts with disorder strength η = 0.12. From left to right: Moiré interference pattern, filtered Moiré interference pattern (OSI > 0.25), smooth layout and the smoothed layout’s amplitude spectrum. Insets show magnified regions. Circles in the amplitude spectrum mark k_c (red) and (blue) (cf. Eq (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 Fig 3.

https://doi.org/10.1371/journal.pcbi.1004602.g005

To characterize the model orientation domain arrangements, we first calculate amplitude spectra for both, unfiltered and smoothed layouts (Fig 5A and 5B), (12) Normalizing and radially averaging yields the so-called marginal amplitude spectrum (Fig 5C and 5D), (13) 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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Pinwheel statistics in the Moiré interference model.

A Column spacing Λ estimated by the wavelet method compared to the Moiré scale Λ_c for different disorder strength. B The pinwheel density ρ as function of disorder strength. Dashed line shows theoretically predicted value . 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 Fig 3.

https://doi.org/10.1371/journal.pcbi.1004602.g006

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 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 (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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Pinwheel statistics of the disordered Moiré interference model fail to match V1 functional architecture.

A Pinwheel density as a function of disorder strength in the statistical connectivity model with noisy hexagonal mosaics. Error bars indicate standard deviation around the mean for 20 model realizations, circles indicate the mean. Green shaded areas indicate the range consistent with the experiments (see Fig 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.

https://doi.org/10.1371/journal.pcbi.1004602.g007

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₁(x), y₂(x)) was chosen from a Gaussian random field with vanishing mean 〈y₁(x)〉 = 〈y₂(x)〉 = 0, fixed standard deviation std(y₁(x)) = std(y₂(x)) = 1 and correlation function with correlation length σ (see Methods) where y₁ and y₂ 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 (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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Impact of spatially correlated positional disorder in hexagonal RGC mosaics on iso-orientation domain layouts.

A Generation of spatially correlated disorder in hexagonal RGC mosaics (see text and Methods). White dots mark ON RGC cells, gray dots marks the perfectly hexagonal ON RGC lattice before distortion. Colors of the arrows and underlay indicate direction of RGC displacements (right colorbar). Magnitude of displacement is indicated by color saturation (bottom colorbar). Note that nearby RGCs move into similar directions and with similar magnitude. Scale bar indicates retinal distances, assuming ON/OFF lattice constant of 170μm. B Larger region of the mosaic shown in A (black square). Colors as in A. Contour lines mark lines of constant det(J(x)) (see text). Scale bar indicates retinal distances. C-E Unfiltered model domain layout (left), thresholded and smoothed layout (middle) its amplitude spectrum (right) for noise correlation length ξ = 5 and different noise amplitudes (η = 0.1 (C), η = 0.2 (D), η = 0.4 (E)). White squares in middle panels indicate wavelength of pattern as measured by wavelet analysis (see Methods). Scale bars indicate cortical distance with parameter choices as in Fig 2 and cortical magnification factor ≈ 1 (see [32]). F Marginal amplitude spectra of smoothed domain layouts for ξ = 5 and different disorder strengths. G Pinwheel density of model layouts as a function of disorder amplitude and correlation length. Red surface indicates experimentally determined value (ρ = 3.14), green surface indicates pinwheel density in the disorder-free case (). Values for ρ ≥ 5 are drawn as a plane at ρ = 5. H As G, but pinwheels per square millimeter . Note that the absolute number of pinwheels is largely constant. All other model parameters as in Fig 2.

https://doi.org/10.1371/journal.pcbi.1004602.g008

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 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² 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: (15) with a, b, β > 0. The filter function was normalized such that (16) With this normalization, one can define the mean column spacing of the resulting layout via with (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⁻¹ 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 Λ². 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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Iso-orientation domain layouts obtained from PIPP RGCs with the statistical connectivity model.

A Generating orientation domain layouts from PIPP RGCs in the statistical connectivity model [83]. Top: Inset of a PIPP RGC mosaic (see Methods). Black (white) dots represent OFF (ON) cells. Middle top: unfiltered layout with RGC mosaics overlaid. Middle bottom: thresholded layout with RGC mosaics overlaid. Bottom: thresholded and smoothed layout (β = 0) with RGC mosaics overlaid. Scale bar indicates retinal distances, assuming PIPP parameters as in [49]. B Left: larger region of the unfiltered layout shown in A (black square). Scale bar indicates retinal distances. Right: normalized amplitude spectrum of unfiltered layout shown on the left. C Thresholded and smoothed layout (top) and corresponding amplitude spectrum (bottom) for filter function parameters (see Eq (15)) β = 0 (left), β = 2 (middle), and β = 10 (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 Λ². 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.

https://doi.org/10.1371/journal.pcbi.1004602.g009

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.

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 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.

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–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.

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 (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 (19) and (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 (21) (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 (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.

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 , (Fig 2B), that may exhibit different lattice constants r and r′, (24) where f = r, r′ is the lattice constant. Describing a rotation of the lattice vectors by the rotation matrix (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 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 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 (26) Evaluating the infinite Gaussian sum (Eq (3)) yields the result for a single sub lattice (either ON or OFF) (27) where (28) (29) (30) (31) (32) (33) and Θ₃ and Θ₄ are the third and fourth Jacobi theta functions [120]. The cortical receptive fields RF_y(x) are obtained by summing the ON and OFF sublattices (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.

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 (35) we refer to the absolute value 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 (36) We calculated the Fourier transform of Eq (34) by transforming Eq (2) and subsequently summing over the two rectangular lattices and in Eq (26). Interchanging summation and integration is valid because all infinite sums are uniformly convergent. The result is (37) where (38) (39) (40) (41) The Fourier transform of cortical receptive fields is given by the sum of the ON and OFF sublattice Fourier transforms (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.

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 will typically look as in Fig 10A. We follow [28, 82] and define the preferred angle as ϑ_pref: = arg(μ)/2, where (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 k_pref = ∣μ∣ [82], referred to as Center-of-Mass Method. More commonly, the circular variance CV(k) [123–125] (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: (45) We refer to this method as CV maximization. Finally, one can use the maximum of the amplitude spectrum (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 (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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. A comparison of different methods for extracting RF parameters from their amplitude spectrum.

A Amplitude spectrum ∣R(k)∣ of a simple cell receptive field calculated with Eq (42) in the manuscript. Circles indicate the preferred spatial frequency as extracted by maximizing the CV (yellow), by the Center-of-Mass Method (brown), and the Maximum Method (pink). B The Tuning curves corresponding to the green, cyan and blue circles in A, normalized relative to their maximum. C The Tuning curves corresponding to the brown, the yellow and the pink circles in A, normalized relative to their minimum. D Circular variance of the various tuning curves calculated via Eq (44) as function of spatial frequency. The brown, the yellow and the pink line correspond to the preferred spatial frequency as extracted by the three methods.

https://doi.org/10.1371/journal.pcbi.1004602.g010

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: (48) with abbreviation . 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 ∣G_y((k)∣² to quadratic order in k: (49) where is the Hessian matrix (50) and . Since G_y(k) = G_y(−k), this Taylor expansion only contains terms of even power in k.

Using the fact that (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, (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] (53) where and 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 (54) with and defined as in Eq (53). The preferred orientation then is arg(μ(y))/2. Rewriting this sum and substituting the respective expressions for and , we obtained (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 (56) is associated with an overall shift of all preferred orientations. Note that ∣u₀∣² = 1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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 ∣k_i∣ 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₀ is given by Eq (56). n.d. means not determined.

https://doi.org/10.1371/journal.pcbi.1004602.t002

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 , where 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 and . 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].