Spatial and color hallucinations in a mathematical model of primary visual cortex

We study a simplified model of the representation of colors in the primate primary cortical visual area V1. The model is described by an initial value problem related to a Hammerstein equation. The solutions to this problem represent the variation of the activity of populations of neurons in V1 as a function of space and color. The two space variables describe the spatial extent of the cortex while the two color variables describe the hue and the saturation represented at every location in the cortex. We prove the well-posedness of the initial value problem. We focus on its stationary, i.e. independent of time, and periodic in space solutions. We show that the model equation is equivariant with respect to the direct product G of the group of the Euclidean transformations of the planar lattice determined by the spatial periodicity and the group of color transformations, isomorphic to O(2), and study the equivariant bifurcations of its stationary solutions when some parameters in the model vary. Their variations may be caused by the consumption of drugs and the bifurcated solutions may represent visual hallucinations in space and color. Some of the bifurcated solutions can be determined by applying the Equivariant Branching Lemma (EBL) by determining the axial subgroups of G . These define bifurcated solutions which are invariant under the action of the corresponding axial subgroup. We compute analytically these solutions and illustrate them as color images. Using advanced methods of numerical bifurcation analysis we then explore the persistence and stability of these solutions when varying some parameters in the model. We conjecture that we can rely on the EBL to predict the existence of patterns that survive in large parameter domains but not to predict their stability. On our way we discover the existence of spatially localized stable patterns through the phenomenon of"snaking".


Introduction
Neural Fields are a useful mathematical formalism for representing the dynamics of cortical areas at a macroscopic level, see the reviews [9,19,24].This formalism has been broadly used to account for the observed activity of the cortical visual area V1.V1 receives inputs from the retina through the LGN.There is a massive feedback from V1 to the LGN.V1 also sends inputs to higher level cortical visual areas such as V2 and V4 and receives feedback signals from them.For a very accessible introduction to the various visual areas, the reader is referred to the book of David Hubel [35], and to [38] for a more recent presentation.
Neural Fields models of V1 are sets of integro-differential equations whose solutions are meant to describe its spatio-temporal activity.The well-posedness of these equations has been studied in depth by various authors [3,27,49,65] with a special attention to the stationary solutions, i.e. those which do not depend upon time.These solutions, also called persistent states, are interesting because they appear to provide good models of memory holding tasks on the time scale of the second [18,28,42].Moreover, they appear to resonate with the fascinating phenomenon of visual hallucinations and their relation with the functional architecture of the visual cortex [10,25] .
With no exception to our knowledge, all previous work in neural fields theory has not taken into account the chromatic aspects of visual perception.In [57], we introduced a neural field model for color perception to explore how the synergy of two antagonistic phenomena, simultaneous contrast and chromatic assimilation, could lead to a "color sensation".
In the present article, we address the questions of how this model can predict visual hallucinations and what are their spatial and chromatic structures.We are guided in this venture by our previous work [64] and make good use of the theory of equivariant bifurcations.
It is structured as follows.In Section 2 we recall the neural field model described in [57] and prove its well-posedness, Section 3 introduces the notion of stationary solutions and their bifurcations.Section 4 is dedicated to the computation of the spectrum of the linear operator in the neural field model.Section 5 describes the symmetries of the model and uses the equivariant branching lemma to predict the type of bifurcations at the primary bifurcation points and the shape of the bifurcated solutions, or planforms.Section 6 shows examples of such planforms.Section 7 goes away from this local analysis and explores numerically, thanks to the development of innovative software, a much larger volume of the set of stationary solutions to the neural field equations.We conclude in Section 8.

The model
In this work we think of V1 as the closure of a regular domain, noted Ω s , s stands for space, of R 2 , in effect the open square (−l /2, l /2)×(−l /2, l /2), where l is a positive number which, for simplicity and without loss of generality, we take equal to 1, except for some of the numerical experiments presented in Section 7.
The visual cortex is organized into hypercolumns, groups of neurons sharing the same receptive field in the retina and coding for specific physical quantities such as edge orientation, spatial frequency, temporal frequency.These signals are mapped from the retina to V1 following an approximately log polar retinotopic transformation (see Remark 21).Unlike in the case of orientation, for which the existence of such hypercolumns in V1 is now well established [58], the anatomical and physiological bases for a functional architecture encoding color are still debated.These bases are most likely connected to the presence of blobs [34,35].Hence, in light of the promising findings made by [13,70] it is reasonable to assume in our work a hypercolumnar organisation of cells tuned to a continuum of colors.Our work also supposes the presence of long-range lateral connections between hypercolumns, in agreement with observations of [41] where horizontal connections tend to link blobs to blobs.Note that this visual information is stored in the cortex in three dimensions, i.e. the cortex has some thickness.We neglect in our work this thickness and consider only its spatial extent.
We now briefly recall the model described in [57].It is based on an opponent representation of colors such as Hering's opponent space [32].In this setting, a color is a pair c := (c 1 , c 2 ) of real numbers which encode the chromaticity of the color 1 .Details are provided in Appendix A. The reader can think of c 1 as encoding the yellow-blue colors and c 2 as encoding the red-green colors in Hering's theory.What is important for us is that the set of chromaticities is symmetric w.r.t the origin, i.e. if c = (c 1 , c 2 ) is a chromaticity, then −c = (−c 1 , −c 2 ) is also a chromaticity, called the opponent chromaticity or color of c.The set of chromaticities is therefore a bounded regular domain of R2 which is symmetric w.r.t the origin, in effect the open disk D(0, c 0 ) centered at the origin and of radius c 0 , where c 0 is a positive number which we take without loss of generality and for convenience equal to 1.

Remark 1.
We implicitly assume that the topology of the chromaticity space is that of the Euclidean plane.This is only a coarse approximation.Note that the problem of defining a metric in color space is still open [6,40,50,61,69].

Remark 2.
The elimination of the semi open radius [0, 1) × {0} from D(0, 1) is practically not important since the functions that we will manipulate, in particular the eigenfunctions of the operator W c , see below, will be smooth and therefore can be continuously extended to the closed disk D(0, 1).
We define to be the bounded domain, in effect the open rectangle of R4 , encapsulating the spatial and chromatic coordinates that will be of interest in the sequel.

Connectivity kernel
Putting all this together, at each point (r, c) = (r 1 , r 2 , ρ, ϕ) of Ω, we consider a neural mass 2 whose average membrane potential is noted V (r, c, t ).It is a function defined on Ω × J , J being an interval of R containing 0. In [57], we assumed that the function V was the solution to an initial value problem3 related to a Hammerstein equation 4 which writes together with the initial condition V (r, c, 0) = V 0 (r, c).This equation describes the time variation of the scalar function V (t ) defined on Ω starting from the initial condition V 0 .At each time t , V (t ) belongs to some functional space, in effect a Hilbert space F , that we describe in the next section.
We now discuss the various elements that appear in this equation.τ is a time constant that defines the speed of the exponential decay toward the initial condition.Without loss of generality we can assume τ = 1.
Υ is a sigmoidal function mapping R to the open interval (0, 1).It is called the activation function, relating the values of the membrane potential V to the neuronal activity a (0 meaning quiet, 1 meaning highly active).It writes ( ε is a parameter that allows us to shift the origin, γ controls the slope of the sigmoid at the (shifted) origin, it is often called the nonlinear gain.I ext is a function representing the input to the neural mass from different brain areas.In the remaining of this paper we take I ext = 0, i.e. we consider that area V 1 is isolated from the rest of the brain.This is clearly an approximation but allows us to do some mathematics and is a first step toward the analysis of the general case.
w is the connectivity kernel 5 .It models the influence of neighboring neural masses at (r , c ) on the neural mass at (r, c) as a (separable in space and color) linear superposition operation where the index s stands for "space" and the index c stands for "chromaticity".If w(r, c, r , c ) is positive (respectively negative) the neural mass at (r , c ) excites (respectively inhibits) the one at (r, c).The product of w s and w c is intended to model the antagonistic effects of color assimilation and contrast which are parts of the class of perceptual phenomena called chromatic interactions [66].
w s is a "classical" two-dimensional "Mexican hat" function, see [9] and Figure 7-Left, which we write as the difference of two circularly symmetric Gaussians centered at 0: where 2 is the usual Euclidean L 2 norm in R 2 .Biology dictates that α s , β s are very small w.r.t 1, the size of Ω s .This indicates that our model only takes into account the neural connections which are local to the visual area V1 and take place within the gray matter while it is known that different visual areas communicate through the fiber bundles (in the biological sense) forming part of the white matter.This would be part of the term I ext in (1) which we have taken to be equal to 0.
If w s (0) > 0, i.e. if µ s > ν s , the neural masses at r such that r − r 2 is small enough excite the neural mass at r , and if α s < β s those sufficiently far away inhibit it.
We furthermore assume i.e. that the spatial excitation and inhibition are balanced.This is both compatible with some biological evidence (balanced networks [20]) and mathematically convenient.
The product w of w s and w c behaves as a function of r, r , c, c as shown in Table 1 adapted from [57].As pointed out in this paper, this is in qualitative agreement with the nonlinear behaviour of color shifts found in [43,44].
Note that by definition of w s and w c the function w is symmetric w.r.t.r and r , and c and c , respectively.

Table 1. Sign of the connectivity kernel
Choice of the appropriate functional space F and well-posedness of (1) Our choice of F is guided by three criteria: (1) The well-posedness of the problem, (2) Its biological relevance, (3) Its suitability for numerical computations.
The choice of a Hilbert space is appealing and a natural choice is F = L 2 (Ω).As argued in [65], this unfortunately allows the membrane potential to be singular since for example the function It is desirable that the average membrane potential stays bounded on the cortex and a way to achieve this is to allow for more spatial and chromatic regularity by assuming that (r, c) → V (r, c) is differentiable almost everywhere.
The choice of the Sobolev space F = W m,2 , is convenient for two reasons.First it is a Hilbert space endowed with the usual inner product: where the multi-index α is a sequence α = (α 1 , α 2 , α 3 , α 4 ) of 4 integers and |α| = 4 i =1 α i , and the symbol D α represents a partial derivative, e.g., The second reason is that, because the boundary of Ω is sufficiently regular (it satisfies the coneproperty [1,2]), F is a commutative Banach algebra for pointwise multiplication [1][Chapter V, Theorem 5.23].This is necessary in the upcoming bifurcation analysis in order to apply the Equivariant Branching Lemma, see Proposition 17, for which we require some smoothness of w.
Remark 5.The reader probably wonders what is the value of m.In order to have the Banach algebra property, we need to have 2m > d , d being the dimension of Ω, i.e. 4. Hence the smallest possible value of m is 3.But in Section 7 we assume that Ω c is one-dimensional making d equal to 3 and the smallest possible value of m is equal to 2 in this case.
We next define the operator W as acting on L 2 (Ω) as follows.Let U ∈ L 2 (Ω), we define Remark 6. Remember that the measure d c is in effect equal to 1  2 d d ϕ .
It is clear that this is well defined and we have the following proposition.
Proposition 7. The operator W maps L 2 (Ω) to F and hence F to F .
Proof.The proof follows from Proposition 2.3 in [65].
We also have the following Proposition about the solutions to (1).
C. R. Mathématique-Draft, 20th October 2021 Proposition 8.For each V 0 ∈ F there exists a unique solution in C (R + , F ) to the following Cauchy problem where the function s is defined by (2).
Proof.The proof is a an immediate consequence of Proposition 2.5 in [65].

Stationary solutions and bifurcations thereof
In this paper, we focus on the stationary (independent of time) solutions to (11), the steadystates.They are important because they are thought to be good models of the memory holding tasks on the timescale of the second as demonstrated experimentally on primates [18,28,42].These solutions may change drastically when some of the parameters such as γ and ε in (2), µ s , ν s , α s , β s in (4), ξ c in (7), and µ c , ν c , α c , β c in (8) vary.In effect, we will concentrate on the first parameter, the nonlinear gain γ, which is important in determining the relation between neuronal activity (a number between 0 and 1) and the associated membrane potential V .It is known that the ingestion of drugs such as LSD and marijuana a) can change this relation and b) can trigger hallucinatory patterns [47,54].It is therefore very much worth our efforts to investigate if, when varying the parameter γ, stationary solutions to (11) do bifurcate and if the bifurcated solutions resemble some of the known visual illusions.
To summarize, we are going to study the bifurcations when γ varies of the solutions to the equation where the operator F is defined by In order to achieve such a task, we need to determine the spectrum of the operator W and the symmetry properties of the operator F with respect to some groups of transformations of Ω.We describe the spectrum of W in Section 4 and the symmetry properties of F in Section 5.
It is natural at this point to introduce the two operators W s and W c acting on L 2 (Ω) associated to the functions w s of (4) and w c of (6) as follows The separability of the function w in the space and color variables is reflected in the notation Similarly the separability of the function w c as the product of w c,m in (7) for the magnitude and w c,a in (8) for the angle implies that where the convolution in the right hand side of the definition of W c,a is a periodic convolution.

The spectrum of W in L 2 (Ω)
The reader can verify that, given the symmetry properties of the functions w c and w s , the operator W is symmetric in L 2 (Ω), i.e. satisfies Another important property of W is that it is compact, i.e. given any bounded sequence (U n ) n≥0 of functions in L 2 (Ω) the sequence (W • U n ) n≥0 contains a converging subsequence for the norm of F .This is a direct consequence of Lemma 2.4 in [65].
A consequence of the fact that W is compact is that its spectrum is compact and at most countable.Moreover each point in the spectrum, except perhaps 0, is isolated.All the non zero elements of the spectrum are eigenvalues and all eigenvalues are real since W is symmetric.
Since w is separable in space and color, the spectrum of W is defined by the spectra of W s and W c : the eigenvalues are the product of those of each operator and the eigenfunctions are separable in space and color, being the products of those of W s and W c .

The spectrum of W s in L 2 (Ω s )
Because we are studying the symmetry-breaking steady-state bifurcations of the solutions to our model, we know from previous work on heat-conduction in fluids [14] or on visual hallucinations [10] that this leads to the formation of spatially periodic patterns.
This brings into play the lattice L generated by the two vectors k 1 and k 2 where (k 1 , k 2 ) is the canonical basis of R 2 and commands that we quotient R 2 by L thus obtaining the 2-torus Ωs = R 2 /L .The dual lattice L * is also generated by the two vectors k 1 and k 2 .
We thus make the following assumption about the solutions to equations (1) which is inspired by biology with the advantage that it simplifies the mathematics: We now work on the Hilbert space L 2 ( Ωs ) of L -periodic functions with the same inner product as before.We note Ws the operator W s restricted to this space i.e. defined by Note that the spatial convolution in ( 3) is now a periodic convolution.Ws is clearly a symmetric compact operator on L 2 ( Ωs ).
It is easy to characterize the spectrum of Ws : Proposition 9.The eigenvalues λ s m,n , m, n ∈ Z, of Ws ∈ L 2 ( Ωs ) are the (real) Fourier coefficients of the L -periodic even function ws : where Because w s is even, so is the sequence of eigenvalues: For (m, n) = (0, 0) the eigenspaces of Ws are of even dimension, larger than or equal to 2. If (m, n) = (0, 0), because of (5), λ s 0,0 = 0 and the dimension of the kernel of Ws is 1.Given an eigenvalue λ s m,n the corresponding eigenspace is generated by the functions: cos 2π〈k p,q , r 〉 and sin 2π〈k p,q , r 〉, for all (p, q) such that λ s p,q = λ s m,n .
Proof.This follows for example from the fact that Ws is diagonalized by the Fourier Transform.
Define the operator τ c acting on the ϕ-periodic functions of period 1 as We have the following Proposition.
Proposition 10.The linear operator W c on L 2 (Ω c ) is symmetric and compact, it commutes with the operator τ c defined by (15).
Proof.Obvious from the definitions.
From this we have Corollary 11.The spectrum of W c is real and at most countable.The eigenfunctions are separable w.r.t., ϕ.
Proof.The first assertion follows from the symmetry and compactness of W c .The second follows from the definitions ( 6)- (8).
It remains to characterize the eigenfunctions and eigenvalues.
Proposition 12.The eigenspaces of W c,a are generated by the functions sin 2πnϕ and cos 2πnϕ, n ∈ N. Hence they are one-dimensional for n = 0, two-dimensional otherwise.The corresponding eigenvalues are the (real) Fourier coefficients of the 1-periodic even function given by (8) of index n.
Proof.This follows immediatly from the fact that the operator W c,a is a 1-periodic convolution and the 1-periodic function Regarding W c,m , we have the following Proposition.
Proposition 13.The eigenspaces of W c,m corresponding to the non zero eigenvalues are onedimensional.Its kernel is reduced to the null function.The eigenvalues λ c n , n ∈ N, are obtained from the countable solutions x n to the transcendental equation from the relation The eigenfunction e c,m λ n corresponding to the eigenvalue λ n = 0 is given by e c,m Proof.We have, according to ( 9) and ( 14) ), so that if u is an eigenfunction corresponding to the eigenvalue λ c we obtain C. R. Mathématique-Draft, 20th October 2021 which are of the form u( ); = ae K + be − K for some constants a and b, real or complex.λ c , a and b are determined by writing that the function u thus defined is indeed an eigenfunction of W c,m corresponding to the eigenvalue λ c .It can be verified that if K ≥ 0, i.e. if λ c ≥ 1 ξ c , there are no solutions.Hence we must assume λ c < 1 ξ c and hence K < 0. The solutions to (19) write A symbolic computation system shows that The two conditions are necessary and sufficient to guarantee that W c,m • u( ) = λ c u( ) for all 0 < < 1.Since a and b are not both equal to 0 (otherwise u = 0 and hence λ c = 0) we must have the eigenvalues are found by solving for x > 0 the equation For each value of x, the corresponding value of λ c is found to be decreases monotonically from +∞ for x = ξ + c to 0 when x → +∞.The corresponding curve therfore intersects the curve representing the π-periodic function x → tan x at a countable number of points x n , n ≥ 0, yielding the eigenvalues λ c n of W c,m .According to (20) and ( 21) the corresponding eigenfunction is

Symmetries of the model and equivariant bifurcations of the solutions
It follows from ( 5) that V = 0 is a solution to (12) and ( 13) for all values of γ.We are interested in the problem of determining how this solution bifurcates when γ increases, while allowing us to change somewhat the value of the threshold ε in (2).
We investigate the properties of the operator F defined by (13) under the action of the Euclidean group E (2) restricted to the lattice L defined at the beginning of Section 4.1 for the spatial part, and the group O(2) for the color part.This group arises from the rotations R c ϕ c acting on the ϕ angle by ϕ → ϕ + ϕ c mod 1 and the reflection τ c : ϕ → −ϕ mod 1.An element ψ of this group is of the form The action of E (2) on the space of L -periodic functions is best understood by considering separately the translations and the rotations.Since the translations of R 2 leave the set of Lperiodic functions invariant and translations in L fix all L -periodic functions, the effective action of the group of translations of R 2 is as the 2-torus T 2 = R 2 /L which is compact.For the rotations, recall that the holohedry H of the lattice L is the largest subgroup of O( 2) that leaves L invariant.In the case of a square lattice H = D 4 is the dihedral group of the symmetries of the square.It follows that the largest subgroup of E (2) that leaves L invariant is the compact semidirect sum E (L ) := D 4 • +T 2 .Therefore, we are interested in the action on F of the compact group G = E (L ) × O(2).We have the following Proposition.Proposition 14.The operator F defined in (13) is equivariant w.r.t. the action of the group G .
Proof.The action of the element g = ( , ψ) ∈ G on F , with ψ given by (22), is defined as follows and, because of Proposition 10, the reader will verify that F is equivariant w.r.t. the action of G , i.e. that we have At this point, branches of planforms are usually obtained by applying the equivariant branching lemma [14,17,30,60] as follows.The kernel V of d F (0, γ p ) is G -invariant.We fix an isotropy subgroup Σ of G (i.e. for which there exists v 0 ∈ V such that g • v 0 = v 0 for all g ∈ Σ) and compute the dimension of the fixed point subspace Fix V (Σ) of V associated with Σ where Fix V (Σ) The equivariant branching lemma states that if Σ is an axial subgroup of G, i.e. such that then generically, there is a branch of steady-state solutions to (1) with symmetries Σ.The genericity condition requires that the eigenvalue that goes through 0 does so with nonzero speed and that some components of the Taylor expansion of F are non zero.We characterize the isotropy subgroups of G and the corresponding fixed point subspaces in the following Proposition.Proof.If we restrict the color group O(2) to SO(2), i.e. if we only consider the action of the color rotations, our symmetry group is the same as for the equivariant Hopf bifurcation with E (L ) symmetry [23,55].We note that the fixed point subspaces have even dimension since SO (2) commutes with E (L ) and all real finite dimensional non trivial representation of SO( 2) are of even dimension.All fixed point subspaces of V for the action of H are therefore of even dimension.
If we now extend SO (2) to O(2), i.e.H to G , given an isotropy subroup Σ of G , either it is an isotropy subgroup of H and we are done, or it is not and the color reflection operator τ c is not reduced to the identity on Fix V (Σ).But τ c is a projector.It has therefore two eigenspaces V ±1 corresponding to its eigenvalues ±1 and V is the direct sum of the corresponding eigenspaces.
so that the projection fall in three classes [36] I: Closed subgroups of SO (2).II: Closed subgroups containing −1: they are of the form Σ ⊕ Z 2 , Σ a subgroup of SO (2).III: Closed subgroups of O(2) which are not a subgroup of SO(2) and do not contain −I .
In Proposition 15 we considered only the subgroups of O(2) of class I and II.As shown in, e.g., [30][Chapter 13, page 131], those subgroups Σ are determined by pairs K ⊂ H of subgroups of SO (2) such that K has index 2 in H , see also [46].If p : SO(2) ⊕ Z 2 → SO(2), p(ψ) = d et (ψ)ψ, is the projection, then H = p(Σ) is isomorphic to Σ and Σ = K ∪ ψK for any −ψ ∈ H \K .It is not difficult to list all such pairs (K , H ) using the results in [45].None of them produces a subgroup of E (L ) × O(2) with a nonzero fixed subspace, basically because the color rotations act on V by multiplications with a magnitude 1 complex exponential.
Since V is G -invariant and G is compact, we can write V as a direct sum of G -irreducible subspaces (subspaces such that their only G -invariant subspaces are 0 or the whole subspace), see e.g.[31][Theorem 1.22]: for some j so that the first step in classifying the planforms associated with a fixed lattice L is to enumerate each G -irreducible subspace of V .Moreover V j is of the form V s j ×V c j where V s j is an eigenspace of W s irreducible under the action of E (L ) and V c j is an eigenspace of W c irreducible under the action of O(2).This has been worked out in the case of E (L ) by several authors including [21].Dropping the upper index j for simplicity, the irreducible representations of V s must be of the form where are no (nontrivial) translations in E (L ) that act trivially on (28).This requirement ensures that we have found the finest lattice, L , that supports the neutral modes (28) [21].
The first one corresponds to the second to where u and v are relatively prime strictly positive integers such that u + v is odd.
With our notations, we have z j e 2πi 〈K j , r 〉 e 2πi nϕ e h,m ( ) + w j e −2πi 〈K j , r 〉 e 2πi nϕ e h,m ( ) which is isomorphic to If we now consider an isotropy subgroup Σ of H , such that dim Fix V (Σ) = 2, Proposition 7.2 in Chapter XVI in [30] asserts that Σ is a twisted group, i.e. there is a pair of subgroups K ⊂ G of E (L ) and a unique homomorphism Θ : G → SO (2) K is the kernel of Θ: K = ker Θ.In [23][Table 16 on page 157 and Table 21 on page 160] the authors compute all such subgroups and their fixed point subspaces.Using Proposition 15 we can obtain the corresponding subgroups of G and their corresponding fixed point subspaces.
In detail D 4 is generated by the π/2 rotation, noted R s π/2 and the reflection, noted τ r 1 through the r 1 -axis.The action of E (L ) on V induces an action on C 2k .The reader will verify that this action is given by and Similarly, for the action of O( 2), we have and Writing shows that the projection of (z, w) on V 1 is 1 2 (z + w, z + w).The results in [23] allow us to determine all the axial subgroups of G in dimension 4 (because of [55]) and only some of them in dimension 8.They are shown in Table 2 where we have noted v d the vector 1  2 (k 1 + k 2 ), e s the identity of D 4 and e c the identity of O(2).We can now apply the equivariant branching lemma .In detail we have the following Proposition.) ) ) ) satisfies the assumptions of the equivariant branching lemma.For each axial subgroup Σ ⊂ G , there exists a unique branch of solutions to F (V, γ n ) = 0 emanating from (0, γ n ) where the symmetry of the solution is Σ.Moreover each of these branches stems from a pitchfork bifurcation.
Proof.The proof is a direct application of [14][Theorem 2.3.2]around γ ≈ γ n .From Section 3 the mapping F is smooth.Its Jacobian at V = 0 is the linear operator where Sig is defined in (2).As W m,2 (Ω) is compactly included in L 2 (Ω) and W ∈ L (L 2 (Ω)) is compact, it follows that W is compact from W m,2 (Ω) to L 2 (Ω).This implies that d F (0, γ n ) is a Fredholm operator with index 0.
From the previous linear analysis, 0 is an isolated eigenvalue of d F (0, γ n ) ∈ L (L 2 (Ω)) with finite multiplicity.
Let us denote the reduced equation by f (v, γ) = 0, v ∈ V : it is given by the Lyapunov-Schmidt method, see e.g.[29].It remains to show that d 2 vγ f (0, 0) = 0 to complete the proof of the generic existence/uniqueness of a bifurcating branch.
We restrict ourselves to the case ε = 0 for which s (0) = 0.The methods described in e.g.[29][Chapter VII, page 295], allow us to find, using the fact that d F is hermitian, that The type of bifurcation that occurs can be determined by applying, e.g., [14][Theorem 2.3.2].For each of the two subgroups G θ of G which appear in Table 2 we verify (using the relation ) is in the normalizer N (G θ ) of G θ since it satisfies g 0 g g −1 0 = g for all generators of G θ glven in Table 2 and acts as −1 on p 1 Fix(G θ ) .
This Theorem allows us to discover a specific class of planforms, i.e. some of those satisfying (24).In general and generically there may exist solutions such that dim(Fix V (Σ)) > 1.Our approach here will not allow us to find these.Note however that the assumption (24) is the most commonly made.Exceptions to this can be found in [11,15,16].

Examples of planforms
Using Table 2 and ( 28), it is easy to write down the analytical expressions of the four types of planforms.In the case where dimV = 4, we have SR(r, , ϕ) = cos 2πr 1 e h,m ( ) cos(2πnϕ), (33) and The first type of planform is called standing rolls (SR) or stripes, the second type is called spots.We show in Figure 2 an example of standing rolls (or stripes) described by equation ( 33) in cortical and retinal coordinates, respectively.The relation between cortical and retinal coordinates is discussed in Remark 21.Similarly we show in Figure 3 an example of spots described by equation (34).2πnϕ), (36) where u and v are relatively prime strictly positive integers such that u + v is odd. Figure 4 shows an example of a planform S2 1,2 and Figure 5 shows an example of a planform S4 3,2 .Remark 18.To determine the color represented at cortical location (r 1 , r 2 ), and produce the images in Figures 2-5, we compute c max (r 1 , r 2 ) = argmax ,ϕ Υ V eq (r 1 , r 2 , , ϕ) ∈ Ω c (corresponding to the maximum activity of the neural mass with spatial coordinates (r 1 , r 2 )), V eq ∈ {SR, S2, S2 u,v , S4 u,v }, and define the Luminance to be the corresponding activity level Υ V eq (r 1 , r 2 , c max (r 1 , r 2 )) .In case there are several argmax, we arbitrarily select one of them.At each cortical location (r 1 , r 2 ) we display the color using the HSL color coordinates, see Appendix A and [37]: Remark 19.In Figures 2-5 the value of n in Proposition 12 is equal to 6.The value of ξ c in Proposition 13 is equal to 2 and x n is the fourth strictly positive root of equation ( 16).This implies that the product e h,m ( ) cos(2πnϕ) has several minima and maxima so that there are several values of ( , ϕ) where the maximum of Υ(SR(r, , ϕ))(respectively of Υ(S2(r, , ϕ)), Υ( S2 u,v (r, , ϕ)) or Υ(S4 u,v (r, , ϕ))) with respect to ( , ϕ) is reached.As in Section 7.1 we have chosen arbitrarily one of them.

Numerical bifurcation analysis
We recall that we are interested in visual hallucinations which are stable solutions of (1) and exist in fairly large regions of the parameter space.The equivariant branching lemma (EBL) provided a set of stationary solutions of (1) for a constant external current I ext = 0.It relies on the Lyapunov-Schmidt reduction (or more generally on the existence of a center manifold) which is local by essence, i.e. valid in a given neighborhood of (0, γ n ) in F × R. The size of this neighborhood quantifies the predictability of the local theory and is bounded by min( Hence, for large cortices, it is vanishingly small.
We thus use numerical bifurcation analysis to assess our theoretical predictions beyond their above domain of validity.We present some numerical results concerning the equilibria of (1) for a constant external current I ext = 0 as functions of different parameters.Because of the volume of data required to explore the full four-dimensional model (2 dimensions in space and 2 in color) we restrict ourselves to a one-dimensional color space by considering a diameter of the unit disc determined by its hue angle ϕ 0 (between 0 and 1 to be consistent with our notations).A point on this diameter is characterized by its polar coordinates (2πϕ 0 , ρ) or (2π(ϕ 0 Note that the problem associated has now the symmetry group (D 4 2).This has been extensively studied [22,56].A key remark is that the additional reflection symmetry has no effect on the bifurcation problems associated with the square lattice 6 .
In the first bifurcation diagram (Section 7.4), we observe that the type of the bifurcated branches and their stability predicted by the EBL are valid outside the neighborhood for the first bifurcation point, of dimension 4. The second bifurcation point, of dimension 8, however provides an example of a branch which becomes stable.Additionally, all stable patterns are alike (stripes) and in agreement with the EBL predictions.
In the second bifurcation diagram (section 7.5), we switched the bifurcation points, forcing the first one to be of dimension 8, in hope to see new stationary states, not like stripes.As before, only the stripes are stable and the EBL predictive power, outside the neighborhood, is not bad.
In search for more interesting hallucinations and based on [65], we changed the criticality of the bifurcation points in the third bifurcation diagram (section 7.6) hoping to stabilize the patterns with the creation of fold bifurcations.The EBL predictions concerning the stability of the branches are challenged very quickly.In passing, this provides an example of stable spots.Additionally, the bistability between the state V = 0 and the branch of stripes suggests snaking branches.These snaking branches are very interesting because they give birth to stable spatially localized visual hallucinations.The last bifurcation diagram shown in Section 7.7 is dedicated to finding snaking branches.
To conclude, it seems that we can rely on the EBL to predict the existence of patterns that survive in large parameter domains but not to predict their stability.From the simulations, it seems that the spatial structure of the solutions does not vary much along the branches and this suggests the possibility to extend the results from [51] to the present setting.If such result were true, it would imply that the EBL is a very valuable tool at elucidating the spatial structure of the visual hallucinations.

Color representation of the equilibria of (1)
The equilibria of (1) are represented by a color image in a way similar to what is described in Section 6.In Figure 6-Left, we show an example of such an equilibrium V eq (r 1 , r 2 , c) where we plot some of its level sets in the three-dimensional space of coordinates (r 1 , r 2 , c).To determine the color represented at cortical location (r 1 , r 2 ), and produce the image in Figure 6-Middle, we compute c(r 1 , r 2 ) = argmax c Υ V eq (r 1 , r 2 , c) ∈ [−1, 1] (corresponding to the maximum activity of the neural mass with spatial coordinates (r 1 , r 2 )) and define the Luminance to be the corresponding activity level a(r 1 , r 2 ) := Υ V eq (r 1 , r 2 , c(r 1 , r 2 )) .In case there are several argmax, we arbitrarily select one of them as in Section 6.At each cortical location (r 1 , r 2 ) we display the color using the HSL color space, see Appendix A and [37]: and obtain the result shown in Figure 6-Middle.ϕ 0 = 1 8 is defined in Appendix A. Remark 20.Note that this way of displaying our results is a severe simplification of our model which says that at every location (r 1 , r 2 ) in the cortex, color is represented by the function c → Υ(V eq (r 1 , r 2 , c)) and not by the three numbers (37), see [40] for an interesting discussion of these issues and much more.
Remark 21.The visual cortex of several species (like monkeys) has the property of being retinotopically organized, e.g.[48,53].That is, there is a one-to-one mapping from retinal coordinates to cortical ones, the mapping in humans being approximately log polar.Thus, we need to apply an inverse log polar transform to the equilibrium shown in Figure 6-Middle in V1 coordinates to obtain the result shown in Figure 6-Right in retinal coordinates.

Connectivity functions
In all subsequent examples, we use the following spatial connectivity function w s (r ) of the form (4) whose parameters µ s , ν s and β s are chosen so that it satisfies (5):

Numerical experiments
All numerical computations were performed in the Julia programming language (version 1.4.2).The bifurcation diagrams were computed using a pseudo-arclength continuation method implemented in the package BifurcationKit.jl[63] with version v0.1.2 .The continuation is based on a Newton-Krylov method to solve (12) with GMRES linear solver [5].The computation of the eigen-elements, to detect the bifurcation points, is based on the Arnoldi algorithm [52].The linear and eigen solvers are both implemented in the package KrylovKit.jl.The nonlinear equations (12) were solved at tolerance 10 −11 in the supremum norm.The bifurcation points were located using a bisection algorithm (on the number of unstable eigenvalues) leading to a precision of 10 −4 on the value of the parameter γ, see (2).Let us describe how ( 13) is implemented.The connectivity kernels were computed with the three-dimensional Fast Fourier Transform (3D FFT) on Graphics Processing Units (GPU) based on the package CUDA.jl(see [7,8]).
In order to compute the bifurcating branches, the reduced equations (see [30] or proof of proposition 17) at the bifurcation point were computed thereby yielding a system of polynomial equations of degree 3 in a number of variables equal to the dimension of the kernel of the Jacobian.The roots of these polynomial equations are then computed and used as guesses for points on the bifurcated branches which are corrected using a Deflated-Krylov-Newton (see [26]) to prevent converging to the trivial solution.This provides an entirely automatic procedure to find the bifurcated branches at a bifurcation point of any dimension.We call this procedure automatic branch switching (aBS).
The whole program runs entirely on GPU, a V100 Nvidia card with 32Gb of RAM.The computations are next to impossible to run without a GPU albeit perhaps on a cluster.In the experiments, we use the values N r 1 = N r 2 = 256, N c = 64.Using a finer discretisation would have helped but we were limited by memory when computing eigen elements.Indeed, computing the branches without stability is a matter of a few minutes.However, we found necessary to use a Krylov Space of size ≈ 100, due to the symmetries, in order to compute the eigenvalues, and this limited the size N r 1 × N r 2 × N c of the discretisation.

First bifurcation diagram
Figure 8 represents the variation of the equilibria (represented by their L 1 norm) as functions of the slope γ of the activation function s, see (2).The trivial equilibrium V = 0 loses stability at a first bifurcation point around γ ≈ 9.78.The dimension of the first primary bifurcation point is four whereas the dimension of the second one is eight.Given the symmetries of the network, it is straightforward to conclude that the first primary bifurcation is a Pitchfork with D 4 symmetry group.
Using the aBS procedure (see Section 7.3), we computed the different bifurcating branches from the first two primary bifurcation points.We know that two branches (at most) bifurcate from the first point: one branch of stable stripes (the blue thick line) and one branch of unstable spots (the blue thin line).At the 8D primary bifurcation point, we found four branches with patterns in agreement with the ones described in [22].
We also computed the bifurcated branches (shown in yellow and brown) from the secondary bifurcation point on the spot branch (in light blue) hoping to find new stable patterns.
This diagram hints at the fact that only the stripes are stable.At the bottom of Figure 8 we display the images corresponding to some of the equilibria shown in the plot at the top of the Figure .The images are embedded in a frame the same color as the curve on which the corresponding equilibrium sits.A black star has been added to the two images corresponding to the two equilibria (also marked with a black star) sitting on the line of equilibria branching out of the line of unstable spots (thin light blue).

Second bifurcation diagram
Next, we modify a bit the connectivity to make the dimension of the first bifurcation point 8 and that of the second one 4.We did this hoping to find more interesting stable patterns than in C. R. Mathématique-Draft, 20th October 2021  3), we compute the different bifurcating branches from the bifurcation points.Again, only the stripes are found to be stable.The bifurcated patterns from the 8D bifurcation point are in agreement with the ones referenced in [22].In fact, the first two bifurcation diagrams collectively show all the possible patterns that can bifurcate from an 8D bifurcation point.

Third bifurcation diagram, changing the criticality
We then change the criticality of the bifurcation points hoping to induce new stable equilibria.To this end, we use the following observation.The normal form of the 4D bifurcation point is [33] 2, and (for example for β 1 , see [62] ) (2) (ε), s 3 = Si g (3) (ε), γ 0 is the value of γ at the bifurcation point and λ 0 is the associated eigenvalue of W .We recall that for a threshold ε = 0, s 2 = 0, s 3 < 0 and thus the Pitchfork is supercritical.When varying the threshold ε, the ratio s  Using this procedure, we obtain the bifurcation diagram shown in Figure 10.We observe that all primary branches bifurcate sub-critically and thus present a saddle-node bifurcation as proved in [65].The first two bifurcation points are 4D and 8D, respectively.The stripes branches (red or light brown) from either bifurcation points become stable after the saddlenode bifurcation.Remarkably, the spots branch (blue) from the first 4D bifurcation point is stable close to its saddle-node bifurcation.Using aBS, we find a connection (in violet) between the spot branch and the stripe branch.
A very similar situation was found in [59] for the Selkov-Schnakenberg model.The fact that there is bistability between the trivial solution and the stripes (or spots) in a gradient system such as (1) (see [62] for a proof) opens the doors for the existence of localized patterns and in particular snaking (see [39] for a review) of localized solutions.
Snaking is interesting in the context of visual hallucinations because it produces stable patterns.These snaking branches usually originate from the bifurcation points on the stripes (resp.spots) branch.We thus computed the secondary branches on the spots branch and show one of them (in green).The green branch connects two bifurcation points on the stripes branch.It consists of localized stripes as seen in Figure 10.Unfortunately, it is neither stable nor does it feature snaking.We come back to this point in the next section.Finally, we also computed a secondary branch (brown) from the branch of spots (primary branch of the 8D bifurcation point).It shows another type of (unstable) localized patterns of spots.

Fourth bifurcation diagram, snaking branch of localized solutions
In this last example, we focus on finding snaking branches for (1).Theoretical analysis (one of the earliest is [12]) shows that snaking occurs in a region which is exponentially small in = γ 1 − γ SN 1 where γ 1 is the value of γ at the 4D bifurcation point and γ SN 1 is that at the location of the saddle-node bifurcation point on the branch of stripes.In [59] the authors study this behaviour numerically.Thus, to observe snaking, we have to increase ε.By performing codim 2 continuation of the saddle-node bifurcation and the first primary bifurcation point in the (γ, ε) plane, we were able to select a better value of the threshold ε than in Figure 10.
We show the results in Figure 11 where we have a snaking branch (green) of solutions arising from a 4D bifurcation point on the branch of stripes (red).We observe that the localization of the pattern is parallel to the stripes.Interestingly, there is another branch (data not shown) of patterns where the localization occurs perpendicular to the stripes and the branch does not snake as explained in [4].

Conclusion
We have studied some aspects of the bifurcations of the solutions to the neural-field equations describing a model of the primary visual area V1.We had proposed a variant of this model in [57] and validated its predictions with some psychovisual data.Here we have focused on the stationary solutions to our equations and their bifurcations which we loosely interpreted as being possible metaphors of visual hallucinations.Visual hallucinations in color are not very well documented in the literature and it is therefore very difficult to compare our model predictions with experimental observations.The closest work to ours is described in [10] and several of the follow-up papers.In these authors' work there was no attempt to model color perception, the focus was on achromatic (i.e.black and white) vision and the interplay between space and edges (borders of objects in the scene) orientations.Like us they studied the bifurcations of the solutions to the neural field equations describing their model and proposed to interpret them as accounting for visual hallucinations.Just like ours, theirs was a rather loose metaphor.
Going into the mathematics behind the model, we have established here and in [57,65] the well-posedness of this class of equations whose solutions are defined on a four-dimensional compact spatio-chromatic space representing part of the cortical organisation of the visual cortex.The group of symmetries that acts on this space arises naturally from its spatial extension (assumed to be a square) and the psychophysics of Hering's color perception theory that emphasizes the symmetry with respect to the origin in chromaticity space (corresponding to what is known as opponent colors).Our model equations are equivariant w.r.t to this spatio-chromatic group and, by identifying some of its axial subgroups we have been able to prove (in part numerically) using the equivariant branching lemma that the simplest stationary solution (which turns out to be 0) bifurcates at pitchforks from which arise branches of stationary solutions enjoying the symmetries of the corresponding axial subgroup.
Going into the numerics, the Julia package developed by the third author has allowed us to explore a much larger range of values of the bifurcation parameter than for example the authors of [10] whose approach, based on the use of normal forms, is not well-suited to the analysis of subcritical bifurcations.It has also allowed us to discover cases of snaking and hence to observe stable localised solutions.To our knowledge this is the first time this has been reported in a model of the primary visual area.It raises the fascinating question of whether such spatio-chromatic patterns can be observed in animal or human perception.
Going further into the numerics, this is to our knowledge the first time that a bifurcation analysis of an infinite dimensional system of integro-differential equations is entirely carried C. R. Mathématique-Draft, 20th October 2021 out on GPUs.By considerably reducing the response time of our computer simulations it has allowed us to explore in much greater depth than what had been possible before the complicated structure of the solutions to our model equations with the result that we will in the future be able to better confront our predictions with perceptual experiments.

Proposition 15 .
The isotropy subgroups Σ of G = E (L ) × O(2) fall in two classes: those which do not contain the color reflection τ c are those of H := E (L ) × SO(2), and those that do contain τ c .In the first case the corresponding fixed subspaces are those of Σ as a subgroup of H .In the second case the fixed subspaces are those of Σ ∩ H projected on V 1 along V −1 , the eigenspaces of τ c corresponding to the eigenvalues 1 and -1, respectively.C. R. Mathématique-Draft, 20th October 2021

Figure 2 .
Figure 2.An example of a planform whose equation is given by (33): Left in cortical coordinates, Right in retinal coordinates.

Figure 3 .
Figure 3.An example of a planform whose equation is given by (34): Left in cortical coordinates, Right in retinal coordinates.

Figure 4 .
Figure 4.An example of a planform whose equation is (35): Left in cortical coordinates, Right in retinal coordinates.

Figure 5 .
Figure 5.An example of a planform whose equation is (36): Left in cortical coordinates, Right in retinal coordinates.

)Figure 6 .
Figure 6.Left: example of level sets of an equilibrium, a stationary solution, of (1).Middle: Color image representation of the equilibrium shown on the left.Right: Same as Middle but in retinal coordinates.

Figure 7 .
Figure 7. Left: Typical shape of the spatial connectivity w s used in the numerical experiments.Right: Heatmap plot of the color connectivity function w c (c, c ) used in the numerical experiments.

Figure 8 .
Figure 8. Parameters ε = 0, Ω s = (−20π, 20π) 2 , Ω c = (−1, 1), β s = 1.9, B = 3π, A = 1.0.Equilibria as functions of the nonlinear gain γ.The stable parts of the branches are indicated with thick lines.The dots indicate bifurcation points.The discretization is 256 × 256 × 64.Only the first two primary bifurcation points are shown.The colours around the images are those of the corresponding branches.The first seven images represent the intersections of the vertical black line with the bifurcated branches.The last two images correspond to the two black stars in the bifurcation diagram.

Figure 9 .
Figure 9. ε = 0.5, Ω s = (−20π, 20π) 2 , Ω c = (−1, 1), β s = 20.9,B = 6π, A = 1.0.Equilibria are shown as functions of the nonlinear gain γ.The stable parts of the branches are indicated with thick lines.The dots indicate bifurcation points.Only the first three primary bifurcation points are shown.As in Figure 8 the colours surrounding the images correspond to those of the branches found at the intersections with the vertical black line.

Figure 10 .
Figure 10.ε = 2.3, Ω s = (−20π, 20π) 2 , Ω c = (−1, 1), β s = 1.9, B = 3π, A = 9.Equilibria as function of the nonlinear gain γ.The stable parts of the branches are indicated with thick lines.The dots indicate bifurcation points.Only the first three primary bifurcation points are shown.The colour surrounding the image is that of its corresponding branch.Their location on the bifurcation diagram is shown with a black star.

s 2 2 takes
on any value, thus we can change the criticality by altering the threshold.Compared to the previous diagram, we use ε = 2.5 and scale the connectivity to restrict the value of γ at the first bifurcation point.

Figure 11 .
Figure 11.ε = 2.9, Ω s = (−20π, 20π) 2 , Ω c = (−1, 1), β s = 1.9, B = 3π, A = 9.Equilibria as functions of the nonlinear gain γ.The stable parts of the branches are indicated with thick lines.The dots indicate bifurcation points.Only the first three primary bifurcation points are shown.The right part of the Figure shows a zoomed version of the diagram in which the snaking is more apparent.Three images are shown at the bottom of the Figure, corresponding to the black stars in the two diagrams above.The colour surrounding the images corresponds to the one of its branch.

Figure 12 .
Figure 12.Stable equilibria in retinal coordinates.The top two images correspond to the first two in Figure 8.The bottom left image corresponds to the spot pattern (second one) in Figure 10.The bottom right image corresponds to the rightmost image in Figure 11.Note that the framing colors match.

Table 2 .
The axial subgroups of G in dimension 4 and two of them in dimension 8, as well as their corresponding fixed spaces.