Coverage, continuity, and visual cortical architecture

Background

The primary visual cortex of many mammals contains a continuous representation of visual space, with a roughly repetitive aperiodic map of orientation preferences superimposed. It was recently found that orientation preference maps (OPMs) obey statistical laws which are apparently invariant among species widely separated in eutherian evolution. Here, we examine whether one of the most prominent models for the optimization of cortical maps, the elastic net (EN) model, can reproduce this common design. The EN model generates representations which optimally trade of stimulus space coverage and map continuity. While this model has been used in numerous studies, no analytical results about the precise layout of the predicted OPMs have been obtained so far.

Results

We present a mathematical approach to analytically calculate the cortical representations predicted by the EN model for the joint mapping of stimulus position and orientation. We find that in all the previously studied regimes, predicted OPM layouts are perfectly periodic. An unbiased search through the EN parameter space identifies a novel regime of aperiodic OPMs with pinwheel densities lower than found in experiments. In an extreme limit, aperiodic OPMs quantitatively resembling experimental observations emerge. Stabilization of these layouts results from strong nonlocal interactions rather than from a coverage-continuity-compromise.

Conclusions

Our results demonstrate that optimization models for stimulus representations dominated by nonlocal suppressive interactions are in principle capable of correctly predicting the common OPM design. They question that visual cortical feature representations can be explained by a coverage-continuity-compromise.

The pattern of orientation columns in the primary visual cortex (V1) of carnivores, primates, and their close relatives are among the most intensely studied structures in the cerebral cortex and a large body of experimental (e.g., [1–13]) and theoretical work (e.g., [14–39]) has been dedicated to uncovering its organization principles and the circuit level mechanisms that underlie its development and operation. Orientation preference maps (OPMs) exhibit a roughly repetitive arrangement of preferred orientations in which adjacent columns preferring the same orientation are separated by a typical distance in the millimeter range [2–5, 10]. This range seems to be set by cortical mechanisms both intrinsic to a particular area [40] but potentially also involving interactions between different cortical regions [41]. The pattern of orientation columns is however not strictly periodic because the precise local arrangement of preferred orientation never exactly repeats. Instead, OPMs appear as organized by a spatially complex aperiodic array of pinwheel centers, around which columns activated by different stimulus orientations are radially arranged like the spokes of a wheel [2–5, 10]. The arrangement of these pinwheel centers, although spatially irregular, is statistically distinct from a pattern of randomly positioned points [38] as well as from patterns of phase singularities in a random pattern of preferred orientations [32, 36, 38, 42] with spatial correlations identical to experimental observations [38, 42]. This suggests that the layout of orientation columns and pinwheels although spatially aperiodic follows a definite system of layout rules. Cortical columns can in principle exhibit almost perfectly repetitive order as exemplified by ocular dominance (OD) bands in the macaque monkey primary visual cortex [43, 44]. It is thus a fundamental question for understanding visual cortical architecture, whether there are layout principles that prohibit a spatially exactly periodic organization of orientation columns and instead enforce complex arrangements of these columns.

Recent comparative data have raised the urgency of answering this question and of dissecting what is constitutive of such complex layout principles. Kaschube et al. [38] quantitatively compared pinwheel arrangements in a large dataset from three species widely separated in the evolution of eutherian mammals. These authors found that the spatial statistics of pinwheels are surprisingly invariant. In particular, the overall pinwheel density and the variability of pinwheel densities in regions from the scale of a single hypercolumn to substantial fractions of the entire primary visual cortex were found to be virtually identical. Characterizing pinwheel layout on the scale of individual hypercolumns, they found the distributions of nearest-neighbor pinwheel distances to be almost indistinguishable. Further supporting common layout rules for orientation columns in carnivores and primates, the spatial configuration of the superficial patch system [45] and the responses to drifting grating stimuli were recently found to be very similar in cat and macaque monkey primary visual cortex [46].

From an evolutionary perspective, the occurrence of quantitatively similar layouts for OPMs in primate tree shrews and carnivorous species appears highly informative. The evolutionary lineages of these taxa diverged more than 65 million years ago during the basal radiation of eutherian mammals [47–49]. According to the fossil record and cladistic reconstructions, their last common ancestors (called the boreo-eutherial ancestors) were small-brained, nocturnal, squirrel-like animals of reduced visual abilities with a telencephalon containing only a minor neocortical fraction [47, 50]. For instance, endocast analysis of a representative stem eutherian from the late cretaceous indicates a total anterior-posterior extent of 4 mm for its entire neocortex [47, 50]. Similarly, the tenrec (Echinops telfari), one of the closest living relatives of the boreoeutherian ancestor [51, 52], has a neocortex of essentially the same size and a visual cortex that totals only 2 mm² [47]. Since the neocortex of early mammals was subdivided into several cortical areas [47] and orientation hypercolumns measure between 0.4 and 1.4 mm² [38], it is difficult to envision ancestral eutherians with a system of orientation columns. In fact, no extant mammal with a visual cortex of such size is known to possess orientation columns [53]. It is therefore conceivable that systems of orientation columns independently evolved in laurasiatheria (such as carnivores) and in euarchonta (such as tree shrews and primates). Because galagos, tree-shrews, and ferrets strongly differ in habitat and ecologically relevant visual behaviors, it is not obvious that the quantitative similarity of pinwheel layout rules in their lineages evolved driven by specific functional selection pressures (see [54] for an extended discussion). Kaschube et al. instead demonstrated that an independent emergence of identical layout rules for pinwheels and orientation columns can be explained by mathematically universal properties of a wide class of models for neural circuit self-organization.

According to the self-organization scenario, the common design would result from developmental constraints robustly imposed by adopting a particular kind of self-organization mechanism for constructing visual cortical circuitry. Even if this scenario is correct, one question still remains: What drove the different lineages to adopt a similar self-organization mechanism? As pointed out above, it is not easy to conceive that this adoption was favored by the specific demands of their particular visual habitats. It is, however, conceivable that general requirements for a versatile and representationally powerful cortical circuit architecture are realized by the common design. If this was true, the evolutionary benefit of meeting these requirements might have driven the adoption of large-scale self-organization and the emergence of the common design over evolutionary times.

The most prominent candidate for such a general requirement is the hypothesis of a coverage-continuity-compromise (e.g., [19, 21, 55, 56]). It states that the columnar organization is shaped to achieve an optimal tradeoff between the coverage of the space of visual stimulus features and the continuity of their cortical representation. On the one hand, each combination of stimulus features should be well represented in a cortical map to avoid 'blindness' to stimuli with particular feature combinations. On the other hand, the wiring cost to establish connections within the map of orientation preference should be kept low. This can be achieved if neurons that are physically close in the cortex tend to have similar stimulus preferences. These two design goals generally compete with each other. The better a cortical representation covers the stimulus space, the more discontinuous it has to be. The tradeoff between the two aspects can be modeled in what is called a dimension reduction framework in which cortical maps are viewed as two-dimensional sheets which fold and twist in a higher-dimensional stimulus space (see Figure 1) to cover it as uniformly as possible while minimizing some measure of continuity [21, 57, 58].

The dimension reduction framework. In the dimension reduction framework, the visual cortex is modeled as a two-dimensional sheet that twists in a higher-dimensional stimulus (or feature) space to cover it as uniformly as possible while minimizing some measure of continuity (left). In this way, it represents a mapping from the cortical surface to the manifold of visual stimulus features such as orientation and retinotopy (right).

Full size image

From prior work, the coverage-continuity-compromise appears to be a promising candidate for a principle to explain visual cortical functional architecture. First, many studies have reported good qualitative agreement between the layout of numerically obtained dimension reducing maps and experimental observations [19, 21, 42, 55–57, 59–71]. Second, geometric relationships between the representations of different visual features such as orientation, spatial frequency, and OD have been reproduced by dimension reduction models [25, 56, 63–65, 67, 68, 72].

Mathematically, the dimension reduction hypothesis implies that the layouts of cortical maps can be understood as optima or near optima (global or local minima) of a free energy functional which penalizes both 'stimulus scotomas' and map discontinuity. Unfortunately, there is currently no dimension reduction model for which the precise layouts of optimal or nearly optimal solutions have analytically been established. To justify the conclusion that the tradeoff between coverage and continuity favors the common rules of OPM design found in experiment, knowledge of optimal dimension-reducing mappings however appears essential.

Precise knowledge of the spatial organization of optimal and nearly optimal mappings is also critical for distinguishing between optimization theories and frozen noise scenarios of visual cortical development. In a frozen noise scenario, essentially random factors such as haphazard wiring [73], the impact of spontaneous activity patterns [74], or an idiosyncratic set of visual experiences [75] determine the emerging pattern of preferred orientations. This pattern is then assumed to be 'frozen' by an unknown mechanism, capable of preventing further modification of preferred orientations by ongoing synaptic turnover and activity-dependent plasticity. Conceptually, a frozen noise scenario is diametrically opposed to any kind of optimization theory. Even if the reorganization of the pattern prior to freezing was to follow a gradient descent with respect to some cost function, the early stopping implies that the layout is neither a local nor a global minimum of this functional. Importantly, the layout of transient states is known to exhibit universal properties that can be completely independent of model details [25, 32]. As a consequence, an infinite set of distinct optimization principles will generate the same spatial structure of transient states. This implies in turn that the frozen transient layout is not specifically shaped by any particular optimization principle. Map layouts will thus in principle only be informative about design or optimization principles of cortical processing architectures if maps are not just frozen transients.

In practice, however, the predictions of frozen noise and optimization theories might be hard to distinguish. Ambiguity between these mutually exclusive theories would result in particular, if the energy landscape of the optimization principle is so 'rugged' that there is essentially a local energy minimum next to any relevant random arrangement. Dimension reduction models are conceptually related to combinatorial optimization problems like the traveling salesman problem (TSP) and many such problems are believed to exhibit rugged energy landscapes [76–78]. It is therefore essential to clarify whether paradigmatic dimension reduction models are characterized by a rugged or a smooth energy landscape. If their energy landscapes were smooth with a small number of well-separated local minima, their predictions would be easy to distinguish from those of a frozen noise scenario.

In this study, we examine the classical example of a dimension reduction model, the elastic network (EN) model. Since the seminal work of Durbin and Mitchison [21], the EN model has widely been used to study visual cortical representations [25, 42, 62–65, 69–72, 79]. The EN model possesses an explicit energy functional which trades off a matching constraint which matches cortical cells to particular stimulus features via Hebbian learning, with a continuity constraint that minimizes Euclidean differences in feature space between neighboring points in the cortex [63]. In two ways, the EN model's explicit variational structure is very appealing. First, coverage and continuity appear as separate terms in the free energy which facilitates the dissection of their relative influences. Second, the free energy allows for the formulation of a gradient descent dynamics. The emergence of cortical selectivity patterns and their convergence toward a minimal energy state in this dynamics might serve as a model for an optimization process taking place in postnatal development.

Following Durbin and Mitchison, we consider the EN model for the joint mapping of two visual features: (i) position in visual space, represented in a retinotopic map (RM) and (ii) line orientation, represented in an OPM. To compute optimal dimension-reducing mappings, we first develop an analytical framework for deriving closed-form expressions for fixed points, local minima, and optima of arbitrary optimization models for the spatial layout of OPMs and RMs in which predicted maps emerge by a supercritical bifurcation as well as for analyzing their stability properties. By applying this framework to different instantiations of the EN model, we systematically disentangle the effects of individual model features on the repertoire of optimal solutions. We start with the simplest possible model version, a fixed uniform retinotopy and an orientation stimulus ensemble with only a single orientation energy and then relax the uniform retinotopy assumption incorporating retinotopic distortions. An analysis for a second widely used orientation stimulus ensemble including also unoriented stimuli is given in Appendix 1. Surprisingly, in all cases, our analysis yields pinwheel-free orientation stripes (OSs) or stereotypical square arrays of pinwheels as local minima or optimal orientation maps of the EN model. Numerical simulations of the EN confirm these findings. They indicate that more complex spatially aperiodic solutions are not dominant and that the energy landscape of the EN model is rather smooth. Our results demonstrate that while aperiodic stationary states exist, they are generally unstable in the considered model versions.

To test whether the EN model is in principle capable of generating complex spatially aperiodic optimal orientation maps, we then perform a comprehensive unbiased search of the EN optima for arbitrary orientation stimulus distributions. We identify two key parameters determining pattern selection: (i) the intracortical interaction range and (ii) the fourth moment of the orientation stimulus distribution function. We derive complete phase diagrams summarizing pattern selection in the EN model for fixed as well as variable retinotopy. Small interaction ranges together with low to intermediate fourth moment values lead to pinwheel-free OSs, rhombic, or hexagonal crystalline orientation map layouts as optimal states. In the regime of large interaction ranges together with higher fourth moment, we find irregular aperiodic OPM layouts with low pinwheel densities as optima. Only in an extreme and previously unconsidered parameter regime of very large interaction ranges and stimulus ensemble distributions with an infinite fourth moment, optimal OPM layouts in the EN model resemble experimentally observed aperiodic pinwheel-rich layouts and quantitatively reproduce the recently described species-invariant pinwheel statistics. Unexpectedly, we find that the stabilization of such layouts is not achieved by an optimal tradeoff between coverage and continuity of a localized population encoding by the maps but results from effectively suppressive long-range intracortical interactions in a spatially distributed representation of localized stimuli.

We conclude our reexamination of the EN model with a comparison between different numerical schemes for the determination of optimal or nearly optimal mappings. For large numbers of stimuli, numerically determined solutions match our analytical predictions, irrespectively of the computational method used.

Model definition and model symmetries

We analyze the EN model for the joint optimization of position and orientation selectivity as originally introduced by Durbin and Mitchison [21]. In this model, the RM is represented by a mapping R(x) = (R₁(x), R₂(x)) which describes the receptive field center position of a neuron at cortical position x. Any RM can be decomposed into an affine transformation x ↦ X from cortical to visual field coordinates, on which a vector-field of retinotopic distortions r(x) is superimposed, i.e.:

with appropriately chosen units for x and R.

The OPM is represented by a second complex-valued scalar field z(x). The pattern of orientation preferences ϑ(x) is encoded by the phase of z(x) via

The absolute value |z(x)| is a measure of the average cortical selectivity at position x. Solving the EN model requires to find pairs of maps {r(x), z(x)} that represent an optimal compromise between stimulus coverage and map continuity. This is achieved by minimizing a free energy functional

in which the functional \mathcal{C} measures the coverage of a stimulus space and the functional \mathcal{R} the continuity of its cortical representation. The stimulus space is defined by an ensemble {S} of idealized point-like stimuli, each described by two features: s _z = |s _z |e^2iθand s _r = (s _x ,s _y ) which specify the orientation θ of the stimulus and its position in visual space s _r (Figure 2b). \mathcal{C} and \mathcal{R} are given by

The EN model. (a) Example OPM (color code) together with a uniform map of visual space (RM) (grid lines). (b) Position s _r = (s _x , s _y ) and orientation θ of a 'pointlike' stimulus. (c) Cortical activity, evoked by the stimulus in b for the model maps in a. Dark regions are activated. Note, that in contrast to SOFM models, the activity pattern does not exhibit a stereotypical Gaussian shape. (d) Modification of orientation preference and retinotopy, caused by the stimulus in b. Orientation preferences prior to stimulus presentation are indicated with grey bars, after stimulus presentation with black bars. Most strongly modified preferences correspond to thick black bars. Modifications of orientation preferences and retinotopy are displayed on an exaggerated scale for illustration purposes.

Full size image

with ∇ = (∂ _x , ∂ _y ) ^T , and η ∈ [0, 1]. The ratios σ ²/η and σ ²/η _r control the relative strength of the coverage term versus the continuity term for OPM and RM, respectively. 〈...〉 _S denotes the average over the ensemble of stimuli.

Minima of the energy functional \mathcal{F} are stable fixed points of the gradient descent dynamics

called the EN dynamics in the following. These dynamics read

where e(x, S, z, r) is the activity-pattern, evoked by a stimulus S = (s _r , s _z ) in a model cortex with retinotopic distortions r(x) and OPM z(x). It is given by

Figure 2 illustrates the general features of the EN dynamics using the example of a single stimulus. Figure 2a shows a model orientation map with a superimposed uniform representation of visual space. A single point-like, oriented stimulus S = (s _r , s _z ) with position s _r = (s _x , s _y ) and orientation θ = 1/2 arg(s _z ) (Figure 2b) evokes a cortical activity pattern e(x, S, z, r) (Figure 2c). The activity-pattern in this example is of roughly Gaussian shape and is centered at the point, where the location s _r of the stimulus is represented in cortical space. However, depending on the model parameters and the stimulus, the cortical activity pattern may assume as well a more complex form (see also 'Discussion' section). According to Equations (3, 4), each stimulus and the evoked activity pattern induce a modification of the orientation map and RM, shown in Figure 2d. Orientation preference in the activated regions is shifted toward the orientation of the stimulus. The representation of visual space in the activated regions is locally contracted toward the position of the stimulus. Modifications of cortical selectivities occur due to randomly chosen stimuli and are set proportional to a very small learning rate. Substantial changes of cortical representations occur slowly through the cumulative effect of a large number of activity patterns and stimuli. These effective changes are described by the two deterministic equations for the rearrangement of cortical selectivities equations (3, 4) which are obtained by stimulus-averaging the modifications due to single activity patterns in the discrete stimulus model [25]. One thus expects that the optimal selectivity patterns and also the way in which cortical selectivities change over time are determined by the statistical properties of the stimulus ensemble. In the following, we assume that the stimulus ensemble satisfies three properties: (i) The stimulus locations s _r are uniformly distributed across visual space. (ii) For the distribution of stimulus orientations, |s _z | and θ are independent. (iii) Orientations θ are distributed uniformly in [0, π].

These conditions are fulfilled by stimulus ensembles used in virtually all prior studies of dimension reduction models for visual cortical architecture (e.g., [19, 21, 25, 64, 65, 71, 72, 80, 81]). They imply several symmetries of the model dynamics equations (3, 4). Due to the first property, the EN dynamics are equivariant under translations

rotations

with 2×2 rotation matrix

and reflections

where Ψ = diag(-1, 1) is the 2×2 reflection matrix. Equivariance means that

with mutatis mutandis the same equations fulfilled by the vector-field F^r [z, r].

As a consequence, patterns that can be converted into one another by translation, rotation, or reflection of the cortical layers represent equivalent solutions of the model equations (3, 4), by construction. Due to the second assumption, the dynamics is also equivariant with respect to shifts in orientation S _ϕ z(x) = e^iϕz(x), i.e.,

Thus, two patterns are also equivalent solutions of the model, if their layout of orientation domains and retinotopic distortions is identical, but the preferred orientations differ everywhere by the same constant angle.

Without loss of generality, we normalize the ensembles of orientation stimuli such that {〈\mid s_z{\mid }^2〉}_{\mathbf{S}}=〈\mid s_z{\mid }^2〉=2 throughout this article. This normalization can always be restored by a rescaling of z(x) (see [25, 69]).

Our formulation of the dimension-reduction problem in the EN model utilizes a continuum description, both for cortical space and the set of visual stimuli. This facilitates mathematical treatment and appears appropriate, given the high number of cortical neurons under one square millimeter of cortical surface (e.g., roughly 70000 in cat V1 [82]). Even an hypothesized neuronal monolayer would consist of more than 20 × 20 neurons per hypercolumn area Λ², constituting a quite dense sampling of the spatial periodicity. Treating the feature space as a continuum implements the concept that the cortical representation has to cover as good as possible the infinite multiplicity of conceivable stimulus feature combinations.

The orientation unselective fixed point

Two stationary solutions of the model can be established from symmetry. The simplest of these is the orientation unselective state with z(x) = 0 and uniform mapping of visual space r(x) = 0. First, by the shift symmetry (Equation (8)), we find that z(x) = 0 is a fixed point of Equation (3). Second, by reflectional and rotational symmetry (Equations (5, 7)), we see that the right-hand side of Equation (4) has to vanish and hence the orientation unselective state with uniform mapping of visual space is a fixed point of Equations (3, 4).

This homogeneous unselective state thus minimizes the EN energy functional if it is a stable solution of Equations (3, 4). The stability can be determined by considering the linearized dynamics of small deviations {r(x), z(x)} around this state. These linearized dynamics read

where (Â) _ij = (x _i - y _i )(x _j - y _j ) -2σ²δ _ij with δ _ij being Kronecker's delta. We first note that the linearized dynamics of retinotopic distortions and orientation preference decouple. Thus, up to linear order and near the homogeneous fixed point, both cortical representation evolve independently and the stability properties of the unselective state can be obtained by a separate examination of the stability properties of both cortical representations.

The eigenfunctions of the linearized retinotopy dynamics L _r [r] can be calculated by Fourier-transforming Equation (10):

where k = |k| and i = 1, 2. A diagonalization of this matrix equation yields the eigenvalues

with corresponding eigenfunctions (in real space)

where k _ϕ = |k|(cos ϕ, sin ϕ) ^T . These eigenfunctions are longitudinal (L) or transversal (T) wave patterns. In the longitudinal wave, the retinotopic distortion vector r(x) lies parallel to k which leads to a 'compression wave' (Figure 3b, left). In the transversal wave pattern (Figure 3b, right), the retinotopic distortion vector is orthogonal to k. We note that the linearized Kohonen model [61] was previously found to exhibit the same set of eigenfunctions [80]. Because both spectra of eigenvalues {\lambda }_T^r, {\lambda }_L^rare smaller than zero for every σ >0, η _r >0, and k >0 (Figure 3a), the uniform retinotopy r(x) = 0 is a stable fixed point of Equation (4) irrespective of parameter choice.

The linearization of the EN model dynamics around the unselective fixed point. (a) Eigenvalue spectra of the linearized retinotopy dynamics for longitudinal mode ({\lambda }_L^r\left(k\right), blue trace) and transversal mode ({\lambda }_T^r\left(k\right), red trace). (b) Longitudinal mode ~{\mathbf{k}}_{\varphi }{e}^{i{\mathbf{k}}_{\varphi }\mathbf{x}}+c.c. (left) and transversal mode ~{\mathbf{k}}_{\varphi +\pi ∕2}{e}^{i{\mathbf{k}}_{\varphi }\mathbf{x}}+c.c. (right). (c) Spectrum of eigenvalues of the linearized OPM dynamics (red trace) for σ < σ*(η). Orange region marks the unstable annulus of Fourier modes (critical circle). (d) Stability regions of the nonselective state in the EN model. The stability border is given by σ*(η) (Equation (12)).

Full size image

The eigenfunctions of the linearized OPM dynamics L _z [z] are Fourier modes ~ e^{i kx} by translational symmetry. By rotational symmetry, their eigenvalues only depend on the wave number k and are given by

(see [25]). This spectrum of eigenvalues is depicted in Figure 3c. For η >0, λ^z (k) has a single maximum at k_c=\frac{1}{\sigma }\sqrt{ln\left(1∕\eta \right)}. For

this maximal eigenvalue r = λ^z (k _c ) is negative. Hence, the unselective state with uniform retinotopy is a stable fixed point of Equations (3, 4) and the only known solution of the EN model in this parameter range. For σ < σ*(η), the maximal eigenvalue r is positive, and the nonselective state is unstable with respect to a band of Fourier modes ~ e^{i kx} with wave numbers around |k| ≈ k _c (see Figure 3c). This annulus of unstable Fourier modes is called the critical circle. The finite wavelength instability [83–85] (or Turing instability [86]) leads to the emergence of a pattern of orientation preference with characteristic spacing Λ = 2π/k _c from the nonselective state on a characteristic timescale τ = 1/r.

One should note that as in other models for the self-organization of orientation columns, e.g., [15, 57], the characteristic spatial scale Λ arises from effective intracortical interactions of 'Mexican-hat' structure (short-range facilitation, longer-ranged suppression). The short-range facilitation in the linearized EN dynamics is represented by the first two terms on the right-hand side of Equation (11). Since σ <1 in the pattern forming regime, the prefactor in front of the first term is positive. Due to the second, Laplacian term, it is favored that neighboring units share selectivity properties, a process mediated by short-range facilitation. Longer-ranged suppression is represented by the convolution term in Equation (11).

Mathematically, this term directly results from the soft-competition in the 'activity-dependent' coverage term of Equation (1). The local facilitation is jointly mediated by coverage (first term) and continuity (second term) contributions.

Figure 3d summarizes the result of the linear stability analysis of the nonselective state. For σ > σ*(η), the orientation unselective state with uniform retinotopy is a minimum of the EN-free energy and also the global minimum. For σ < σ*(η), this state represents a maximum of the energy functional and the minima must thus exhibit a space-dependent pattern of orientation selectivities.

Orientation stripes

Within the potentially infinite set of orientation selective fixed points of the model, one class of solutions can be established from symmetry: {r(x) = 0, z(x) = A₀e^{i kx}}. In these pinwheel-free states, orientation preference is constant along one axis in cortex (perpendicular to the vector k), and each orientation is represented in equal proportion (see Figure 4a). Retinotopy is perfectly uniform. Although this state may appear too simple to be biologically relevant, we will see that it plays a fundamental role in the state space of the EN model. It is therefore useful to establish its existence and basic characteristics. The existence of OS solutions follows directly from the model's symmetries (Equations (5) to (9)). Computing

Exact and approximate orientation selective fixed points of OPM optimization models. (a) Pinwheel-free OS pattern. Diagram shows the position of the wave vector in Fourier space. (b) rPWC with four nonzero wave vectors. (c) Essentially complex planforms (ECPs). The index n indicates the number of nonzero wave vectors. The index i enumerates nonequivalent configurations of wave vectors with the same n, starting with i = 0 for the most anisotropic planform. For n = 3, 5, and 15, there are 2, 4, and 612 different ECPs, respectively. OPM layouts become more irregular with increasing n.

Full size image

demonstrates that F^z [e^{i kx}, 0] is proportional to e^{i kx}. This establishes that the subspace of functions ~ e^{i kx} is invariant under the dynamics given by Equation (3). For A₀ = 0, we recover the trivial fixed point of the EN dynamics by construction, as shown above. This means that within this subspace A₀ = 0 is either a minimum or a maximum of the EN energy functional (Equation (1)). Furthermore, for A₀ → ∞ the EN energy tends to infinity. If the trivial fixed point is unstable, it corresponds to a maximum of the EN energy functional. Therefore, there must exist at least one minimum with A₀ ≠ 0 in the subspace of functions ~ e^{i kx} which then corresponds to a stationary state of the EN dynamics.

Regarding the dynamics of retinotopic deviations, the model's symmetries equations can be invoked to show that for the state {0, A₀e^{i kx}}, the right-hand side of Equation (4) has to be constant in space:

If this constant was nonzero the RM would drift with constant velocity. This, however, is impossible in a variational dynamics such that this constant must vanish. The OS solution (Figure 4a) is to the best of our knowledge the only exact nontrivial stationary solution of Equations (3, 4) that can be established without any approximations.

Doubly periodic and quasi-periodic solutions

In the EN model as considered in this study, the maps of visual space and orientation preference are jointly optimized to trade off coverage and continuity leading to mutual interactions between the two cortical representations. These mutual interactions vanish in the rigid retinotopy limit η _r → ∞ and the perfectly uniform retinotopy becomes an optimal solution for arbitrary orientation column layout z(x). As it is not clear how essential the mutual interactions with position specificity are in shaping the optimal orientation column layout, we continue our investigation of solution classes by considering global minima of optimization models with fixed uniform retinotopy. The mutual interactions will be taken into account in a subsequent step.

In the rigid retinotopy limit, minima of the energy functional are stable stationary states of the dynamics of the OPM (Equation (3)) with r(x) = 0. To compute orientation selective stationary solutions of this OPM dynamics, we employ that in the vicinity of a supercritical bifurcation where the nonselective fixed point z(x) = 0 becomes unstable, the entire set of nontrivial fixed points is determined by the third-order terms of the Volterra series representation of the operator F^z [z, 0] [35, 84, 85, 87]. The symmetries given by Equations (5) to (9) restrict the general form of such a third-order approximation for any model of OPM optimization to

where the cubic operator N_3^z is written in trilinear form, i.e.,

In particular, all even terms in the Volterra Series representation of F^z [z, 0] vanish due to the Shift-Symmetry (Equations (8, 9)). Explicit analytic computation of the cubic nonlinearities for the EN model is cumbersome but not difficult (see 'Methods' section) and yields a sum

The individual nonlinear operators N_3^j are with one exception nonlocal convolution-type operators and are given in the 'Methods' section (Equation (38)), together with a detailed description of their derivation.

Only the coefficients a _j depend on the properties of the ensemble of oriented stimuli.

To calculate the fixed points of Equation (13), we use a perturbative method called weakly nonlinear analysis that enables us to analytically examine the structure and stability of inhomogeneous stationary solutions in the vicinity of a finite-wavelength instability. Here, we examine the stability of so-called planforms [83–85]. Planforms are patterns that are composed of a finite number of Fourier components, such as

for a pattern of orientation columns. With the above planform ansatz, we neglect any spatial dependency of the amplitudes A _j (t) for example due to long-wave deformations for the sake of simplicity and analytical tractability. When the dynamics is close to a finiite wavelength instability, the essential Fourier components of the emerging pattern are located on the critical circle |k _j | = k _c . The dynamic equations for the amplitudes of these Fourier components are called amplitude equations. For a discrete number of N Fourier components of z(x) whose wave vectors lie equally spaced on the critical circle, the most general system of amplitude equations compatible with the model's symmetries (Equations (5) to (9)) has the form [35, 87]

with r >0. Here, g _ij and f _ij are the real-valued coupling coefficients between the amplitudes A _i and A _j . They depend on the differences between indices |i - j| and are entirely determined by the nonlinearity N_3^z\left[z,z,\bar{z}\right] in Equation (13). If the wave vectors k _i = (cos α _i , sin α _i )k _c are parameterized by the angles α _i , then the coefficients g _ij and f _ij are functions only of the angle α = |α _i - α _j | between the wave vectors k _i and k _j . One can thus obtain the coupling coefficients from two continuous functions g(α) and f(α) that can be obtained from the nonlinearity N_3^z\left[z,z,\bar{z}\right] (see 'Methods' section for details). In the following, these functions are called angle-dependent interaction functions. The amplitude equations are variational if and only if g _ij and f _ij are real-valued. In this case they can be derived through

from an energy

If the coefficients g _ij and f _ij are derived from Equation (1), the energy U _A for a given planform solution corresponds to the energy density of the EN energy functional considering only terms up to fourth-order in z(x).

The amplitude equations (15) enable to calculate an infinite set of orientation selective fixed points. For the above OS solution with one nonzero wave vector z(x) = A₀e^{i kx}, the amplitude equations predict the so far undetermined amplitude

and its energy

Since g _ii >0, this shows that OS stationary solutions only exist for r >0, i.e., in the symmetry breaking regime. As for all following fixed-points, U_OS specifies the energy difference to the homogeneous unselective state z(x) = 0.

A second class of stationary solutions can be found with the ansatz

with amplitudes A _j = |A _j |e^iϕ _j and ∠(k₁, k₂) = α >0. By inserting this ansatz into Equation (15) and assuming uniform amplitude \mid A_1\mid =\mid A_2\mid =\mid A_2\mid =\mid A_4\mid =\mathcal{A}, we obtain

The phase relations of the four amplitudes are given by

These solutions describe a regular rhombic lattice of pinwheels and are therefore called rhombic pinwheel crystals (rPWCs) in the following. Three phases can be chosen arbitrarily according to the two above conditions, e.g., ϕ₀, Δ₀ = ϕ₁ - ϕ₃ and Δ₁ = ϕ₂ - ϕ₄. For an rPWC parameterized by these phases, Δ₀ shifts the absolute positions of the pinwheels in x-direction, Δ₁ shifts the absolute positions of the pinwheels in y-direction, and ϕ₀ shifts all the preferred orientations by a constant angle. The energy of an rPWC solution is

An example of such a solution is depicted in Figure 4b. We note that rPWCs have been previously found in several other models for OPM development [27, 31, 37, 39, 88]. The pinwheel density ρ of an rPWC, i.e., the number of pinwheels in an area of size Λ², is equal to ρ = 4 sin α [54]. The angle α which minimizes the energy U_rPWC can be computed by maximizing the function

in the denominator of Equation (20).

The two solution classes discussed so far, namely OS and rPWCs, exhibit one prominent feature, absent in experimentally observed cortical OPMs, namely perfect spatial periodicity. Many cortical maps including OPMs do not resemble a crystal-like grid of repeating units. Rather the maps are characterized by roughly repetitive but aperiodic spatial arrangement of feature preferences (e.g., [5, 10]). This does not imply that the precise layout of columns is arbitrary. It rather means that the rules of column design cannot be exhaustively characterized by mapping a 'representative' hypercolumn.

Previous studies of abstract models of OPM development introduced the family of so-called essentially complex planforms (ECPs) as stationary solutions of Equation (15). This solution class encompasses a large variety of realistic quasi-periodic OPM layouts and is therefore a good candidate solution class for models of OPM layouts. In addition, Kaschube et al. [38] demonstrated that models in which these are optimal solutions can reproduce all essential features of the common OPM design in ferret, tree-shrew, and galago. An n-ECP solution can be written as

with n = N/2 wave vectors k _j = k _c (cos(πj/n), sin(πj/n)) distributed equidistantly on the upper half of the critical circle, complex amplitudes A _j and binary variables l _j = ±1 determining whether the mode with wave vector k _j or -k _j is active (nonzero). Because these planforms cannot realize a real-valued function they are called essentially complex [35]. For an n-ECP, the third term on the right-hand side of Equation (15) vanishes and the amplitude equations for the active modes A _i reduce to a system of Landau equations

where g _ij is the n × n-coupling matrix for the active modes. Consequently, the stationary amplitudes obey

The energy of an n-ECP is given by

We note that this energy in general depends on the configuration of active modes, given by the l _j 's, and therefore planforms with the same number of active modes may not be energetically degenerate.

Families of n-ECP solutions are depicted in Figure 4c. The 1-ECP corresponds to the pinwheel-free OS pattern discussed above. For fixed n ≥ 3, there are multiple planforms not related by symmetry operations which considerably differ in their spatial layouts. For n ≥ 4, the patterns are spatially quasi-periodic, and are a generalization of the so-called Newell-Pomeau turbulent crystal [89, 90]. For n ≥ 10, their layouts resemble experimentally observed OPMs. Different n-ECPs however differ considerably in their pinwheel density. Planforms whose nonzero wave vectors are distributed isotropically on the critical circle typically have a high pinwheel density (see Figure 4c, n = 15 lower right). Anisotropic planforms generally contain considerably fewer pinwheels (see Figure 4c, n = 15 lower left). All large n-ECPs, however, exhibit a complex quasi-periodic spatial layout and a nonzero density of pinwheels.

In order to demonstrate that a certain planform is an optimal solution of an optimization model for OPM layouts in which patterns emerge via a supercritical bifurcation, we not only have to show that it is a stationary solution of the amplitude equations but have to analyze its stability properties with respect to the gradient descent dynamics as well as its energy compared to all other candidate solutions.

Many stability properties can be characterized by examining the amplitude equations (15). In principle, the stability range of an n-ECPs may be bounded by two different instability mechanisms: (i) an intrinsic instability by which stationary solutions with n active modes decay into ones with lower n. (ii) an extrinsic instability by which stationary solutions with a 'too low' number of modes are unstable to the growth of additional active modes. These instabilities can constrain the range of stable n to a small finite set around a typical n [35, 87]. A mathematical evaluation of both criteria leads to precise conditions for extrinsic and intrinsic stability of a planform (see 'Methods' section). In the following, a planform is said to be stable, if it is both extrinsically and intrinsically stable. A planform is said to be an optimum (or optimal solution) if it is stable and possesses the minimal energy among all other stationary planform solutions.

Taken together, this amplitude equation approach enables to analytically compute the fixed points and optima of arbitrary optimization models for visual cortical map layout in which the functional architecture is completely specified by the pattern of orientation columns z(x) and emerges via a supercritical bifurcation. Via a third-order expansion of the energy functional together with weakly nonlinear analysis, the otherwise analytically intractable partial integro-differential equation for OPM layouts reduces to a much simpler system of ordinary differential equations, the amplitude equations. Using these, several families of solutions, OSs, rPWCs, and essentially complex planforms, can be systematically evaluated and comprehensively compared to identify sets of unstable, stable and optimal, i.e., lowest energy fixed points. As already mentioned, the above approach is suitable for arbitrary optimization models for visual cortical map layout in which the functional architecture is completely specified by the pattern of orientation columns z(x) which in the EN model is fulfilled in the rigid retinotopy limit. We now start by considering the EN optimal solutions in this limit and subsequently generalize this approach to models in which the visual cortical architecture is jointly specified by maps of orientation and position preference that are matched to one another.

Representing an ensemble of 'bar'-stimuli

We start our investigation of optimal dimension-reducing mappings in the EN model using the simplest and most frequently used orientation stimulus ensemble, the distribution with s _z -values uniformly arranged on a ring with radius r_{s_z}=\sqrt{2}[57, 64–66, 91]. We call this stimulus ensemble the circular stimulus ensemble in the following. According to the linear stability analysis of the nonselective fixed point, at the point of instability, we choose σ = σ*(η) such that the linearization given in Equation (11) is completely characterized by the continuity parameter η. Equivalent to specifying η is to fix the ratio of activation range σ and column spacing Λ

as a more intuitive parameter. This ratio measures the effective interaction-range relative to the expected spacing of the orientation preference pattern. In abstract optimization models for OPM development a similar quantity has been demonstrated to be a crucial determinant of pattern selection [35, 87]. We note, however, that due to the logarithmic dependence of σ/Λ on η, a slight variation of the effective interaction range may correspond to a variation of the continuity parameter η over several orders of magnitude. In order to investigate the stability of stationary planform solutions in the EN model with a circular orientation stimulus ensemble, we have to determine the angle-dependent interaction functions g(α) and f(α). For the coefficients a _j in Equation (14) we obtain

The angle-dependent interaction functions of the EN model with a circular orientation stimulus ensemble are then given by

These functions are depicted in Figure 5 for two different values of the interaction range σ/Λ. We note that both functions are positive for all σ/Λ which is a sufficient condition for a supercritical bifurcation from the homogeneous nonselective state in the EN model.

Angle-dependent interaction functions for the EN model with fixed retinotopy and circular orientation stimulus ensemble. (a, b) g(α) and f(α) for σ/Λ = 0.1 (a) and σ/Λ = 0.2 (b).

Full size image

Finally, by minimizing the function s(α) in Equation (21), we find that the angle α which minimizes the energy of the rPWC fixed-point is α = π/2. This corresponds to a square array of pinwheels (sPWC). Due to the orthogonal arrangement oblique and cardinal orientation columns and the maximized pinwheel density of ρ = 4, the square array of pinwheels has the maximal coverage among all rPWC solutions.

Optimal solutions close to the pattern formation threshold

We first tested for the stability of pinwheel-free OS solutions and the sPWCs, by analytical evaluation of the criteria for intrinsic and extrinsic stability (see 'Methods' section). We found both, OSs and sPWCs, to be intrinsically and extrinsically stable for all σ/Λ. Next, we tested for the stability of n-ECP solutions with 2 ≤ n ≤ 20. We found all n-ECP configurations with 2 ≤ n ≤ 20 to be intrinsically unstable for all σ/Λ. Hence, none of these planforms represent optimal solutions of the EN model with a circular stimulus ensemble, while both OSs and sPWC are always local minima of the energy functional.

By evaluating the energy assigned to the sPWC (Equation 20) and the OS pattern (Equation 18), we next identified two different regimes: (i) For short interaction range σ/Λ ≲ 0.122 the sPWC possesses minimal energy and is therefore the predicted global minimum. (ii) For σ/Λ ≳ 0.122 the OS pattern is optimal.

Figure 6a shows the resulting simple phase diagram. sPWCs and OSs are separated by a phase border at σ/Λ ≈ 0.122. We numerically confirmed these analytical predictions by extensive simulations of Equation (3) with r(x) = 0 and the circular stimulus ensemble (see 'Methods' section for details). Figure 6b,c shows snapshots of a representative simulation with short interaction range (r = 0.1, σ/Λ = 0.1 (η = 0.67)) (see also Additional file 1). After the phase of initial pattern emergence (symmetry breaking), the OPM layout rapidly approaches a square array of pinwheels, the analytically predicted optimum (Figure 6d). Pinwheel density time courses (see 'Methods' section) display a rapid convergence to a value close to the predicted density of 4 (Figure 6e). Figure 6f shows the stationary mean squared amplitudes of the pattern obtained for different values of the control parameter r (black circles). For small control parameters, the pattern amplitude is perfectly predicted by Equation (19) (solid green line). Figure 6g,h shows snapshots of a typical simulation with longer interaction range (r = 0.1, σ/Λ = 0.15 (η = 0.41)) (see also Additional file 2). After the emergence of an OPM with numerous pinwheels, pinwheels undergo pairwise annihilation as previously described for various models of OPM development and optimization [25, 27, 35]. The OP pattern converges to a pinwheel-free stripe pattern, which is the analytically computed optimal solution in this parameter regime (Figure 6i). Pinwheel densities decay toward zero over the time course of the simulations (Figure 6j). Also in this parameter regime, the mean squared amplitude of the pattern is well-predicted by Equation (17) for small r (Figure 6k).

Optimal solutions of the EN model with a circular orientation stimulus ensemble [57, 64–66, 91]and fixed representation of visual space. (a) At criticality, the phase space of this model is parameterized by either the continuity parameter η (blue labels) or the interaction range σ/Λ (red labels, see text). (b, c) OPMs (b) and their power spectra (c) in a simulation of Equation (3) with r(x) = 0 and r = 0.1, σ/Λ = 0.1 (η = 0.67) and circular stimulus ensemble (see also Additional file 1). (d) Analytically predicted optimum for σ/Λ ≲ 0.122 (quadratic pinwheel crystal). (e) Pinwheel density time courses for four different simulations (parameters as in b; gray traces, individual realizations; black trace, simulation in b; red trace, mean value). (f) Mean squared amplitude of the stationary pattern, obtained in simulations (parameters as in b) for different values of the control parameter r (black circles) and analytically predicted value (solid green line). (g, h) OPMs (g) and their power spectra (h) in a simulation of Equation (3) with r(x) = 0 and σ/Λ = 0.15 (η = 0.41) (other parameters as in b, see also Additional file 2). (i) Analytically predicted optimum for σ/Λ ≳ 0.122 (orientation stripes). (j) Pinwheel density time courses for four different simulations (parameters as in g; gray traces, individual realizations; black trace, simulation in g; red trace, mean value). (k) Mean squared amplitude of the stationary pattern, obtained in simulations (parameters as in g) for different values of the control parameter r (black circles) and analytically predicted value (solid green line).

Full size image

In summary, the phase diagram of the EN model with a circular stimulus ensemble close to threshold is divided into two regions: (i) for a small interaction range (large continuity parameter) a square array of pinwheels is the optimal dimension-reducing mapping and (ii) for a larger interaction range (small continuity parameter) OSs are the optimal dimension-reducing mapping. Both states are stable throughout the entire parameter range. All other planforms, in particular quasi-periodic n-ECPs are unstable. At first sight, this structure of the EN phase diagram may appear rather counterintuitive. A solution with many pinwheel-defects is energetically favored over a solution with no defects in a regime with large continuity parameter where discontinuity should be strongly penalized in the EN energy term. However, a large continuity parameter at pattern formation threshold inevitably leads to a short interaction range σ compared to the characteristic spacing Λ (see Equation (24)). In such a regime, the gain in coverage by representing many orientation stimuli in a small area spanning the typical interaction range, e.g., with a pinwheel, is very high. Our results show that the gain in coverage by a spatially regular positioning of pinwheels outweighs the accompanied loss in continuity above a certain value of the continuity parameter.

EN dynamics far from pattern formation threshold

Close to pattern formation threshold, we found only two stable solutions, namely OSs and sPWCs. Neither of the two exhibits the characteristic aperiodic and pinwheel-rich organization of experimentally observed OPMs. Furthermore, the pinwheel densities of both solutions (ρ = 0 for OSs and ρ = 4 for sPWCs) differ considerably from experimentally observed values [38] around 3.14. One way toward more realistic stable stationary states might be to increase the distance from pattern formation threshold. In fact, further away from threshold, our perturbative calculations may fail to correctly predict optimal solutions of the model due to the increasing influence of higher order terms in the Volterra series expansion of the right-hand side in Equation (3).

To asses this possibility, we simulated Equation (3) with r(x) = 0 and a circular stimulus ensemble for very large values of the control parameter r. Figure 7 displays snapshots of such a simulation for r = 0.8 as well as their pinwheel density time courses for two different values of σ/Λ. Pinwheel annihilation in the case of large σ/Λ is less rapid than close to threshold (Figure 7a,b). The OPM nevertheless converges toward a layout with rather low pinwheel density with pinwheel-free stripe-like domains of different directions joined by domains with essentially rhombic crystalline pinwheel arrangement. The linear zones increase their size over the time course of the simulations, eventually leading to stripe-patterns for large simulation times. For smaller interaction ranges σ/Λ, the OPM layout rapidly converges toward a crystal-like rhombic arrangement of pinwheels, however containing several dislocations (Figure 24a in Appendix 1) [84]. Dislocations are defects of roll or square patterns, where two rolls or squares merge into one, thus increasing the local wavelength of the pattern [83, 85]. Nevertheless, for all simulations, the pinwheel density rapidly reaches a value close to 4 (Figure 7c) and the square arrangement of pinwheels is readily recognizable. Both features, the dislocations in the rhombic patterns and domain walls in the stripe patterns, have been frequently observed in pattern-forming systems far from threshold [84, 85].

Numerical analysis of the EN dynamics with circular orientation stimulus ensemble and fixed representation of visual space far from pattern formation threshold. (a) OPMs and their power spectra in a simulation of Equation (3) with r(x) = 0, r = 0.8, σ/Λ = 0.3 (η = 0.028) and circular orientation stimulus ensemble. Pinwheel density time courses for four different simulations (parameters as in a; gray traces, individual realizations; black trace, simulation in a; red trace, mean value) (c, d) OPMs and their power spectra in a simulation of Equation (3) with r(x) = 0, r = 0.8, σ/Λ = 0.12 (η = 0.57) and circular orientation stimulus ensemble. (d) Pinwheel density time courses for four different simulations (parameters as in c; gray traces, individual realizations; black trace, simulation in c; red trace, mean value).

Full size image

In summary, the behavior of the EN dynamics with circular stimulus ensemble far from pattern formation threshold agrees very well with our analytical predictions close to threshold. Again, orientation stripes and square pinwheel crystals are identified as the only stationary solutions. Aperiodic and pinwheel-rich patterns which resemble experimentally OPM layouts were not observed.

Taking retinotopic distortions into account

So far, we have examined the optimal solutions of the EN model for the simplest and most widely used orientation stimulus ensemble. Somewhat unexpected from previous reports, the optimal states in this case do not exhibit the irregular structure of experimentally observed orientation maps. Our treatment however differs from previous approaches in that the mapping of visual space so far was assumed to be undistorted and fixed, i.e., r(x) = 0. We recall that in their seminal publication, Durbin and Mitchison [21] in particular demonstrated interesting correlations between the map of orientation preference and the map of visual space. These correlations suggest a strong coupling between the two that may completely alter the model's dynamics and optimal solutions.

It is thus essential to clarify whether the behavior of the EN model observed above changes or persists if we relax the simplifying assumption of undistorted retinotopy and allow for retinotopic distortions. By analyzing the complete EN model dynamics (Equations (3, 4)), we study the EN model exactly as originally introduced by Durbin and Mitchison [21].

We again employ the fact that in the vicinity of a supercritical bifurcation where the nonorientation selective state becomes unstable, the entire set of nontrivial fixed points of Equations (3, 4) is determined by the third-order terms of the Volterra series representation of the nonlinear operators F^z [z, r] and F^r [z, r]. The model symmetries equations (5) to (9) restrict the general form of the leading order terms for any model for the joint optimization of OPM and RM to

Because the uniform retinotopy is linearly stable, retinotopic distortions are exclusively induced by a coupling of the RM to the OPM via the quadratic vector-valued operator {\mathbf{Q}}^r\left[z,\bar{z}\right]. These retinotopic distortions will in turn alter the dynamics of the OPM via the quadratic complex-valued operator Q^z [r, z]. Close to the point of pattern onset (r ≪ 1), the timescale of OPM development, τ = 1/r, becomes arbitrarily large and retinotopic deviations evolve on a much shorter timescale. This separation of timescales allows for an adiabatic elimination of the variable r(x), assuming it to always be at the equilibrium point of Equation (27):

We remark that as for all finite wave numbers k >0, the operator L _r [r] is indeed invertible when excluding global translations in the set of possible perturbations of the trivial fixed point. From Equation (28), the coupled dynamics of OPM and RM is thus reduced to a third-order effective dynamics of the OPM:

The nonlinearity N_3^r\left[z,z,\bar{z}\right] accounts for the coupling between OPM and RM. Its explicit analytical calculation for the EN model is rather involved and yields a sum

The individual nonlinear operators N_r^j are nonlinear convolution-type operators and are presented in the 'Methods' section together with a detailed description of their derivation. Importantly, it turns out that the coefficients a_r^j are completely independent of the orientation stimulus ensemble.

The adiabatic elimination of the retinotopic distortions results in an equation for the OPM (Equation (29)) which has the same structure as Equation (13), the only difference being an additional cubic nonlinearity. Due to this similarity, its stationary solutions can be determined by the same methods as presented for the case of a fixed retinotopy. Again, via weakly nonlinear analysis we obtain amplitude equations of the form Equation (15). The nonlinear coefficients g _ij and f _ij are determined from the angle-dependent interaction functions g(α) and f(α). For the operator N_3^r\left[z,z,\bar{z}\right], these functions are given by

verifying that, N_3^r\left[z,z,\bar{z}\right] is independent of the orientation stimulus ensemble. Besides the interaction range σ/Λ the continuity parameter η _r ∈ [0, ∞] for the RM appears as an additional parameter in the angle-dependent interaction function. Hence, the phase diagram of the EN model will acquire one additional dimension when retinotopic distortions are taken into account. We note, that in the limit η _r → ∞, the functions g _r (α) and f _r (α) tend to zero and as expected one recovers the results presented above for fixed uniform retinotopy. The functions g _r (α) and f _r (α) are depicted in Figure 8 for various interaction ranges σ/Λ and retinotopic continuity parameters η _r .

Angle-dependent interaction functions for the coupling between OPM and RM in the EN model. (a, b) g _r (α) and f _r (α) for η _r = 0.005 and σ/Λ = 0.3 (a) and 0.1 (b). (c, d) g _r (α) and f _r (α) for η _r = 0.05 and σ/Λ = 0.3 (c) and 0.1 (d). (e, f) g _r (α) and f _r (α) for η _r = 0.5 and σ/Λ = 0.3 (e) and 0.1 (f).

Full size image

Coupled essentially complex n-planforms

In the previous section, we found that by an adiabatic elimination of the retinotopic distortions in the dynamics equations (26, 27) the system of partial integro-differential equations can be reduced to a single equation for the OPM. In this case, the stationary solutions of the OPM dynamics are again planforms composed of a discrete set of Fourier modes

with |k| = k _c . However, each of these stationary planform OPM solutions induces a specific pattern of retinotopic distortions by Equation (28). The joint mapping {X + r(x), z(x)} is then an approximate stationary solution of Equations (26, 27) and will be termed coupled planform solution in the following. In contrast to other models for the joint mapping of orientation and visual space (e.g., [31, 33, 92]), the coupling between the representation of visual space and orientation in the EN model is not induced by model symmetries but a mere consequence of the joint optimization of OPM and RM that requires them to be matched to one another.

For planforms given by Equation (30), it is possible to analytically evaluate Equation (28) and compute the associated retinotopic distortions r(x). After a somewhat lengthy calculation (see 'Methods' section), one obtains

with Δ _jk = k _j - k _k and {\lambda }_L^r\left(k\right)=-k^2\left({\eta }_r+e^{-{\sigma }^2{k}^2}{\sigma }^2\right). These retinotopic distortions represent superpositions of longitudinal modes (see Figure 3b). Hence, coupled planform stationary solutions of the EN dynamics do not contain any transversal mode components. According to Equation (31), the pinwheel-free coupled 1-ECP state has the functional form {r(x) = 0, z(x) = A₀e^{i kx}}. This means that the OS solution does not induce any deviations from the perfect retinotopy as shown previously from symmetry. This is not the case for the square pinwheel crystal (sPWC)

the second important solution for undistorted retinotopy. Inserting this ansatz into Equation (31) and neglecting terms of order \mathcal{O}\left({\left(e^{-k_c^2{\sigma }^2}\right)}^2\right) or higher, we obtain

These retinotopic distortions are a superposition of one longitudinal mode in x-direction and one in y-direction, both with doubled wave number ~ 2k _c . The doubled wave number implies that the form of retinotopic distortions is independent of the topological charge of the pinwheels. Importantly, the gradient of the retinotopic mapping R(x) = X + r_sPWC(x) is reduced at all pinwheel locations. The coupled sPWC is therefore in two ways a high coverage mapping as expected. First, the representations of cardinal and oblique stimuli (real and imaginary part of z(x)) are orthogonal to each other. Second, the regions of highest gradient in the orientation map correspond to low gradient regions in the RM.

In Figure 9, the family of coupled n-ECPs is displayed, showing simultaneously the distortions of the RM and the OPM. Retinotopic distortions are generally weaker for anisotropic n-ECPs and stronger for isotropic n-ECPs. However, for all stationary solutions the regions of high gradient in the orientation map coincide with low gradient regions (the folds of the grid) in the RM. This is precisely what is generally expected from a dimension-reducing mapping [21, 62, 63, 91]. In the following section, we will investigate which of these solutions become optimal depending on the two parameters σ/Λ and η _r that parameterize the model.

Coupled n -ECPs as dimension-reducing solutions of the EN model. Coupled n-ECP are displayed in visual space showing simultaneously the distortion of the RM and the OPM (σ/Λ = 0.3 (η = 0.028), η _r = η, circular stimulus ensemble). The distorted grid represents a the cortical square array of cells. Each grid intersection is at the receptive field center of the corresponding cell. Preferred stimulus orientations are color-coded as in Figure 2a. As in Figure 4, n and i enumerate the number of nonzero wave vectors and nonequivalent configurations of wave vectors with the same n, respectively. The coupled 1-ECP is a pinwheel-free stripe pattern without retinotopic distortion. Only the most anisotropic and the most isotropic coupled n-ECPs are shown for each n. Note that for all ECPs, high gradients within the orientation mapping coincide with low gradients of the retinotopic mapping and vice versa. Retinotopic distortions are displayed on a fivefold magnified scale for visualization purposes.

Full size image

The impact of retinotopic distortions

According to our analysis, at criticality, the nontrivial stable fixed points of the EN dynamics are determined by the continuity parameter η ∈ (0, 1) for the OPM or, equivalently, the ratio \sigma ∕\Lambda =\frac{1}{2\pi }\sqrt{log\left(1∕\eta \right)} and the continuity parameter η _r for the mapping of visual space. We first tested for the stability of pinwheel-free orientation stripe (OS) solutions and rPWC solutions of Equation (15), with coupling matrices g _ij and f _ij as obtained from the nonlinearities in Equation (29). The angle which minimizes the energy U_rPWC (Equation (20)) is not affected by the coupling between retinotopic and OPM and is thus again α = π/4. By numerical evaluation of the criteria for intrinsic and extrinsic stability, we found both, OSs and sPWCs, to be intrinsically and extrinsically stable for all σ/Λ and η _r .

Next, we tested for the stability of coupled n-ECP solutions for 2 ≤ n ≤ 20. We found all coupled n-ECP configurations with n ≥ 2 to be intrinsically unstable for all σ/Λ and η _r . Evaluating the energy assigned to sPWCs and OSs, we identified two different regimes: (i) for shorter interaction range σ/Λ the sPWC is the minimal energy state and (ii) for larger interaction range σ/Λ the optimum is an OS pattern as indicated by the phase diagram in Figure 10a. The retinotopic continuity parameter has little influence on the energy of the two fixed points. The phase border separating stripes from rhombs runs almost parallel to the η _r -axis. We numerically confirmed these analytical predictions by extensive simulations of Equation (3, 4) (see 'Methods' section for details). Figure 10c shows snapshots of a representative simulation with small interaction range (r = 0.1, σ/Λ = 0.1 (η = 0.67), η _r = η). After the initial symmetry breaking phase, the OPM layout rapidly converges toward a crystalline array of pinwheels, the predicted optimum in this parameter regime (Figure 10c). Retinotopic deviations are barely visible. Figure 10b displays pinwheel density time courses for four such simulations. Note that in one simulation, the pinwheel density drops to almost zero. In this simulation, the OP pattern converges to a stripe-like layout. This is in line with the finding of bistability of rhombs and stripes in all parameter regimes. Although the sPWC represents the global minimum in the simulated parameter regime, OSs are also a stable fixed point and, depending on the initial conditions, may arise as the final state of a fraction of the simulations. In the two simulations with pinwheel densities around 3.4, patterns at later simulation stages consist of different domains of rhombic pinwheel lattices with α < π/2.

Phase diagram of the EN model with variable retinotopy for a circular stimulus ensemble [57, 64–66, 91]. (a) Regions of the η _r -σ/Λ-plane in which n-ECPs or rPWCs have minimal energy. (b) Pinwheel density time courses for four different simulations of Equations (3, 4) with r = 0.1, σ/Λ = 0.13 (η = 0.51), η _r = η (grey traces, individual realizations; red trace, mean value; black trace, realization shown in c). (c) OPMs (upper row), their power spectra (middle row), and RMs (lower row) obtained in a simulation of Equations (3, 4); parameters as in b. (d) Pinwheel density time courses for four different simulations of Equations (3, 4) with r = 0.1, σ/Λ = 0.3 (η = 0.03), η _r = η (grey traces, individual realizations; red trace, mean value; black trace, realization shown in e). (e) OPMs (upper row), their power spectra (middle row), and RMs (lower row) in a simulation of Equations (3, 4); parameters as in d.

Full size image

Figure 10d,e shows the corresponding analysis with parameters for larger interaction range r = 0.1, σ/Λ = 0.15 (η = 0.41), η _r = η. Here after initial pinwheel creation, pinwheels typically annihilate pairwisely and the OPM converges to an essentially pinwheel-free stripe pattern, the predicted optimal solution in this parameter regime (Figure 10e). Retinotopic deviations are slightly larger. The behavior of the EN model for the joint optimization of RM and OPM thus appears very similar compared to the fixed retinotopy case. Perhaps surprisingly, the coupling of both feature maps has little effect on the stability properties of the fixed points and the resulting optimal solutions.

As in the previous case, the structure of the phase diagram in Figure 10a appears somewhat counterintuitive. A high coverage and pinwheel-rich solution is the optimum in a regime with large OPM continuity parameter where discontinuities in the OPM such as pinwheels should be strongly penalized. A pinwheel-free solution with low coverage and high continuity is the optimum in a regime with small continuity parameter. As explained above, a large OPM continuity parameter at pattern formation threshold implies a small interaction range σ/Λ (see Equation (24)). In such a regime, the gain in coverage by representing many orientation stimuli in a small area spanning the typical interaction range, e.g., with a pinwheel, is very high. Apparently this gain in coverage by a regular positioning of pinwheels outweighs the accompanied loss in continuity for very large OPM continuity parameters. This counterintuitive interplay between coverage and continuity thus seems to be almost independent of the choice of retinotopic continuity parameters.

The circular orientation stimulus ensemble contains only stimuli with a fixed and finite 'orientation energy' or elongation |s _z |. This raises the question of whether the simple nature of the circular stimulus ensemble might restrain the dynamics of the EN model. The EN dynamics are expected to depend on the characteristics of the activity patterns evoked by the stimuli and these will be more diverse and complex with ensembles containing a greater diversity of stimuli. Therefore, we repeated the above analysis of the EN model for a richer stimulus ensemble where orientation stimuli are uniformly distributed on the disk {s _z , |s _z | ≤ 2}, a choice adopted by a subset of previous studies, e.g., [19, 25, 81]. In particular, this ensemble contains unoriented stimuli with |s _z | = 0. Intuitively, the presence of these unoriented stimuli might be expected to change the role of pinwheels in the optimal OPM layout. Pinwheels' population activity is untuned for orientation. Pinwheel centers may therefore acquire a key role for the representation of unoriented stimuli. Nevertheless, we found the behavior of the EN model when considered with this richer stimulus ensemble to be virtually indistinguishable from the circular stimulus ensemble. Details of the derivations, phase diagrams and numerically obtained solutions are given in Appendix 1.

Are there stimulus ensembles for which realistic, aperiodic maps are optimal?

So far, we have presented a comprehensive analysis of optimal dimension-reducing mappings of the EN model for two widely used orientation stimulus distributions (previous sections and Appendix 1). In both cases, optima were either regular crystalline pinwheel lattices or pinwheel-free orientation stripes. These results might indicate that the EN model for the joint optimization of OPM and RM is per se incapable of reproducing the structure of OPMs as found in the visual cortex. Drawing such a conclusion is suggested in view of the apparent insensitivity of the model's optima to the choice of stimulus ensemble. The two stimulus ensembles considered so far however do not exhaust the infinite space of stimulus distributions that are admissible in principle. From the viewpoint of 'biological plausibility' it is certainly not obvious that one should strive to examine stimulus distributions very different from these, as long as the guiding hypothesis is that the functional architecture of the primary visual cortex optimizes the joint representation of the classical elementary stimulus features. If, however, stimulus ensembles were to exist, for which optimal EN mappings truly resemble the biological architecture, their characteristics may reveal essential ingredients of alternative optimization models for visual cortical architecture.

Adopting this perspective raises the technical question of whether an unbiased search of the infinite space of stimulus ensembles only constrained by the model's symmetries (Equations (5) to (9)) is possible. To answer this question, we examined whether the amplitude equations (15) can be obtained for an arbitrary orientation stimulus distribution. Fortunately, we found that the coefficients of the amplitude equations are completely determined by the finite set of moments of order less than 5 of the distributions. The approach developed so far can thus be used to comprehensively examine the nature of EN optima resulting for any stimulus distribution with finite fourth-order moment. While such a study does not completely exhaust the infinite space of all eligible distributions, it appears to only exclude ensembles with really exceptional properties. These are probability distributions with diverging fourth moment, i.e., ensembles that exhibit a heavy tail of essentially 'infinite' orientation energy stimuli.

Since the coupling between OPM and RM did not have a large impact in the case of the two classical stimulus ensembles, we start the search through the space of orientation stimulus ensembles by considering the EN model with fixed retinotopy r(x) = 0. The coefficients a _i for the nonlinear operators N_i^3\left[z,z,\bar{z}\right] in Equation (14) for arbitrary stimulus ensembles are given by

The corresponding angle-dependent interaction functions are given by (see 'Methods' section)

Again, without loss of generality, we set 〈|s_z|²〉 = 2. At criticality, both functions are parameterized by the continuity parameter η ∈ (0, 1) for the OPM or, equivalently, the interaction range \sigma ∕\Lambda =\frac{1}{2\pi }\sqrt{log\left(1∕\eta \right)} and the fourth moment 〈|s_z|⁴〉 of the orientation stimulus ensemble. The fourth moment, is a measure of the peakedness of a stimulus distribution. High values generally indicate a strongly peaked distribution with a large fraction of nonoriented stimuli (|s _z |⁴ ≈ 0), together with a large fraction of high orientation energy stimuli (|s _z |⁴ large).

The dependence of g(α) on the fourth moment of the orientation stimulus distribution and f(α) suggests that different stimulus distributions may indeed lead to different optimal dimension-reducing mappings.

The circular stimulus ensemble possesses the minimal possible fourth moment, with 〈|s_z|⁴〉 = (〈|s_z|²〉)² = 4. The fourth moment of the uniform stimulus ensemble is 〈|s_z|⁴〉 = 16/3. The angle-dependent interaction functions for both ensembles (Equation (25), Figure 22 in Appendix 1) are recovered, when inserting these values into Equation (33).

To simplify notation in the following, we define

as the parameter characterizing an orientation stimulus distribution. This parameter ranges from zero for the circular stimulus ensemble to infinity for ensembles with diverging fourth moments. Figure 11 displays the angle-dependent interaction functions for different values of σ/Λ and s₄. In all parameter regimes, g(α) and f(α) are larger than zero. The amplitude dynamics are therefore guaranteed to converge to a stable stationary fixed point and the bifurcation from the nonselective fixed point in the EN model is predicted to be supercritical in general.

Angle-dependent interaction functions for the EN model with fixed retinotopy for different fourth-moment values of the orientation stimulus distribution and effective interaction-widths. (a, b) g(α) and f(α) for s₄ = 8 and σ/Λ = 0.1 (a) and σ/Λ = 0.3 (b). (c, d) g(α) and f(α) for s₄ = 20 and σ/Λ = 0.1 (c) and σ/Λ = 0.3 (d). (e, f) g(α) and f(α) for s₄ = 100 and σ/Λ = 0.1 (e) and σ/Λ = 0.3 (f).

Full size image