# Coverage, Continuity and Visual Cortical Architecture

## Abstract

###### Abstract

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

##### Results

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

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

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

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

## Introduction

The pattern of orientation columns in the primary visual cortex (V1) of carnivores, primates and their close relatives are among the most intensely-studied structures in the cerebral cortex and a large body of experimental (e.g. [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]) and theoretical work (e.g. [14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39]) has been dedicated to uncovering its organization principles and the circuit level mechanisms that underlie its development and operation.

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

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

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

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

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

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

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

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

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

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

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

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

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

## Results and Discussion

### Model definition and model symmetries

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

with appropriately chosen units for and .

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

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

(1) |

in which the functional measures the coverage of a stimulus space and the functional the continuity of its cortical representation. The stimulus space is defined by an ensemble {} of idealized pointlike stimuli, each described by two features: and which specify the orientation of the stimulus and its position in visual space (Fig. 2b). and are given by

with , and . The ratios and control the relative strength of the coverage term versus the continuity term for OPM and RM, respectively. denotes the average over the ensemble of stimuli.

Minima of the energy functional are stable fixed points of the gradient descent dynamics

(2) | |||||

called the EN dynamics in the following. These dynamics read

(3) | |||||

(4) |

where is the activity-pattern, evoked by a stimulus in a model cortex with retinotopic distortions and OPM . It is given by

Figure 2 illustrates the general features of the EN dynamics using the example of a single stimulus. Fig. 2a shows a model orientation map with a superimposed uniform representation of visual space. A single pointlike, oriented stimulus with position and orientation (Fig. 2b) evokes a cortical activity pattern (Fig. 2c). The activity-pattern in this example is of roughly Gaussian shape and is centered at the point, where the location 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). According to Eqs. (3, 4), each stimulus and the evoked activity pattern induce a modification of the orientation map and retinotopic map, shown in Fig. 2d. Orientation preference in the activated regions is shifted towards the orientation of the stimulus. The representation of visual space in the activated regions is locally contracted towards 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 Eqs. (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 are uniformly distributed across visual space. (ii) For the distribution of stimulus orientations, and are independent. (iii) Orientations are distributed uniformly in .

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

rotations

with 22 rotation matrix

and reflections

where is the 22 reflection matrix. Equivariance means that

(5) | |||||

(6) | |||||

(7) |

with mutatis mutandis the same equations fulfilled by the vector-field .

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 Eqs. (3, 4), by construction. Due to the second assumption, the dynamics is also equivariant with respect to shifts in orientation , i.e.

(8) | |||||

(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 throughout this paper. This normalization can always be restored by a rescaling of (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 [83]). Even an hypothesized neuronal mono-layer would consist of more than 20x20 neurons per hypercolumn area , constituting a quite dense sampling of the spatial periodicity. Treating the feature space as a continuum implements the concept that the cortical representation has to cover as good as possible the infinite multiplicity of conceivable stimulus feature combinations.

### The orientation unselective fixed point

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

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

(10) | |||||

(11) |

where with 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 can be calculated by Fourier-transforming Eq. (10):

where and . A diagonalization of this matrix equation yields the eigenvalues

with corresponding eigenfunctions (in real space)

where . These eigenfunctions are longitudinal (L) or transversal (T) wave patterns. In the longitudinal wave, the retinotopic distortion vector lies parallel to which leads to a “compression wave” (Fig. 3b, left). In the transversal wave pattern (Fig. 3b, right), the retinotopic distortion vector is orthogonal to . We note that the linearized Kohonen model [61] was previously found to exhibit the same set of eigenfunctions [81]. Because both spectra of eigenvalues are smaller than zero for every , , and (Fig. 3a), the uniform retinotopy is a stable fixed point of Eq. (4) irrespective of parameter choice.

The eigenfunctions of the linearized OPM dynamics are Fourier modes by translational symmetry. By rotational symmetry, their eigenvalues only depend on the wave number and are given by

(see [25]). This spectrum of eigenvalues is depicted in Fig. 3c. For , has a single maximum at . For

(12) |

this maximal eigenvalue is negative. Hence, the unselective state with uniform retinotopy is a stable fixed point of Eqs. (3, 4) and the only known solution of the EN model in this parameter range.

For , the maximal eigenvalue is positive, and the nonselective state is unstable with respect to a band of Fourier modes with wave numbers around (see Fig. 3c). This annulus of unstable Fourier modes is called the critical circle. The finite wavelength instability [84, 85, 86] (or Turing instability [87]) leads to the emergence of a pattern of orientation preference with characteristic spacing from the nonselective state on a characteristic timescale .

One should note that as in other models for the self-organization of orientation columns, e.g. [57, 15], 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 Eq. (11). Since 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 Eq. (11). Mathematically, this term directly results from the soft-competition in the “activity-dependent” coverage term of Eq. (1). The local facilitation is jointly mediated by coverage (first term) and continuity (second term) contributions.

Fig. 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: . In these pinwheel-free states, orientation preference is constant along one axis in cortex (perpendicular to the vector ), and each orientation is represented in equal proportion (see Fig. 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 orientation stripe (OS) solutions follows directly from the model’s symmetries Eqs. (5-9). Computing

demonstrates that is proportional to . This establishes that the subspace of functions is invariant under the dynamics Eq. (3). For , we recover the trivial fixed point of the EN dynamics by construction, as shown above. This means that within this subspace is either a minimum or a maximum of the EN energy functional Eq. (1). Furthermore, for 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 in the subspace of functions which then corresponds to a stationary state of the EN dynamics.

Regarding the dynamics of retinotopic deviations, the model’s symmetries Eqs. (5-9) can be invoked to show that for the state , the right-hand side of Eq. (4) has to be constant in space:

If this constant was non-zero the retinotopic map would drift with constant velocity. This, however, is impossible in a variational dynamics such that this constant must vanish. The orientation stripe solution (Fig. 4a) is to the best of our knowledge the only exact nontrivial stationary solution of Eqs. (3, 4) that can be established without any approximations.

### Doubly periodic and quasiperiodic 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 and the perfectly uniform retinotopy becomes an optimal solution for arbitrary orientation column layout . 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 orientation preference map Eq. (3) with . 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 becomes unstable, the entire set of nontrivial fixed points is determined by the third-order terms of the Volterra series representation of the operator [85, 35, 88, 86]. The symmetries Eqs. (5-9) restrict the general form of such a third-order approximation for any model of OPM optimization to

(13) |

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

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

(14) |

The individual nonlinear operators are with one exception nonlocal convolution-type operators and are given in the Methods section (Eq. (38)), together with a detailed description of their derivation. Only the coefficients depend on the properties of the ensemble of oriented stimuli.

To calculate the fixed points of Eq. (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 [85, 84, 86]. Planforms are patterns that are composed of a finite number of Fourier components, such as

for a pattern of orientation columns. With the above planform ansatz, we neglect any spatial dependency of the amplitudes for example due to long-wave deformations for the sake of simplicity and analytical tractability. When the dynamics is close to a finite wavelength instability, the essential Fourier components of the emerging pattern are located on the critical circle . The dynamic equations for the amplitudes of these Fourier components are called amplitude equations. For a discrete number of Fourier components of whose wave vectors lie equally spaced on the critical circle, the most general system of amplitude equations compatible with the symmetries Eqs. (5-9) has the form [88, 35]

(15) |

with . Here, and are the real-valued coupling coefficients between the amplitudes and . They depend on the differences between indices and are entirely determined by the nonlinearity in Eq. (13). If the wave vectors are parametrized by the angles , then the coefficients and are functions only of the angle between the wave vectors and . One can thus obtain the coupling coefficients from two continuous functions and that can be obtained from the nonlinearity (see Methods for details). In the following, these functions are called angle-dependent interaction functions. The amplitude equations Eq. (15) are variational if and only if and are real-valued. In this case they can be derived via

from an energy

(16) |

If the coefficients and are derived from Eq. (1), the energy for a given planform solution corresponds to the energy density of the EN energy functional considering only terms up to fourth order in .

The amplitude equations Eq. (15) enable to calculate an infinite set of orientation selective fixed points. For the above orientation stripe solution with one nonzero wave vector , the amplitude equations predict the so far undetermined amplitude

(17) |

and its energy

(18) |

Since , this shows that OS stationary solutions only exist for , i.e. in the symmetry breaking regime. As for all following fixed-points, specifies the energy difference to the homogeneous unselective state .

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

with amplitudes and . By inserting this ansatz into Eq. (15) and assuming uniform amplitude , we obtain

(19) |

The phase relations of the four amplitudes are given by

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

(20) |

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

(21) |

in the denominator of Eq. (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 Eq. (15). This solution class encompasses a large variety of realistic quasiperiodic OPM layouts and is therefore a good candidate solution class for models of OPM layout. In addition, Kaschube et al. 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 [38].

An n-ECP solution can be written as

with wave vectors distributed equidistantly on the upper half of the critical circle, complex amplitudes and binary variables determining whether the mode with wave vector or 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 Eq. (15) vanishes and the amplitude equations for the active modes reduce to a system of Landau equations

where is the -coupling matrix for the active modes. Consequently, the stationary amplitudes obey

(22) |

The energy of an n-ECP is given by

(23) |

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

Families of n-ECP solutions are depicted in Fig. 4c. The 1-ECP corresponds to the pinwheel-free OS pattern discussed above. For fixed , there are multiple planforms not related by symmetry operations which considerably differ in their spatial layouts. For , the patterns are spatially quasiperiodic, and are a generalization of the so-called Newell-Pomeau turbulent crystal [90, 91]. For , 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 Fig. 4c, n=15 lower right). Anisotropic planforms generally contain considerably fewer pinwheels (see Fig. 4c, n=15 lower left). All large n-ECPs, however, exhibit a complex quasiperiodic 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 active modes decay into ones with lower . (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 to a small finite set around a typical [88, 35]. A mathematical evaluation of both criteria leads to precise conditions for extrinsic and intrinsic stability of a planform (Methods). 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 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, orientation stripes, rhombic pinwheel crystals 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 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 -values uniformly arranged on a ring with radius [92, 57, 66, 64, 65]. 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 Eq. (11) is completely characterized by the continuity parameter . Equivalent to specifying is to fix the ratio of activation range and column spacing

(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 [88, 35]. 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 and . For the coefficients in Eq. (14) we obtain

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

(25) | |||||

These functions are depicted in Fig. 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 in Eq. (21), we find that the angle which minimizes the energy of rPWC fixed-point is . 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 , 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 (Methods). We found both, OS and sPWCs, to be intrinsically and extrinsically stable for all . Next, we tested for the stability of n-ECP solutions with . We found all n-ECP configurations with 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 OS and sPWC are always local minima of the energy functional.

By evaluating the energy assigned to the sPWC (Eq. (20)) and OS (Eq. (18)), we next identified two different regimes: (i) For short interaction range the sPWC possesses minimal energy and is therefore the predicted global minimum. (ii) For the OS is optimal. Figure 6a shows the resulting simple phase diagram. The sPWC and OS are separated by a phase border at . We numerically confirmed these analytical predictions by extensive simulations of Eq. (3) with and the circular stimulus ensemble (see Methods for details). Fig. 6b-c show snapshots of a representative simulation with short interaction range (, ). After the phase of initial pattern emergence (symmetry breaking), the OPM layout rapidly approaches a square array of pinwheels, the analytically predicted optimum (Fig. 6d). Pinwheel density time courses (see Methods) display a rapid convergence to a value close to the predicted density of 4 (Fig. 6e). Fig. 6f shows the stationary mean squared amplitudes of the pattern obtained for different values of the control parameter (black circles). For small control parameters, the pattern amplitude is perfectly predicted by Eq. (19) (solid green line). Fig. 6g-h show snapshots of a typical simulation with longer interaction range (, ). 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 (Fig. 6i). Pinwheel densities decay towards zero over the time course of the simulations (Fig. 6j). Also in this parameter regime, the mean squared amplitude of the pattern is well-predicted Eq. (17) for small . (Fig. 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) orientation stripes are the optimal dimension-reducing mapping. Both states are stable throughout the entire parameter range. All other planforms, in particular quasiperiodic 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 Eq. (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 OS and sPWCs. Neither of the two exhibits the characteristic aperiodic and pinwheel-rich organization of experimentally observed orientation preference maps. Furthermore, the pinwheel densities of both solutions ( for OS and for sPWCs) differ considerably from experimentally observed values [38] around . One way towards 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 Eq. (3).

To asses this possibility, we simulated Eq. (3) with and a circular stimulus ensemble for very large values of the control parameter . Fig. 7 displays snapshots of such a simulation for 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 (Fig. 7a, b). The OPM nevertheless converges towards 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 towards a crystal-like rhombic arrangement of pinwheels, however containing several dislocations (Fig. 24a) [85]. Dislocations are defects of roll or square patterns, where two rolls or squares merge into one, thus increasing the local wavelength of the pattern [84, 86]. Nevertheless, for all simulations, the pinwheel density rapidly reaches a value close to 4 (Fig. 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 [85, 86].

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 quadratic 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. . We recall that in their seminal publication, Durbin and Mitchison in particular demonstrated interesting correlations between the map of orientation preference and the map of visual space [21]. 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 system of equations Eqs. (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 non-orientation selective state becomes unstable, the entire set of nontrivial fixed points of Eqs. (3, 4) is determined by the third-order terms of the Volterra series representation of the nonlinear operators and . The model symmetries Eqs. (5-9) restrict the general form of the leading order terms for any model for the joint optimization of OPM and RM to

(26) | |||||

(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 . These retinotopic distortions will in turn alter the dynamics of the OPM via the quadratic complex-valued operator . Close to the point of pattern onset (), the timescale of OPM development, , becomes arbitrarily large and retinotopic deviations evolve on a much shorter timescale. This separation of timescales allows for an adiabatic elimination of the variable , assuming it to always be at the equilibrium point of Eq. (27):

(28) |

We remark that as for all finite wave numbers , the operator is indeed invertible when excluding global translations in the set of possible perturbations of the trivial fixed point. Via Eq. (28), the coupled dynamics of OPM and RM is thus reduced to a third-order effective dynamics of the OPM:

(29) | |||||

The nonlinearity accounts for the coupling between OPM and RM. Its explicit analytical calculation for the EN model is rather involved and yields a sum

The individual nonlinear operators 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 are completely independent of the orientation stimulus ensemble.

The adiabatic elimination of the retinotopic distortions results in an equation for the OPM (Eq. (29)) which has the same structure as Eq. (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 Eq. (15). The nonlinear coefficients and are determined from the angle-dependent interaction functions and . For the operator , these functions are given by

verifying that , is independent of the orientation stimulus ensemble. Besides the interaction range the continuity parameter 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 , the functions and tend to zero and as expected one recovers the results presented above for fixed uniform retinotopy. The functions and are depicted in Fig. 8 for various interaction ranges and retinotopic continuity parameters .