Coordinated Optimization of Visual Cortical Maps (I) Symmetry-based Analysis

In the primary visual cortex of primates and carnivores, functional architecture can be characterized by maps of various stimulus features such as orientation preference (OP), ocular dominance (OD), and spatial frequency. It is a long-standing question in theoretical neuroscience whether the observed maps should be interpreted as optima of a specific energy functional that summarizes the design principles of cortical functional architecture. A rigorous evaluation of this optimization hypothesis is particularly demanded by recent evidence that the functional architecture of orientation columns precisely follows species invariant quantitative laws. Because it would be desirable to infer the form of such an optimization principle from the biological data, the optimization approach to explain cortical functional architecture raises the following questions: i) What are the genuine ground states of candidate energy functionals and how can they be calculated with precision and rigor? ii) How do differences in candidate optimization principles impact on the predicted map structure and conversely what can be learned about a hypothetical underlying optimization principle from observations on map structure? iii) Is there a way to analyze the coordinated organization of cortical maps predicted by optimization principles in general? To answer these questions we developed a general dynamical systems approach to the combined optimization of visual cortical maps of OP and another scalar feature such as OD or spatial frequency preference. From basic symmetry assumptions we obtain a comprehensive phenomenological classification of possible inter-map coupling energies and examine representative examples. We show that each individual coupling energy leads to a different class of OP solutions with different correlations among the maps such that inferences about the optimization principle from map layout appear viable. We systematically assess whether quantitative laws resembling experimental observations can result from the coordinated optimization of orientation columns with other feature maps.


Introduction
In case of the low order inter-map coupling energies strong inter-map coupling leads to a suppression of OP selectivity. This suppressive effect can be avoided by restricting coupling strengths. One aim of this article is to test different optimization principles and potentially rule out some optimization principles. When comparing our results from different optimization principles to biological data such parameter tuning reduces the practicability. In this supporting information we complement our study using the high order inter-map coupling energies. We show that in this case a suppression of OP selectivity cannot occur.
We derive coupled amplitude equations which, however, involve several mathematical assumptions. A systematic treatment as it is shown in the main article would imply that low order and higher order intermap coupling energies are in general non-zero. Low order energy terms would enter at third order in the expansion and higher order corrections could potentially alter the stability properties. In addition, higher order inter-map coupling energies can affect the stability of patterns. In the following, we assume all low order inter-map coupling energies to be zero and that contributions entering the amplitude equations at higher orders can be neglected. The obtained results are confirmed numerically in part (II) of this study.

Coupled amplitude equations: Higher order terms
We studied the coupled Swift-Hohenberg equations [z, z, z, o, o, o, o] o, o, z, z, z, z] , with the higher order inter-map coupling energies using weakly nonlinear analysis. We study Eq. (1) close to the pattern forming bifurcation where r z and r o are small. We therefore expand both control parameters in powers of the small expansion parameter µ r z = µr z1 + µ 2 r z2 + µ 3 r z3 + . . .
Close to the bifurcation the fields are small and thus nonlinearities are weak. We therefore expand both fields as We further introduced a common slow timescale T = r z t and insert the expansions in Eq.
We consider amplitude equations up to seventh order as this is the order where the nonlinearity of the higher order coupling energy enters first. For Eq. (5) and Eq. (6) to be fulfilled each individual order in µ has to be zero. At linear order in µ we get the two homogeneous equationŝ Thus z 1 and o 1 are elements of the kernel ofL 0 z andL 0 o . Both kernels contain linear combinations of modes with a wavevector on the critical circle i.e.
with the complex amplitudes A (1) j = B j e ıψj and ⃗ k j = k c,z (cos(jπ/n), sin(jπ/n)), ⃗ k ′ j = k c,o (cos(jπ/n), sin(jπ/n)). In view of the hexagonal or stripe layout of the OD pattern shown in Fig. 1, n = 3 is an appropriate choice. In the following sections we assume k c,o = k c,z = k c i.e. the Fourier components of the emerging pattern are located on a common circle. To account for species differences we also analyzed models with detuned OP and OD wavelengths in part (II) of this study.
At second order in µ we getL As z 1 and o 1 are elements of the kernel r z1 = r o1 = 0. At third order, when applying the solvability condition (see Methods), we get We insert the leading order fields Eq. (8) and obtain the amplitude equations These uncoupled amplitude equations obtain corrections at higher order. There are fifth order, seventh order and even higher order corrections to the uncoupled amplitude equations. In addition, at seventh order enters the nonlinearity of the higher order inter-map coupling energies. The amplitude equations up to seventh order are thus derived from and corresponding equations for the fields o 1 , o 3 , and o 5 . The field z 1 is given in Eq. (8) and its amplitudes A (1) and B (1) are determined at third order. The field z 3 contains contributions from modes off the critical Its amplitude A (3) are determined at fifth order. The field z 5 also contains contributions from modes off the critical circle z 5,of f and on the critical circle i.e.
Its amplitude A (5) are determined at seventh order. This leads to a series of amplitude equations which are solved order by order. We set r z2 = r z and r o2 = r o and rescale to the fast time. This leads to We can combine the amplitude equations up to seventh order by introducing the amplitudes A j = j . This leads to the amplitude equations For simplicity we have written only the simplest inter-map coupling terms. Depending on the configuration of active modes additional contributions may enter the amplitude equations. In addition, for the product-type coupling energy, there are coupling terms which contain the constant δ, see Methods. In case of A ≪ B ≪ 1 the inter-map coupling terms in dynamics of the modes B are small. In this limit the dynamics of the modes B decouples from the modes A and we can use the uncoupled OD dynamics, see Methods. In the following, we use the effective inter-map coupling strength ϵB 4 (and τ B 4 ).

Optima of particular optimization principles: Higher order coupling terms
In this article we demonstrated that the low order coupling terms can lead to a complete suppression of OP selectivity i.e. vanishing magnitude of the order parameter |z|. As the coupling terms are effectively linear they not only influence pattern selection but also whether there is a pattern at all. This is in general not the case for higher order coupling energies using the amplitude equations Eq. (15). In this case the coupling is an effective cubic interaction term and complete selectivity suppression is impossible.
Moreover, as in the low-order energy case, we could identify the limit r z ≪ r o in which the backreaction onto the OD map formally becomes negligible. The potential is of the form where δ i,j denotes the Kronecker delta and the uncoupled contributions Amplitude equations can be derived from the potential by ∂ t A i = −δV /δA i . We have not written terms involving the modes A j − or B j . The complete amplitude equations involving all modes and the corresponding coupling coefficients are given in Text S2. As for the low order coupling energies terms involving the constant δ depend only on the coupling coefficient of the product-type energy τ . In the fol-lowing we specify the amplitude equations for negligible backreaction where B = B hex , B = B st , or B = 0.

Product-type energy
First, we studied the higher order product-type inter-map coupling energy U = τ o 4 |z| 4 . As for the lower order version of this coupling energy the shift δ(γ) explicitly enters the amplitude equations resulting in a rather complex parameter dependence, see Eq. (68) in the Methods section.
The equation for the mode A 3 is given by interchanging the modes A 2 and A 3 in Eq. (18). The equations for the modes A i − are given by interchanging the modes A i and A i − and interchanging the modes B i and B i .
In this case, at low inter-map coupling the OP stripes given by with ϕ 1 − ϕ 1 − = 2ψ 1 + π run parallel to the OD stripes. Their stationary amplitudes are given by The parameter dependence of these stripe solutions is shown in Fig. S1A.
At large inter-map coupling the attractor states of the OP map consist of a stripe pattern containing only two preferred orientations, namely ϑ = ϕ 1 and ϑ = ϕ 1 + π/2. The zero contour lines of the OD map are along the maximum amplitude of orientation preference minimizing the energy term.
In addition there are rhombic solutions which exist also in the uncoupled case, see Fig. S1B. However, these rhombic solutions are energetically not favored compared to stripe solutions, see Fig. S1C. The inclusion of the inter-map coupling makes these rhombic solution even more stripe-like.
In case of a OD constant solution the amplitude equations read with g ii = 1 + δ 4 τ , g ij = 2 + 2δ 4 τ and f ij = 2 + 2δ 4 τ . Inter-map coupling thus leads to a renormalization of the uncoupled interaction terms. Stationary solutions are stripes with the amplitude and rhombic solutions with the stationary phases ϕ 1 + ϕ 1 − − ϕ 2 − ϕ 2 − = π and the stationary amplitudes In the case of OD hexagons we identify, in addition to stripe-like and rhombic solutions, uniform solutions When solving the amplitude equations numerically we have seen that the phase relations vary with the inter-map coupling strength τ for non-uniform solutions. But for the uniform solution the phase relations are independent of the inter-map coupling strength. We use the ansatz for uniform solutions where δ i,j is the Kronecker delta and ∆ a constant parameter. This leads to the stationarity condition Four types of stationary solutions exist namely the ∆ = 0, ∆ = π, which we already observed in case of the low order energies, and the solutions which depends on B and δ and thus on the bias γ. The course of Eq. (27) as a function of γ is shown in Fig. S2B. Stationary amplitudes for these solutions are given by We study the stability properties of OP stripe-like, rhombic and uniform solutions using linear stability analysis. The eigenvalues of the stability matrix , see also Text S3, are calculated numerically. Linear

Bifurcation diagram
For increasing inter-map coupling strength the amplitudes of the OP solutions are shown in Fig. S1.
In case of inter-map coupling strength dependent stationary phases, stationary solutions are calculated numerically using a Newton method and initial conditions close to these solutions. We followed the unstable solutions (dashed lines in Fig. S1) until this method did not converge anymore. Not shown are solutions which are unstable in general. The parameter dependence of OP solutions when interacting with OD stripes is shown in Fig. S1A,B. Similar to the low order variant of this coupling energy the amplitude of the stripes pattern A 1 is suppressed while the amplitude of the opposite mode A 1 − grows.
Finally both amplitudes collapse, leading to an orientation scotoma solution. In contrast to the low order variant this stripe pattern is stable for arbitrary large inter-map coupling. In case of OP rhombic solutions inter-map coupling transforms this solution by reducing the amplitudes A 2 = A 2 − while increasing the amplitudes A 3 = A 3 − . Without OD bias this solution is then transformed into the orientation scotoma stripe pattern, similar to the low order variant of this energy. In contrast to the low order energy, for non-zero bias the amplitudes A 2 and A 3 stay small but non-zero.
The parameter dependence of OP solutions when interacting with OD hexagons is shown in Fig. S1C,D.
For a small OD bias (γ = γ * ) OP rhombic solutions decay into OP stripe-like patterns. These stripe-like patterns stay stable also for large-inter map coupling. In case of a larger OD bias (γ = 3γ * ), both the OP stripe and the OP rhombic solutions decay into the uniform PWC solution. Thus for small bias there is a bistability between stripe-like and uniform PWC solutions while for larger OD bias the uniform PWC solution is the only stable solution. The potential of OP stripe and OP rhombic solutions is shown in

Phase diagram
The stability properties of all stationary solutions are summarized in the phase diagram Fig. S2. Compared to the gradient-type interaction energy we cannot scale out the dependence on r o . The phase diagram is thus plotted for r o = 0.2. We rescale the inter-map coupling strength as τ B 4 where B is the stationary amplitude of the OD hexagons. In the regime of stable OD stripes there is a transition from OP stripes towards the orientation scotoma stripe solution. In the regime of stable OD hexagons there is a transition from OP stripes towards PWC solutions (red line). The stability border of PWC solutions is strongly OD bias dependent and has a peak at γ ≈ 2γ * . For small OD bias γ the uniform solution Eq. (27) is stable. With increasing bias there is a smooth transition of this solution until at γ = γ c the d = π uniform solution becomes stable. In Fig. S2C the stability border γ c between the two types of uniform solutions is plotted as a function of r o . We observe that there is only a weak dependence on the control parameter and γ c ≈ 2γ * . Figure S3 illustrates the uniform solutions Eq. (27) for different values of the OD bias γ. For small bias, the OP pattern has six pinwheels per unit cell. Two of them are located at OD maxima while one is located at an OD minimum. The remaining three pinwheels are located near the OD border. With increasing bias, these three pinwheels are pushed further away from the OD border, being attracted to the OD maxima. With further increasing bias three shifted pinwheels merge with the one at the OD maximum building a single charge 1 pinwheel centered on a contralateral peak. The remaining two pinwheels are located at an ispi and contra peak, respectively. Note, compared to the Braitenberg PWC of the ∆ = 0 uniform solution the charge 1 pinwheel here is located at the contralateral OD peak. Finally, the charge 1 pinwheels split up again into four pinwheels. With increasing bias the solution more and more resembles the Ipsi-center PWC (∆ = π solution) which is stable also in the lower order version of the coupling energy. Finally, at γ/γ * ≈ 2 the Ipsi-center PWC becomes stable and fixed for γ > 2γ * . The distribution of preferred orientations for different values of the bias γ is shown in Fig. S3E,F, reflecting the symmetry of each pattern. The distribution of intersection angles is shown in Fig. S3G. Remarkably, all solutions show a tendency towards perpendicular intersection angles. This tendency is more pronounced with increasing OD bias. At about γ/γ * ≈ 1.9 parallel intersection angles are completely absent and at γ/γ * ≈ 2 there are exclusively perpendicular intersection angles.

Gradient-type energy U = ϵ |∇o·∇z| 4
Finally, we examine the higher order version of the gradient-type inter-map coupling. The interaction terms are independent of the OD shift δ. In this case the coupling strength can be rescaled as βB 4 and is therefore independent of the bias γ. The bias in this case only determines the stability of OD stripes, hexagons or the constant solution.

Stationary solutions and their stability
As for its lower order pendant a coupling to OD stripes is relatively easy to analyze. The energetic ground state corresponds to OP stripes with the direction perpendicular to the OD stripes for which U = 0.
In addition, there are rhombic solutions with the stationary amplitudes r z /(5 + 80ϵB 4 ). In case the OD map is a constant, o(x) = δ, the gradient-type inter-map coupling leaves the OP unaffected. As for its lower order pendant the stationary states are therefore OP stripes running in an arbitrary direction and the uncoupled rhombic solutions.
In case of OD hexagons we identified three types of non-uniform solutions. Besides stripe-like solutions of z(x) with one dominant mode we find rPWCs A j = A j − = (A, a, A) with a ≪ A and distorted We calculated the stability properties of all stationary solutions by linear stability analysis considering perturbations of the amplitudes A j → A + a j , A j − → A + a j − and of the phases ϕ j → ϕ j + φ j , This leads to a perturbation matrix M . In general amplitude and phase perturbations do not decouple. We therefore calculate the eigenvalues of the perturbation M matrix numerically. It turns out that for this type of coupling energy only the uniform solutions with ∆ = ± arccos ( 5

13
) are stable while the ∆ = 0 and ∆ = π solutions are unstable in general.

Bifurcation diagram
For increasing inter-map coupling strength the amplitudes of the OP stripe and OP rhombic solutions are shown in Fig. S4A. In case of inter-map coupling strength dependent stationary phases, stationary solutions are calculated numerically using a Newton method and initial conditions close to these solutions.
We followed the unstable solutions (dashed lines in Fig. S4) until this method did not converge anymore.
Not shown are solutions which are unstable in general. In case of stable OD hexagons there is a transition from rPWC (blue) towards distorted rPWC (green). The distorted rPWCs then decay into the hPWC (red). In case of OP stripes (black dashed lines) inter-map coupling leads to a slight suppression of the dominant mode and a growth of the remaining modes. This growth saturates at small amplitudes and thus the solution stays stripe-like. This stripe-like solution remains stable for arbitrary large inter-map coupling. Therefore there is a bistability between hPWC solutions and stripe-like solutions for large inter-map coupling.
The stability borders for the rPWC and distorted rPWC solutions were obtained by calculating their bifurcation diagram numerically from the amplitude equations, see Text S2. With increasing map coupling we observe a transition from a rPWC towards a distorted rPWC at ϵB 4 ≈ 0.033 (blue dashed line in Above ϵB 4 ≈ 0.12 the hPWC is energetically preferred compared to stripe-like solutions (red dashed line in in Fig. S5) and thus corresponds to the energetic ground state for large inter-map coupling, see

Phase diagram
We calculated the phase diagram of the coupled system in the limit r z ≪ r o , shown in Fig. S5. The phase diagram contains the stability borders of the uncoupled OD solutions γ * , γ * 2 , γ * 3 , γ * 4 . They correspond to vertical lines, as they are independent of the inter-map coupling in the limit r z ≪ r o . At γ = γ * hexagons become stable. Stripe solutions become unstable at γ = γ * 2 . At γ = γ * 3 the homogeneous solution becomes stable while at γ = γ * 4 hexagons loose their stability. In units γ/γ * the borders γ * 2 , γ * 3 , γ * 4 vary slightly with r o , see Figure 9, and are drawn here for r o = 0.2. We rescale the inter-map coupling strength as ϵB 4 where B is the stationary amplitude of the OD hexagons. The stability borders of OP solutions are then horizontal lines. For γ < γ * or for γ > γ * 4 pinwheel free orientation stripes are dynamically selected. For γ * < γ < γ * 4 and above a critical effective coupling strength ϵB 4 ≈ 0.042 hPWC solutions are stable and become the energetic ground state of Eq. (16) above ϵB 4 ≈ 0.117. Below ϵB 4 ≈ 0.065, rPWC solutions are stable leading to a bistability region between rPWC and hPWC solutions. We find in this region that rhombic solutions transform into distorted rhombic solutions above an effective coupling strength of ϵB 4 ≈ 0.033.

Interaction induced pinwheel crystals
First, we studied the spatial layout of the rhombic solutions which is illustrated in Fig. S6. The rPWC solutions are symmetric under rotation by 180 degree. The rhombic solution has 4 pinwheels per unit cell and its pinwheel density is thus ρ = 4 cos(π/6) ≈ 3.5. One may expect that the energy term Eq. (2) favors pinwheels to co-localize with OD extrema. In case of the rhombic layout there is only one pinwheel at an OD extremum while the other three pinwheels are located at OD saddle-points which are also energetically favorable positions with respect to U . The orientation selectivity |z(x)| for the rPWC is shown in Fig. S6B. The pattern of selectivity is arranged in small patches of highly selective regions.
The hexagonal layout of the two stable uniform solutions is shown in Fig. S7. The ∆ = ± arccos(5/13) solutions have six pinwheels per unit cell. Their pinwheel density is therefore ρ = 6 cos π/6 ≈ 5.2. Three pinwheels of the same topological charge are located at the extrema of the OD map. Two of these are located at the OD maximum while one is located at the OD minimum. The remaining three pinwheels are not at an OD extremum but near the OD border. The distance to the OD border depends on the OD bias, see Fig. S7D. For a small bias (γ ≈ γ * ) these three pinwheels are close to the OD borders and with increasing bias the OD border moves away from the pinwheels. The pinwheel in the center of the OP hexagon is at the contralateral OD peak. Because these pinwheels organize most of the map while the others essentially only match one OP hexagon to its neighbors we refer to this pinwheel crystal as the Contra-center pinwheel crystal. Note, that some pinwheels are located at the vertices of the hexagonal pattern. Pinwheels located between these vertices (on the edge) are not in the middle of this edge. Solutions with ∆ = ± arccos(5/13) are therefore not symmetric under a rotation by 60 degree but symmetric under a rotation by 120 degree. Therefore the solution ∆ = + arccos(5/13) cannot be transformed into the solution ∆ = − arccos(5/13) by a rotation of the OD and OP pattern by 180 degrees. This symmetry is also reflected by the distribution of preferred orientations, see

Summary
We derived amplitude equations and analyzed ground states of the higher order inter-map coupling energies. We calculated local and global optima and derived corresponding phase diagrams. A main difference between phase diagrams for low order and high order coupling energies consists in the collapse of orientation selectivity above a critical coupling strength that occurs only in the low order models. In contrast, for the high order versions, orientation selectivity is preserved for arbitrarily strong inter-map coupling. In order to neglect the backreaction on the dynamics of the modes B we assumed A ≪ B ≪ 1.
decrease in B cannot be compensated by another parameter (as it would be r z in case of the low order inter-map coupling energies). For a finite B even higher order corrections to the amplitude equations than those presented here can thus become significant. Such terms are we neglected in the present treatment.
In part (II) of this study we numerically confirm our main results for the higher order inter-map coupling energies.
From a practical point of view, the analyzed phase diagrams and pattern properties indicate that the higher order gradient-type coupling energy is the simplest and most convenient choice for constructing models that reflect the correlations of map layouts in the visual cortex. For this coupling, intersection angle statistics are reproduced well, pinwheels can be stabilized, and pattern collapse cannot occur.         1,3,4). B U = τ |z| 4 o 4 , using stationary amplitudes from Fig. (S1). D Positions of rPWCs move continuously (pinwheel 5,6).