Coverage, continuity, and visual cortical architecture

Abstract

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.

Introduction

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

Figure 1

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

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.

Results and discussion

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

$$R\left(x\right)=X+r\left(x\right)$$

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

$$\vartheta \left(x\right)=\frac{1}{2}arg\left(z\left(x\right)\right).$$

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

$$\mathcal{F}={\sigma }^2\mathcal{C}+\mathcal{R}$$

(1)

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

Figure 2

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.

$$\begin{array}{ll}\hfill \mathcal{C}\left[z,\mathbf{r}\right]& =-{〈ln\int d^{2}y{e}^{-\left(\mid s_z-z{(\mathbf{y})\mid }^2+\mid {\mathbf{s}}_r-\mathbf{X}-\mathbf{r}\left(\mathbf{y}\right){\mid }^2\right)∕2{\sigma }^2}〉}_{\mathbf{S}}\\ \hfill \mathcal{R}\left[z,\mathbf{r}\right]& =\frac{1}{2}\int d^{2}y\eta \left|\right|\nabla z\left(\mathbf{y}\right)|{|}^2+{\eta }_r\sum _{j=1}^2||\nabla r_j\left(\mathbf{y}\right)|{|}^2,\end{array}$$

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

$$\begin{array}{c}{\partial }_{t}z\left(\mathbf{x}\right)=-2\frac{\delta \mathcal{F}\left[z,\mathbf{r}\right]}{\delta \bar{z}\left(\mathbf{x}\right)}\equiv F^z\left[z,\mathbf{r}\right]\left(\mathbf{x}\right)\\ {\partial }_t\mathbf{r}\left(\mathbf{x}\right)=-\frac{\delta \mathcal{F}\left[z,\mathbf{r}\right]}{\delta \mathbf{r}\left(\mathbf{x}\right)}\equiv {\mathbf{F}}^r\left[z,\mathbf{r}\right]\left(\mathbf{x}\right)\end{array}$$

(2)

called the EN dynamics in the following. These dynamics read

$${\partial }_{t}z\left(\mathbf{x}\right)={〈\left[s_z-z\left(\mathbf{x}\right)\right]e\left(\mathbf{x},\mathbf{S},z,\mathbf{r}\right)〉}_{\mathbf{S}}+\eta \Delta z\left(\mathbf{x}\right)$$

(3)

$${\partial }_t\mathbf{r}\left(\mathbf{x}\right)={〈\left[{\mathbf{s}}_r-\mathbf{X}-\mathbf{r}\left(\mathbf{x}\right)\right]e\left(\mathbf{x},\mathbf{S},z,\mathbf{r}\right)〉}_{\mathbf{S}}+{\eta }_r\Delta \mathbf{r}\left(\mathbf{x}\right),$$

(4)

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

$$e\left(\mathbf{x},\dots \right)=\frac{e^{-\left(\mid {\mathbf{s}}_r-\mathbf{X}-\mathbf{r}{(\mathbf{x})\mid }^2\right)∕2{\sigma }^2}{e}^{-\left(\mid s_z-z{(\mathbf{x})\mid }^2\right)∕2{\sigma }^2}}{\int d^{2}y{e}^{-\left(\mid {\mathbf{s}}_r-\mathbf{X}-\mathbf{r}{(\mathbf{y})\mid }^2\right)∕2{\sigma }^2}{e}^{-\left(\mid s_z-z{(\mathbf{y})\mid }^2\right)∕2{\sigma }^2}}.$$

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

$$\begin{array}{c}{\widehat{T}}_{\mathbf{y}}z\left(\mathbf{x}\right)=z\left(\mathbf{x}+\mathbf{y}\right)\\ {\widehat{T}}_{\mathbf{y}}\mathbf{r}\left(\mathbf{x}\right)=\mathbf{r}\left(\mathbf{x}+\mathbf{y}\right),\end{array}$$

rotations

$$\begin{array}{c}{\widehat{R}}_{\beta }z\left(\mathbf{x}\right)=e^{2i\beta }z\left({\mathrm{\Omega }}_{-\beta }\mathbf{x}\right)\\ {\widehat{R}}_{\beta }\mathbf{r}\left(\mathbf{x}\right)={\mathrm{\Omega }}_{\beta }\mathbf{r}\left({\mathrm{\Omega }}_{-\beta }\mathbf{x}\right)\end{array}$$

with 2×2 rotation matrix

$${\mathrm{\Omega }}_{\beta }=\left(\begin{array}{cc}\hfill cos\beta \hfill & \hfill -sin\beta \hfill \\ \hfill sin\beta \hfill & \hfill cos\beta \hfill \end{array}\right),$$

and reflections

$$\begin{array}{c}\widehat{P}z\left(\mathbf{x}\right)=\bar{z}\left(\Psi \mathbf{x}\right)\\ \widehat{P}\mathbf{r}\left(\mathbf{x}\right)=\Psi \mathbf{r}\left(\Psi \mathbf{x}\right),\end{array}$$

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

$${\widehat{T}}_{\mathbf{y}}{F}^z\left[z,\mathbf{r}\right]=F^z\left[{\widehat{T}}_{y}z,{\widehat{T}}_{\mathbf{y}}\mathbf{r}\right]$$

(5)

$${\widehat{R}}_{\beta }{F}^z\left[z,\mathbf{r}\right]=F^z\left[{\widehat{R}}_{\beta }z,{\widehat{R}}_{\beta }\mathbf{r}\right]$$

(6)

$$\widehat{P}{F}^z\left[z,\mathbf{r}\right]=F^z\left[\widehat{P}z,\widehat{P}\mathbf{r}\right],$$

(7)

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

$$e^{i\varphi }{F}^z\left[z,\mathbf{r}\right]=F^z\left[e^{i\varphi }z,\mathbf{r}\right]$$

(8)

$${\mathbf{F}}^r\left[z,\mathbf{r}\right]={\mathbf{F}}^r\left[e^{i\varphi }z,\mathbf{r}\right].$$

(9)

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

$${\partial }_t\mathbf{r}\left(\mathbf{x}\right)\simeq {\mathbf{L}}_r\left[\mathbf{r}\right]=\frac{1}{16\pi {\sigma }^4}\int d^{2}y{e}^{-\frac{{\left(\mathbf{x}-\mathbf{y}\right)}^2}{4{\sigma }^2}}\widehat{\mathbf{A}}\mathbf{r}\left(\mathbf{y}\right)+{\eta }_r△\mathbf{r}\left(\mathbf{x}\right)$$

(10)

$${\partial }_{t}z\left(\mathbf{x}\right)\simeq L_z\left[z\right]=\left(\frac{1}{{\sigma }^2}-1\right)z\left(\mathbf{x}\right)+\eta \Delta z\left(\mathbf{x}\right)-\frac{1}{4\pi {\sigma }^4}\int d^{2}y{e}^{-\frac{{\left(\mathbf{x}-\mathbf{y}\right)}^2}{4{\sigma }^2}}z\left(\mathbf{y}\right),$$

(11)

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

$${\partial }_t{\tilde{r}}_i\left(\mathbf{k}\right)=-\sum _{j=1}^2\left({\sigma }^2{e}^{-k^2{\sigma }^2}{k}_i{k}_j+{\eta }_r{k}^2{\delta }_{ij}\right){\tilde{r}}_j\left(\mathbf{k}\right),$$

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

$${\lambda }_L^r=-k^2\left({\eta }_r+e^{-{\sigma }^2{k}^2}{\sigma }^2\right),{\lambda }_T^r=-{\eta }_r{k}^2$$

with corresponding eigenfunctions (in real space)

$$\begin{array}{c}{\mathbf{r}}_L\left(\mathbf{x}\right)={\mathbf{k}}_{\varphi }{e}^{i{\mathbf{k}}_{\varphi }\mathbf{x}}+\mathsf{\text{c}}.\mathsf{\text{c}}.\\ {\mathbf{r}}_T\left(\mathbf{x}\right)={\mathbf{k}}_{\varphi +\pi ∕2}{e}^{i{\mathbf{k}}_{\varphi }\mathbf{x}}+\mathsf{\text{c}}.\mathsf{\text{c}}.,\end{array}$$

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^r$are 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.

Figure 3

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

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

$${\lambda }^z\left(k\right)=-1+\frac{1}{{\sigma }^2}\left(1-e^{-k^2{\sigma }^2}\right)-\eta k^2$$

(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

$$\sigma >{\sigma }^{*}\left(\eta \right)=\sqrt{1+\eta ln\eta -\eta },$$

(12)

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

Figure 4

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.

$$T_{\mathbf{y}}\left[F^z\left[e^{i\mathbf{k}x},0\right]\right]=F^z\left[T_y\left[e^{i\mathbf{k}x}\right],T_y\left[0\right]\right]=F^z\left[e^{i\mathbf{k}y}{e}^{i\mathbf{k}x},0\right]=e^{i\mathbf{k}y}{F}^z\left[e^{i\mathbf{k}x},0\right]$$

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:

$$T_{\mathbf{y}}\left[{\mathbf{F}}^r\left[e^{i\mathbf{k}x},0\right]\right]={\mathbf{F}}^r\left[T_y\left[e^{i\mathbf{k}x}\right],T_y\left[0\right]\right]={\mathbf{F}}^r\left[e^{i\mathbf{k}y}{e}^{i\mathbf{k}x},0\right]={\mathbf{F}}^r\left[e^{i\mathbf{k}x},0\right]$$

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

$${\partial }_{t}z\left(\mathbf{x}\right)\approx L_z\left[z\right]+N_3^z\left[z,z,\bar{z}\right],$$

(13)

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

$$N_3^z\left[\sum _j{\alpha }_j{z}_j,\sum _k{\beta }_k{z}_k,\sum _l{\gamma }_l{\bar{z}}_l\right]=\sum _{j,k,l}{\alpha }_j{\beta }_k{\gamma }_l{N}_3^z\left[z_j,z_k,{\bar{z}}_l\right].$$

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

$$N_3^z\left[z,z,\bar{z}\right]=\sum _{j=1}^{11}{a}_j{N}_3^j\left[z,z,\bar{z}\right].$$

(14)

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

$$z\left(\mathbf{x}\right)=\sum _j{A}_j\left(t\right)e^{i{\mathbf{k}}_j\mathbf{x}}$$

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]

$${Ȧ}_i=r{A}_i-A_i\sum _{j=1}^N{g}_{ij}\mid A_j{\mid }^2-{Ā}_{i^{-}}\sum _{j=1}^N{f}_{ij}{A}_j{A}_{j^{-}},$$

(15)

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

$${Ȧ}_j\left(t\right)=-\frac{\partial U_A}{\partial {Ā}_j}$$

from an energy

$$U_A=-r\sum _{i=1}^N\mid A_i{\mid }^2+\frac{1}{2}\sum _{i,j=1}^N{g}_{ij}\mid A_i{\mid }^2\mid A_j{\mid }^2+\frac{1}{2}\sum _{i,j=1}^N{f}_{ij}{Ā}_i{Ā}_{i^{-}}{A}_j{A}_{j^{-}}.$$

(16)

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

$$\mid A_0{\mid }^2=\frac{r}{g_{ii}}$$

(17)

and its energy

$$U_{\mathsf{\text{OS}}}=-\frac{r}{2g_{ii}}.$$

(18)

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

$$z\left(\mathbf{x}\right)=A_1{e}^{i{\mathbf{k}}_1\mathbf{x}}+A_2{e}^{i{\mathbf{k}}_2\mathbf{x}}+A_3{e}^{-i{\mathbf{k}}_1\mathbf{x}}+A_4{e}^{-i{\mathbf{k}}_2\mathbf{x}}$$

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

$${\mathcal{A}}^2=\frac{r}{g_{00}+g_{0\pi }+g_{0\alpha }+g_{0\pi -\alpha }-2f_{0\alpha }}.$$

(19)

The phase relations of the four amplitudes are given by

$$\begin{array}{c}{\varphi }_1+{\varphi }_3={\varphi }_0\\ {\varphi }_2+{\varphi }_4={\varphi }_0+\pi .\end{array}$$

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

$$U_{\mathsf{\text{rPWC}}}=-\frac{2r}{g_{00}+g_{0\pi }+g_{0\alpha }+g_{0\pi -\alpha }-2f_{0\alpha }}.$$

(20)

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

$$s\left(\alpha \right)=g_{0\alpha }+g_{0\pi -\alpha }-2f_{0\alpha }$$

(21)

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

$$z\left(\mathbf{x}\right)=\sum _{j=1}^n{A}_j{e}^{i{l}_j{\mathbf{k}}_j\mathbf{x}},$$

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

$${Ȧ}_i=r{A}_i-A_i\sum _{j=1}^n{g}_{ij}\mid A_j{\mid }^2,$$

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

$$\mid A_i{\mid }^2=r\sum _{j=1}^n{\left({\mathbf{g}}^{-1}\right)}_{ij}.$$

(22)

The energy of an n-ECP is given by

$$U_{\mathsf{\text{ECP}}}=-\frac{r}{2}\sum _{i,j}{\left({\mathbf{g}}^{-1}\right)}_{ij}.$$

(23)

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 Λ

$$\sigma ∕\Lambda =\frac{1}{2\pi }\sqrt{log\left(1∕\eta \right)}$$

(24)

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

$$\begin{array}{ccc}\hfill a_1=\frac{1}{4{\sigma }^6}-\frac{1}{{\sigma }^4}+\frac{1}{2{\sigma }^2}\hfill & \hfill a_2=\frac{1}{4\pi {\sigma }^6}-\frac{1}{8\pi {\sigma }^8}\hfill & \hfill a_3=-\frac{1}{16\pi {\sigma }^8}+\frac{1}{8\pi {\sigma }^6}\hfill \\ \hfill a_4=-\frac{1}{8\pi {\sigma }^8}+\frac{1}{4\pi {\sigma }^6}-\frac{1}{8\pi {\sigma }^4}\hfill & \hfill a_5=-\frac{1}{16\pi {\sigma }^8}\hfill & \hfill a_6=\frac{1}{8\pi {\sigma }^6}-\frac{1}{16\pi {\sigma }^8}\hfill \\ \hfill a_7=\frac{1}{12{\pi }^2{\sigma }^{10}}-\frac{1}{12{\pi }^2{\sigma }^8}\hfill & \hfill a_8=\frac{1}{24{\pi }^2{\sigma }^{10}}\hfill & \hfill a_9=-\frac{3}{64{\pi }^3{\sigma }^{12}}\hfill \\ \hfill a_{10}=\frac{1}{12{\pi }^2{\sigma }^{10}}-\frac{1}{12{\pi }^2{\sigma }^8}\hfill & \hfill a_{11}=\frac{1}{24{\pi }^2{\sigma }^{10}}.\hfill \end{array}$$

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

$$\begin{array}{c}g\left(\alpha \right)=\frac{1}{{\sigma }^4}\left(1-2e^{-k_c^2{\sigma }^2}-e^{2k_c^2{\sigma }^2\left(cos\alpha -1\right)}\left(1-2e^{-k_c^2{\sigma }^{2}cos\alpha }\right)\right)\\ +\frac{1}{2{\sigma }^2}\left(e^{2k_c^2{\sigma }^2\left(cos\alpha -1\right)}-1\right)+\frac{8}{{\sigma }^6}{e}^{-2k_c^2{\sigma }^2}{sinh}^4\left(1∕2k_c^2{\sigma }^{2}cos\alpha \right)\\ f\left(\alpha \right)=\frac{1}{{\sigma }^4}\left(1-e^{-2k_c^2{\sigma }^2}\left(cosh\left(2k_c^2{\sigma }^{2}cos\alpha \right)+2cosh\left(k_c^2{\sigma }^{2}cos\alpha \right)\right)+2e^{-k_c^2{\sigma }^2}\right)\\ +\frac{1}{2{\sigma }^2}\left(e^{-2k_c^2{\sigma }^2}cosh\left(2k_c^2{\sigma }^{2}cos\alpha \right)-1\right)+\frac{4}{{\sigma }^6}{e}^{-2k_c^2{\sigma }^2}{sinh}^4\left(1∕2k_c^2{\sigma }^{2}cos\alpha \right).\end{array}$$

(25)

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.

Figure 5

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

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

Figure 6

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

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

Figure 7

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

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

$${\partial }_{t}z\left(\mathbf{x}\right)=L_z\left[z\right]+Q^z\left[\mathbf{r},z\right]+N_3^z\left[z,z,\bar{z}\right]+\cdots$$

(26)

$${\partial }_t\mathbf{r}\left(\mathbf{x}\right)={\mathbf{L}}_r\left[\mathbf{r}\right]+{\mathbf{Q}}^r\left[z,\bar{z}\right]+\cdots .$$

(27)

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

$$\mathbf{r}\left(\mathbf{x}\right)=-{\mathbf{L}}_r^{-1}\left[{\mathbf{Q}}^r\left[z,\bar{z}\right]\right].$$

(28)

We remark that as ${\lambda }_{T∕L}^r\left(k\right)<0$ 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:

$$\begin{array}{c}{\partial }_{t}z\left(\mathbf{x}\right)\approx L_z\left[z\right]\underset{N_3^r\left[z,z,\bar{z}\right]}{\underbrace{-Q^z\left[{\mathbf{L}}_r^{-1}\left[{\mathbf{Q}}^r\left[z,\bar{z}\right]\right],z\right]}}+N_3^z\left[z,z,\bar{z}\right]\\ =L_z\left[z\right]+N_3^r\left[z,z,\bar{z}\right]+N_3^z\left[z,z,\bar{z}\right].\end{array}$$

(29)

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

$$N_3^r\left[z,z,\bar{z}\right]=\sum _{j=1}^{12}{a}_r^j{N}_r^j\left[z,z,\bar{z}\right].$$

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

$$\begin{array}{c}{g}_r\left(\alpha \right)=\frac{{\left(\left(1-{\sigma }^2-2e^{-k_c^2{\sigma }^2}\right)e^{2k_c^2{\sigma }^2\left(cos\alpha -1\right)}+e^{-k_c^2{\sigma }^2}\right)}^2}{2{\sigma }^4\left({\eta }_r+{\sigma }^2{e}^{-2k_c^2{\sigma }^2\left(cos\alpha -1\right)}\right)}\\ f_r\left(\alpha \right)=\frac{1}{2}\left(g_r\left(\alpha \right)+g_r\left(\alpha +\pi \right)\right),\end{array}$$

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 .

Figure 8

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

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

$$z\left(\mathbf{x}\right)=\sum _j^N{A}_j{e}^{i{\mathbf{k}}_j\mathbf{x}},$$

(30)

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

$$\begin{array}{c}\mathbf{r}\left(\mathbf{x}\right)=-\sum _{k=1,j<k}^n\frac{{\Delta }_{jk}}{{\lambda }_L^r\left(\mid {\Delta }_{jk}\mid \right)}\left(\frac{1}{{\sigma }^2}{\left(e^{-k_c^2{\sigma }^2∕2}-e^{-{\Delta }_{jk}^2{\sigma }^2∕2}\right)}^2-e^{-{\Delta }_{jk}^2{\sigma }^2}\right)\\ *\left(\Im \left(A_j{Ā}_k\right)cos\left({\Delta }_{jk}\mathbf{x}\right)+\Re \left(A_j{Ā}_k\right)sin\left({\Delta }_{jk}\mathbf{x}\right)\right),\end{array}$$

(31)

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)

$$z_{\mathsf{\text{sPWC}}}\left(\mathbf{x}\right)\propto sin\left(k_c{x}_1\right)+isin\left(k_c{x}_2\right),$$

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

$${\mathbf{r}}_{s\mathsf{\text{PWC}}}\left(\mathbf{x}\right)\propto \frac{e^{-k_c^2{\sigma }^2}}{{\sigma }^2{\lambda }_2^r\left(2k_c\right)}\left(\begin{array}{c}\hfill k_{c}sin\left(2k_c{x}_1\right)\hfill \\ \hfill k_{c}sin\left(2k_c{x}_2\right)\hfill \end{array}\right).$$

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.

Figure 9

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.

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.

Figure 10

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.

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

$$\begin{array}{ccc}\hfill a_1=\frac{〈\mid s_z{\mid }^4〉}{16{\sigma }^6}-\frac{〈\mid s_z{\mid }^2〉}{2{\sigma }^4}+\frac{1}{2{\sigma }^2}\hfill & \hfill a_2=\frac{〈\mid s_z{\mid }^2〉}{8\pi {\sigma }^6}-\frac{〈\mid s_z{\mid }^4〉}{32\pi {\sigma }^8}\hfill & \hfill a_3=-\frac{〈\mid s_z{\mid }^4〉}{64\pi {\sigma }^8}+\frac{〈\mid s_z{\mid }^2〉}{16\pi {\sigma }^6}\hfill \\ \hfill a_4=-\frac{〈\mid s_z{\mid }^4〉}{32\pi {\sigma }^8}+\frac{〈\mid s_z{\mid }^2〉}{8\pi {\sigma }^6}-\frac{1}{8\pi {\sigma }^4}\hfill & \hfill a_5=-\frac{〈\mid s_z{\mid }^4〉}{64\pi {\sigma }^8}\hfill & \hfill a_6=\frac{〈\mid s_z{\mid }^2〉}{16\pi {\sigma }^6}-\frac{〈\mid s_z{\mid }^4〉}{64\pi {\sigma }^8}\hfill \\ \hfill a_7=\frac{〈\mid s_z{\mid }^4〉}{48{\pi }^2{\sigma }^{10}}-\frac{〈\mid s_z{\mid }^2〉}{24{\pi }^2{\sigma }^8}\hfill & \hfill a_8=\frac{〈\mid s_z{\mid }^4〉}{96{\pi }^2{\sigma }^{10}}\hfill & \hfill a_9=-\frac{3〈\mid s_z{\mid }^4〉}{256{\pi }^3{\sigma }^{12}}\hfill \\ \hfill a_{10}=\frac{〈\mid s_z{\mid }^4〉}{48{\pi }^2{\sigma }^{10}}-\frac{〈\mid s_z{\mid }^2〉}{24{\pi }^2{\sigma }^8}\hfill & \hfill a_{11}=\frac{〈\mid s_z{\mid }^4〉}{96{\pi }^2{\sigma }^{10}}.\hfill \end{array}$$

(32)

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

$$\begin{array}{c}g\left(\alpha \right)=\frac{〈\mid s_z{\mid }^2〉}{2{\sigma }^4}\left(1-2e^{-k_c^2{\sigma }^2}-e^{2k_c^2{\sigma }^2\left(cos\alpha -1\right)}\left(1-2e^{-k_c^2{\sigma }^{2}cos\alpha }\right)\right)\\ +\frac{1}{2{\sigma }^2}\left(e^{2k_c^2{\sigma }^2\left(cos\alpha -1\right)}-1\right)+\frac{2〈\mid s_z{\mid }^4〉}{{\sigma }^6}{e}^{-2k_c^2{\sigma }^2}{sinh}^4\left(1∕2k_c^2{\sigma }^{2}cos\alpha \right)\\ f\left(\alpha \right)=\frac{〈\mid s_z{\mid }^2〉}{2{\sigma }^4}\left(1-e^{-2k_c^2{\sigma }^2}\left(cosh\left(2k_c^2{\sigma }^{2}cos\alpha \right)+2cosh\left(k_c^2{\sigma }^{2}cos\alpha \right)\right)+2e^{-k_c^2{\sigma }^2}\right)\\ +\frac{1}{2{\sigma }^2}\left(e^{-2k_c^2{\sigma }^2}cosh\left(2k_c^2{\sigma }^{2}cos\alpha \right)-1\right)+\frac{〈\mid s_z{\mid }^4〉}{{\sigma }^6}{e}^{-2k_c^2{\sigma }^2}{sinh}^4\left(1∕2k_c^2{\sigma }^{2}cos\alpha \right).\end{array}$$

(33)

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

$$s_4=〈\mid s_z{\mid }^4〉-{〈\mid s_z{\mid }^2〉}^2=〈\mid s_z{\mid }^4〉-4$$

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.

Figure 11

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

By evaluating the energy assigned to the rPWC and n-ECPs, we investigated the structure of the two-dimensional phase space of the EN model with an arbitrary orientation stimulus distribution. First, it is not difficult to show that the angle α which minimizes the energy U_rPWC (Equation (20)) of an rPWC is α = π/4 for all σ/Λ and s₄. Hence, a square lattice of pinwheels (sPWC) is in all parameter regimes energetically favored over any other rhombic lattice configuration of pinwheels. Figure 12 displays the phase diagram of the EN model with an arbitrary orientation stimulus distribution. For orientation stimulus distributions with small fourth moments, optimal mappings consist of either parallel pinwheel-free stripes or quadratic pinwheel crystals. These distributions include the circular and the uniform stimulus ensembles with s₄ = 0 and s₄ = 4/3. Above a certain value of the fourth moment around s₄ = 6, n-ECPs with n >2 become optimal mappings. For a short interaction range σ/Λ, hexagonal pinwheel crystals dominate the phase diagram in a large region of parameter space. With increasing interaction range, we observe a sequence of phase transitions by which higher n-ECPs become optimal. For n >3, these optima are spatially aperiodic. In all parameter regimes, we found that the n-ECP with the most anisotropic mode configuration (Figure 4c, left column) is the energetically favored state for n >3. Pinwheel densities of these planforms are indicated in Figure 12 and are typically smaller than 2.0. We note that this is well below experimentally observed pinwheel density values [38]. Optimal mappings of orientation preference for finite fourth moment in the EN model are thus either orientation stripes, periodic arrays of pinwheels (hexagonal, square), or aperiodic pinwheel arrangements with low pinwheel density.

Figure 12

Stripe-like, crystalline, and quasi-crystalline cortical representations as optimal solutions to the mapping of orientation preference with fixed uniform retinotopy in the EN model. The graph shows the regions of the s₄-σ/Λ-plane in which n-ECPs or sPWCs have minimal energy. For n ≥ 3, pinwheel densities of the energetically favored n-ECP configuration are indicated.

We numerically tested these analytical predictions by simulations of Equation (3) (r(x) = 0) with two additional stimulus ensembles with s₄ = 6 and s₄ = 8 (see 'Methods' section). Figure 13a shows snapshots of a simulation with (r = 0.1, σ/Λ = 0.2 (η = 0.2)) and s₄ = 6 (see also Additional file 3). After the initial phase of pattern emergence, the OPM layout converges toward an arrangement of fractured stripes which resembles the 2-ECP state (Figure 13a, most right), the optimum predicted in this regime. In the power spectra, two distinct peaks of the active modes are clearly visible in the final stages of the simulation (Figure 13a, lower row). The 2-ECP state is exotic in the sense that it is the only n-ECP containing line defects and thus the pinwheel density is not a well-defined quantity. This explains the pronounced numerical variability in the measured pinwheel densities in simulations during the convergence toward a 2-ECP state (Figure 13b).

Figure 13

Approaching crystalline n -ECP optima in the EN model with fixed retinotopy. (a) OPMs (upper row) and their power spectra (lower row) in a simulation of Equation (3) with r(x) = 0, r = 0.1, σ/Λ = 0.2 and s₄ = 6 (see also Additional file 3). The predicted optimum is the 2-ECP (black frame). (b) 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) OPMs (upper row) and their power spectra (lower row) in a simulation of Equation (3) with r(x) = 0, r = 0.1, σ/Λ = 0.3 and s₄ = 8 (see also Additional file 4). The predicted optimum is the anisotropic 3-ECP (black frame). (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).

Figure 13c shows snapshots of a simulation with (r = 0.1, σ/Λ = 0.2 (η = 0.2)) and s₄ = 8, Gaussian stimulus ensemble) (see also Additional file 4). After the initial phase of pattern emergence, the OPM layout converges toward a regular hexagonal arrangement of pinwheels which resembles the anisotropic 3-ECP (Figure 13c, far right), the optimum predicted in this regime. In the power spectra, three distinct peaks forming an angle of 60 degrees are clearly visible in the later stages of the simulation (Figure 13c, lower row). Pinwheel densities in the simulations consistently approach the theoretically predicted value of 2 cos(π/6) ≃ 1.73 (Figure 13d).

Permutation symmetric limit

In the previous section, we uncovered a parameter regime for the EN model in which optimal solutions are spatially aperiodic. This can be viewed as a first step toward realistic optimal solutions. In the identified regime, however, among the family of n-ECPs only those with pinwheel densities well below experimentally observed values [38] are energetically favored (see Figure 12). In this respect, the repertoire of aperiodic optima of the EN model differs from previously considered abstract variational models for OPM development [35, 36, 38, 39]. In these models, an energetic degeneracy of aperiodic states with low and high pinwheel densities has been found which leads to a pinwheel statistics of the repertoire of optimal solutions that quantitatively reproduces experimental observations [38, 93]. What is the reason for this difference between the two models? In [35], the energetic degeneracy of aperiodic states with low and high pinwheel densities was derived from a so-called permutation symmetry

$$N_3^z\left[u,v,w\right]=N_3^z\left[w,u,v\right],$$

(34)

of the cubic nonlinearities of the model. It can be easily seen, that the cubic nonlinearities obtained in the third order expansion of the EN model do not exhibit this permutation symmetry (see 'Methods' section). As shown by Reichl [94], the absence of permutation symmetry can lead to a selection of a subrange of pinwheel densities in the repertoire of optima of OPM models. Depending on the degree of permutation symmetry breaking, the family of optima of such models, albeit encompassing aperiodic OPM layouts, may consist of layouts with either unrealistically low or high pinwheel densities. Furthermore, for very strong permutation symmetry breaking, stationary solutions from solution classes other than the n-ECPs and rPWCs with low or high pinwheel densities may become optima of models for OPM development. In order to determine a regime in which the EN model optima quantitatively resemble experimentally observed OPM layouts, it is therefore important to quantify the degree of permutation symmetry breaking in the EN model and to examine whether permutation symmetric limits exist. As shown in the 'Methods' section, any cubic nonlinearity $N_3^z\left[z,z,\bar{z}\right]$ that obeys Equation (34) has a corresponding angle-dependent interaction function g(α) which is π-periodic. Therefore, we examine the degree of permutation symmetry breaking in the EN model by comparing the angle-dependent interaction function g(α) of its third order expansion (see Equation 33 and Figure 11) to the π-periodic function g _pm (α) = 1/2 (g(α) + g(α + π)). This 'permutation-symmetrized' part of the angle-dependent interaction function of the EN model for general orientation stimulus ensembles reads

$$\begin{array}{ll}\hfill g_{pm}\left(\alpha \right)& =\frac{2〈\mid s_z{\mid }^4〉}{{\sigma }^6}{e}^{-2k_c^2{\sigma }^2}{sinh}^4\left(1∕2k_c^2{\sigma }^{2}cos\alpha \right)\\ -\frac{〈\mid s_z{\mid }^2〉}{2{\sigma }^4}{e}^{-2k_c^2{\sigma }^2}\left(\left(cosh\left(2k_c^2{\sigma }^{2}cos\alpha \right)-2cosh\left(k_c^2{\sigma }^{2}cos\alpha \right)\right)-2e^{k_c^2{\sigma }^2}-e^{2k_c^2{\sigma }^2}\right)\\ +\frac{1}{2{\sigma }^2}\left(1+e^{-2k_c^2{\sigma }^2}cosh\left(2k_c^2{\sigma }^{2}cos\alpha \right)\right).\end{array}$$

(35)

A comparison between g _pm (α) and g(α) is depicted in Figure 14a-d. It shows that essentially insensitive to the interaction range σ/Λ, at large values of the fourth moment original and permutation symmetrized angle-dependent interaction functions converge. We quantified the degree of permutation symmetry breaking with the parameter

Figure 14

Quantifying permutation symmetry breaking in the EN model. (a-d) g(α) (red traces) and the 'permutation symmetrized' function g _pm (α) = 1/2(g(α) + g(α + π)) (blue traces, see Equation (35)) for σ/Λ = 0.1 and 0.3 and s₄ = 6 and 100. (e) Permutation symmetry parameter d (Equation 36) in the EN model with fixed retinotopy. Permutation symmetry breaking is largest for σ/Λ ≈ 0.25 and small s₄. In the limit s₄ → ∞, permutation symmetry is restored.

$$d=\frac{\parallel g-g_{pm}{\parallel }_2}{\parallel g{\parallel }_2}\mathsf{\text{sgn}}\left(g\left(0\right)-g\left(\pi \right)\right).$$

(36)

This parameter is zero in the case of a permutation symmetric cubic nonlinearity. In the case of a g-function completely antisymmetric around α = π/2, the parameter is either plus or minus one, depending on whether the maximum of g _pm is at zero or π. If d is smaller than zero, low pinwheel densities are expected to be energetically favored and vice versa. The values of d in parameter space is depicted in Figure 14e. It is smaller than zero in the entire phase space, implying a tendency for low pinwheel density optimal states, in agreement with the phase diagram in Figure 12. Permutation symmetry breaking is largest for σ/Λ around 0.25 and small fourth moment values of the orientation stimulus distribution. It decays to zero for large fourth moments proportionally to 1/s₄ as can be seen by inserting Equations (33) and (35) into Equation (36). In the infinite fourth moment limit s₄ → ∞, the cubic nonlinearities of the third-order expansion of the EN model become permutation symmetric.

In this case, the EN model is parameterized by only one parameter, the effective intracortical interaction range σ/Λ and we obtain a rather simple phase diagram (Figure 15). Optimal solutions are n-ECPs for increasing σ/Λ and we observe a sequence of phase transitions toward a higher number of active modes and therefore more complex spatially aperiodic OPM layouts. Importantly, for a subregion in the phase diagram with given number of active modes, all possible n-ECP mode configurations are energetically degenerate. It is precisely this degeneracy that has been previously shown to result in a pinwheel statistics of the repertoire of aperiodic optima which quantitatively agrees with experimental observations [38]. Therefore, our unbiased search in fact identified a regime, namely a very large effective interaction range and infinite fourth moment of the orientation stimulus ensemble, in which the EN model formally predicts which quantitatively reproduce the experimentally observed V1 architecture.

Figure 15

Phase diagram of the EN model with fixed retinotopy in the permutation symmetric limit s ₄ → ∞. The graphs show the regions on the σ/Λ-axis (lower axis) and the corresponding η-axis (upper axis), where n-ECPs or sPWCs have minimal energy. High n-ECPs (n ≳ 10) exhibit universal pinwheel statistics. Note however the extremely small η-values for large σ/Λ.

Unexpectedly, however, this regime coincides with the limit of applicability of our approach. Permutation symmetry is exactly obtained by approaching stimulus distribution with diverging fourth moment for which the amplitude equations may become meaningless. We would generally expect that the EN for very large but finite fourth moment can closely resemble a permutation symmetric model. However, to consolidate the relevance of this regime, it appears crucial to establish the robustness of the limiting behavior to inclusion of retinotopic distortions.

Optimal solutions of the EN model with variable retinotopy and arbitrary orientation stimulus ensembles

In the EN model for the joint mapping of visual space and orientation preferences, the angle-dependent interaction functions depend on four parameters: η, σ, the fourth moment 〈|s _z |⁴〉of the stimulus ensemble and η _r . By setting σ = σ*(η), we are left with three free parameters at criticality. Therefore, a three-dimensional phase diagram now completely describes pattern selection in the EN model. For better visualization, in Figure 16, we show representative cross sections through this three-dimensional parameter space for fixed η _r . First, we note the strong similarity between the phase diagram for fixed retinotopy (Figure 12) and the cross sections through the phase diagrams for the joint mappings shown in Figure 16. This expresses the fact that retinotopic mapping and OPM are only weakly coupled or mathematically, g _r (α) ≪ g(α) in all parameter regimes (see Appendix 2). Again, for distributions with small fourth moment, optimal mappings consist of either pinwheel-free orientation stripes or sPWCs. Above a certain fourth moment value around s₄ = 6, higher coupled n-ECPs are optimal. For small interaction range σ/Λ, hexagonal pinwheel crystals (coupled 3-ECPs) represent optimal mappings in a large fraction of parameter space. With increasing σ/Λ, we observe a sequence of phase transitions by which higher n-ECPs become optimal. Anisotropic planforms at the lower end of the spectrum of pinwheel densities are always energetically favored over high pinwheel density layouts. The only difference between the cross sections is that the region covered by sPWCs increases for decreasing η _r . The phase diagram for large η _r = 10 is virtually indistinguishable from the phase diagram in Figure 16.

Figure 16

Stripe-like, crystalline, and quasi-crystalline cortical representations as optimal solutions to the joint mapping problem of visual space and orientation preference in the EN. (a-d) Phase diagrams for the joint mapping of visual space and orientation preference in the EN near criticality for η _r = 0 (a), η _r = 0.01 (b), η _r = 0.1 (c), and η _r = 10 (d). The graphs show the regions of the s₄-σ/Λ-plane in which coupled n-ECPs or sPWCs have minimal energy. For n ≥ 3, pinwheel densities of the energetically favored n-ECP configuration are indicated. Note the strong similarity between the phase diagrams and the phase diagrams in the fixed retinotopy case (Figure 12).

Optimal mappings of orientation preference are thus either orientation stripes, periodic arrays of pinwheels (hexagonal, quadratic) or quasi-periodic pinwheel arrays with low pinwheel density. Retinotopic distortions lead to lower gradients of the retinotopic mapping at high gradient regions of the OPM. This is in line with some of the experimental evidence [55, 95] but contradicts others [96].

Most importantly, we note that the results on permutation symmetry breaking in the fixed retinotopy case are not altered by allowing for retinotopic distortions. Since g _r (α) does not depend on the fourth moment of the orientation stimulus distribution, non-permutation symmetric terms decay as 1/s₄ for large s₄.

Hence, in the limit s₄ → ∞, permutation symmetry is restored and we recover the phase diagram in Figure 15 also for the EN model with variable retinotopy independent of η _r . As the energy contribution of retinotopic deviations r(x) becomes negligible in the infinite fourth moment limit, the optima are then simply the corresponding coupled n-ECPs and these states are energetically degenerate for fixed n. For very large effective interaction range and infinite fourth moment of the orientation stimulus ensemble, the EN model with variable retinotopy is able to quantitatively reproduce the experimentally observed pinwheel statistics in OPMs. It furthermore predicts reduced gradients of the visual space mapping at high gradient regions of the OPM.

Finite stimulus samples and discrete stimulus ensembles

Our reexamination of the EN model for the joint optimization of position and orientation selectivity has been so far carried out without addressing the apparently fundamental discrepancy between our results and the large majority of previous reports. Since the seminal publication of Durbin and Mitchison [21], numerous studies have used the EN model to simulate the development of visual cortical maps or to examine the structure of optimal mappings by numerical simulation [58, 62–65, 79, 97, 98]. These studies have either used the circular or the uniform orientation stimulus ensemble for which, to the best of our knowledge, the only two nontrivial stationary solutions are square pinwheel crystals or orientation stripes. Furthermore, we found that the gradient descent dynamics seems to readily converge to the respective minima of the EN free energy. This indicates that other local minima and more complex intrinsically aperiodic states are not dominant in this model. In fact, we found that all aperiodic stationary solutions we could perturbative calculate analytically are unstable and thus represent hyperbolic saddle points and not local minima. As these stable solutions barely resemble experimentally observed OPMs, it is not obvious how the EN model in all of these studies could appear as a model well suited to describe the complex layout of real cortical orientation maps. Prior studies however often used computational methods different from our fixed parameter steepest descent simulations.

Two alternative approaches have been used predominantly to study dimension reducing mappings for cortical representations. These methods have been applied to both the EN model and the other widely used dimension reduction model, the self-organizing feature map (SOFM), originally introduced by Kohonen [59]. The simplest way to compute mappings from a high dimensional feature space onto the two-dimensional model cortex is by iterating the following procedure for a large number of randomly chosen stimuli (e.g., [56, 57, 66–68, 99, 100]): (i) Stimuli are chosen one at a time randomly from the complete feature space. (ii) The activation function for a particular stimulus is computed. In the case of the EN model, this activation function can acquire a rather complex form with multiple peaks (see 'Discussion' section). In the case of an SOFM, this activation function is a 2D-Gaussian. (iii) The preferred features of the cortical grid points are updated according to a discretized version of Equations (3, 4) or the corresponding equations for the SOFM model. Typically, this procedure is repeated on the order of 10⁶ times. The resulting layout is then assumed to at least approximately solve the dimension reduction problem. In many studies, small stimulus sets have been chosen presumably for computational efficiency and not assuming specifically that the cortex is optimized for a discrete finite set of stimuli. In [21] for instance, a set of 216 stimuli was used, that was likely already at the limit of computing power available at this time.

In a more refined approach, the EN model as well as Kohonen's SOFM model have been trained with a finite set of stimuli (typically with on the order of 10³ to 10⁴) and the final layout of the model map has been obtained by deterministic annealing [101], i.e., by gradually reducing the numerical value of σ in a numerical minimization procedure for the energy functional $\mathcal{F}$ at each value of σ (see e.g., [21, 64, 65, 79] and see 'Methods' section). In such simulations, often nonperiodic boundary conditions were used. One might suspect in particular the second approach to converge to OPM layouts deviating from our results. It is conceivable in principle, that deterministic annealing might track stationary solutions across parameter space that are systematically missed by both, our continuum limit analytical calculations as well as our descent numerical simulations.

To assess the potential biases of the different approaches, we implemented (i) finite stimulus sampling in our gradient descent simulations and (ii) studied the results of deterministic annealing simulations varying both the size of the stimulus set as well as the type of boundary conditions applied.

We simulated Equations (3, 4) with finite sets of stimuli of different sizes (see 'Methods' section), drawn from the circular stimulus ensemble. Following [21, 25], η was set to a small value (η = 0.025) such that the optimal configuration for the joint mapping of visual space and orientation preference is the coupled 1-ECP (see Figure 10), i.e., a pattern of parallel orientation stripes without any retinotopic distortion (see Figure 9). Figure 17 displays representative simulations for stimulus sets of size N = 216 (as used in [21]) (a), N = 10⁵ (b), N = 10⁶ (c) stimuli. Simulation time t is measured in units of the intrinsic time scale τ (see 'Methods' section). For N = 216 stimuli, RM and OPM quickly reach an apparently stationary configuration with a large number of pinwheels at around t = 20τ. Power is distributed roughly isotropically around the origin of Fourier space (k = 0). The stable OPM lacks a typical length scale and, expressing the same fact, the power spectrum lacks the characteristic ring of enhanced Fourier amplitude. Retinotopic distortions are fairly pronounced. Both obtained maps resemble the configurations reported in [21].

Figure 17

Development of OPM and retinotopic distortions in EN simulations with fixed stimulus sets of different sizes. (a) OPMs (left), their power spectra (middle) and RMs (right) for t = 10τ (upper row) and t = 300τ (lower row) obtained in simulations with fixed stimulus set (η = 0.028, σ/Λ = 0.3, s4 = 4/3, 216 stimuli). (b) 10⁵ stimuli (all other parameters as in a). (c) 10⁶ stimuli (all other parameters as in a). Large stripe-like OP domains are generated via pairwise pinwheel annihilation for large simulation times. Retinotopic distortions are fairly weak.(d) Pinwheel density time course for EN simulations with fixed stimulus sets of different sizes, including the simulations from a to c (red, green, blue traces 216, 10⁵, 10⁶ stimuli) (all other parameters as in a). Dashed lines represent individual simulations, solid lines an average over four simulations. Note, that the pinwheel density rapidly decays below 2.0 in both cases, and in particular for 10⁶ stimuli, the OPM pattern acquires large stripe-like regions.

For N = 10⁵ stimuli, we find that OPMs exhibit a characteristic scale (see dark shaded ring in the power spectrum) and a dynamic rearrangement of the maps persists at least until t = 200τ. Stripe-like OP domains are rapidly generated via pairwise pinwheel annihilation for t >10τ. Retinotopic distortions are fairly weak. For N = 10⁶ stimuli, again OPMs exhibit a characteristic scale (see dark shaded ring in the power spectrum) and the map dynamics persists beyond t = 200τ. A larger fraction of the pinwheels annihilate pairwisely compared to N = 10⁵ stimuli, leading progressively to a pattern with large stripe-like domains. Retinotopic distortions are fairly weak. For both cases with massive stimulus sampling (N = 10⁵, N = 10⁶), the pinwheel density rapidly drops below the range observed in tree shrews, galagos and ferrets and than further decreases during subsequent map rearrangement. In summary, the more stimuli are chosen for the optimization procedure, the less pinwheels are preserved in the pattern of orientation preference and the more the resulting map resembles the analytically obtained optimal solution. Deterministic annealing approaches which change parameters of the energy functional during the computational minimization process differ more fundamentally from our gradient descent simulations than the iterative schemes used with fixed parameters. Studies using deterministic annealing in addition frequently used nonperiodic boundary conditions (e.g. [64, 65, 79]). To study all potential sources of deviating results, we implemented deterministic annealing for the EN energy function (see 'Methods' section, Equation (46)) for periodic boundaries, nonperiodic boundary conditions as well as random and grid-like finite stimulus ensembles (see 'Methods' section). We closely follow the refined methods used in [64, 65, 79] and performed deterministic annealing simulations for the EN model with retinotopic distortions and stimuli drawn from the circular stimulus ensemble.

Figures 18a and 19a display representative simulations for random stimulus sets of size N = 10³, N = 10⁴ and N = 10⁵ for periodic boundary conditions (Figures 18a) and nonperiodic boundary conditions (Figures 19a). Furthermore depicted are the pinwheel densities of stationary solutions as well as their energies, relative to the energy of a pinwheel-free stripe solution (see 'Methods' section) for different annealing rates ξ (Figures 18b-d and 19b-d). Figures 18e-g and 19e-g additionally show the statistics of nearest neighbor (NN) pinwheel distances as well as the SD of the pinwheel densities for randomly selected subregions in the OPM as introduced in [38], averaged over four simulations with N = 10⁵. To facilitate comparison, we superimposed fits to the experimentally observed statistics [38] for orientation maps in tree shrews, ferrets and galagos.

Figure 18

The EN model with periodic boundary conditions, solved with deterministic annealing. (a) OPMs (left) and RMs (right) for N = 10³ (upper row), N = 10⁴ (middle row) and N = 10⁵ (lower row) random stimuli and periodic boundary conditions (annealing rate χ = 0.999, see 'Methods' section). β is the continuity parameter in the conventional definition of the EN model (see 'Methods' section, Equation 46) and is scaled, such that a comparable number of columns is emerging in the simulations for each size of the stimulus set. (b) Pinwheel densities of EN solutions for different numbers of stimuli, χ = 0.999. (c) Pinwheel densities of EN solutions for 10⁵ stimuli and different annealing rates. (d) Energies of solutions for 10⁵ stimuli, relative to the energy of a pinwheel-free stripe solution (see 'Methods' section) for different annealing rates. (b-d) Crosses mark individual simulations, red line indicates average values. (e, f) Statistics of nearest neighbor pinwheel distances for pinwheels of (e) arbitrary and (f) opposite and equal charge for 10⁵ random stimuli and periodic boundary conditions, averaged over four simulations (red curves). Black curves represent fits to the experimental data from [38]. (g) Standard deviations (SD) of pinwheel densities estimated from randomly selected regions in the OPM. Black dashed curve indicates SD for a two-dimensional Poisson process of equal density.

Figure 19

The EN model with nonperiodic boundary conditions, solved with deterministic annealing. (a-g) As Figure 18, but for nonperiodic boundary conditions.

When annealing with periodic boundary conditions, the maps found with deterministic annealing essentially resemble our gradient descent dynamics simulations. The larger the set of stimuli, the more stripe-like are the OPMs obtained (Figure 18a,b). Furthermore, the more carefully we annealed, the lower the pinwheel density of the obtained layouts (Figure 18c). For N = 10⁵, the pinwheel density averaged over four simulations with annealing rate 0.999 was ρ = 2.04 As expected, the energy of the final layouts decreased with slower annealing rates (Figure 18d). However, when starting from random initial conditions, the energy of the final layouts found was always higher compared to the energy of a pinwheel-free stripe solution (see 'Methods' section for details), which is the predicted optimum for the circular stimulus ensemble. NN-pinwheel distance histograms are concentrated around half the typical column spacing and in particular pinwheel pairs with short distances are lacking completely (Figure 18e,f). For nonperiodic boundary conditions and random stimuli, we found that retinotopic distortions are more pronounced than for periodic boundary conditions. They however decreased with increasing number of stimuli. For large the stimulus numbers, we observed stripe-like orientation preference domains which are interspersed with lattice-like pinwheel arrangements (see Figure 19c), lower row, upper left corner of the OPM). For N = 10⁵, the pinwheel density averaged over four simulations with annealing rate 0.999 was ρ = 2.71.

Similarly to the results for periodic boundary conditions, short distance pinwheel pairs occur less frequently than in the experimentally observed maps, indicating an increased regularity in the pinwheel distances compared to real OPMs (Figure 19e,f). This regularity is further indicated by a smaller exponent of the SD compared to the Poisson process (Figure 19g). The perfect stripe-like solution is not the optimum for nonperiodic boundaries. The energy of the map layouts found with very slow annealing rates is slightly lower than the energy of the pinwheel-free OPM layout (Figure 19d). We note that the layout of the OPM at the boundaries does not differ substantially from the layout inside the simulated domain, suggesting that boundary effects affect the entire simulated domain for the relatively small region treated. Finally, we performed simulations with grid-like stimulus patterns as e.g., used in [64, 65]. These simulations displayed a strong tendency toward rhombic pinwheel arrangements, i.e., the second stable stationary solution found for the circular stimulus ensemble. We refer to Appendix 3 for further details. In summary, our results for the discrete EN model with deterministic annealing largely agree with the analytical results. Irrespective of the numerical methodology, the emerging map structure for large numbers of stimuli is confined to the states predicted by our analytical treatment of the continuum formulation of the EN. This behavior is expected because the energies underlying the deterministic annealing and the steepest descent simulations are mathematically equivalent (see 'Methods' section). In any kind of deterministic annealing simulation we tested, resulting patterns were patchworks of the two fundamental stable solutions identified by the analytical treatment: pinwheel free stripes and square lattices of pinwheels. Such patchworks are spatially more complicated than perfect stripes or crystals. Nevertheless, they qualitatively differ in numerous respects from the experimentally observed spatial arrangements (see Figures 18, 19, and 28 in Appendix 3). How the fundamental stable solutions are stitched together somewhat differs between the different kinds of simulations. For instance, using a grid-like stimulus ensemble with nonperiodic boundary conditions apparently energetically favors the rPWC compared to the pinwheel-free stripe regions (see Figure 27 in Appendix 3). In summary, while some of the patterns obtained by deterministic annealing might be called 'good-looking' maps, all of them substantially deviate from the characteristics of experimentally observed pinwheel arrangements.

We conclude that the differences between our results and those of previous studies are most likely due to the small finite stimulus samples used largely for reasons of computational tractability. Deterministic annealing using stimulus samples that fill the feature space converges to the same types of patterns found by perturbation theory. We further conclude that our methods do not systematically miss biologically relevant local minima of the classical EN energy function.

Discussion

Summary

In this study, we examined the solutions of what is perhaps the most prominent optimization model for the spatial layout of orientation and RMs in the primary visual cortex, the EN model. We presented an analytical framework that enables us to derive closed-form expressions for hyperbolic fixed points, local and global minima, and to analyze their stability properties for arbitrary optimization models for the spatial layout of OPMs and RMs. Using this framework, we systematically reexamined previously used instantiations of the EN model, dissecting the impact of stimulus ensembles and of interactions between the two maps on optimal map layouts. To our surprise, the analysis yielded virtually identical results for all of these model instantiations that substantially deviate from previous numerical reports. Pinwheel-free orientation stripes and crystalline square lattices of pinwheels are the only optimal dimension-reducing OPM layouts of the EN model. Both states are generally stable but exchange their roles as optima and local minima at a phase border. Numerical simulations of the EN gradient descent dynamics as well as simulations utilizing deterministic annealing confirmed our analytical results. For both processes, the initially spatially irregular layouts rapidly decayed into a patchwork of stripe-like or crystal-like local regions that then became globally more coherent on longer timescales. Pinwheel-free solutions were approached after an initial phase of pattern emergence by pairwise pinwheel annihilation. Crystalline configurations were reached by the generation of additional pinwheels and pinwheel annihilation together with a coordinated rearrangement toward a square lattice. These results indicate that layouts which represent an optimal compromise of coverage and continuity for retinotopy and orientation do not per se reproduce the spatially aperiodic and complex structure of orientation maps in the visual cortex.

To clarify whether the EN model is in principle capable of reproducing the biological observations, we performed an unbiased comprehensive inspection of EN optima for arbitrary stimulus distributions possessing finite fourth moments. This analysis identified two key parameters determining pattern selection: (i) the effective intracortical interaction range and (ii) the fourth moment of the orientation stimulus distribution. We derived complete phase diagrams summarizing pattern selection in the EN model for fixed as well as variable retinotopy. Small interaction ranges together with low fourth moment values lead to either pinwheel-free orientation stripes, rhombic or hexagonal crystalline orientation map layouts as optimal states. Large interaction ranges together with orientation stimulus distributions with high fourth moment values lead to the stabilization of irregular aperiodic OPM layouts. These solutions belong to a class of solutions previously called n-ECPs. This solution class encompasses a large variety of OPM layouts and has been identified as optimal solutions of abstract variational models of OPM development [35]. We showed that in the EN model due to a lack of a so-called permutation symmetry, among this family of solutions, states with low pinwheel densities are selected as global minima. In the extreme and previously unexplored parameter regime of very large effective interaction ranges and stimulus ensemble distributions with infinite fourth moment, permutation symmetry is restored and spatially aperiodic OPM layouts with higher pinwheel density are included in the repertoire of optimal solutions. Only in this limit, the repertoire of optima reproduces the recently described species-insensitive OPM design [38] and quantitatively matches experimentally observed orientation map layouts. None of these findings depend on whether the EN model is considered with variable or fixed retinotopy.

Comparison to previous studies

It is an important and long-standing question, whether the structure of cortical maps of variables such as stimulus orientation or receptive field position can be explained by a simple general principle. The concept of dimension reduction is a prominent candidate for such a principle (see e.g., [58, 102] for reviews) and the qualitative agreement between experimental data and previous numerical results from dimension reduction models [21, 42, 60, 62–66, 68, 98, 102–104] can be viewed as evidence in favor of the dimension reduction hypothesis. Yet comprehensive analytical investigations of dimension reduction problems and in particular the determination of their optimal and nearly optimal solutions have been impeded by the mathematical complexity of these problems. For the EN algorithm applied to the TSP, previous analytical results established the unselective fixed point above the first bifurcation point as well as the parameters at which this solution becomes unstable [105]. Subsequent work extended these results to the EN model for cortical map formation. The periodicity of solutions depending on the model parameters has been obtained by computing the eigenvalues of the Hessian matrix of the energy function [63, 97, 106]. Hoffsümmer et al. [72] confirmed these results, and computed the periodicity of the emerging patterns in the continuous EN model formulation by linear stability analysis of the EN gradient descent dynamics as used in this study. Our results extend these findings and for the first time provide analytical expressions for the precise layout of optimal and nearly optimal dimension-reducing maps.

In the light of the qualitative agreement between experimental data and numerical solutions of the EN model previously described, it is perhaps our most surprising result that the model's optimal dimension-reducing maps are regular periodic crystalline structures or pinwheel-free stripe patterns in large regions of parameter space. In particular, the species-insensitive pinwheel statistics observed experimentally [38] are not exhibited by optimal solutions of the classical EN in any of the previously considered parameter regimes.

Our comparison of different numerical approaches indicates that the differences to previous studies are mainly attributable to differences in the sampling of the stimulus manifold in the numerical optimization procedures. In their seminal publication, Durbin and Mitchison used sets of 216 stimuli from the circular stimulus ensemble and applied a Gauss-Seidel procedure to obtain stationary configurations [21]. A similar procedure was used in [104]. Quite frequently, the number of stimuli used for optimization is of the same order of magnitude as the number of model neurons or centroids in feature space. This provides a relatively sparse sampling of the stimulus manifold [63–65]. Finite stimulus sampling effects are expected to worsen when feature spaces of higher dimension are considered.

The choice of small stimulus sets in previous dimension reduction studies was imposed mainly by the limitations of computing power. Using a parallelized implementation of the Cholesky method for deterministic annealing [62–65] on a multicore architecture with 2 TB working memory, we explored the dependence of the obtained near optimal solutions on the sampling of the feature space manifold over two orders of magnitude. We find that, the more stimuli are sampled, the closer the numerically obtained configurations resemble our analytical predictions. Our results on the classical EN model with deterministic annealing suggest that in the limit of large stimulus numbers, one would perfectly recover our analytical results both for periodic conditions or nonperiodic boundary conditions with realistic system sizes. This dense stimulus sampling limit is also readily visible in our reproduction of the original Durbin and Mitchison sampling and the modification of the predicted map structure with stimulus number (Figure 17). The finding that computational limitations prevented Durbin and Mitchison from obtaining the genuine predictions of their dimension reduction model should not be viewed as diminishing the importance of their contribution. The dimension reduction approach has played a unique and extremely productive role in guiding the conceptualization of cortical functional architecture. It has established an abstract view on cortical representations without which most of our current theoretical knowledge about candidate theories for cortical architectures could not have been obtained.

Our results about optimal states of the EN for the circular and uniform stimulus ensembles however agree with some prior work. In [25], the gradient descent dynamics of the EN model was used as a model for the emergence and refinement of cortical maps during development. Simulated visual stimulus features included retinotopy, orientation and eye dominance. The numerical procedures were similar to the one developed in this study. Parameters were chosen such that s₄ = 4/3 and σ/Λ ≈ 0.366. This study found that an initially large number of pinwheels decayed via pairwise annihilation of pinwheels with opposite topological charge. Our analysis predicts a stripe-like OP pattern as optimal solution in this regime, both in the case of a fixed uniform retinotopy as well as with variable retinotopy. In our simulations, this state is reached after an initial phase of symmetry breaking with the generation of numerous pinwheels via pairwise pinwheel annihilation. Our analytical and numerical results thus confirm, explain, and generalize these previous findings.

The previous results also indicated that the inclusion of eye dominance in the EN model slightly slows down but does not stop the pinwheel annihilation process (see [25], Figure 3). This raises the possibility that the main features of our analysis of optimal solutions for the EN model may persist when additional feature dimensions are taken into account. Reichl et al. in fact observed that models with interacting OPM and OD maps (ODMs) exhibit a transition from pinwheel-free stripes to periodic pinwheel crystals similar to the transitions found in the EN [37] and demonstrated that this transition is a general feature of models with interacting OPM and ODMs [107]. A rigorous characterization of map structures predicted by the simultaneous optimization of multiple periodic feature representations such as orientation preference and OD constitutes an important goal for future studies. The recent study by Reichl et al. [37] suggests that this issue can successfully be approached using concepts from the nonlinear dynamics of pattern formation. Finally, one recent study used the continuous formulation of the EN model to investigate the impact of postnatal cortical growth on the formation of OD columns in cat visual cortex [69]. Consistent with our results, this study also observed perfectly regular stripe-like patterns as stationary states in gradient descent simulations. The dynamics of the convergence of the ODMs toward the stripes was modified by including cortical growth into the model. However, as soon as growth terminated, simulated ODM layouts readily converged toward regular stripes. How cortical growth interacts with the formation of orientation columns is currently not understood and represents a further interesting topic for future studies.

Geometric relationship between retinotopic distortions and OPMs

Experimental results on the geometric relationships between the map of visual space and the map of orientation preference are ambiguous. Optical imaging experiments in cat V1 suggested a systematic covariation of inhomogeneities in the RM with singularities in the pattern of orientation columns in optical imaging experiments [96]. Regions of high gradient in the map of visual space preferentially appeared to overlap with regions of high gradient of the OPM. In ferret, however, it has been reported that high gradient regions of the map of visual space correspond to regions of low gradient in the OPM [67]. In tree shrew V1, no local relationships between the mapping of stimulus orientation and position seem to exist and the map of visual space appears to be ordered up to very fine scales [108]. In line with this, single unit recordings in cat area 17 revealed no correlation between receptive-field position scatter and orientation scatter across local cell ensembles [109, 110].

Our analysis of the EN model shows that its optimal states exhibit a negative correlation between the rates of change of orientation preference and retinotopic position, similar to what has been observed in the ferret [67]. This is expected from the principle of dimension reduction and in agreement with the original numerical results by Durbin and Mitchison [21]. However, both in simulations of the gradient descent dynamics and in deterministic annealing simulations with periodic boundary conditions as well as in analytically obtained optimal solutions, deviations from a perfectly uniform mapping of visual space are surprisingly weak (see Figures 9, 10, 18, and 25 in Appendix 1).

Deterministic annealing simulations with open nonperiodic boundary conditions showed a substantially increased magnitude of retinotopic distortions. This raises the possibility that different behaviors observed in different experiments might be at least partially related to the influence of boundary effects. The influence of boundary effects is expected to decline into the interior of an area, in particular for large areas as V1 (see [111]). In the bulk of V1, we thus expect only a weak coupling of orientation map and retinotopic distortions according to the EN model. In this regime, the predictions from models with reduced rotational symmetry (so-called Shift-Twist symmetry [112]) about the coupling between retinotopic distortions and OPMs [33] appear to be more promising than the weak effects resulting from the coverage-continuity-compromise. Consistent with the measurements of Das and Gilbert [96], such models predict small but significant positive correlations between the rates of change of orientation preference and retinotopic position [33]. Moreover, the form of the retinotopic distortions in such models is predicted to differ for pinwheels with positive and negative topological charge [92]. This interesting prediction of OPM models with Shift-Twist symmetry deserves to be tested by measuring the receptive field center positions around the two types of pinwheels with single cell resolution [12].

Aperiodic OPMs reflect long-range intracortical suppression

Our unbiased search through the space of stimulus ensembles with finite fourth moment revealed the existence of spatially aperiodic optimal solutions in the EN. It is important to realize that the selection of these solutions is not easily viewed as resulting from an optimal compromise between coverage and continuity. In fact, the continuity parameter in the respective parameter regime is so small that solutions essentially maximize coverage (see Figures 12, 15, and 16). Instead, this phenomenon reflects a different key factor in the stabilization of pinwheel-rich aperiodic layouts, namely the dominance of long-ranged and effectively suppressive interactions. This is illustrated in Figure 20 which depicts different forms of stimulus-evoked activity patterns in the EN model. For a short-range interaction (Figure 20a), the activity evoked by low as well as high orientation energy stimuli is an almost Gaussian activity peak located near the stimulus position. The peak is shallow for low (left) and sharp for high 'orientation energy' (right). In the corresponding parameter regime, square pinwheel crystals are the optimal solution of the EN. For a longer range of interaction where aperiodic OPM layouts are the optimal states, the activity evoked by a single point-like stimulus is qualitatively different. Here, the activity pattern is extended and spans several hypercolumns (Figure 20b). It is weakly modulated for low orientation energy stimuli (left) and consists of several distinct peaks for high orientation energy stimuli (right). In this regime, neurons at a distance of several columns compete for activity through the normalization term in the EN which leads to a nonlocal and effectively suppressive intracortical interaction.

Figure 20

Different patterns of evoked activity for different effective ranges of intracortical interaction in the EN model. The component 〈s _z |z(x)〉 of the orientation map z(x) in the direction of the stimulus s _z is plotted as a meshed 3D graph in a 6Λ × 6Λ patch. Color code and height of the projection below indicate the strength of activation. The stimulus is presented in the center of the displayed cortical subregion. (a) Evoked activity patterns e(x, S, z, r) for small interaction range σ/Λ = 0.1 and weakly oriented stimulus with s _z = 0.01 (left) and strongly oriented stimulus with s _z = 4 (right). rPWCs (see upper right) are optimal in this regime. (b) Evoked activity patterns e(x, S, z, r) for large interaction range σ/Λ = 0.9 and weakly oriented stimulus with s _z = 0.01 (left) and strongly oriented stimulus with s _z = 4 (right). Spatially aperiodic 8-planforms (see upper right) are optimal in this regime. A uniform retinotopy was assumed in all cases for simplicity.

It is presumably not a mere coincidence that recent studies of abstract variational models of OPM development [35, 38, 93] mathematically identified this type of interaction as a key mechanism for stabilizing realistic OPM layouts. It has been shown that all models for OPM development that share the basic symmetries (i) translational symmetry (ii) rotational symmetry (iii) shift symmetry and (iv) permutation symmetry and in addition are dominated by long-range suppressive interactions, form a universality class that generates maps exhibiting a universal and realistic pinwheel statistics. In such models, suppressive long-range interactions are key to stabilizing irregular arrangement of pinwheels, which otherwise largely disappear or crystalize during optimization. We have stressed that the EN model as considered here obeys the symmetries (i) to (iii). In the limit of infinite orientation stimulus ensemble fourth moment, permutation symmetry (iv) is restored. The EN can thus be tuned into the above universality class by sending the orientation stimulus distribution fourth moment to infinity and choosing an exponentially small continuity parameter to realize effective long-range coupling. Indeed, the phase diagrams for abstract variational models of OPM development [35] and those of the EN model found here are structurally very similar. In both cases, a rather large orientation stripe phase is complemented by a cascade of phase transitions toward more complex, aperiodic and pinwheel-rich OPM layouts induced by long-range suppressive interactions. Using abstract variational models, it has been shown recently that the stabilization of regular crystalline pinwheel layouts can alternatively be achieved by a strong coupling between the map of orientation and the map of eye dominance [37, 107]. The structure of the phase diagrams of such models however appears fundamentally different from the structure of the EN phase diagrams.

The parameter regime in which the EN model's optimal solutions exhibit the experimentally observed pinwheel statistics is not at all intuitive and in our opinion questions the conventional interpretation of the EN model for the formation of cortical feature maps. Firstly, the extremely small continuity parameter questions the fundamental role of a tradeoff between coverage and continuity. We note that such a parameter regime is currently not accessible to numerical simulations. In addition, an apparently fundamental property for any adequate model for OPM optimization or development, namely a Turing-type finite wavelength instability of the unselective state [32], is lost in the limit η → 0. At first sight the infinite fourth value required may appear reminiscent of the power-law distributions for orientation energy found in the statistics of natural images [113, 114]. However, as visualized in Figure 20b, the essential property of the EN model in the infinite fourth moment regime is the occurrence of patterns of activity spatially extended beyond a single hypercolumn representing spatially localized point-like stimuli. These activity patterns mediate the long-range interactions between distant orientation columns which in turn cause the stability of realistic pinwheel-rich aperiodic OPM layouts. It is obvious that spatially extended stimuli provide a much more plausible and realistic source of extended activity patterns in models for visual cortical development (for an extended discussion see [54]). Optimization models for cortical maps based on the representation of more complex spatially extended visual stimuli, such as natural scenes, rather than a model based on point-like stimuli with extreme statistics would then be a more appropriate basis for understanding visual cortical functional architecture.

Comparison to the SOFM model

Several alternatives to the EN model have been proposed as optimization approaches that can account for the structure of visual cortical maps. One prominent alternative dimension reduction model is the so-called SOFM, originally introduced by Kohonen [59]. It is widely believed that this model, albeit lacking an exact energy functional [115], implements a competition between coverage and continuity very similar to the EN model [56, 57, 66, 115]. The SOFM has been reported to reproduce many of the experimentally observed geometric properties of visual cortical feature maps (e.g., [56, 57, 61, 66–68]). The numerical procedures used in all of these studies were either the deterministic annealing procedure or the nonrecurring application of a stimulus set without systematic assessment of pattern convergence. An analysis of the nontrivial stationary states of a dynamical systems formulation of the SOFM model is currently lacking. The main difference between the SOFM model and the EN model is that the activation function by definition has the form of a stereotypical Gaussian and competition is incorporated by a hard winner-takes-all mechanism. As a consequence, it is not obvious that a long-range suppressive interaction regime can be realized in this model. According to our analysis, one would thus expect orientation stripes and rPWCs as nontrivial stationary states of the SOFM model. In a very recent study of the SOFM algorithm that used a numerical procedure similar to the gradient descent simulations developed in this article, both pinwheel annihilation and rhombic pinwheel crystallization have been observed [116]. In addition, one study that examined the SOFM model for orientation and retinotopy found a fast convergence to pinwheel-free stripe-like solutions for a wide parameter range [25]. In view of these results, it seems worthwhile to also reexamine the SOFM model with respect to its stationary states.

Rugged or smooth energy landscape

As for many optimization problems in biology, the optimization of visual cortical functional architecture has been considered a problem characterized by a rugged energy landscape [76]. In case of the EN model the expectation of a rugged energy landscape at first sight seems quite plausible. Originally, the elastic network algorithm was invented as a fast analogue method to approximately solve NP-hard problems in combinatorial optimization such as the TSP [77, 78]. In the TSP, the stimulus positions correspond to the locations of cities a salesman has to visit on the shortest possible tour. In problems such as the TSP, the energy functionals to be minimized are known to possess many local minima and the global minimization of these functionals generally represents an extremely difficult problem [78]. Our analysis reveals that the trade-off between coverage and continuity for the mapping of a continuous feature space manifold leads to a much simpler structure of the energy landscape. This is also indicated by the fact that almost all of our gradient descent dynamics simulations readily converged to the predicted global minimum of the energy functional. Figure 21 illustrates the smooth structure of the EN energy landscape close to pattern formation threshold for different model parameters for a one dimensional path through the state space. In this landscape, the small set of stable planforms correspond to local minima of the EN energy functional, and unstable planforms to saddle points in the energy landscape. The optimal states correspond to global minima. Note that along the depicted state space path, unstable stationary solutions may appear as local minima if the unstable directions along which the energy decreases are orthogonal to the path.

Figure 21

Illustration of the EN energy landscape close to pattern formation threshold. The variation of the energy between states of ideal OS, the 2-ECP state, sPWC and possible mode configurations for 8-ECPs is shown for the case that (a) the OS state has the lowest energy (s₄ = 0, σ/Λ = 0.3), (b) the sPWC state has the lowest energy (s₄ = 0, σ/Λ = 0.1), and (c) the most anisotropic 8-ECP has lowest energy (s₄ = 100, σ/Λ = 0.8). The energy values between the state are computed from a state obtained by linear interpolation between two neighboring states on the x-axis. Note that not all local minima in a-c correspond to a stable fixed point of the amplitude dynamics (see text).

What is the origin of this qualitative difference in the shape of the energy landscapes? In the traveling salesman problem, the finite repertoire of possible tours consists of all permutations of the N cities that the salesman has to visit. By self-organized competition between the aim to visit all cities and the aim to minimize the path length, the elastic network algorithm converges to a specific ordering of the cities that eventually yields a very short tour. Most likely, the qualitative difference to the EN model for visual cortical map architecture originates from the transition from a finite number of cities to a continuum. When the elastic network algorithm is considered with an ensemble of cities (or stimuli) distributed according to a continuous probability density function, there is no discrete repertoire of tours. Both, the repertoire of tours as well as the path through the landscape of cities or equivalently the space of visual stimulus features are determined by self-organization. The first is generated by the symmetry breaking mechanism that leads to the instability of the homogeneous state. The second corresponds to the selection of one of the many nontrivial stable steady states.

An interesting property of the EN model dynamics that can be inferred from the energy landscape depicted in Figure 21 is the type of competition between two stable stationary states, where both are present in the system with a wall or a domain boundary between them. The motion of the wall or domain boundary is predicted to proceed in the direction that increases the fraction of the pattern with lower energy. An example of such competition can be seen in Figure 23g in Appendix. At t = 100τ, a small domain with an sPWC state is present. The area of this region is gradually reduced over the time course of the simulation until the pinwheel-free optimal state is reached.

Are simple OPM layouts an artifact of model simplicity?

The perfectly periodic types of stationary solutions (stripes, crystals) that appear to dominate the classical EN model for retinotopy and orientation have been found in other models of visual cortical layouts that are relatively abstract. One might therefore suspect that they represent a mere artifact of model simplicity. One conceptually appealing approach where perfectly periodic layouts have been found is wiring-length minimization [27]. According to this hypothesis, the structure of an OPM can be understood by minimizing the total length of dendritic and axonal processes. Maps obtained by stopping minimization of wire length exhibited qualitatively realistic layouts (see Figure 6 in [27]). Complete optimization, however, leads to either stripe-like pinwheel-free patterns or rPWCs, identical to the ones obtained in our investigation of optimal solutions of the EN model [27]. Similarly, stripe-like and rhombic optima have been found in several abstract vector-field approaches for OPM development [31, 33, 117].

It is ruled out by two observations, that the crystalline and perfectly periodic optima observed in all four optimization models, the EN model, the SOFM, the wiring-length minimization model, and the vector-field models are a mere artifact of the abstract order parameter field description of cortical selectivity patterns that is common to these approaches. Firstly, equally simplistic order parameter models for OPM development with long-range interactions have been shown to reproduce spatially irregular map layouts [35, 38, 93]. The occurrence of periodic optimal solutions is thus not a necessity in this model class. Secondly, pinwheel crystallization has also been observed in detailed network models for the development of OPMs, notably in the first ever model for the self-organization of orientation selectivity by von der Malsburg in 1973 [14, 118]. Thus, on the one hand the phenomenon of pinwheel crystallization is thus not restricted to simple order parameter models and on the other hand abstract and mathematically relatively simple models can exhibit complex and biologically realistic optimal solutions.

Map rearrangement and layout optimization

Irrespective of the optimization principle invoked to describe the structure of visual cortical maps, several common features of the resulting dynamics have been observed. The dynamics of optimization models usually starts with a phase of pattern emergence, where selectivity to visual features arises from an initially homogeneous unselective or weakly selective state. As we and others have shown, feature maps in these models continue to evolve after single cell selectivities reach mature levels. In fact, the phase of initial pattern emergence is typically followed by a prolonged phase of rearrangement of selectivities and preferences until a stable configuration is reached that represents a genuine optimum. This is not an exceptional type of dynamics but rather constitutes the generic expectation for a spatially extended system [83, 84].

What drives the second phase of map rearrangement? The initial emergence of feature selectivity is predominantly a local process in which merely neighboring units interact with each other to roughly match their selectivities. In the resulting spatial layout, selectivities are therefore far from being optimally arranged in space with respect to the global organization of selectivities on larger scales. Depending on the interactions incorporated in the model, local matching processes may (i) effectively propagate through space optimizing the pattern over gradually increasing spatial scales or (ii) distant sites may start to directly interact with each other to guide a rearrangement toward a globally optimized pattern after their initial emergence.

An illustrative example is provided by the emergence of pinwheel-sparse orientation stripes. Qualitatively, it is easy to see that a pattern of orientation stripes satisfies the continuity constraint very well. In a stripe pattern, preferred orientations are constant along one direction in space, realizing the absolute minimum of the orientation gradient in this direction. Reaching such a configuration obviously requires to select the preferred orientation at widely separated sites (along the stripe axis) to be identical. Because initially such sites develop independent preferred orientations, the optimized column layout can only emerge through a secondary rearrangement process. If the dominant low energy state has low pinwheel density, the later phase is governed by pinwheel motion and pairwise pinwheel annihilation. If this state is pinwheel-rich, e.g., a pinwheel crystal or an aperiodic pinwheel-rich state, both pinwheel annihilation and pinwheel creation together with a coordinated rearrangement of pinwheels are expected to occur.

The local, essentially random processes during the initial emergence of a fist pattern are in principle incapable of directly generating an optimized layout. In fact, it has been established that this initial so-called symmetry breaking phase will in general produce a random arrangement of selectivities of model-insensitive statistics [25, 32, 36]. The occurrence of some form of secondary reorganization is thus a qualitative prediction of any optimization model, provided that the optimal map is not seeded by an innate mechanism. The results presented in this study and many reports demonstrate that Hebbian plasticity is capable and often expected to achieve such rearrangements.

In gradient descent dynamics simulations of the EN model for retinotopy and stimulus orientation with conventional stimulus ensembles, pinwheel densities were found to be strongly time-dependent after the initial column formation (see e.g. Figures 6, 7). In particular, the timescale for the establishment of full orientation selectivity and the time needed for either annihilation of a substantial fraction of pinwheels or their crystallization into periodic pinwheel crystals are in the same range of tens of tau. A similar time dependence of pinwheel density has also been observed in other models for OPM development with periodic optima [35, 38]. Pinwheel annihilation in the EN can be slightly slowed down by additional features such as retinotopy Figures 10 and 17 or OD [25, 37] but not by orders of magnitude. For this reason, signatures of the periodic optima of a developmental dynamics become visible at rather early simulation stages. Long-term minimization is apparently not essential to express the main layout features of the global minimum.

Because the main features of the dominant optimal solutions become apparent immediately after orientation selectivity saturates it appears not easy to reproduce the species-independent map layout in models with periodic crystalline optima by pattern freezing. In our simulations to match even only the pinwheel density, a very precise timing of the freezing point would be required. There is currently no evidence for such a freezing mechanism in early development. In cats and ferrets, cortical maps for OD, orientation or direction arise on a timescale between hours and a few days (e.g., [13, 119, 120]). The underlying circuits can be rapidly modified, e.g. by deprivation experiments, even on the timescale of hours [121, 122] weeks after full selectivity has been established. Recently, evidence for long-term visual cortical circuit reorganization after the emergence of feature selectivity during normal development has emerged in diverse systems. In mouse, for example, activity-dependent changes induced by normal visual experience during the critical period, i.e., long after the primary emergence of orientation selectivity, have been shown to gradually match eye-specific inputs in the cortex [123]. Specifically, the data from mouse indicates that preferred orientations in the two eyes initially often emerge unmatched and subsequently change toward one binocularly matched orientation preference. Because preferred orientations in the two eyes initially are statistically independent, this suggests that neurons can rotate their orientation preferences up to at least 45° during postnatal development. This is reminiscent of pairing experiments in kitten visual cortex in which Frégnac and coworkers induced neurons to changed their preferred orientation by up to 90° after pairing of a visual stimulus with intracortical stimulation [124, 125] (see also [126]). Also in the cat, visual cortical orientation columns in visual areas V1 and V2 have been found to undergo rearrangement during the late phase of the critical period [41]. In this process, columns in mutually connected regions of areas V1 and V2 or in retinotopically matched regions in left and right hemisphere areas become progressively better matched in size. In the same species, a systematic reorganization of OD columns during postnatal development has been observed [69]. Essential features of this columnar rearrangement are reproduced by the EN model for OD patterns simulated in a growing domain.

In view of these observations, it seems unlikely that aperiodic orientation maps in the visual cortex represent frozen transient states of a developmental dynamics whose attracting layouts are pinwheel crystals or pinwheel free states. In fact, models for the activity-dependent development of OPMs with aperiodic optima predict only subtle changes of the OPM layout during the convergence after the establishment of selectivity [35, 38]. This might also explain the apparent stability of cortical maps during normal development over short periods [119]. Further studies of the long-term rearrangement and stabilization of cortical functional architecture are needed to exhaustively characterize such processes. Given the fundamental role of map reorganization for any optimization theory of visual cortical development, chronic imaging experiments tracking the spatial arrangement of feature selectivities in individual animals beyond the emergence of selectivity and through later developmental stages are expected to be highly informative about fundamental principles of visual cortical optimization.

Conclusions

Together with recent progress on the quantitative characterization of cortical functional architecture [38, 69, 93], this study lays the foundation for a mathematically rigorous and biologically informative search for optimization principles that successfully explain the architecture of columnar contour representations in the primary visual cortex. A mathematically controlled and quantitatively precise determination of the predictions of candidate optimization principles is demanded by accumulating evidence indicating that geometrical features of visual cortical representations are biologically laid down with a precision in the range of a few percent [38, 127, 128]. Such data is expected to substantially reduce the range of candidate optimization principles that are consistent with biological observations. In particular, for the principle that cortical orientation maps are designed to optimally compromise stimulus coverage and feature continuity, our analysis demonstrates that the classical EN model for orientation preference and retinotopy essentially fails at explaining the biologically observed architecture. Our finding that the EN model exhibits biologically realistic optima only in a limit in which point-like stimuli are represented by complex spatially extended activity patterns corroborates that large-scale interactions are essential for the stabilization of OPM layouts with realistic geometry [35, 39, 87, 93]. In the light of these results, principles for the optimal representation of entire visual scenes by extended cortical activity patterns appear as promising candidates for future studies (see also [54]). In fact, there is recent evidence that visual cortical activity becomes progressively better matched to the statistics of natural stimuli but not to simplistic artificial stimulus ensembles [129]. We expect the methods developed here to facilitate a comprehensive characterization of such candidate principles.

Methods

Expansion of EN equation

In order to analytically calculate the approximate optimal dimension-reducing mappings in the EN model with fixed retinotopy, an expansion of the nonlinear EN OPM dynamics (Equation (3)) up to third-order around the unselective fixed point has to be derived. This expansion is briefly sketched in the following. Equation (3) with r(x) = 0 is of the form

$${\partial }_{t}z\left(\mathbf{x},t\right)={\mathcal{N}}_{\mathbf{x}}\left[z\right]+\eta \Delta z\left(\mathbf{x},t\right),$$

where ${\mathcal{N}}_{\mathbf{x}}\left[z\right]$ is a nonlinear functional of z(·, t), parameterized by the position x. Clearly, the diffusion term contains no nonlinear terms in z(·, t) and therefore third order terms of the dynamics ∂ _t z(x, t) exclusively stem from third order terms of the Volterra series expansion of the functional ${\mathcal{N}}_{\mathbf{x}}\left[z\right]$ around the fixed point z(x, t)≡ 0. By the Shift symmetry (Equation (8)), only third-order contributions of the form $N_3\left[z,z,\bar{z}\right]$ are allowed, i.e.,

$$N_3\left[z,z,\bar{z}\right]=\frac{1}{2}\iiint d^{2}y{d}^{2}w{d}^{2}v{\left.\frac{{\delta }^3{\mathcal{N}}_{\mathbf{x}}\left[z\right]}{\delta z\left(\mathbf{y}\right)\delta z\left(\mathbf{w}\right)\delta \bar{z}\left(\mathbf{v}\right)}\right|}_{z\equiv 0}z\left(\mathbf{y}\right)z\left(\mathbf{w}\right)\bar{z}\left(\mathbf{v}\right).$$

Collecting all the terms yields

$$N_3\left[z,z,\bar{z}\right]=\sum _{j=1}^{11}{a}_j{N}_3^j\left[z\right],$$

(37)

where

$$\begin{array}{ll}\hfill N_3^1\left[z\right]& ={\left|z\left(\mathbf{x}\right)\right|}^{2}z\left(\mathbf{x}\right)\\ \hfill N_3^2\left[z\right]& ={\left|z\left(\mathbf{x}\right)\right|}^2\int d^{2}y{K}_2\left(\mathbf{y}-\mathbf{x}\right)z\left(\mathbf{y}\right)\\ \hfill N_3^3\left[z\right]& =z{\left(\mathbf{x}\right)}^2\int d^{2}y{K}_2\left(\mathbf{y}-\mathbf{x}\right)\bar{z}\left(\mathbf{y}\right)\\ \hfill N_3^4\left[z\right]& =z\left(\mathbf{x}\right)\int d^{2}y{K}_2\left(\mathbf{y}-\mathbf{x}\right){\left|z\left(\mathbf{y}\right)\right|}^2\\ \hfill N_3^5\left[z\right]& =\bar{z}\left(\mathbf{x}\right)\int d^{2}y{K}_2\left(\mathbf{y}-\mathbf{x}\right)z{\left(\mathbf{y}\right)}^2\\ \hfill N_3^6\left[z\right]& =\int d^{2}y{K}_2\left(\mathbf{y}-\mathbf{x}\right){\left|z\left(\mathbf{y}\right)\right|}^{2}z\left(\mathbf{y}\right)\\ \hfill N_3^7\left[z\right]& =z\left(\mathbf{x}\right)\iint d^{2}y{d}^{2}w{K}_3\left(\mathbf{y}-\mathbf{x},\mathbf{w}-\mathbf{x},\mathbf{y}-\mathbf{w}\right)\bar{z}\left(\mathbf{w}\right)z\left(\mathbf{y}\right)\\ \hfill N_3^8\left[z\right]& =\bar{z}\left(\mathbf{x}\right)\iint d^{2}y{d}^{2}w{K}_3\left(\mathbf{y}-\mathbf{x},\mathbf{w}-\mathbf{x},\mathbf{y}-\mathbf{w}\right)z\left(\mathbf{w}\right)z\left(\mathbf{y}\right)\\ \hfill N_3^9\left[z\right]& =\iiint d^{2}y{d}^{2}w{d}^{2}v{K}_4\left(\mathbf{y}-\mathbf{x},\mathbf{w}-\mathbf{x},\mathbf{v}-\mathbf{x},\mathbf{y}-\mathbf{w},\mathbf{v}-\mathbf{w},\mathbf{y}-\mathbf{v}\right)\bar{z}\left(\mathbf{v}\right)z\left(\mathbf{w}\right)z\left(\mathbf{y}\right)\\ \hfill N_3^{10}\left[z\right]& =\iint d^{2}y{d}^{2}w{K}_3\left(\mathbf{y}-\mathbf{x},\mathbf{w}-\mathbf{x},\mathbf{y}-\mathbf{w}\right){\left|z\left(\mathbf{w}\right)\right|}^{2}z\left(\mathbf{y}\right)\\ \hfill N_3^{11}\left[z\right]& =\iint d^{2}y{d}^{2}w{K}_3\left(\mathbf{y}-\mathbf{x},\mathbf{w}-\mathbf{x},\mathbf{y}-\mathbf{w}\right)z{\left(\mathbf{w}\right)}^2\bar{z}\left(\mathbf{y}\right)\end{array}$$

(38)

and

$$\begin{array}{c}{K}_2\left(\mathbf{x}\right)=e^{-{\mathbf{x}}^2∕\left(4{\sigma }^2\right)}\\ K_3\left({\mathbf{x}}_1,{\mathbf{x}}_2,{\mathbf{x}}_3\right)=e^{-\left({\mathbf{x}}_1^2+{\mathbf{x}}_2^2+{\mathbf{x}}_3^2\right)∕\left(6{\sigma }^2\right)}\\ K_4\left({\mathbf{x}}_1,{\mathbf{x}}_2,{\mathbf{x}}_3,{\mathbf{x}}_4,{\mathbf{x}}_5,{\mathbf{x}}_6\right)=e^{-\left({\mathbf{x}}_1^2+{\mathbf{x}}_2^2+{\mathbf{x}}_3^2+{\mathbf{x}}_4^2+{\mathbf{x}}_5^2+{\mathbf{x}}_6^2\right)∕\left(8{\sigma }^2\right)}.\end{array}$$

The coefficients a _j for various orientation stimulus ensembles are given in 'Results' section.

Adiabatic elimination of r(x, t)

In order to analytically calculate the approximate optimal dimension-reducing mappings in the EN model with variable retinotopy, an expansion of the nonlinear EN retinotopy and orientation map dynamics (Equations (3, 4)) up to third-order around the nonselective fixed point has to be derived and retinotopic distortions have to be adiabatically eliminated. Both of these calculations are briefly sketched in the following. Equation (3) is of the form

$${\partial }_{t}z\left(\mathbf{x},t\right)={\mathcal{N}}_{\mathbf{x}}\left[z,\mathbf{r}\right]+\eta \Delta z\left(\mathbf{x},t\right),$$

where ${\mathcal{N}}_{\mathbf{x}}\left[z,\mathbf{r}\right]$ is a nonlinear functional of z(·, t) and r(·, t), parameterized by the position x. The diffusion term contains no nonlinear terms in z(·, t) and therefore third order terms of the dynamics of z(x, t) exclusively stem from third-order terms of the Volterra series expansion of the functional ${\mathcal{N}}_{\mathbf{x}}\left[z,\mathbf{r}\right]$ around the fixed point {z(x, t)≡ 0, r(x, t)≡ 0}. By the shift symmetry (Equation (8)), only terms in form of a cubic operator $N_3\left[z,z,\bar{z}\right]$ and a quadratic operator Q^z [r, z] are allowed when expanding up to third order. $N_3\left[z,z,\bar{z}\right]$ is given in Equation (37). Q[z, r] can be calculated via

$$Q^z\left[\mathbf{r},z\right]=\iint d^{2}y{d}^{2}w\left(\frac{{\delta }^2}{\delta z\left(\mathbf{y}\right)\delta r_1\left(\mathbf{w}\right)}{\left.{\mathcal{N}}_{\mathbf{x}}\left[z,\mathbf{r}\right]\right|}_{z\equiv 0,\mathbf{r}\equiv 0}{r}_1\left(\mathbf{w}\right)+\frac{{\delta }^2}{\delta z\left(\mathbf{y}\right)\delta r_2\left(\mathbf{w}\right)}{\left.{\mathcal{N}}_{\mathbf{x}}\left[z,\mathbf{r}\right]\right|}_{z\equiv 0,\mathbf{r}\equiv 0}{r}_2\left(\mathbf{w}\right)\right)z\left(\mathbf{y}\right)$$

and this yields

$$\begin{array}{ll}\hfill Q^z\left[\mathbf{r},z\right]& =\frac{\left(〈\mid s_z{\mid }^2〉-2{\sigma }^2\right)}{16\pi {\sigma }^6}z\left(\mathbf{x}\right)\int d^{2}y〈\mathbf{r}\left(\mathbf{y}\right),{\mathbf{K}}_2^r\left(\mathbf{y}-\mathbf{x}\right)〉\\ -\frac{〈\mid s_z{\mid }^2〉}{16\pi {\sigma }^6}\int d^{2}y〈\mathbf{r}\left(\mathbf{x}\right),z\left(\mathbf{y}\right){\mathbf{K}}_2^r\left(\mathbf{y}-\mathbf{x}\right)〉+\frac{〈\mid s_z{\mid }^2〉}{16\pi {\sigma }^6}\int d^{2}y〈\mathbf{r}\left(\mathbf{y}\right),z\left(\mathbf{y}\right){\mathbf{K}}_2^r\left(\mathbf{y}-\mathbf{x}\right)〉\\ -\frac{〈\mid s_z{\mid }^2〉}{36{\pi }^2{\sigma }^8}\iint d^{2}y{d}^{2}wz\left(\mathbf{y}\right)〈\mathbf{r}\left(\mathbf{w}\right),{\mathbf{K}}_3^r\left(\mathbf{y}-\mathbf{x},\mathbf{w}-\mathbf{x},\mathbf{y}-\mathbf{w}\right)〉,\end{array}$$

where 〈·, ·〉 denotes the scalar product between two vectors and

$${\mathbf{K}}_2^r\left(\mathbf{x}\right)=e^{-{\mathbf{x}}^2∕4{\sigma }^2}\mathbf{x}$$

(39)

$${\mathbf{K}}_3^r\left({\mathbf{x}}_1,{\mathbf{x}}_2,{\mathbf{x}}_3\right)=e^{-\frac{{\mathbf{x}}_1^2+{\mathbf{x}}_2^2+{\mathbf{x}}_3^2}{6{\sigma }^2}}\left[{\mathbf{x}}_1+{\mathbf{x}}_3\right].$$

(40)

In complete analogy, by expanding the right hand side of the dynamical equation for the retinotopic distortions (Equation (4)) up to second-order, the vector-valued quadratic operator ${\mathbf{Q}}^r\left[z,\bar{z}\right]$ can be obtained as

$$\begin{array}{ll}\hfill {\mathbf{Q}}^r\left[z,\bar{z}\right]& =-\Re \left(\frac{〈\mid s_z{\mid }^2〉}{16\pi {\sigma }^6}\bar{z}\left(\mathbf{x}\right)\int d^{2}y{\mathbf{K}}_2^r\left(\mathbf{y}-\mathbf{x}\right)z\left(\mathbf{y}\right)\right)\\ +\frac{2{\sigma }^2-〈\mid s_z{\mid }^2〉}{32\pi {\sigma }^6}\int d^{2}y{\mathbf{K}}_2^r\left(\mathbf{y}-\mathbf{x}\right)\mid z\left(\mathbf{y}\right){\mid }^2\\ +\frac{〈\mid s_z{\mid }^2〉}{72{\pi }^2{\sigma }^8}\iint d^{2}y{d}^{2}w{\mathbf{K}}_3^r\left(\mathbf{y}-\mathbf{x},\mathbf{y}-\mathbf{w},\mathbf{w}-\mathbf{x}\right)z\left(\mathbf{w}\right)\bar{z}\left(\mathbf{y}\right).\end{array}$$

(41)

Inserting $\mathbf{r}\left(\mathbf{x}\right)=-{\mathbf{L}}_r^{-1}\left[{\mathbf{Q}}^r\left[z,\bar{z}\right]\right]$ into Q^z [r, z] and using the linearity of ${\mathbf{L}}_r^{-1}$ as well as the bilinearity of both, ${\mathbf{Q}}^r\left[z,\bar{z}\right]$ and Q^z [r, z], yields a sum $N_3^r\left[z,z,\bar{z}\right]={\sum }_{j=1}^{12}{a}_r^j{N}_r^j$, with

$$\begin{array}{ll}\hfill N_r^1& =z\left(\mathbf{x}\right)\int d^{2}y〈{\mathbf{L}}_r^{-1}\left[\Re \left(\bar{z}\left(\mathbf{y}\right)\int d^{2}w{K}_2\left(\mathbf{w}-\mathbf{y}\right)z\left(\mathbf{w}\right)\right)\right],{\mathbf{K}}_2^r\left(\mathbf{y}-\mathbf{x}\right)〉\\ \hfill N_r^2& =z\left(\mathbf{x}\right)\int d^{2}y〈{\mathbf{L}}_r^{-1}\left[\int d^{2}w{K}_2\left(\mathbf{w}-\mathbf{y}\right)\mid z\left(\mathbf{w}\right){\mid }^2\right],{\mathbf{K}}_2^r\left(\mathbf{y}-\mathbf{x}\right)〉\\ \hfill N_r^3& =z\left(\mathbf{x}\right)\int d^{2}y〈{\mathbf{L}}_r^{-1}\left[\iint d^{2}w{d}^{2}v{\mathbf{K}}_3^r\left(\mathbf{w}-\mathbf{y},\mathbf{w}-\mathbf{v},\mathbf{v}-\mathbf{y}\right)z\left(\mathbf{w}\right)\bar{z}\left(\mathbf{v}\right)\right],{\mathbf{K}}_2^r\left(\mathbf{y}-\mathbf{x}\right)〉\\ \hfill N_r^4& =〈{\mathbf{L}}_r^{-1}\left[\Re \left(\bar{z}\left(\mathbf{x}\right)\int d^{2}y{\mathbf{K}}_2^r\left(\mathbf{y}-\mathbf{x}\right)z\left(\mathbf{y}\right)\right)\right],\int d^{2}yz\left(\mathbf{y}\right){\mathbf{K}}_2^r\left(\mathbf{y}-\mathbf{x}\right)〉\\ \hfill N_r^5& =〈{\mathbf{L}}_r^{-1}\left[\int d^{2}y{\mathbf{K}}_2^r\left(\mathbf{y}-\mathbf{x}\right)\mid z\left(\mathbf{y}\right){\mid }^2\right],\int d^{2}yz\left(\mathbf{y}\right){\mathbf{K}}_2^r\left(\mathbf{y}-\mathbf{x}\right)〉\\ \hfill N_r^6& =〈{\mathbf{L}}_r^{-1}\left[\iint d^{2}y{d}^{2}w{\mathbf{K}}_3^r\left(\mathbf{y}-\mathbf{x},\mathbf{y}-\mathbf{w},\mathbf{w}-\mathbf{x}\right)z\left(\mathbf{w}\right)\bar{z}\left(\mathbf{y}\right)\right],\int d^{2}yz\left(\mathbf{y}\right){\mathbf{K}}_2^r\left(\mathbf{y}-\mathbf{x}\right)〉\\ \hfill N_r^7& =\int d^{2}yz\left(\mathbf{y}\right)〈{\mathbf{L}}_r^{-1}\left[\Re \left(\bar{z}\left(\mathbf{y}\right)\int d^{2}w{K}_2\left(\mathbf{w}-\mathbf{y}\right)z\left(\mathbf{w}\right)\right)\right],{\mathbf{K}}_2^r\left(\mathbf{y}-\mathbf{x}\right)〉\\ \hfill N_r^8& =\int d^{2}yz\left(\mathbf{y}\right)〈{\mathbf{L}}_r^{-1}\left[\int d^{2}w{\mathbf{K}}_2^r\left(\mathbf{w}-\mathbf{y}\right)\mid z\left(\mathbf{w}\right){\mid }^2\right],{\mathbf{K}}_2^r\left(\mathbf{y}-\mathbf{x}\right)〉\\ \hfill N_r^9& =\int d^{2}yz\left(\mathbf{y}\right)〈{\mathbf{L}}_r^{-1}\left[\iint d^{2}v{d}^{2}w{\mathbf{K}}_3^r\left(\mathbf{w}-\mathbf{y},\mathbf{w}-\mathbf{v},\mathbf{v}-\mathbf{y}\right)z\left(\mathbf{w}\right)\bar{z}\left(\mathbf{y}\right)\right],{\mathbf{K}}_2^r\left(\mathbf{y}-\mathbf{x}\right)〉\\ \hfill N_r^{10}& =\iint d^{2}y{d}^{2}wz\left(\mathbf{y}\right)〈{\mathbf{L}}_r^{-1}\left[\Re \left(\bar{z}\left(\mathbf{w}\right)\int d^{2}v{\mathbf{K}}_2^r\left(\mathbf{v}-\mathbf{w}\right)z\left(\mathbf{v}\right)\right)\right],{\mathbf{K}}_3^r\left(\mathbf{y}-\mathbf{x},\mathbf{y}-\mathbf{w},\mathbf{w}-\mathbf{x}\right)〉\\ \hfill N_r^{11}& =\iint d^{2}y{d}^{2}wz\left(\mathbf{y}\right)〈{\mathbf{L}}_r^{-1}\left[\int d^{2}v{\mathbf{K}}_2^r\left(\mathbf{v}-\mathbf{w}\right)\mid z\left(\mathbf{v}\right){\mid }^2\right],{\mathbf{K}}_3^r\left(\mathbf{y}-\mathbf{x},\mathbf{y}-\mathbf{w},\mathbf{w}-\mathbf{x}\right)〉\\ \hfill N_r^{12}& =\iint d^{2}y{d}^{2}wz\left(\mathbf{y}\right)〈{\mathbf{L}}_r^{-1}\left[\iint d^{2}v{d}^{2}u{\mathbf{K}}_3^r\left(\mathbf{v}-\mathbf{w},\mathbf{v}-\mathbf{u},\mathbf{u}-\mathbf{w}\right)z\left(\mathbf{v}\right)\bar{z}\left(\mathbf{u}\right)\right],{\mathbf{K}}_3^r\left(\mathbf{y}-\mathbf{x},\mathbf{y}-\mathbf{w},\mathbf{w}-\mathbf{y}\right)〉.\end{array}$$

The coefficients $a_r^j$ are given by

$$\begin{array}{ll}\hfill a_r^1& =-\frac{\left(〈\mid s_z{\mid }^2〉-2{\sigma }^2\right)〈\mid s_z{\mid }^2〉}{256{\pi }^2{\sigma }^{12}}=-\frac{\left(1-{\sigma }^2\right)}{64{\pi }^2{\sigma }^{12}}\\ \hfill a_r^2& =\frac{{\left(〈\mid s_z{\mid }^2〉-2{\sigma }^2\right)}^2}{512{\pi }^2{\sigma }^{12}}=\frac{{\left(1-{\sigma }^2\right)}^2}{128{\pi }^2{\sigma }^{12}}\\ \hfill a_r^3& =\frac{\left(〈\mid s_z{\mid }^2〉-2{\sigma }^2\right)〈\mid s_z{\mid }^2〉}{1152{\pi }^3{\sigma }^{14}}=\frac{\left(1-{\sigma }^2\right)}{288{\pi }^3{\sigma }^{14}}\\ \hfill a_r^4& =-\frac{{〈\mid s_z{\mid }^2〉}^2}{256{\pi }^2{\sigma }^{12}}=-\frac{1}{64{\pi }^2{\sigma }^{12}}\\ \hfill a_r^5& =-\frac{\left(〈\mid s_z{\mid }^2〉-2{\sigma }^2\right)〈\mid s_z{\mid }^2〉}{512{\pi }^2{\sigma }^{12}}=-\frac{\left(1-{\sigma }^2\right)}{128{\pi }^2{\sigma }^{12}}\\ \hfill a_r^6& =\frac{{〈\mid s_z{\mid }^2〉}^2}{1152{\pi }^3{\sigma }^{14}}=\frac{1}{288{\pi }^3{\sigma }^{14}}\\ \hfill a_r^7& =-\frac{{〈\mid s_z{\mid }^2〉}^2}{256{\pi }^2{\sigma }^{12}}=-\frac{1}{64{\pi }^2{\sigma }^{12}}\\ \hfill a_r^8& =-\frac{\left(〈\mid s_z{\mid }^2〉-2{\sigma }^2\right)〈\mid s_z{\mid }^2〉}{512{\pi }^2{\sigma }^{12}}=-\frac{\left(1-{\sigma }^2\right)}{128{\pi }^2{\sigma }^{12}}\\ \hfill a_r^9& =\frac{{〈\mid s_z{\mid }^2〉}^2}{1152{\pi }^3{\sigma }^{14}}=\frac{1}{288{\pi }^3{\sigma }^{14}}\\ \hfill a_r^{10}& =\frac{{〈\mid s_z{\mid }^2〉}^2}{576{\pi }^3{\sigma }^{14}}=\frac{1}{144{\pi }^3{\sigma }^{14}}\\ \hfill a_r^{11}& =\frac{\left(〈\mid s_z{\mid }^2〉-2{\sigma }^2\right)〈\mid s_z{\mid }^2〉}{1152{\pi }^3{\sigma }^{14}}=\frac{\left(1-{\sigma }^2\right)}{288{\pi }^3{\sigma }^{14}}\\ \hfill a_r^{12}& =-\frac{{〈\mid s_z{\mid }^2〉}^2}{2592{\pi }^4{\sigma }^{16}}=-\frac{1}{648{\pi }^4{\sigma }^{16}},\end{array}$$

where the second equal sign is valid for 〈|s _z |²〉 = 2.

Amplitude equations from $N_3^z\left[z,z,\bar{z}\right]$

We catalog the numerous stationary solutions of Equation (13) following [35], by considering planforms

$$z\left(\mathbf{x},t\right)=\sum _{j=1}^N{A}_j\left(t\right)e^{i{\mathbf{k}}_j\mathbf{x}}$$

with an even number N of modes with amplitudes A _j and k _j = k _c (cos(2πj/N), sin(2πj/N )). In the vicinity of a finite wavelength instability, where the nonselective state z(x) = 0 becomes unstable with respect to a band of Fourier modes around a finite wave number k _c -by symmetry, the dynamics of the amplitudes A _j at threshold has the form

$${Ȧ}_i=A_i-A_i\sum _{j=1}^N{g}_{ij}\mid A_j{\mid }^2-\sum _{j=1}^N{f}_{ij}{A}_j{A}_{j^{-}}{Ā}_{i^{-}},$$

(42)

where j^- denotes the index of the mode antiparallel to the mode j, ${\mathbf{k}}_j=-{\mathbf{k}}_{j^{-}}$, and the coefficients $g_{ij}=\left(1-\frac{1}{2}{\delta }_{ij}\right)g\left(\mid {\alpha }_i-{\alpha }_j\mid \right)$ and $f_{ij}=\left(1-{\delta }_{ij}-{\delta }_{i^{-}j}\right)f\left(\mid {\alpha }_i-{\alpha }_j\mid \right)$ only depend on the angle |α _i - α _j | between mode i and j. The angle-dependent interaction functions g(α) and f(α) are obtained from Equation (13) by a multi scale expansion [35, 83, 84, 87] as

$$\begin{array}{c}g\left(\alpha \right)=-e^{-i{\mathbf{k}}_0\mathbf{x}}\left[N_3^z\left(e^{i{\mathbf{k}}_0\mathbf{x}},e^{i\mathbf{h}\left(\alpha \right)\mathbf{x}},e^{-i\mathbf{h}\left(\alpha \right)\mathbf{x}}\right)\right.\\ +\left.N_3^z\left(e^{i\mathbf{h}\left(\alpha \right)\mathbf{x}},e^{-i\mathbf{h}\left(\alpha \right)\mathbf{x}},e^{i{\mathbf{k}}_0\mathbf{x}}\right)\right]\end{array}$$

(43)

$$\begin{array}{c}f\left(\alpha \right)=-\frac{1}{2}{e}^{-i{\mathbf{k}}_0\mathbf{x}}\left[N_3^z\left(e^{i\mathbf{h}\left(\alpha \right)\mathbf{x}},e^{-i\mathbf{h}\left(\alpha \right)\mathbf{x}},e^{i{\mathbf{k}}_0\mathbf{x}}\right)\right.\\ +\left.N_3^z\left(,e^{-i\mathbf{h}\left(\alpha \right)\mathbf{x}},e^{i\mathbf{h}\left(\alpha \right)\mathbf{x}},e^{i{\mathbf{k}}_0\mathbf{x}}\right)\right],\end{array}$$

(44)

where k₀ = k _c (1, 0) and h(α) = k _c (cos α, sin α). f(α) is π-periodic, since the right hand side of Equation (44) is invariant with respect to the transformation h(α) → h(α + π) = -h(α). g(α) is 2π-periodic in general. If, however, the nonlinearity is permutation symmetric (Equation (34)) it can be seen from Equation (43) that g(α) is π-periodic as well.

Stability of stationary planform solutions

To determine the stability of fixed points of the amplitude equations (42), the eigenvalues of their stability matrices have to be determined. In general, for any fixed point A = A⁰ of the dynamical system $\dot{\mathbf{A}}=\mathbf{F}\left(\mathbf{A}\right)$ with complex-valued A and F, we have to compute the eigenvalues of the Hermitian 2N × 2N matrix

$$M={\left(\begin{array}{cc}\frac{\partial F}{\partial A}& \frac{\partial F}{\partial \overline{A}}\\ \frac{\partial \overline{F}}{\partial A}& \frac{\partial F}{\partial A}\end{array}\right)|}_{A=A_0}.$$

For the system of amplitude equations, we obtain

$$\begin{array}{ll}\hfill \frac{\partial F_i}{\partial A_k}& =r{\delta }_{ik}-{\delta }_{ik}\left(\sum _j^N{g}_{ij}\mid A_j{\mid }^2\right)-A_i{g}_{ik}{Ā}_k-{Ā}_{i^{-}}{f}_{ik}\left(A_{k^{-}}+A_k\right)\\ \hfill \frac{\partial F_i}{\partial {Ā}_k}& =-A_i{g}_{ik}{A}_k-{\delta }_{i^{-}k}\left(\sum _j^N{f}_{ij}{A}_j{A}_{j^{-}}\right).\end{array}$$

Stability of a solution, or more precise intrinsic stability is given, if all eigenvalues of M are negative definite. Extrinsic stability is given, if the growth of additional Fourier modes is suppressed. To test whether a planform solution is extrinsically stable, we introduce a test mode T such that

$$z\left(\mathbf{x}\right)=T{e}^{i{\mathbf{k}}_{\beta }\mathbf{x}}+\sum _j^N{A}_j{e}^{i{\mathbf{k}}_j\mathbf{x}},$$

with k _β = (cos β, sin β)k _c . We insert this ansatz into Equation 15 and obtain

$${\partial }_{t}T=rT-\sum _j^{N}g\left(\beta -{\beta }_j\right)\mathsf{\text{∣}}{A}_j{\mathsf{\text{∣}}}^{2}T+\mathcal{O}\left(T^2\right)$$

as the dynamics of the test mode T, where g(β) is the angle-dependent interaction function corresponding to $N_3\left[z,z,\bar{z}\right]$. For the solution T = 0 to be stable, we therefore obtain the condition

$$r-\sum _j^{N}g\left(\beta -{\beta }_j\right)\mathsf{\text{∣}}{A}_j{\mathsf{\text{∣}}}^2<0,\forall \alpha \in \left[0,2\pi \right],$$

where we assumed ${\mathbf{k}}_{\beta }\ne {\mathbf{k}}_j,{\mathbf{k}}_j^{-}$. These conditions for intrinsic and extrinsic stability were numerically evaluated to study the stability of n-ECPs and rPWCs.

Coupled essentially complex planforms

In 'Results' section, we presented a closed form expression for the retinotopic distortions associated via Equation (28) with stationary planform solutions of Equation (29). Here, we sketch how to explicitly calculate this representation. We start with the ansatz

$$z\left(\mathbf{x}\right)=\sum _j^n{A}_j{e}^{i{\mathbf{k}}_j\mathbf{x}}\mid {\mathbf{k}}_j\mid =k_c$$

(45)

for the OPM z(x). Note that this general ansatz includes essentially complex planforms as well as rPWCs. To simplify notation, we denote the individual terms in Equation (41)

$$\begin{array}{ll}\hfill {\mathbf{Q}}_1\left[z,\bar{z}\right]& =-\Re \left(\frac{〈\mid s_z{\mid }^2〉}{16\pi {\sigma }^6}\bar{z}\left(\mathbf{x}\right)\int d^{2}y{\mathbf{K}}_2^r\left(\mathbf{y}-\mathbf{x}\right)z\left(\mathbf{y}\right)\right)\\ \hfill {\mathbf{Q}}_2\left[z,\bar{z}\right]& =\frac{2{\sigma }^2-〈\mid s_z{\mid }^2〉}{32\pi {\sigma }^6}\int d^{2}y{\mathbf{K}}_2^r\left(\mathbf{y}-\mathbf{x}\right)\mid z\left(\mathbf{y}\right){\mid }^2\\ \hfill {\mathbf{Q}}_3\left[z,\bar{z}\right]& =\frac{〈\mid s_z{\mid }^2〉}{72{\pi }^2{\sigma }^8}\iint d^{2}y{d}^{2}w{\mathbf{K}}_3^r\left(\mathbf{y}-\mathbf{x},\mathbf{y}-\mathbf{w},\mathbf{w}-\mathbf{x}\right)z\left(\mathbf{w}\right)\bar{z}\left(\mathbf{y}\right).\end{array}$$

Each of the ${\mathbf{Q}}_i\left[z,\bar{z}\right]$, i = 1, 2, 3 can be evaluated for the planform ansatz in Equation (45) and we obtain

$$\begin{array}{c}{Q}_1\left[{\displaystyle \sum _j^n{A}_j}{e}^{i{k}_{j}x},{\displaystyle \sum _k^n{\overline{A}}_k}{e}^{-i{k}_{k}x}\right]=\frac{〈\mid s_z{\mid }^2〉e^{-k_c^2{\sigma }^2}}{2{\sigma }^2}{\displaystyle \sum _{k,j<k}^n\{\Re (A_j{\overline{A}}_k)(k_j-k_k)\sin ((k_j-k_k)x)}\\ +\Im (A_j{\overline{A}}_k)(k_j-k_k)\cos ((k_j-k_k)x)\}\\ Q_2\left[{\displaystyle \sum _j^n{A}_j}{e}^{i{k}_{j}x},{\displaystyle \sum _k^n{\overline{A}}_k}{e}^{-i{k}_{k}x}\right]=-\frac{(2{\sigma }^2-〈\mid s_z{\mid }^2〉)}{2{\sigma }^2}{\displaystyle \sum _{j<k}^n{e}^{-{(k_j-k_k)}^2{\sigma }^2}}(k_j-k_k)\{\Re (A_j{\overline{A}}_k)\sin ((k_j-k_k)x)\\ +\Im (A_j{\overline{A}}_k)\cos ((k_j-k_k)x)\}\\ Q_3\left[{\displaystyle \sum _j^n{A}_j}{e}^{i{k}_{j}x},{\displaystyle \sum _k^n{\overline{A}}_k}{e}^{-i{k}_{k}x}\right]=-\frac{〈\mid s_z{\mid }^2〉e^{-\frac{k_c^2{\sigma }^2}{2}}}{{\sigma }^2}{\displaystyle \sum _{j<k}^n(k_j-k_k)e^{-\frac{{(k_j-k_k)}^2{\sigma }^2}{2}}}\{\Re (A_j{\overline{A}}_k)\sin ((k_j-k_k)x)\\ +\Im (A_j{\overline{A}}_k)\cos ((k_j-k_k)x)\}.\end{array}$$

All resulting terms are proportional to either (k _i - k _j ) sin((k _i - k _j )x) or (k _i - k _j ) cos((k _i - k _j )x) i ≠ j. These functions are longitudinal modes (see Figure 3b) which have been identified as eigenfunctions of the linearized dynamics of retinotopic distortions L _r [r] with eigenvalue

$${\lambda }_L^r\left(\mid {\mathbf{k}}_i-{\mathbf{k}}_j\mid \right)=-\mid {\mathbf{k}}_i-{\mathbf{k}}_j{\mid }^2\left({\eta }_r+e^{-{\sigma }^2\mid {\mathbf{k}}_i-{\mathbf{k}}_j{\mid }^2}{\sigma }^2\right)$$

Hence, they are eigenfunctions of ${\mathbf{L}}_r^{-1}\left[\mathbf{r}\right]$ with eigenvalue $1/{\lambda }_L^r\left(\mid {\mathbf{k}}_i-{\mathbf{k}}_j\mid \right)$. Using this when inserting in Equation (28) and setting 〈|s _z |²〉 = 2, we obtain expression (31) for the retinotopic distortions belonging to an arbitrary planform.

Phase diagrams

To compute the regions of minimal energy shown in Figures 6, 10, 12, 15, and 16 as well as Figures 23, 25 in Appendix 1, we first computed the fixed points of Equation (42) at each point in parameter space. For n-ECPs, we constructed the coupling matrix g in Equation (22) for all mode configurations not related by any combination of the symmetry operations: (i) Translation: $A_j\to A_j{e}^{i{\mathbf{k}}_j\mathbf{y}}$, (ii) Rotation: A _j → A _{j +Δ j} , (iii) Parity: $A_j\to {Ā}_{{\left(N-j\right)}^{-}}$. By Equation (22), we then computed the absolute values of the corresponding amplitudes. If ${\sum }_{j=1}^n{\left({\mathbf{g}}^{-1}\right)}_{ij}\ge 0$ for all i, a valid n-ECP fixed point of Equation 42 was identified. Its energy was then computed by Equation (23). For orientation stripes and rPWCs, the derived analytical expressions for their energy (Equations (18, 20)) were evaluated. To analyze the stability of the fixed points, the conditions for intrinsic and extrinsic stability (see above) were numerically evaluated.

Numerical procedures-gradient descent simulations

To test our analytical calculations and explore their range of validity, we simulated Equations (3) and (4) on a 64 × 64 grid with periodic boundary conditions. Simulated systems were spatially discretized with typically 8 grid points per expected column spacing Λ_max of the orientation preference pattern (see 'Results' section) to achieve sufficient resolution in space. Test simulations with finer discretization (16 and 32 grid points per Λ_max) did not lead to different results. Progression of time was measured in units of the intrinsic timescale τ (see 'Results' section) of the pattern formation process. The integration time step δt was bounded by the relevant decay time constant of the Laplacian in Equation (3) around k _c and by the intrinsic timescale τ of the system, using $\delta t=min\left\{1/\left(20\eta \underset{c}{\overset{2}{k}}\right),\tau /10\right\}$. This ensured good approximation to the temporally continuous changes of the patterns. We used an Adams-Bashforth scheme for the first terms on the respective r.h.s. of Equations (3, 4). The second terms (diffusion) were treated by spectral integration, exhibiting unconditional numerical stability. The stimulus positions s _r were chosen to be uniformly distributed in retinal coordinates. The stimulus averages in Equations (3, 4) were approximated by choosing a random representative sample of N _s stimuli at each integration time step, with

$$N_s=max\left\{{10}^5,\frac{N_0{\Gamma }^2}{{\left({\epsilon }_s\right)}^n}\frac{\delta t}{\tau }\right\},$$

where n corresponds to the dimensions of the feature space in addition to the two retinal positions (in our case, n = 2), Γ² = (L/Λ_max)² the squared aspect ratio of the simulated system in units of Λ², ε _s the resolution in feature space, N₀ the number of stimuli we required to approximate the cumulative effect of the ensemble of stimuli within each feature space voxel ε^{n +2} . With N₀ = 100 and ε _s = 0.05, we ensured a high signal-to-noise ratio for all the simulations. Typical values for N _s were between 2.5 × 10⁵ and 4 × 10⁶. All simulations were initialized with z(x, t = 0) = 10^-6e^{i 2πξ} (^x) and r(x, t = 0) = 0, where the ξ(x) are independent identically distributed random numbers uniformly in [0, 1]. Different realizations were obtained by using different stimulus samples.

Stimuli were drawn from different distributions, each with 〈|s _z |²〉 = 2. We considered (i) stimuli uniformly distributed on a ring with $\mid s_z{\mid }^2=\sqrt{2}$ (circular stimulus ensemble), (ii) stimuli uniformly distributed within a circle {s _z , |s _z | ≤ 2} (uniform stimulus ensemble), and (iii) a Gaussian stimuli ensemble with ${\rho }_{s_z}=1/\left(2\pi \right)exp\left(-\mid s_z{\mid }^2/2\right)$. In addition, we considered mixtures of a high-fourth moment Pearson type VII distribution and the circular stimulus ensemble. The Pearson distribution is given by

$${\rho }_{s_z}=\frac{1}{\alpha B\left(m-\frac{1}{2},\frac{1}{2}\right)}{\left[1+\frac{\mid s_z{\mid }^2}{{\alpha }^2}\right]}^{-m},$$

where B(·,·) is the Beta function [130] and $\alpha =\sqrt{2m-3}$, and $m=\frac{5}{2}+\frac{12}{\gamma -12}$ such that 〈|s _z |²〉 = 2, 〈|s _z |⁴〉 = γ or equivalently s₄ = γ - 4.

In addition to simulations in which independent sets of stimuli were used for evaluating the stimulus average in Equations (3, 4) for every time step, we also performed simulations in which the same fixed set of N stimuli was used (see 'Results' section). To determine the time step δt in these simulations, we first calculated

$$N_{\tau }=\frac{N_0{\Gamma }^2}{{\left({\epsilon }_s\right)}^n}$$

(parameters as in regular simulations) which yields the number of stimuli presented to the model in one intrinsic time unit τ in regular simulations. To subject the network to the same number of stimuli per intrinsic time scale τ in fixed stimulus set simulations, the integration time step δt was in this case chosen as

$$\delta t=min\left\{\frac{N}{N_s}\tau ,\frac{1}{20\eta \underset{c}{\overset{2}{k}}},\frac{\tau }{10}\right\}.$$

For small N, this resulted in very small integration steps. For very large N, time steps were identical to the regular simulations. Different realizations were obtained by different but fixed stimulus sets. In all simulations with fixed stimulus sets, stimuli were drawn from the circular stimulus ensemble. All other numerical methods were chosen as in regular simulations.

Numerical procedures - solving the EN model with deterministic annealing

A large body of previous study has solved the EN models for various aspects of visual cortical architecture for discrete fixed sets of stimuli and using deterministic annealing. We therefore also used deterministic annealing with fixed discrete sets of stimuli to solve the EN model for the most frequently used stimulus distribution. This allowed us to better compare our analytical and numerical results based on the gradient descent dynamics for a continuous stimulus with prior results. For the discrete deterministic annealing approach, cortical maps are described by a collection of M centroids ${\left\{{\mathbf{y}}_m\right\}}_{m=1}^M\subset {ℝ}^d$ that can be represented as a D × M matrix Y = (y₁, . . . , y _M ). Maps forming a compromise of coverage and continuity are obtained for ${\left\{{\mathbf{x}}_n\right\}}_{n=1}^N\subset {ℝ}^d$ represented as a D × N matrix X = (x₁, . . . , x _M ). In our case, d = 4. The trade-off between coverage and continuity is formulated by the energy function

$$E\left(\mathbf{Y},\sigma \right)=-\alpha \sigma \sum _{n=1}^{N}log\sum _{m=1}^M{e}^{-\frac{1}{2}{∥\frac{{\mathbf{x}}_n-{\mathbf{y}}_m}{\sigma }∥}^2}+\frac{\beta }{2}\mathsf{\text{tr}}\left({\mathbf{Y}}^T\mathbf{Y}S\right).$$

(46)

The matrix S determines the topology of the network as well as the boundary conditions and is typically derived from a discretized derivative based on a finite difference scheme or stencil (see below). For large N and M, the energy function in Equation (46) is equivalent to the energy functional of the continuum formulation given in Equation (1) for β = ηN and S implementing the discretized Laplacian operator in two dimensions.

Following [62–65, 71] we minimized the EN energy function (Equation (46)) by an iterative deterministic annealing algorithm, starting with a minimization for large σ and tracking this minimum to a small value of σ. As in [4], we reduced σ from σ_init = 0.2 to the point at which the amplitude of the orientation maps saturate (σ ≈ 0.03), following σ = σ_init × χ^j where j counts the annealing step. This choice tracks stationary solutions of the EN to parameters that are very far from threshold. For high precision tracking of solutions, we used an annealing rate of χ = 0.999.

At each value of σ, setting the gradient of Equation (46) to zero yields a nonlinear system of equations

$$\mathbf{Y}\mathbf{A}=\mathbf{X}\mathbf{W}\mathsf{\text{with}}\mathbf{A}=\mathbf{G}+\sigma \beta \left(\frac{\mathbf{S}+{\mathbf{S}}^T}{2}\right),$$

(47)

Here, the N × M-matrix W is given by

$$w_{nm}=\frac{e^{-\frac{1}{2}{∥\frac{{\mathbf{x}}_n-{\mathbf{y}}_m}{\sigma }∥}^2}}{{\sum }_{m^{\prime }=1}^M{e}^{-\frac{1}{2}{∥\frac{{\mathbf{x}}_n-{\mathbf{y}}_{m^{\prime }}}{\sigma }∥}^2}}$$

and the M × M-matrix G is

$$g_{ij}={\delta }_{ij}\sum _{n=1}^N{w}_{ni}.$$

A is a symmetric positive-definite M × M matrix. The M × M matrix A is symmetric and positive-definite. Since both G and W depend on Y, this equation is nonlinear in Y and has to be solved iteratively. Following [62–65, 71], we solved Equation (47) at each value of σ and for each iteration directly via Cholesky-factorization.

We implemented periodic and nonperiodic boundary conditions by appropriate choice of the matrix S. S must be positive (semi)definite for the energy to be bounded from below. We used the canonical 2D Laplacian stencil of order 2, to construct the M × M matrix

$$S=\left(\begin{array}{ccccccccccccc}\hfill -4+2a\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 1\hfill & \hfill -4+a\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill -4+a\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 1\hfill & \hfill \cdots \hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill \\ \hfill 1\hfill & \hfill -4+2a\hfill & \hfill 1\hfill \\ \hfill 1\hfill & \hfill -4+a\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill \cdots \hfill & \hfill 1\hfill & \hfill -4\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill -4\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 1\hfill \\ \hfill ⋮\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill -4+a\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill -4+a\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 1\hfill & \hfill -4+2a\hfill \end{array}\right),$$

Here, a = 0 for periodic boundary conditions and a = 1 for nonperiodic boundary conditions. In Appendix 3, we also present simulation results for a fourth derivative stencil, in which S² was used for the continuity term. We used random stimulus positions and orientations as well as stimuli arranged on a grid in feature space. For random stimuli, positions were drawn from a uniform distribution in [0, 1] × [0, 1]. Orientations s _z were drawn from the circular stimulus ensemble, with |s _z | = 0.08 as in [65]. Stimuli from grid-like ensembles were distributed evenly-spaced in [0, 1] × [0, 1] and contained 2 ^k evenly space orientations with |s _z | = 0.08.

To compute the energy of pinwheel-free configuration, we initiated simulations with a stripe-like orientation preference pattern with the same typical spacing as the observed orientation maps and annealed from σ = 0.035 to σ = 0.03.

To enable comparison between simulations with different numbers of stimuli, we scaled the continuity parameter proportionally to N such that the equivalent η was approximately constant. The simulated domain then contained a comparable number of hyper columns for all stimulus numbers.

Pinwheel density from simulations

Pinwheels locations in models were identified by the crossings of the zero contour lines of real and imaginary parts of the orientation map. Estimation of local column spacing Λ(x) was done using the wavelet analysis introduced in [127, 128]. In short, an overcomplete basis of complex Morlet wavelets at various scales and orientations was compared to the OPM pattern at each location. Λ(x) was estimated by the scale of the best matching wavelet. The mean column spacing 〈Λ(x)〉 _x of a given map was then calculated from the local column spacing by spatial averaging. For details we refer to [38, 127, 128]. Given N _pw pinwheels in a simulated cortical area of size L², we defined the pinwheel density [25, 38, 102]

$$\rho =N_{pw}\frac{{〈\Lambda 〉}_x^2}{L^2}.$$

The pinwheel density ρ is a dimensionless quantity and depends only on the layout of orientation columns. The pinwheel density defined in this way is large for patchy and small for more band-like columnar layouts.

Appendix 1

The impact of nonoriented stimuli

The main text of this article contains a complete analysis of optimal dimension-reducing mappings of the EN model with a circular ensemble of orientation stimuli. These optima are simple regular orientation stripes or square pinwheel crystals. The circular orientation stimulus ensemble, however, 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 realm of complex dynamics in 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 also examined the EN model in detail for a richer ensemble of stimuli. In this ensemble, called a uniform stimulus ensemble in the following, orientation stimuli are uniformly distributed on the disk {s _z , |s _z | ≤ 2}, a choice common to many previous studies, e.g., [19, 25, 81]. The uniform ensemble in particular contains unoriented stimuli with |s _z | = 0. Intuitively, the presence of these unoriented stimuli might be expected to fundamentally change the importance 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. As such an effect should be independent of retinotopic distortions and to aid comparison with our previous results, we will again start with a fixed uniform retinotopy r(x) = 0.

The linear stability properties of the unselective fixed point are independent of the ensemble of orientation stimuli (〈|s _z |²〉 = 2 throughout this article). The coefficients in Equation (14), however, of course depend on the fourth moment of the stimulus distribution. Inserting 〈|s _z |⁴〉 = 16/3 into Equation (32), we obtain

$$\begin{array}{ccc}\hfill a_1=\frac{1}{3{\sigma }^6}-\frac{1}{{\sigma }^4}+\frac{1}{2{\sigma }^2}\hfill & \hfill a_2=\frac{1}{4\pi {\sigma }^6}-\frac{1}{6\pi {\sigma }^8}\hfill & \hfill a_3=-\frac{1}{12\pi {\sigma }^8}+\frac{1}{8\pi {\sigma }^6}\hfill \\ \hfill a_4=-\frac{1}{6\pi {\sigma }^8}+\frac{1}{4\pi {\sigma }^6}-\frac{1}{8\pi {\sigma }^4}\hfill & \hfill a_5=-\frac{1}{12\pi {\sigma }^8}\hfill & \hfill a_6=\frac{1}{8\pi {\sigma }^6}-\frac{1}{12\pi {\sigma }^8}\hfill \\ \hfill a_7=\frac{1}{9{\pi }^2{\sigma }^{10}}-\frac{1}{12{\pi }^2{\sigma }^8}\hfill & \hfill a_8=\frac{1}{18{\pi }^2{\sigma }^{10}}\hfill & \hfill a_9=-\frac{1}{16{\pi }^3{\sigma }^{12}}\hfill \\ \hfill a_{10}=\frac{1}{9{\pi }^2{\sigma }^{10}}-\frac{1}{12{\pi }^2{\sigma }^8}\hfill & \hfill a_{11}=\frac{1}{18{\pi }^2{\sigma }^{10}}.\hfill \end{array}$$

The angle-dependent interaction functions are then given by

$$\begin{array}{ll}\hfill g\left(\alpha \right)& =\frac{1}{{\sigma }^4}\left(1-2e^{-k_c^2{\sigma }^2}-e^{2k_c^2{\sigma }^2\left(cos\alpha -1\right)}\left(1-2e^{-k_c^2{\sigma }^{2}cos\alpha }\right)\right)\\ +\frac{1}{2{\sigma }^2}\left(e^{2k_c^2{\sigma }^2\left(cos\alpha -1\right)}-1\right)+\frac{32}{3{\sigma }^6}{e}^{-2k_c^2{\sigma }^2}{sinh}^4\left(1/2k_c^2{\sigma }^{2}cos\alpha \right)\\ \hfill f\left(\alpha \right)& =\frac{1}{{\sigma }^4}\left(1-e^{-2k_c^2{\sigma }^2}\left(cosh\left(2k_c^2{\sigma }^{2}cos\alpha \right)+2cosh\left(k_c^2{\sigma }^{2}cos\alpha \right)\right)+2e^{-k_c^2{\sigma }^2}\right)\\ +\frac{1}{2{\sigma }^2}\left(e^{-2k_c^2{\sigma }^2}cosh\left(2k_c^2{\sigma }^{2}cos\alpha \right)-1\right)+\frac{16}{3{\sigma }^6}{e}^{-2k_c^2{\sigma }^2}{sinh}^4\left(1/2k_c^2{\sigma }^{2}cos\alpha \right).\end{array}$$

Both functions are depicted in Figure 22 for two different values of the effective intracortical interaction range σ/Λ. They qualitatively resemble the functions depicted in Figure 5. Figure 23 displays the phase diagram of the EN model with uniform stimulus ensemble. As summarized in the main part of this article, it is almost identical to that obtained for the circular stimulus ensemble (Figure 6). Two different optimal states are found, square pinwheel crystals (sPWCs) and orientation stripes (OSs) separated by a phase border at σ/Λ ≃ 0.15. Both fixed points are stable for all σ/Λ. Figure 23b-k demonstrates, that these analytical results are confirmed by direct numerical simulations of Equation (3) with r(x) = 0. As for the circular stimulus ensemble, we also tested the stability of stationary n-ECP solutions with 2 ≤ n ≤ 20 by numerical evaluation of the criteria for intrinsic and extrinsic stability (see 'Methods' section). We find all n-ECPs with 2 ≤ n ≤ 20 intrinsically unstable for all interaction ranges σ/Λ. The simple phase space structure furthermore apparently remains unchanged if we consider the model far from pattern formation threshold as shown in Figure 24. Simulations bear a close resemblance to the simulations with circular orientation stimulus ensemble (Figure 7). Either convergence to sPWC-like patterns or patterns with large orientation stripe domains is observed. Again, pinwheel annihilation in the case of large σ/Λ is less rapid than close to threshold (Figure 24a,b). The linear pinwheel-free zones increase their size over the time course of the simulations, eventually leading to a stripe pattern. For smaller interaction ranges σ/Λ, the OPM layout rapidly converges toward a crystal-like rhombic arrangement of pinwheels with dislocations and pinwheel density close to 4.

Figure 22

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

Figure 23

Optimal solutions of the EN model for uniform stimulus ensemble 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 effective 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, r = 0.1, σ/Λ = 0.12 (η = 0.57) and uniform stimulus ensemble. (d) Analytically predicted optimum for σ/Λ ≲ 0.15 (rPWC). (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 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) obtained in a simulation of Equation (3) with r(x) = 0, r = 0.1, σ/Λ = 0.15 (η = 0.41) and uniform stimulus ensemble. (i) Analytically predicted optimum for σ/Λ ≳ 0.15 (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 in simulations (parameters as in g) for different values of the control parameter r (black circles) and analytically predicted value (solid green line).

Figure 24

Numerical analysis of the EN dynamics with uniform orientation stimulus ensemble and fixed representation of visual space far from pattern formation threshold. (a) OPMs and their power spectra in a representative simulation of Equation (3) with r(x) = 0, r = 0.8, σ/Λ = 0.3 (η = 0.028) and uniform stimulus ensemble. (b) 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) OPMs and their power spectra in a representative simulation of Equation (3) with r(x) = 0 and σ/Λ = 0.12 (η = 0.57), other parameters as in a. (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).

Figure 25 shows that taking retinotopic distortions into account yields an almost identical picture compared to the circular stimulus ensemble. For small interaction range σ/Λ, the analytically predicted optimum is a quadratic pinwheel crystal with pinwheel density ρ = 4. For larger σ/Λ, the analytically predicted optimum is an orientation stripe pattern with pinwheel density ρ = 0. Our results are confirmed by direct simulations of Equations (3, 4) (Figure 25b-e). The simulation results are virtually indistinguishable from the circular stimulus ensemble.

Figure 25

Phase diagram of the EN model for the joint mapping of visual space and orientation with a uniform orientation stimulus ensemble. (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) 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.

All together, the EN dynamic given by Equations (3, 4) and in particular the set of ground states of the EN model and their stability regions appear almost identical when considering either a circular or a uniform orientation stimulus ensemble. We found two different optima depending on the parameter regime, orientation stripes for larger interaction ranges and quadratic pinwheel crystals for shorter interaction ranges. In addition, the EN dynamics appears to be unchanged by the presence of unoriented stimuli.

Appendix 2

Strength of retinotopic coupling

In our manuscript, we have shown that retinotopic distortions only have a weak influence on the optima of the EN model as well as its dynamics (see Figures 10 and 12). Here, we quantify the influence of retinotopic distortions on the pattern formation process by comparing the angle-dependent interaction function for retinotopic coupling g _r (α) with angle-dependent interaction function of the EN model with fixed retinotopy. We use the ratio

$$c=\frac{\left|\right|g_r\left(\cdot \right)|{|}_2}{\left|\right|g\left(\cdot \right)|{|}_2}$$

as a measure to quantify the influence of retinotopic distortions. || · ||₂ denotes the 2-Norm,

$$\left|\right|f\left(\cdot \right)|{|}_2={\int }_0^{2\pi }{f}^2\left(\alpha \right)d\alpha .$$

If c is larger than one, g _r (α) dominates the total interaction function g _r (α) + g(α) and retinotopic distortions may strongly influence the layout and stability of solutions of the EN model. On the other hand, if c is small, the solutions and their stability properties are expected to not change substantially when including variable retinotopy into the EN model. Figure 26 displays the parameter c in the s₄-σ/Λ-plane for the EN model at threshold for two different conditions, η = η _r and η _r = 0. In the latter case, retinotopic distortions are expected to have the strongest impact. However, in both cases, c ≪ 1, in almost all of the parameter space, implying little influence of retinotopic deviations. Only for small σ/Λ and small s₄, c is larger than one. As shown in Figure 12, this leads to slight deformations of the stability regions for rhombs, and stripes in this region of parameter space but does not result in novel optimal solutions.

Figure 26

Strength of coupling between orientation map and retinotopy in the EN model. (a) Ratio of ||g _r (·)||₂/||g(·)||₂ in the s₄-σ/Λ-plane for the EN model at threshold and η = η _r . (b) As a, but for η _r = 0, i.e., strongest coupling. Note the logarithmic scaling of the colormap.

Appendix 3

Grid-like stimulus ensembles

References [64, 65]) performed simulations with stimuli distributed in regular intervals in feature space, called grid-like ensemble. For comparison, we also performed deterministic annealing simulations with grid-like stimulus sets of varying size with nonperiodic boundary conditions (see 'Methods' section). For these grid-like stimulus patterns, a competition between stripes and rhombs is observed (Figure 27a). Notably, these are the only two stable states identified by our analysis for the circular stimulus ensemble. For nonperiodic boundary conditions, rhombic pinwheel arrangements seem energetically favored for grid-like stimuli, almost independently of the size of the stimulus set. The average pinwheel density for N = 100 × 100 × 8 stimuli was ρ = 3.4 (Figure 27b). As expected from the predominantly rhombic arrangement of pinwheels, NN-pinwheel distances concentrate around half the typical column spacing and pinwheel pairs at short distances are not observed (Figure 27c). With these features, the maps obtained substantially differ from the experimentally observed pinwheel statistics [38].

The discrete EN model with fourth derivative

In previous studies of the EN model, alternative definitions of the continuity term in the EN model have been explored [64]. A general continuity term for the spatially continuous formulation of the EN for OPM and retinotopy is a linear differential operator which will suppress the emergence of high-frequencies during the EN dynamics. A finite-wavelength instability is expected in this case, although the precise expressions for the critical σ and the typical wavelength will differ. As linear terms do not enter in the higher-order derivatives of the EN functional, changing the continuity term is not expected to alter the stability results obtained in this study.

To numerically test this expected robustness of our results for the EN model with discrete fixed sets of stimuli (see Figures 18 and 19), we also performed simulations using deterministic annealing with a fourth derivative stencil (see 'Methods' section). Figure 28 illustrates that the results almost perfectly match the ones for the second-order derivative, considered in the main part of this article (Figures 18, 19 and Figure 27).

Figure 27

The EN model with deterministic annealing and stimuli, distributed on a grid in feature space. (a) OPMs (left) and RMs (right) for N = 20 × 20 × 8 (upper row), N = 50 × 50 × 8 (middle row) and N = 100 × 100 × 8 (lower row) stimuli and nonperiodic boundary conditions (annealing rate χ = 0.999). β is the continuity parameter in the conventional definition of the EN model (see 'Methods' section, Equation (46)) and is scaled, such that a comparable number of columns emerges in all simulation for each N. (b) Pinwheel densities of EN solutions for different numbers of stimuli (annealing rate χ = 0.999). Crosses mark individual simulations, red line indicates average values. (c) Statistics of nearest neighbor pinwheel distances for pinwheels of (upper left) arbitrary and (upper right) opposite and equal charge for 100×100×8 stimuli and nonperiodic boundary conditions, averaged over four simulations (red curves). Black curves represent fits to the experimental data from [38]. Lower left: SD of pinwheel densities estimated from randomly selected regions in the OPM. Black dashed curve indicates SD for a two-dimensional Poisson process of equal density.

Figure 28

The EN model with deterministic annealing and fourth derivative stencil. (a) OPMs (left) and RMs (right) for N = 103 (upper row), N = 10⁴ (middle row) and N = 10⁵ (lower row) stimuli, non-periodic boundary conditions, and annealing rate χ = 0.999. β is the continuity parameter in the conventional definition of the EN model (see 'Methods' section, Equation (46)) and is scaled, such that a comparable number of columns is emerging in the simulations for each stimulus set. (b) Pinwheel densities (upper left) of EN solutions, SD of pinwheel densities estimated from randomly selected regions in the solutions (upper right). Crosses mark individual simulations, red line indicates average values. Black dashed curve indicates SD for a two-dimensional Poisson process of equal density. Statistics of nearest neighbor pinwheel distances for pinwheels of arbitrary (lower left) and (lower right) opposite and equal charge for 10⁵ random stimuli and periodic boundary conditions, averaged over four simulations (red curves). Black curves represent fits to the experimental data from [38]. (c) As a but for nonperiodic boundary conditions. (d) As b, but for non-periodic boundary conditions. (e) OPMs (left) and RMs (right) for N = 20 × 20 × 8 (upper row), N = 50 × 50 × 8 (middle row) and N = 100 × 100 × 8 (lower row) stimuli, nonperiodic boundary conditions, annealing rates χ = 0.999. (f) As b, but for nonperiodic boundary conditions and grid-like stimuli.

When annealing with periodic boundary conditions, the solutions very much resemble our gradient descent dynamics simulations with Laplacian term. The larger the set of stimuli, the more stripe-like are the OPMs obtained (Figure 28a) and consequently pinwheel densities decrease (Figure 28b, upper right). The exponent for the SD is considerably lower than for the Poisson process (Figure 28b, upper right).Typical NN-pinwheel distances concentrate around half the typical column spacing and in particular pinwheel pairs with short distances lack completely (Figure 28b, lower left and right).

For nonperiodic boundary conditions and random stimuli, we found that retinotopic distortions are much more pronounced. They however decreased with increasing number of stimuli. For large stimulus numbers, we observed stripe-like orientation preference domains which are interspersed with lattice-like pinwheel arrangements (see Figure 28c), lower row, upper left corner of the OPM). Similarly to the periodic boundary conditions, short distance pinwheel pairs occur much less frequently than in the experimentally observed maps, indicating an increased regularity in the pinwheel arrangements compared to realistic OPMs (Figure 28d, lower left and right). This regularity also manifests itself in a smaller exponent of the SD compared to the Poisson process (Figure 28d).

Simulations with grid-like stimulus as, e.g., used in [64, 65] displayed a strong tendency toward rhombic pinwheel arrangements analogous to the second derivative case (Figure 27e,f)