On the Effects on Cortical Spontaneous Activity of the Symmetries of the Network of Pinwheels in Visual Area V1

This paper challenges and extends earlier seminal work. We consider the problem of describing mathematically the spontaneous activity of V1 by combining several important experimental observations including (1) the organization of the visual cortex into a spatially periodic network of hypercolumns structured around pinwheels, (2) the difference between short-range and long-range intracortical connections, the first ones being rather isotropic and producing naturally doubly periodic patterns by Turing mechanisms, the second one being patchy, and (3) the fact that the Turing patterns spontaneously produced by the short-range connections and the network of pinwheels have similar periods. By analyzing the PO maps, we are able to classify all possible singular points (the pinwheels) as having symmetries described by a small subset of the wallpaper groups. We then propose a description of the spontaneous activity of V1 using a classical voltage-based neural field model that features isotropic short-range connectivities modulated by non-isotropic long-range connectivities. A key observation is that, with only short-range connections and because the problem has full translational invariance in this case, a spontaneous doubly periodic pattern generates a 2-torus in a suitable functional space which persists as a flow-invariant manifold under small perturbations, for example when turning on the long-range connections. Through a complete analysis of the symmetries of the resulting neural field equation and motivated by a numerical investigation of the bifurcations of their solutions, we conclude that the branches of solutions which are stable over an extended range of parameters are those that correspond to patterns with an hexagonal (or nearly hexagonal) symmetry. The question of which patterns persist when turning on the long-range connections is answered by (1) analyzing the remaining symmetries on the perturbed torus and (2) combining this information with the Poincaré–Hopf theorem. We have developed a numerical implementation of the theory that has allowed us to produce the predicted patterns of activities, the planforms. In particular we generalize the contoured and non-contoured planforms predicted by previous authors.


Introduction
The primary area (V1) of the visual cortex is one of the first locations targeted by connections from the thalamus which relays (and processes) inputs from the retina. In some mammals like primates, cats or ferrets (see [1,2,3]), this cortical area is very precisely organized in modules, called cortical hypercolumns, which process visual attributes (like local orientation, spatial frequency) of different points in the visual field. Note that there is almost no experimental evidence of hypercolumns from histology (see [4]). Most of our knowledge derives from functional evidence, i.e. when an external stimulus is applied. In this work, we focus on the processing of the local orientation of visual stimuli which is reflected by the ability of some neurons to fire only when a drifting grating of specific (called preferred) orientation is presented in their receptive field. The distribution, on the cortical sheet, of the preferred orientation of these particular neurons (see [5,6]) defines a Preferred Orientation map (hereafter called the PO map) which has a near lattice structure and which is continuous except at particular points called pinwheels.
It has been argued by several authors (Ermentrout [7], Bressloff and Cowan [8,9]) that V1 is to some extent structured like a crystal: it can be approximated by a plane where the pinwheels are arranged in a doubly periodic lattice and the main features of cortical activity in V1 can be interpreted in this framework. This idealization naturally introduces symmetries in the problem, which makes deeper analysis accessible. As long as the pinwheels are nearly arranged on a periodic lattice we can expect that the main conclusions of our analysis will remain valid.
There is also experimental evidence (see below) that the spatial distribution of connections emanating from one neuron in V1 differs according to whether the connections are local (within one hypercolumn) or long-range (between different hypercolumns). Local connections are considered to be isotropic. In the absence of long-range connections the network (or field) of local connections would be Euclidean-invariant, that is, invariant under rigid displacements and reflections in the plane, a property that is transmitted to the model equations [7]. On the other hand, long-range connections are subject to the constraint of respecting the symmetries of the PO map, thereby reducing the full Euclidean group symmetries to a crystallographic subgroup associated with the lattice of pinwheels. Moreover, experimental observations suggest that the strength of long-range connections is significantly weaker than that of local connections, which allows for treating the long-range connections as a perturbation of the local ones.
However, as shown initially by Ermentrout and Cowan [7], spontaneous activity of V1 without long-range connections also produces naturally doubly periodic steady patterns by the Turing mechanism [10]. This fact has been at the origin of a series of papers where hallucinatory patterns (seen by patients under the influence of drugs or other stimuli of the visual cortex which are not coming from the thalamus) have been explained as Turing patterns bifurcating in V1 [8,11]. It was assumed in these papers that the lattice of pinwheels can be itself approximated as a continuum in the plane, that is, every point in the plane is a pinwheel. Then the full system of connections, local and long-range, keeps translation invariance: any planar translation applied to a pattern yields another one, up to periodicity. Therefore a pattern is not an isolated solution but rather generates a manifold of solutions under the action of translations, which is called an orbit under the action of the group of translations. This orbit, thanks to the periodicity of the pattern, can be identified with a 2-torus in a suitable functional space. Moreover, as long as this torus is a normally hyperbolic manifold (it means here that steady-state patterns are hyperbolic along directions transversal to the orbit), it persists as a flow-invariant manifold under small perturbations.
When the lattice of pinwheels is discrete, long-range connections reduce the translation group R 2 to a discrete subgroup isomorphic to Z 2 . Our aim in this paper is to study the effect of switching on such long-range connections on the tori of solutions of the system with local connections only. As in the above cited papers we assume that the strength of long-range connections is small compared to that of local connections. The introduction of the long-range connections destroys part of the symmetries but not all of them. The perturbed torus features these remaining symmetries. We classify possible dynamics on this invariant manifold by applying topological methods (the Poincaré-Hopf theorem) together with the symmetries. Numerical (direct) simulations of the equations then allow one to determine which phase portrait is actually observed among those which have been theoretically identified.
Our motivation for this work is to understand how introducing a discrete lattice of pinwheels would modify the states and dynamics of spontaneous activity in V1, and we expect this could have some consequences on the theory of hallucinations of Bressloff and Golubitsky et al. in [8].
The problem of the effect of long-range connections with discrete translation invariance on the Turing patterns which bifurcate when these connections are inactive has been addressed by Bressloff [12] (see also Baker and Cowan [13]) in the following context. These authors assumed that the periods of the patterns and of the lattice of pinwheels were not correlated. An analysis inspired by methods initially introduced for hydrodynamical systems allowed them to build Ginzburg-Landau equations, in order to describe the slow modulation of the bifurcating patterns under the effect of spatial periodic forcing due to the long-range connections.
There are, however, some experimental observations which suggest that in fact, due to synaptic plasticity, the periods of the pinwheel lattice and of the Turing patterns are close to each other [14,15]. It is therefore not unrealistic to assume that these periods are in fact equal, which we do in this study. Our approach differs on another point: we do not assume the system to be close to bifurcation from homogeneous state. Therefore the validity range of our results goes well beyond bifurcation analysis, and allows a rigorous mathematical treatment in the spirit of Lauterbach and Roberts [16] about forced symmetry breaking of group orbits of equilibria with spherical symmetry.
The way in which neurons in V1 acquire orientation selectivity is still controversial as the connections from the thalamus provide a very small percentage of the inputs: 95 % of these inputs are made of recurrent connections, i.e. intracortical connections (see [17]). Despite this experimental evidence there are still two extreme attitudes to account for these facts: either we assume that the recurrent connections provide much of the input to each neuron, or that each neuron mostly follows the thalamic input, discarding the recurrent connections; see [18]. Our approach lies in between.
We use the following standard neural field representation for the neural activity V of the neuron located at x in the connected domain Ω of R 2 representing V1: where S is the sigmoidal function with range (0, 1) and T /σ is the threshold. The slope is σ at x = T /σ . The domain Ω is usually taken as the infinite plane. This is not an unrealistic approximation in the analysis, because the number of pinwheels in V1 is rather large, typically several thousands [8]. However, in numerical computations one has to choose a bounded domain and the simplest choice is a square with periodic boundary conditions: opposite sides are identified. The connectivity function J represents the intracortical connections between neurons and I thalamus the thalamic input. Following [12], we decompose the connectivity function as follows: where J loc models the local connections and J LR (x, y) models the long-range connections. The small factor LR in front of J LR accounts for the fact mentioned above that the strength of the long-range connections is significantly weaker than that of the short-range connections [19]. From [4,15] we can assume that the local connections are isotropic (rotation invariant), and we further assume for simplicity that they are also homogeneous (translation invariant). Hence J loc only depends on the distance between x and y: x − y is the usual Euclidean norm in R 2 . The form of the function J LR is discussed in Sect. 3.1. The external (thalamic) input function I thalamus will be assumed equal to 0 throughout the paper in order to consider spontaneous activity only.

Remark 1 For all numerical experiments, we choose
Indeed, from [4,15], we know that the local connectivity is homogeneous in the cat and in the monkey. Note that we assume that the excitatory (respectively, inhibitory) connections are modeled as Gaussian with width σ loc (respectively, 2σ loc ) [15]. The The connectivity is tuned such that the first bifurcation point for the nonlinear gain σ is a supercritical D 4 -pitchfork (respectively, D 6 -pitchfork) bifurcation with wave-vector k = 1 for which the spots are stable, while the stripes are not. The size of the square cortex is 8 × 2π and the numerical mesh size is 1024 2 (top) and 3 × 512 2 (bottom), the threshold is T = 0.1 and σ loc = π · 0.395. We only used the 15 eigenvalues with the largest real parts to find the bifurcation points. See Appendix B for details concerning the numerical methods constant a = 4e −σ 2 loc /2 is tuned such that the first bifurcation of the homogeneous state has a wave-vector of norm 1; see [20] for more details and for bifurcation diagrams. σ loc controls the stability of stripes vs. spots.
Let us have a look at the bifurcation diagram in Fig. 1, top. It shows bifurcated patterns in a square Ω with periodic boundary conditions when only local connections are turned on. The bifurcation parameter in abscissa is the slope σ of the sigmoid function S. It shows, as expected, two primary bifurcated branches: one of stripes and one of "spots", which correspond to a periodic pattern with square symmetry. Near bifurcation other types of solutions branch off these primary branches. Observe that the stripes are always unstable and that the spots are stable in a small interval, up to the first secondary bifurcation. On the other hand two secondary branches of doubly periodic states bifurcate sub-critically and recover stability after bending back in the increasing σ direction. These states are then stable on large intervals of values of σ . In their stability domain they look pretty much like hexagonal patterns. In fact, they approximate the exact hexagonal patterns that would bifurcate on a domain with hexagonal symmetry (and periodic boundary conditions). If the size of Ω is increased, these two branches come closer to each other and bifurcate closer from the primary bifurcation point. In the limit of a domain with infinite size, they merge into one branch, which bifurcates at the same point as the spots and stripes, and they correspond exactly to the hexagonal patterns which are well known from Turing theory of pattern formation. We compare this diagram to the case where the domain Ω is an hexagon with periodic boundary conditions; see Fig. 1, bottom. We were able to find numerically the branches of spots (transcritical) and stripes (pitchfork). We did not compute the secondary branches because our point was to show that the hexagons (spots) are stable over an extended region.
Since we are interested in stable states, we can consider in the analysis, instead of these "imperfect hexagons", the exact hexagonal patterns that will show up in an hexagonal lattice. This will make the analysis simpler and richer at the same time, and the results will still be relevant for the secondary "approximate" hexagonal patterns observed in Fig. 1, top. Despite the fact that square patterns (spots in Fig. 1, top) are stable only in a small region near bifurcation, we shall discuss both square and hexagonal patterns for completeness.
Our primary motivation here is to understand how the long-range connections with discrete lattice symmetry affect the dynamics of the network. This is made possible by using the fact that these connections act as a perturbation [19] of the dynamics generated by the local connections. Looking at the spontaneous activity, i.e. without thalamic input, is motivated by two reasons. The first reason is that it has been argued that the hallucinating patterns generated by some drugs can be explained by the spontaneous activity of networks [7,8] similar to the one studied here. The second reason is that if we were to include an input with amplitude that we vary, it would be the same as working with a fixed amplitude but varying the slope of the sigmoid function S [21]. Hence, looking at the case of a thalamic input can be thought of as a deformation of the case considered here.
The paper is organized as follows. We first describe in Sect. 2 the Turing patterns that appear when the long-range connections are switched off. In Sect. 3 we compute the network symmetries that are induced by the long-range connections. We then study in Sect. 4 the perturbation of the Turing patterns, before drawing some conclusions in Sect. 5.

Turing Patterns in the Unperturbed Case
In this section, we formulate our assumptions on the symmetries of the unperturbed solutions, 1 written V 0 (x), that will be analyzed upon switching on the long-range connections. Since the computation of Turing patterns has been extensively documented in the literature (see for example [22,23] for a review and [7] for neural fields), we will not cover it here and only briefly recall the main results. We have seen in the introduction that, when LR = 0, the system is invariant under any translation in the periodic domain Ω. Adjusting Ω to be a square or an hexagon, branches of steady patterns with square or hexagonal symmetry bifurcate from the homogeneous state as the slope σ of the sigmoid function S reaches a threshold value. These patterns are stable in some range of values of σ as shown in Fig. 1. Stability here means "orbital stability". Indeed any translation applied to such a state V 0 gives another solution, thanks to the translational invariance. The group of translations acting in Ω (which has periodic boundaries) is the torus R 2 /Z 2 , which, moreover, acts faithfully on any V 0 with square or hexagonal pattern. Therefore the set (called translation group orbit) of translated states from V 0 is a torus manifold. Moreover, this torus is an attractor if V 0 is orbitally stable.
A central hypothesis of the current study is the assumption that the spatial period of V 0 matches the spatial period of the PO map. The experimental data that support this assumption were already mentioned in the introduction. For example, we will not study the approximate "hexagonal" patterns (in the square cortex as in Fig. 1) because they do not satisfy our assumption. Also, the stripes 2 patterns will not be studied here. Indeed, they have less symmetries and their study would only require minor modifications of the present argument.

Network Symmetries
For reasons which have been explained in the introduction we are interested in the perturbation of Turing patterns with square or hexagonal symmetries. Such patterns have maximal lattice symmetry, meaning that they are invariant under the dihedral group of rotations/reflections D k with k = 4 (squares) or 6 (hexagons), and of course they are invariant under the discrete translations which define the lattice group. It is therefore convenient to suppose that the domain Ω representing the primary visual cortex is a square or a hexagon with periodic boundary conditions (opposite sides are identified). Let e 1 , e 2 be two vectors generating the lattice and Ω 0 be a fundamental domain centered at the origin O of Ω (see Fig. 2). We choose e 1 and e 2 as unit vectors, making an angle of π/2 in the square case, and of π/3 in the hexagonal case. Then being an arbitrary length scale. By construction there exists an integer N such that, by periodicity, There are (2N + 1) 2 copies of Ω 0 in Ω. We further assume for convenience that = 2π . The number N will be adjusted to match the number of elementary cells produced in Ω by the bifurcation of patterns when only local connections are active. The largest group of isometries leaving the lattice spanned by Ω 0 invariant is D k Z 2 2N (k = 4 in the square lattice case and k = 6 in the hexagonal lattice case). Let us now describe how the pinwheels are distributed in Ω. Pinwheels are the singular points of the PO map at which no orientation is preferred. Together with the orientation map they introduce an anisotropy in the lattice, reducing the symmetry group D k Z 2 2N of the lattice to a crystallographic subgroup (see [24,25], which is called, in the two-dimensional case, a wallpaper group. 3 See [27] for a description of these groups and for the shortened crystallographic notation, which is used in the following. It is well known that there are 17 such groups, up to isomorphism, see Fig. 3 for an illustration of these groups and the corresponding tilings. However, there are restrictions on the number of possible patterns of pinwheels. First, we only consider square and hexagonal lattices. Second, in order to be biologically plausible the patterns built from the orientation map must be continuous except at the pinwheels. We have determined by inspection that over all possibilities, only those five patterns highlighted in Fig. 3 are compatible with these constraints. Among those five we have selected the two which are shown in Fig. 2, i.e. pmm and p3m1, because they best agree with the experimental data. Pinwheels are the points (black/white dots in the pictures) at which all colors meet. Each color in Fig. 2 defines an orientation, which we identify with an angle between −π/2 and π/2. In the case (a) of Fig. 2 the tiling is invariant by reflection across the thick black lines (which also give invariance by rotations of angle π around their intersection points). In the Fig. 3 The 17 planar tilings. Each elementary cell represents an hypercolumn with a pinwheel at its center. The colored edges represent a coarse preferred orientation domain. The names of the five biologically plausible tilings are highlighted in red and a more accurate representation is shown on the right. The names of the tilings (and the associated wallpaper groups) are the shortened crystallographic notations; see [27]. pmm, cmm, and p3m1 are those that best agree with the experimental data case (b) the thick black lines are also axes of reflection. Observe that this gives 3-fold rotational symmetries around their points of intersection. The corresponding group is named p3m1. In the following we shall refer to these two PO maps by the names of the associated wallpaper groups.

Remark 2
We centered the hypercolumns around the pinwheels and used a uniform preferred orientation density such that the proportion of neurons coding for each orientation is the same for each orientation. This is done in order to avoid any preference of the network for a particular orientation. Note that pinwheels come in pairs called counterclockwise and clockwise: at a counterclockwise (respectively, clockwise) pinwheel, the orientation map shows increasing orientations when moving counterclockwise (respectively, clockwise) around it. Finally notice that the fundamental domain is either centered at a pinwheel ( Fig. 2(a)), or it is not ( Fig. 2(b)).
To be rigorous the set of orientations should be identified with the projective line P 1 (R). We can further identify P 1 (R) with the interval (− π 2 , π 2 ]. Let P be the set of pinwheels in Ω, then the PO map is a map Note that P 1 (R) is diffeomorphic to the circle S 1 (−π, π] mod 2π through the map (− π 2 , π 2 ] ϕ → 2ϕ. This allows one to naturally define an action of the circle group From this, we can define an action of S 1 on the map θ as follows. Let φ 0 = π/2 in the square case and φ 0 = 2π/3 in the hexagonal case and let R p φ be the rotation of angle φ centered on a pinwheel p. Then, with the PO maps of Fig. 2, where = 1 if the pinwheel is counterclockwise and −1 otherwise. Similarly, let K be the reflection w.r.t. a line passing through a pinwheel and parallel to the vector e 1 (horizontal axis). For the PO maps of Fig. 2, This defines an action of the reflection K on the map θ . In order to transfer these relations as (possible) symmetries of the network, we use the action of D k Z 2 2N on L 2 (Ω, R): This allows one to write (4) as (R Remark 3 Note that we could add an arbitrary value θ 0 ∈ (− π 2 , π 2 ] to the orientation map and still obtain an orientation map. Except when θ 0 = 0 or ±π/2, this would have the effect of breaking the reflection symmetry of θ , because applying the new reflection would not preserve the lattice. In the case of the tree shrew, the cortical coordinates (here assumed to be equal to visual field coordinates) must be such that the zero level of the PO map is parallel to the x 1 -axis, because of the anisotropy of the long-range connections [28]. In this case, the reflection K acts indeed on the network. On the other hand, for the macaque, any value of θ 0 is relevant because the long-range connections are approximately isotropic [29]. Hence in this case the assumption that K acts on the network is somewhat artificial. Nevertheless the presence of a reflection symmetry has strong consequences on the dynamics and we shall subsequently consider both the reflection and the non-reflection cases.
We now draw an easy consequence of Eq. (4), which will play a major role in the following analysis. Lemma 3.1 Let αe 1 , α > 0, be the vector between two closest pinwheels of different types in the direction e 1 . α is the smallest distance between them. Then (4) implies that Proof Let us write R p cc φ 0 (respectively, R p c φ 0 ) for the rotation of angle φ 0 around a counterclockwise (respectively, clockwise) pinwheel. According to (4) . This gives t = α(e 1 + e 2 ) and T t · θ = θ + φ 0 .

Model and Symmetries of the Long-Range Connections
In macaques, the anisotropy of horizontal connections follows from the anisotropy of the visual field representation in V1, i.e. it is not correlated to a feature like orientation [29]: the patchiness of the connections is hence isotropic. It seems reasonable to assume that this patchiness comes from the connections between populations of neurons which share similar preferred orientations [28]. Our long-range connectivity model reads where J 0 (χ, x) = e −[(1−χ) 2 x 2 1 +x 2 2 ]/2σ 2 LR . When χ = 0 the connectivity is isotropic, while when χ = 1 it is "the most" anisotropic. G σ θ is a centered Gaussian 4 with variance σ θ ≈ 35 • . It produces inhomogeneous patchy connections.

Remark 4
The following model of long-range connections was introduced by Bressloff [12]: where H is the Heaviside step function and P 0 is some constant. It differs mainly 5 by the first factor, which enforces homogeneity of the connections. It is not clear to us that this factor in (8) enforces connections between neurons with similar preferred orientations. Establishing which one of the two factors proposed in (7) and (8) best agrees with observations requires one to perform more experiments.
We now turn to the examination of the invariance properties of the long-range connections with respect to the symmetries of the lattice of pinwheels and PO map. It is important to emphasize that the symmetries of the long-range connections differ from the symmetry groups pmm or p3m1 of the PO maps as we shall see below. The function , y) for all x, y. This implies that the equations are equivariant with respect to the action defined in the previous section for the transformation γ . We consider successively the cases of square and hexagonal lattices. However, there are some general features which are shared by both types of networks and we start by stating them.
Proof From (4), the factor G σ θ (θ (x) − θ(y)) is clearly unaffected by the rotations. Now thanks to (4) and the expression R p cc Note that in the square case the above result still holds for = −1 because 2φ = π in this case and J 0 is an even function. Proof The factor G σ θ (θ (x)−θ(y)) is clearly unaffected by the reflection. Now thanks to (5),

Lemma 3.4 Let t be the vector as in Lemma
Proof We can rewrite Lemma 3.1 as T −1 t · θ = θ − φ 0 . Hence, we have γ · θ = θ − φ 0 2 . We now look at the effect on the long-range connections. The factor G σ θ (θ (x) − θ(y)) is clearly unaffected by the transformation γ . Now, thanks to the previous relation: Applying γ to the previous equation gives the γ invariance.
We will use these lemmas to derive the generators of the symmetry group in the square and hexagonal cases. We will also compute the specific crystallographic group that they generate.

The Square Case
We are now in a position to derive the symmetry group of the network equations. Note that the case of the Bressloff connectivity (8) leads to a different symmetry group. We start with the generators. Proposition 3. 5 We write R p c φ 0 the rotation of angle π 2 centered on a clockwise pinwheel. For the (pmm) PO map in Fig. 2(a), the symmetry group associated with the connectivity (7) in the case χ > 0, LR = 0 is: Finally, the subgroup of translation symmetries is the lattice L π(e 1 + e 2 ), π(e 1 − e 2 ) .
Proof This is a direct consequences of the lemmas in the previous section. There is, however, a simplification in the square case, namely the translation T π(e 1 +e 2 ) is a symmetry. Indeed, using Lemma 3.1, it is straightforward to show that J LR (T π(e 1 +e 2 ) x, T π(e 1 +e 2 ) y) = J LR (x, y). Note that the translation T π(e 1 +e 2 ) commutes with T 2πe 1 and R p c φ 0 . Finally, let us show that T 2πe 1 is generated by the group elements listed in the lemma. Indeed, T π(e 1 −e 2 ) = (R p c φ 0 ) −1 T π(e 1 +e 2 ) and T 2πe 1 = T π(e 1 −e 2 ) + T π(e 1 +e 2 ) .
To prove that the lattice of translations of the symmetry group is L[ π 2 (e 1 + e 2 ), π 2 (e 1 − e 2 )], we use the relation which is a translation. Using the rotations with axis a = 0 and b = πe 1 allows one to conclude.
The previous result gives the generators of the group. These groups are very well known as crystallographic groups in the literature (see [24,25,27] for an introduction). More precisely, we find that the symmetry group Γ of the equations, in the case χ, LR = 0, is the crystallographic group P4M if θ 0 ∈ π 2 Z and P4 otherwise. It is interesting to note that D 2 is the point group associated to the pmm PO map, whereas the point group of the network equations can be D 4 .

The Hexagonal Case
The main difference with the square case is that the clockwise rotation is not a network symmetry because of the anisotropic function J 0 in (7) when χ > 0. Hence, only half of the pinwheels are centers of rotation for the network equations (1), namely the counterclockwise pinwheels. A direct consequence of the lemmas in the previous section is the following. Proposition 3. 6 For the (p3m1) PO map in Fig. 2(b), the symmetry group Γ associated with the connectivity (7) in the case χ > 0, LR = 0 is:

Finally, the subgroup of translation symmetries is the lattice
Proof The only thing left to prove is about the lattice of translations. We can look at (see Lemma A.1) As for the square lattice, we can identify the group in the hexagonal case. The symmetry group Γ of the equations, in the case χ, LR = 0, is the crystallographic group p3m1 if θ 0 ∈ π 2 Z and P3 otherwise.

Forced Symmetry Breaking of Patterns
We now study the perturbation of the torus of translated states from V 0 when it is an attractor. It is well known that in this case, a torus, flow-invariant manifold persists when the equations are perturbed, as long as the perturbations are small (and smooth) enough [30]. Our aim in this section is to determine how many steady-states do actually persist when the system is perturbed by turning on long-range connections LR = 0, and which phase portrait is observed on the perturbed torus.
Our method is as follows. We first analyze the remaining symmetries on the perturbed torus when LR = 0. This allows one to assert the persistence of some steadystates (equilibria) which stand for points of maximal isotropy for the action of Γ on the perturbed torus. Moreover, when these isotropy subgroups contain the m-fold rotation group C m with m ≥ 3, these equilibria are foci (attractive or repulsive) or nodes (the case when the Jacobian matrix has a double real eigenvalue). Now, the topology of the torus is an important constraint for the distribution of equilibria of saddle and other types. This follows from the Poincaré-Hopf theorem, which can be found in an abundance of literature and textbooks, and which can be stated as follows: In our case V is the torus, and hence χ = 0. Therefore there are an equal number of equilibria of saddle type and of non-saddle type. These two pieces informations (symmetries and Poincaré-Hopf theorem) greatly help to classify the possible phase portraits on the torus. The idea of analyzing the dynamics on a group orbit of equilibria under symmetry-breaking perturbations was introduced by [16]. It was applied in a theoretical setting to the case when the group orbit is a torus of square patterns [31] or hexagonal patterns [32] with the aim of showing the existence of robust heteroclinic cycles under certain conditions. Our aim here is different and we only focus on the cases which correspond to our model. We next consider the square and hexagonal cases and we list the simplest possible phase portraits, that is, with the minimal number of equilibria, assuming that these are the pictures which will naturally arise in our model, unless further degeneracies are assumed. Then direct numerical simulations of the dynamics on the perturbed torus allow us to fix the actual dynamics which is induced on the invariant tori when switching on the long-range connections. Then direct numerical simulations of the dynamics on the perturbed torus allow us to select the dynamics, among the ones we predicted, which is induced on the invariant tori when switching on the long-range connections.

On the Perturbed Torus
As the perturbed torus T LR is diffeomorphic to the unperturbed torus T 0 , we can work in the coordinates of T 0 for which the group action expression is known (see below). In the numerical experiments, the diffeomorphism is not known and thus, as we use T 0 coordinates, the projection of the patterns on this coordinate system is only approximative (see Fig. 4, for example). To fix ideas, it is useful to be a bit more explicit. When LR = 0, we write the unperturbed torus as where V 0 is a stationary solution for LR = 0. We assume that the torus solution is invariant by rotation and reflection which is equivalent to assuming V 0 invariant by rotation and reflection. Under this assumption, the actions of rotations/reflections on the torus satisfy This follows from R o T t · V 0 = T R o t V 0 . It implies that the action of a rotation R a of axis a is given by In our model, the lattice of translations symmetries of the torus matches the one of the PO map (see Sect. 2). Hence, V 0 has the pinwheel periodicity: If we identify V 0 (· − t) and t, we can further simplify the study of the perturbed torus by decomposing t as follows: The assumption (10) implies that φ i ∈ [0, 2π).

Remark 5
We cannot apply directly our method to the branch of stripes in Fig. 1, nor to the branch of hexagonal patterns, because the unperturbed torus generated by these patterns is not invariant by rotations.

Square Case
Using the decomposition t = φ 1 e 1 + φ 2 e 2 , one finds: We collect the main results concerning the dynamics on the square in the following proposition. It is the backbone for determining the possible phase portraits on the perturbed torus.
which are centers of rotation.
Proof It is easy to prove (12). Fixed point subspaces are flow invariant, this implies that Fix( (R o ) 2 ) is composed of stationary solutions. We also note that Fix R o = (0, 0), (π, π) .
From γ being affine and (0, π) ∈ Fix(γ ), we find that dF (0, π) commute with dγ = dR o seen as a map on the torus. This allows one to conclude that (0, π) is a node/focus, and also (π, 0). Being fixed points of rotations, the four nodes/foci are center of rotation symmetry. Assume now that there are a finite number of zeros (φ i ) i=1,...,n on T LR which are all non-degenerate. Thanks to Theorem 4.1 n i=1 sign det dF (φ i ) = 0.

This gives
which implies the existence of at least four saddles.
A convenient way to find the foci is to look at the fundamental domain in Fig. 2(a). These foci corresponds to the pinwheels in the fundamental domain.
We have seen that the minimal configuration, under the hypothesis that all equilibria on the perturbed torus are hyperbolic, is that of eight equilibrium points with four foci and four saddles. How does the associated phase portrait look like on the torus? The answer crucially depends upon the presence of the reflection symmetry (case when θ 0 = 0). In this case the axes of reflection symmetry go through the equilibria. Therefore the foci are necessarily of node type: the eigenvalues of the Jacobian are double and real. The axes of reflection symmetry are invariant by the flow, which constrains the phase portrait to look like the one shown in Fig. 4. On the other hand when there is no reflection symmetry the foci are typically "true" foci: the eigenvalues of In order to observe limit cycles numerically, we had to change the connectivity. Indeed, if we use the prefactor G σ θ = cos in (7), the imaginary part of the eigenvalues of the equilibria coming from the breaking of the reflection symmetry by choosing θ 0 / ∈ π 2 Z is tiny: at least 3 orders of magnitude smaller than the real part. In effect, even if we break the reflection symmetry, we observe a dynamics similar to the one in Fig. 4. To have larger imaginary parts, we connect neurons with opposite preferred orientation by choosing the prefactor G σ θ = sin in (7). We further choose the connectivity with largest imaginary part among connectivities in Appendix B. Note that despite varying almost all parameters, we only observed the two situations as in Figs. 4 and 5 (up to a time reversal of the absence of limit cycles), as if naturally, the network equations (1) produced the simplest possibilities for all parameters that we investigated.

Remark 6
We would like to mention that great care was taken to code the equivariance and that numerically, the error on symmetries was around 10 −16 for the 2-norm of an arbitrary vector of 2-norm around 37. The numerical errors on the equivariance relations comes mainly from the PO map θ (see Eqs. (4)-(6)). Therefore, we computed the PO map by first building its fundamental domain by rotating/reflecting a basic cell and then padding this fundamental domain to cover Ω. This numerical accuracy of the equivariance allows one to check the predicted values of the stationary points with great accuracy using a Newton algorithm. In particular, we find numerically in Fig. 4 that the points (R o ) k · ( π 2 , 0), k = 0, . . . , 3 are indeed saddle points.

Hexagonal Case
We now consider the hexagonal lattice. This case is different from the square lattice because it seems from Proposition 3.6 that only counterclockwise pinwheels are center of rotations. However, it turns out that γ as in Lemma 3.4 is a rotation on the torus, hence yielding the other pinwheels as center of rotations. We prove the next proposition using a different method from the one used in the proof of Proposition 4.2. Using a decomposition t = φ 1 e 1 + φ 2 e 2 , one finds Proof In the fundamental domain, only counterclockwise pinwheels lead to a rotational symmetry. The center of rotation is then a node/focus point. In particular, we find the following node/foci points (see Sect. 3.1.2 for a definition) We now look at the symmetry γ ≡ T −1 2π 3 (e 1 +e 2 ) R p c defined in Proposition 3.6, where the axis of rotation is p c = 2π 3 (2e 1 +e 2 ). On the hexagonal torus, we find that γ = R o , which can be seen by writing the equations in the basis e 1 , e 2 . It yields Fix(γ ) = ( π 3 , π 3 )Z. Hence, these points are equilibria of node/foci type. It gives three additional nodes/foci.
We can use these centers of rotation to rotate each node/foci in order to find other equilibria. Using Lemma A.2, we have 3 (e 1 + 2e 2 ) and u = 4π 3 e 2 . This yields the additional centers of rotations (hence node/foci): It follows that there are (at least) nine foci. Assuming that there is a finite number of zeros on T LR which are all non-degenerate, Theorem 4.1 implies that there are as many saddles as foci.
We take the opportunity to show how the different nodes/foci found in Proposition 3.6 can be interpreted in cortical coordinates, as shown in Fig. 6. Briefly, we add the cortical activity as a semi-transparent overlay (transparent is when the cortical activity is high) to the PO map θ . We then plot small edges or hexagonal patches depending on which a subset of the hypercolumn is activated.
As in the square case, we can deduce from these results the possible phase portraits when assuming that the equilibria on the perturbed torus are all hyperbolic and that Fig. 7 Sketch of the possible phase portraits on the perturbed torus. One case (left) is shown, assuming that a periodic orbit surrounds the focus at the center. Another case (right) is shown in the case of reflection symmetry. The presence or absence of periodic orbits is essentially related to the sign of the real part of the eigenvalues of the Jacobian at the central focus there are exactly 18 of them, nine foci and nine saddles. The situation is slightly more complicated than in the square case, but it is not difficult to show that a typical phase portrait looks like one of the diagrams shown in Fig. 7. It is numerically checked that this "minimal" situation indeed occurs.
As in the square case, we had to use G σ θ = sin in order to see foci with nonvanishing rotation number. Compared to the square case, we had more difficulty to keep small enough errors in the equivariance relations. This numerical error is 3 orders of magnitude bigger than for the square lattice 6 and given that we need to numerically solve (1) for a very long time, the errors of the symmetries seem to build up. Nevertheless, we were still able to produce simulations, corresponding to one of the possible diagrams (see Fig. 7). Except for the points ( 2π 3 , 2π 3 )Z, we find the stationary points predicted by Proposition 3.6 using a Newton algorithm. Around the points ( 2π 3 , 2π 3 )Z, the Newton algorithm does not converge: this seems to be caused by the presence of saddle points (also predicted in Proposition 3.6 which produces the small kinks on the trajectories around the red points (see Fig. 8)).
In the numerical simulation displayed in Fig. 8, no periodic orbit was found. The simulation seems to correspond to the predicted dynamics shown in Fig. 7, middle. On the other hand, we observe the remarkable fact that the simulations also lead to simplest possible scenarios with the smallest number of stationary points and a periodic solution.

Discussion
In this work, we have extended the seminal and very influential work [8] on cortical hallucinations after the original work of Ermentrout and Cowan [7]. Our idea is to assume a discrete lattice of pinwheels (as opposed to a continuous one) and to describe the cortical activity in the laminar zone which is located in between pinwheels. This is important, for example, if we want to confront our predictions to optical imaging experiments. This approach has led us to shed some new light onto two outstanding questions. First, what are the possible lattices of pinwheels? They turn out to be closely related to the wallpaper groups. Second, what is the simplest spontaneous dynamics generated by these networks? They turn out to be determined by the perturbation of invariant tori.
The first question is natural but, to the best of our knowledge, has never been addressed theoretically, despite the fact that it allows one to apply the equivariant theory of dynamical systems in a very similar way to [8] albeit in a more biologically plausible setting. In [33] the authors describe a mechanism that allows them to describe the probability of observing a network of pinwheels with a given density, but because their equations are driven by white noise they lose all symmetries.
The second question is more subtle and in effect stems from numerical work where we found it very difficult to implement the ideas of [8,12], at least for a square cortex. The computation of the bifurcation diagram in Fig. 1 led us to study the perturbation of solutions that were not close to bifurcation points (see the brown lines in Fig. 1), indeed this accounts for the most probable dynamics. As such, the bifurcation diagram shown in Fig. 1 is an indication of the difficulty to apply the theory of cortical hallucinations developed in [8] where it was assumed that such stereotyped cortical patterns could be explained by adjusting the network parameters close to bifurcations. Indeed, in such a setting, the validity of the theory shrinks rapidly with the size of the cortex and, as far as we know, for all practical purposes it is difficult to use it to account for observations. The problem is the same as in the work described in [12], where the author studies the perturbation of a system close to a bifurcation with a spatial forcing close to resonance. Secondary bifurcations might seriously restrict the validity of the approach.
Let us examine the consequences of our investigations onto our current understanding of the functioning of V1.
The first consequence can be drawn from Fig. 1; the hexagonal case is the "only" robust one from a modeling point of view. Indeed, even with the square lattice, the branches which are stable over an extended range of parameters, have a near hexagonal symmetry. Hence, there is a mismatch between the solution approximate sym- The second consequence is that, within the class of pinwheel lattices that we consider, those displaying the reflection symmetry are non-generic (see Sect. 3). This induces the presence of foci and (possibly) limit cycles as described in Propositions 4.2 and 4.3. However, the imaginary part of the eigenvalues at the foci is very small for the excitatory long-range connections that are biologically plausible [28,29] and this makes the observation of limit cycles difficult. 7 To increase the imaginary part in simulations and observe limit cycles (see Fig. 5), we artificially use a connectivity that connects neurons with orthogonal preferred orientations e.g. G σ θ = sin. This is different from [8] where all bifurcated patterns are stationary. Note that [11] reports time-periodic states in model similar to [8] but with additional symmetry.
The last consequence of our work is that it generalizes both theories described in [7] and [8]. Indeed, the first theory allows the description of non-contoured cortical spontaneous activations like in Fig. 9(a)-(b). The second theory generalizes the first one by allowing the prediction of contoured patterns such as those in Fig. 9(c) as well as the non-contoured ones. Our theory describes contoured and non-contoured patterns as well as a mixture of such activation patterns as depicted in Fig. 6(c).
Our work can be extended in several ways. First, there is a need for a bifurcation study of the dynamics on the torus w.r.t. the parameters LR , χ , etc. Indeed, we have not been very quantitative concerning the role of these parameters in shaping the dynamics. There are also lattices that have not been investigated (cmm, p4m and p6m in Fig. 2) and it would be interesting to study the kind of planforms they can produce. We have not looked at the case where the stripes are stable in the D 4 -pitchfork. That would be a minor modification of the present work but it would show another type of planforms. More generally, we have not discussed the cases where the unperturbed solution is a circle rather than a torus. Hence, it is desirable to classify the different planforms that can be produced from the unperturbed invariant manifolds for the different lattices. Another extension concerns the study of cases where the lattice of symmetry of the unperturbed torus T 0 and the pinwheel lattice differ. The perturbation from the long-range connections would act as a periodic forcing on the unperturbed torus. Some models are available [12,13], but we believe that they lack an important component: synaptic plasticity. A relatively simple extension would be to consider the effect of synaptic/propagation delays [34,35]. Synaptic delays [36] will not affect the unperturbed torus but are likely to increase the imaginary part of the eigenvalues which we found to be nonzero albeit small when the reflection symmetry is broken. Other very exciting extensions concern the modeling of the spatial frequency tuning and the ocular dominance domains. It would be very interesting to re-visit some recent work on binocular rivalry [37,38] in the light of the conclusions presented in this paper.