Model definition and model symmetries

with appropriately chosen units for x and R.

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

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.

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

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

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

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.

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.

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

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.

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.

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.

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

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.

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.

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

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.

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

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

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

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.

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 .

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.

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

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