### 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*_{1}(**x**), *R*_{2}(**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

\begin{array}{ll}\hfill \mathcal{C}\left[z,\mathbf{r}\right]& =-{\u3008ln\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)\u22152{\sigma}^{2}}\u3009}_{\mathbf{S}}\phantom{\rule{2em}{0ex}}\\ \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},\phantom{\rule{2em}{0ex}}\end{array}

with ∇ = (*∂*_{
x
} , *∂*_{
y
} ) ^{T} , and *η* ∈ [0, 1]. The ratios *σ* ^{2}/*η* and *σ* ^{2}/*η*_{
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)\phantom{\rule{1em}{0ex}}=\phantom{\rule{1em}{0ex}}-2\frac{\delta \mathcal{F}\left[z,\mathbf{r}\right]}{\delta \stackrel{\u0304}{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)\phantom{\rule{1em}{0ex}}=\phantom{\rule{1em}{0ex}}-\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)\phantom{\rule{1em}{0ex}}=\phantom{\rule{1em}{0ex}}{\u3008\left[{s}_{z}-z\left(\mathbf{x}\right)\right]e\left(\mathbf{x},\mathbf{S},z,\mathbf{r}\right)\u3009}_{\mathbf{S}}+\eta \Delta z\left(\mathbf{x}\right)

(3)

{\partial}_{t}\mathbf{r}\left(\mathbf{x}\right)\phantom{\rule{1em}{0ex}}=\phantom{\rule{1em}{0ex}}{\u3008\left[{\mathbf{s}}_{r}-\mathbf{X}-\mathbf{r}\left(\mathbf{x}\right)\right]e\left(\mathbf{x},\mathbf{S},z,\mathbf{r}\right)\u3009}_{\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)\phantom{\rule{1em}{0ex}}=\phantom{\rule{1em}{0ex}}\frac{{e}^{-\left(\mid {\mathbf{s}}_{r}-\mathbf{X}-\mathbf{r}{(\mathbf{x})\mid}^{2}\right)\u22152{\sigma}^{2}}{e}^{-\left(\mid {s}_{z}-z{(\mathbf{x})\mid}^{2}\right)\u22152{\sigma}^{2}}}{\int {d}^{2}y{e}^{-\left(\mid {\mathbf{s}}_{r}-\mathbf{X}-\mathbf{r}{(\mathbf{y})\mid}^{2}\right)\u22152{\sigma}^{2}}{e}^{-\left(\mid {s}_{z}-z{(\mathbf{y})\mid}^{2}\right)\u22152{\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)\phantom{\rule{1em}{0ex}}=\phantom{\rule{1em}{0ex}}z\left(\mathbf{x}+\mathbf{y}\right)\\ {\widehat{T}}_{\mathbf{y}}\mathbf{r}\left(\mathbf{x}\right)\phantom{\rule{1em}{0ex}}=\phantom{\rule{1em}{0ex}}\mathbf{r}\left(\mathbf{x}+\mathbf{y}\right),\end{array}

rotations

\begin{array}{c}{\widehat{R}}_{\beta}z\left(\mathbf{x}\right)\phantom{\rule{1em}{0ex}}=\phantom{\rule{1em}{0ex}}{e}^{2i\beta}z\left({\mathrm{\Omega}}_{-\beta}\mathbf{x}\right)\\ {\widehat{R}}_{\beta}\mathbf{r}\left(\mathbf{x}\right)\phantom{\rule{1em}{0ex}}=\phantom{\rule{1em}{0ex}}{\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)\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}\stackrel{\u0304}{z}\left(\Psi \mathbf{x}\right)\\ \widehat{P}\mathbf{r}\left(\mathbf{x}\right)\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}\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]\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}{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]\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}{F}^{z}\left[{\widehat{R}}_{\beta}z,{\widehat{R}}_{\beta}\mathbf{r}\right]

(6)

\widehat{P}{F}^{z}\left[z,\mathbf{r}\right]\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}{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]\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}{F}^{z}\left[{e}^{i\varphi}z,\mathbf{r}\right]

(8)

{\mathbf{F}}^{r}\left[z,\mathbf{r}\right]\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}{\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 {\u3008\mid {s}_{z}{\mid}^{2}\u3009}_{\mathbf{S}}=\u3008\mid {s}_{z}{\mid}^{2}\u3009=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 Λ^{2}, 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)\phantom{\rule{1em}{0ex}}\simeq \phantom{\rule{1em}{0ex}}{\mathbf{L}}_{r}\left[\mathbf{r}\right]=\frac{1}{16\pi {\sigma}^{4}}\int {d}^{2}y\phantom{\rule{0.3em}{0ex}}{e}^{-\frac{{\left(\mathbf{x}-\mathbf{y}\right)}^{2}}{4{\sigma}^{2}}}\widehat{\mathbf{A}}\phantom{\rule{0.3em}{0ex}}\mathbf{r}\left(\mathbf{y}\right)+{\eta}_{r}\u25b3\mathbf{r}\left(\mathbf{x}\right)

(10)

{\partial}_{t}z\left(\mathbf{x}\right)\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}\simeq \phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}{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*σ*^{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}{\stackrel{\u0303}{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){\stackrel{\u0303}{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)\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}{\mathbf{k}}_{\varphi}{e}^{i{\mathbf{k}}_{\varphi}\mathbf{x}}+\mathsf{\text{c}}.\mathsf{\text{c}}.\\ {\mathbf{r}}_{T}\left(\mathbf{x}\right)\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}{\mathbf{k}}_{\varphi +\pi \u22152}{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.

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\u2215\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*_{0}*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

{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} = 0, we recover the trivial fixed point of the EN dynamics by construction, as shown above. This means that within this subspace *A*_{0} = 0 is either a minimum or a maximum of the EN energy functional (Equation (1)). Furthermore, for *A*_{0} → ∞ 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} ≠ 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*_{0}*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,\stackrel{\u0304}{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}{\stackrel{\u0304}{z}}_{l}\right]=\sum _{j,k,l}{\alpha}_{j}{\beta}_{k}{\gamma}_{l}{N}_{3}^{z}\left[{z}_{j},{z}_{k},{\stackrel{\u0304}{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,\stackrel{\u0304}{z}\right]=\sum _{j=1}^{11}{a}_{j}{N}_{3}^{j}\left[z,z,\stackrel{\u0304}{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]

{\u0226}_{i}=r{A}_{i}-{A}_{i}\sum _{j=1}^{N}{g}_{ij}\mid {A}_{j}{\mid}^{2}-{\u0100}_{{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,\stackrel{\u0304}{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,\stackrel{\u0304}{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

{\u0226}_{j}\left(t\right)=-\frac{\partial {U}_{A}}{\partial {\u0100}_{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}{\u0100}_{i}{\u0100}_{{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*_{0}*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}{2{g}_{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**_{1}, **k**_{2}) = *α >*0. By inserting this ansatz into Equation (15) and assuming uniform amplitude \mid {A}_{1}\mid \phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\mid {A}_{2}\mid \phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\mid {A}_{2}\mid \phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\mid {A}_{4}\mid \phantom{\rule{2.77695pt}{0ex}}=\mathcal{A}, we obtain

{\mathcal{A}}^{2}=\frac{r}{{g}_{00}+{g}_{0\pi}+{g}_{0\alpha}+{g}_{0\pi -\alpha}-2{f}_{0\alpha}}.

(19)

The phase relations of the four amplitudes are given by

\begin{array}{c}{\varphi}_{1}+{\varphi}_{3}\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}{\varphi}_{0}\\ {\varphi}_{2}+{\varphi}_{4}\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}{\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., *ϕ*_{0}, Δ_{0} = *ϕ*_{1} - *ϕ*_{3} and Δ_{1} = *ϕ*_{2} - *ϕ*_{4}. For an rPWC parameterized by these phases, Δ_{0} shifts the absolute positions of the pinwheels in *x*-direction, Δ_{1} shifts the absolute positions of the pinwheels in *y*-direction, and *ϕ*_{0} 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}-2{f}_{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 Λ^{2}, 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}-2{f}_{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

{\u0226}_{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 \u2215\Lambda =\frac{1}{2\pi}\sqrt{log\left(1\u2215\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}}\phantom{\rule{1em}{0ex}}\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}}\phantom{\rule{1em}{0ex}}\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)\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}\frac{1}{{\sigma}^{4}}\left(1-2{e}^{-{k}_{c}^{2}{\sigma}^{2}}-{e}^{2{k}_{c}^{2}{\sigma}^{2}\left(cos\alpha -1\right)}\left(1-2{e}^{-{k}_{c}^{2}{\sigma}^{2}cos\alpha}\right)\right)\\ \phantom{\rule{1em}{0ex}}+\frac{1}{2{\sigma}^{2}}\left({e}^{2{k}_{c}^{2}{\sigma}^{2}\left(cos\alpha -1\right)}-1\right)+\frac{8}{{\sigma}^{6}}{e}^{-2{k}_{c}^{2}{\sigma}^{2}}{sinh}^{4}\left(1\u22152{k}_{c}^{2}{\sigma}^{2}cos\alpha \right)\\ f\left(\alpha \right)\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}\frac{1}{{\sigma}^{4}}\left(1-{e}^{-2{k}_{c}^{2}{\sigma}^{2}}\left(cosh\left(2{k}_{c}^{2}{\sigma}^{2}cos\alpha \right)+2cosh\left({k}_{c}^{2}{\sigma}^{2}cos\alpha \right)\right)+2{e}^{-{k}_{c}^{2}{\sigma}^{2}}\right)\\ \phantom{\rule{1em}{0ex}}+\frac{1}{2{\sigma}^{2}}\left({e}^{-2{k}_{c}^{2}{\sigma}^{2}}cosh\left(2{k}_{c}^{2}{\sigma}^{2}cos\alpha \right)-1\right)+\frac{4}{{\sigma}^{6}}{e}^{-2{k}_{c}^{2}{\sigma}^{2}}{sinh}^{4}\left(1\u22152{k}_{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.

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

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)\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}{L}_{z}\left[z\right]+{Q}^{z}\left[\mathbf{r},z\right]+{N}_{3}^{z}\left[z,z,\stackrel{\u0304}{z}\right]+\cdots

(26)

{\partial}_{t}\mathbf{r}\left(\mathbf{x}\right)\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}{\mathbf{L}}_{r}\left[\mathbf{r}\right]+{\mathbf{Q}}^{r}\left[z,\stackrel{\u0304}{z}\right]+\cdots \phantom{\rule{0.3em}{0ex}}.

(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,\stackrel{\u0304}{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,\stackrel{\u0304}{z}\right]\right].

(28)

We remark that as {\lambda}_{T\u2215L}^{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)\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}\approx \phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}{L}_{z}\left[z\right]\underset{{N}_{3}^{r}\left[z,z,\stackrel{\u0304}{z}\right]}{\underset{\u23df}{-{Q}^{z}\left[{\mathbf{L}}_{r}^{-1}\left[{\mathbf{Q}}^{r}\left[z,\stackrel{\u0304}{z}\right]\right],z\right]}}+{N}_{3}^{z}\left[z,z,\stackrel{\u0304}{z}\right]\\ =\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}{L}_{z}\left[z\right]+{N}_{3}^{r}\left[z,z,\stackrel{\u0304}{z}\right]+{N}_{3}^{z}\left[z,z,\stackrel{\u0304}{z}\right].\end{array}

(29)

The nonlinearity {N}_{3}^{r}\left[z,z,\stackrel{\u0304}{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,\stackrel{\u0304}{z}\right]=\sum _{j=1}^{12}{a}_{r}^{j}{N}_{r}^{j}\left[z,z,\stackrel{\u0304}{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,\stackrel{\u0304}{z}\right], these functions are given by

\begin{array}{c}{g}_{r}\left(\alpha \right)\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}\frac{{\left(\left(1-{\sigma}^{2}-2{e}^{-{k}_{c}^{2}{\sigma}^{2}}\right){e}^{2{k}_{c}^{2}{\sigma}^{2}\left(cos\phantom{\rule{2.77695pt}{0ex}}\alpha -1\right)}+{e}^{-{k}_{c}^{2}{\sigma}^{2}}\right)}^{2}}{2{\sigma}^{4}\left({\eta}_{r}+{\sigma}^{2}{e}^{-2{k}_{c}^{2}{\sigma}^{2}\left(cos\phantom{\rule{2.77695pt}{0ex}}\alpha -1\right)}\right)}\\ {f}_{r}\left(\alpha \right)\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}\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,\stackrel{\u0304}{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)\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}-\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}\u22152}-{e}^{-{\Delta}_{jk}^{2}{\sigma}^{2}\u22152}\right)}^{2}-{e}^{-{\Delta}_{jk}^{2}{\sigma}^{2}}\right)\\ *\left(\Im \left({A}_{j}{\u0100}_{k}\right)cos\left({\Delta}_{jk}\mathbf{x}\right)+\Re \left({A}_{j}{\u0100}_{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*_{0}*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(2{k}_{c}\right)}\left(\begin{array}{c}\hfill {k}_{c}sin\left(2{k}_{c}{x}_{1}\right)\hfill \\ \hfill {k}_{c}sin\left(2{k}_{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 ~ 2*k*_{
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 \u2215\Lambda =\frac{1}{2\pi}\sqrt{log\left(1\u2215\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.

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,\stackrel{\u0304}{z}\right] in Equation (14) for arbitrary stimulus ensembles are given by

\begin{array}{ccc}\hfill {a}_{1}=\frac{\u3008\mid {s}_{z}{\mid}^{4}\u3009}{16{\sigma}^{6}}-\frac{\u3008\mid {s}_{z}{\mid}^{2}\u3009}{2{\sigma}^{4}}+\frac{1}{2{\sigma}^{2}}\hfill & \hfill {a}_{2}=\frac{\u3008\mid {s}_{z}{\mid}^{2}\u3009}{8\pi {\sigma}^{6}}-\frac{\u3008\mid {s}_{z}{\mid}^{4}\u3009}{32\pi {\sigma}^{8}}\phantom{\rule{1em}{0ex}}\hfill & \hfill {a}_{3}=-\frac{\u3008\mid {s}_{z}{\mid}^{4}\u3009}{64\pi {\sigma}^{8}}+\frac{\u3008\mid {s}_{z}{\mid}^{2}\u3009}{16\pi {\sigma}^{6}}\hfill \\ \hfill {a}_{4}=-\frac{\u3008\mid {s}_{z}{\mid}^{4}\u3009}{32\pi {\sigma}^{8}}+\frac{\u3008\mid {s}_{z}{\mid}^{2}\u3009}{8\pi {\sigma}^{6}}-\frac{1}{8\pi {\sigma}^{4}}\phantom{\rule{1em}{0ex}}\hfill & \hfill {a}_{5}=-\frac{\u3008\mid {s}_{z}{\mid}^{4}\u3009}{64\pi {\sigma}^{8}}\hfill & \hfill {a}_{6}=\frac{\u3008\mid {s}_{z}{\mid}^{2}\u3009}{16\pi {\sigma}^{6}}-\frac{\u3008\mid {s}_{z}{\mid}^{4}\u3009}{64\pi {\sigma}^{8}}\hfill \\ \hfill {a}_{7}=\frac{\u3008\mid {s}_{z}{\mid}^{4}\u3009}{48{\pi}^{2}{\sigma}^{10}}-\frac{\u3008\mid {s}_{z}{\mid}^{2}\u3009}{24{\pi}^{2}{\sigma}^{8}}\hfill & \hfill {a}_{8}=\frac{\u3008\mid {s}_{z}{\mid}^{4}\u3009}{96{\pi}^{2}{\sigma}^{10}}\hfill & \hfill {a}_{9}=-\frac{3\u3008\mid {s}_{z}{\mid}^{4}\u3009}{256{\pi}^{3}{\sigma}^{12}}\hfill \\ \hfill {a}_{10}=\frac{\u3008\mid {s}_{z}{\mid}^{4}\u3009}{48{\pi}^{2}{\sigma}^{10}}-\frac{\u3008\mid {s}_{z}{\mid}^{2}\u3009}{24{\pi}^{2}{\sigma}^{8}}\hfill & \hfill {a}_{11}=\frac{\u3008\mid {s}_{z}{\mid}^{4}\u3009}{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)\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}\frac{\u3008\mid {s}_{z}{\mid}^{2}\u3009}{2{\sigma}^{4}}\left(1-2{e}^{-{k}_{c}^{2}{\sigma}^{2}}-{e}^{2{k}_{c}^{2}{\sigma}^{2}\left(cos\alpha -1\right)}\left(1-2{e}^{-{k}_{c}^{2}{\sigma}^{2}cos\alpha}\right)\right)\\ \phantom{\rule{1em}{0ex}}+\frac{1}{2{\sigma}^{2}}\left({e}^{2{k}_{c}^{2}{\sigma}^{2}\left(cos\alpha -1\right)}-1\right)+\frac{2\u3008\mid {s}_{z}{\mid}^{4}\u3009}{{\sigma}^{6}}{e}^{-2{k}_{c}^{2}{\sigma}^{2}}{sinh}^{4}\left(1\u22152{k}_{c}^{2}{\sigma}^{2}cos\alpha \right)\\ f\left(\alpha \right)\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}\frac{\u3008\mid {s}_{z}{\mid}^{2}\u3009}{2{\sigma}^{4}}\left(1-{e}^{-2{k}_{c}^{2}{\sigma}^{2}}\left(cosh\left(2{k}_{c}^{2}{\sigma}^{2}cos\alpha \right)+2cosh\left({k}_{c}^{2}{\sigma}^{2}cos\alpha \right)\right)+2{e}^{-{k}_{c}^{2}{\sigma}^{2}}\right)\\ \phantom{\rule{1em}{0ex}}+\frac{1}{2{\sigma}^{2}}\left({e}^{-2{k}_{c}^{2}{\sigma}^{2}}cosh\left(2{k}_{c}^{2}{\sigma}^{2}cos\alpha \right)-1\right)+\frac{\u3008\mid {s}_{z}{\mid}^{4}\u3009}{{\sigma}^{6}}{e}^{-2{k}_{c}^{2}{\sigma}^{2}}{sinh}^{4}\left(1\u22152{k}_{c}^{2}{\sigma}^{2}cos\alpha \right).\end{array}

(33)

Again, without loss of generality, we set 〈|*s*_{z}|^{2}〉 = 2. At criticality, both functions are parameterized by the continuity parameter *η* ∈ (0, 1) for the OPM or, equivalently, the interaction range \sigma \u2215\Lambda =\frac{1}{2\pi}\sqrt{log\left(1\u2215\eta \right)} and the fourth moment 〈|*s*_{z}|^{4}〉 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
} |^{4} ≈ 0), together with a large fraction of high orientation energy stimuli (|*s*_{
z
} |^{4} 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}|^{4}〉 = (〈|*s*_{z}|^{2}〉)^{2} = 4. The fourth moment of the uniform stimulus ensemble is 〈|*s*_{z}|^{4}〉 = 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}=\u3008\mid {s}_{z}{\mid}^{4}\u3009-{\u3008\mid {s}_{z}{\mid}^{2}\u3009}^{2}=\u3008\mid {s}_{z}{\mid}^{4}\u3009-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*_{4}. 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.

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*_{4}. 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*_{4} = 0 and *s*_{4} = 4/3. Above a certain value of the fourth moment around *s*_{4} = 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.

We numerically tested these analytical predictions by simulations of Equation (3) (**r**(**x**) = **0**) with two additional stimulus ensembles with *s*_{4} = 6 and *s*_{4} = 8 (see 'Methods' section). Figure 13a shows snapshots of a simulation with (*r* = 0.1, *σ*/Λ = 0.2 (*η* = 0.2)) and *s*_{4} = 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 13c shows snapshots of a simulation with (*r* = 0.1, *σ*/Λ = 0.2 (*η* = 0.2)) and *s*_{4} = 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,\stackrel{\u0304}{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)\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}& =\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}\frac{2\u3008\mid {s}_{z}{\mid}^{4}\u3009}{{\sigma}^{6}}{e}^{-2{k}_{c}^{2}{\sigma}^{2}}{sinh}^{4}\left(1\u22152{k}_{c}^{2}{\sigma}^{2}cos\alpha \right)\phantom{\rule{2em}{0ex}}\\ \phantom{\rule{1em}{0ex}}-\frac{\u3008\mid {s}_{z}{\mid}^{2}\u3009}{2{\sigma}^{4}}{e}^{-2{k}_{c}^{2}{\sigma}^{2}}\left(\left(cosh\left(2{k}_{c}^{2}{\sigma}^{2}cos\alpha \right)-2cosh\left({k}_{c}^{2}{\sigma}^{2}cos\alpha \right)\right)-2{e}^{{k}_{c}^{2}{\sigma}^{2}}-{e}^{2{k}_{c}^{2}{\sigma}^{2}}\right)\phantom{\rule{2em}{0ex}}\\ \phantom{\rule{1em}{0ex}}+\frac{1}{2{\sigma}^{2}}\left(1+{e}^{-2{k}_{c}^{2}{\sigma}^{2}}cosh\left(2{k}_{c}^{2}{\sigma}^{2}cos\alpha \right)\right).\phantom{\rule{2em}{0ex}}\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

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*_{4} as can be seen by inserting Equations (33) and (35) into Equation (36). In the infinite fourth moment limit *s*_{4} → ∞, 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.

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
} |^{4}〉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*_{4} = 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.

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*_{4} for large *s*_{4}.

Hence, in the limit *s*_{4} → ∞, 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^{6} 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^{3} to 10^{4}) 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^{5} (b), *N* = 10^{6} (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].

For *N* = 10^{5} 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^{6} 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^{5} 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^{5}, *N* = 10^{6}), 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^{3}, *N* = 10^{4} and *N* = 10^{5} 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^{5}. To facilitate comparison, we superimposed fits to the experimentally observed statistics [38] for orientation maps in tree shrews, ferrets and galagos.

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^{5}, 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^{5}, 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.