O<EMAIL_ADDRESS>A<EMAIL_ADDRESS>1 R.
<EMAIL_ADDRESS>
# Spatial and color hallucinations in a mathematical model of primary visual
cortex

Olivier D. Faugeras 0000-0001-8363-1550  UniversitÃ© CÃ´te dâAzur, Inria,
2004, route des Lucioles - BP 93,
06902 Sophia Antipolis Cedex, France  [ ,  Anna Song  Department of
Mathematics, Imperial College London, London, UK  Haematopoietic Stem Cell
Laboratory, The Francis Crick Institute, London, UK  [ and  Romain Veltz
0000-0003-4653-1475  [

###### Abstract.

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 $\mathcal{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 $\mathcal{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".

###### keywords:

Integro-partial differential equations, Initial Value Problem, Equivariant
Bifurcations, Snaking, Equivariant Branching Lemma, Neural Fields, Color
Perception, Visual Hallucinations, Numerical Bifurcation Analysis, Julia
programming

###### 1991 Mathematics Subject Classification:

37L10,35B06,35B32,35G25,37M20,45B05,45K05,46E35,47G20,47H30,
65J15,65R0,9208,9210,92B20

##  1\. Introduction

Neural Fields are a useful mathematical formalism for representing the
dynamics of cortical areas at a macroscopic level, see the reviews [24, 9,
19]. 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, 65, 49] 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 [25, 10] .

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.

##  2\. The model

In this work we think of V1 as the closure of a regular domain, noted
$\Omega_{s}$, $s$ stands for space, of $\mathbb{R}^{2}$, in effect the open
square $(-l/2,l/2)\times(-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 7.2). 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 [70, 13] 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 color111We do not consider the achromatic component of a
color, except for display, see Section 7.1, and restrict ourselves to the
chromatic components.. 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 $\mathbb{R}^{2}$ 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 2.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 [69, 40, 50, 6, 61].

For technical reasons, we consider only the open disk minus a radius, i.e.
$D(0,1)\backslash\\{[0,1)\times\\{0\\}\\}$. We use polar coordinates to
parametrize $D(0,1)$. A chromaticity is therefore represented by a pair
$(\rho,\varphi)$ with $\rho\in(0,1)$ and $\varphi\in(0,1)$, keeping in mind
that the actual polar coordinates are $(\rho,2\pi\varphi)$.

######  Remark 2.2.

The elimination of the semi open radius $[0,1)\times\\{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
$\bar{D}(0,1)$.

The perceptual interpretation of $\rho$ is the saturation of the color $c$
while $\varphi$ is its hue. The opponent chromaticity of $(\rho,\varphi)$ is
therefore $(\rho,\varphi+\frac{1}{2}\mod 1)$. The set $\Omega_{c}$ of
chromaticities is therefore the open square $(0,1)\times(0,1)$.

We define

$\Omega=\Omega_{s}\times\Omega_{c}$

to be the bounded domain, in effect the open rectangle of $\mathbb{R}^{4}$,
encapsulating the spatial and chromatic coordinates that will be of interest
in the sequel.

###  2.1. Connectivity kernel

Putting all this together, at each point $(r,c)=(r_{1},r_{2},\rho,\varphi)$ of
$\Omega$, we consider a neural mass222A neural mass is an aggregate of neurons
in which the spatial dimension is ignored, see e.g. [67]. whose average
membrane potential is noted $V(r,c,t)$. It is a function defined on
$\Omega\times\mathcal{J}$, $\mathcal{J}$ being an interval of $\mathbb{R}$
containing 0. In [57], we assumed that the function $V$ was the solution to an
initial value problem333More precisely we used instead of $V$ the ”activity”
variable $a=\Upsilon(V)$, see the definition (2) of the sigmoid function
below. The two formulations are equivalent in the case $I_{\rm ext}=0$
considered here. related to a Hammerstein equation444In mathematical
neuroscience this integro-differential equation is often called a Wilson-Cowan
equation [68] which writes

$\tau\frac{\partial V}{\partial
t}(r,c,t)=-V(r,c,t)+w\star\Upsilon\left(V(r,c,t)\right)+I_{\rm ext}(r,c,t),$
(1)

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 $\Omega$ starting from the initial condition $V_{0}$. At each time
$t$, $V(t)$ belongs to some functional space, in effect a Hilbert space
$\mathcal{F}$, that we describe in the next section.

We now discuss the various elements that appear in this equation.

$\tau$ is a time constant that defines the speed of the exponential decay
toward the initial condition. Without loss of generality we can assume
$\tau=1$.

$\Upsilon$ is a sigmoidal function mapping $\mathbb{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

$\Upsilon(x)={\rm Sig}(\gamma x-\varepsilon)\quad\text{ where }\quad{\rm
Sig}(x)=\frac{1}{1+e^{-x}}.$ (2)

$\varepsilon$ is a parameter that allows us to shift the origin, $\gamma$
controls the slope of the sigmoid at the (shifted) origin, it is often called
the nonlinear gain.

$I_{\rm ext}$ is a function representing the input to the neural mass from
different brain areas. In the remaining of this paper we take $I_{\rm ext}=0$,
i.e. we consider that area $V1$ 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 kernel555It is only qualitatively the same as the one
in [57].. It models the influence of neighboring neural masses at
$(r^{\prime},c^{\prime})$ on the neural mass at $(r,c)$ as a (separable in
space and color) linear superposition operation

$w\star
s\left(V(r,\rho,\varphi,t)\right)=\int_{\Omega_{s}\times\Omega_{c}}w(r,\rho,\varphi,r^{\prime},\rho^{\prime},\varphi^{\prime})\Upsilon\left(V(r^{\prime},\rho^{\prime},\varphi^{\prime},t)\right)\,\rho^{\prime}d\rho^{\prime}\,d\varphi^{\prime}\,dr^{\prime},$
(3)

with

$w(r,c,r^{\prime},c^{\prime})=w_{s}(r-r^{\prime})w_{c}(\rho,\rho^{\prime},\varphi-\varphi^{\prime}),$

where the index $s$ stands for “space” and the index $c$ stands for
“chromaticity”.

If $w(r,c,r^{\prime},c^{\prime})$ is positive (respectively negative) the
neural mass at $(r^{\prime},c^{\prime})$ 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:

$w_{s}(r)=\mu_{s}e^{-\|r\|_{2}^{2}/2\alpha_{s}^{2}}-\nu_{s}e^{-\|r\|_{2}^{2}/2\beta_{s}^{2}},$
(4)

where $\|\ \|_{2}$ is the usual Euclidean $\mathbb{L}_{2}$ norm in
$\mathbb{R}^{2}$. Biology dictates that $\alpha_{s}$, $\beta_{s}$ are very
small w.r.t $1$, the size of $\Omega_{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_{\rm ext}$
in (1) which we have taken to be equal to 0.

If $w_{s}(0)>0$, i.e. if $\mu_{s}>\nu_{s}$, the neural masses at $r^{\prime}$
such that $\|r-r^{\prime}\|_{2}$ is small enough excite the neural mass at
$r$, and if $\alpha_{s}<\beta_{s}$ those sufficiently far away inhibit it.

We furthermore assume

$\int_{\Omega_{s}}w_{s}(r)\,dr=0,$ (5)

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 function $w_{c}$ is separable in polar coordinates,

$w_{c}(\rho,\rho^{\prime},\varphi)=w_{c,m}(\rho,\rho^{\prime})w_{c,a}(\varphi),$
(6)

with

$w_{c,m}(\rho,\rho^{\prime})=e^{-\xi_{c}|\rho^{2}-\rho^{\prime 2}|},$ (7)

and

$w_{c,a}(\varphi)=\mu_{c}e^{-2\pi\alpha_{c}|(\varphi-\frac{1}{2}\mod
1)-\frac{1}{2}|}-\nu_{c}e^{-2\pi\beta_{c}|(\varphi\mod 1)-\frac{1}{2}|}$ (8)

The parameters $\xi_{c}$, $\mu_{c}$, $\nu_{c}$, $\alpha_{c}$ and $\beta_{c}$
are positive and such that $w_{c,a}(c,c^{\prime})$ is positive if $c^{\prime}$
is close to $c$ and negative if $c^{\prime}$ is close to $-c$, i.e. in terms
of polar coordinates if $\rho^{\prime}$ is close to $\rho$ and either
$\varphi^{\prime}$ is close to $\varphi$ or $\varphi^{\prime}$ is close to
$\varphi+\frac{1}{2}\mod 1$.

Figure 1. Plot of $w_{c,a}(\varphi)$ defined in (8) for the values of the
parameters $\alpha_{c}=0.3$, $\beta_{c}=0.4$, $\mu_{c}=0.6$ and
$\nu_{c}=0.69$.

######  Remark 2.3.

Note that the function $w_{c,a}:\mathbb{R}\to\mathbb{R}$ is 1-periodic and
even.

######  Remark 2.4.

For technical convenience we perform the change of variable
$\rho^{2}\to\varrho$ so that $w_{c,m}$ simplifies to

$w_{c,m}(\varrho,\varrho^{\prime})=e^{-\xi_{c}|\varrho-\varrho^{\prime}|},$
(9)

and the measure $\rho\,d\rho$ on $(0,1)$ is equal to $\frac{1}{2}d\varrho$. A
chromaticity is therefore represented by a pair $(\varrho,\varphi)$,
$\varrho,\,\varphi\in(0,1)$.

The product $w$ of $w_{s}$ and $w_{c}$ behaves as a function of
$r,\,r^{\prime},\,c,\,c^{\prime}$ as shown in Table 1 adapted from [57].

Table 1. Sign of the connectivity kernel $w$ $w_{s}(r-r^{\prime})w_{c}(\varrho,\varrho^{\prime},\varphi-\varphi^{\prime})$ |  $c^{\prime}$ close to $c$ |  $c^{\prime}$ close to $-c$
---|---|---
$r^{\prime}$ close to $r$ | $>0$ | $<0$
$r^{\prime}$ far from $r$ | $<0$ | $>0$

As pointed out in this paper, this is in qualitative agreement with the
nonlinear behaviour of color shifts found in [44, 43].

Note that by definition of $w_{s}$ and $w_{c}$ the function $w$ is symmetric
w.r.t. $r$ and $r^{\prime}$, and $c$ and $c^{\prime}$, respectively.

###  2.2. Choice of the appropriate functional space $\mathcal{F}$ and well-
posedness of (1)

Our choice of $\mathcal{F}$ is guided by three criteria:

  1. (1)

The well-posedness of the problem,

  2. (2)

Its biological relevance,

  3. (3)

Its suitability for numerical computations.

The choice of a Hilbert space is appealing and a natural choice is
$\mathcal{F}=\mathbb{L}^{2}(\Omega)$. As argued in [65], this unfortunately
allows the membrane potential to be singular since for example the function
$(r,c)\to\|r-(\frac{1}{2},\frac{1}{2})\|_{2}^{-1/2}\
\in\mathbb{L}^{2}(\Omega)$. 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)\to V(r,c)$
is differentiable almost everywhere.

The choice of the Sobolev space $\mathcal{F}=W^{m,2}$, is convenient for two
reasons. First it is a Hilbert space endowed with the usual inner product:

$\langle V_{1},V_{2}\rangle_{\mathcal{F}}=\sum_{|\alpha|=0}^{m}\langle
D^{\alpha}V_{1},\,D^{\alpha}V_{2}\rangle_{\mathbb{L}^{2}(\Omega)}\quad\forall
V_{1},\,V_{2}\in\mathcal{F}$ (10)

where the multi-index $\alpha$ is a sequence
$\alpha=(\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4})$ of 4 integers and
$|\alpha|=\sum_{i=1}^{4}\alpha_{i}$, and the symbol $D^{\alpha}$ represents a
partial derivative, e.g.,

$D^{\alpha}V=\frac{\partial^{\alpha_{1}+\cdots+\alpha_{4}}}{\partial
r_{1}^{\alpha_{1}}\partial
r_{2}^{\alpha_{2}}\partial\varrho^{\alpha_{3}}\partial\varphi^{\alpha_{4}}}V$

The second reason is that, because the boundary of $\Omega$ is sufficiently
regular (it satisfies the cone-property [1, 2]), $\mathcal{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 5.6, for which we require some
smoothness of $w$.

######  Remark 2.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
$\Omega$, i.e. 4. Hence the smallest possible value of $m$ is 3. But in
Section 7 we assume that $\Omega_{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 $\mathbb{L}^{2}(\Omega)$ as
follows. Let $U\in L^{2}(\Omega)$, we define

$W\cdot U(r,c):=w\star
U(r,c)=\int_{\Omega}w_{s}(r-r^{\prime})w_{c}(c,c^{\prime})U(r^{\prime},c^{\prime})\,dr^{\prime}\,dc^{\prime}$

######  Remark 2.6.

Remember that the measure $dc^{\prime}$ is in effect equal to
$\frac{1}{2}d\varrho^{\prime}\,d\varphi^{\prime}$.

It is clear that this is well defined and we have the following proposition.

######  Proposition 2.7.

The operator $W$ maps $\mathbb{L}^{2}(\Omega)$ to $\mathcal{F}$ and hence
$\mathcal{F}$ to $\mathcal{F}$.

######  Proof 2.8.

The proof follows from Proposition 2.3 in [65].

We also have the following Proposition about the solutions to (1).

######  Proposition 2.9.

For each $V_{0}\in\mathcal{F}$ there exists a unique solution in
$\mathcal{C}(\mathbb{R}^{+},\mathcal{F})$ to the following Cauchy problem

$\left\\{\begin{array}[]{ccl}\frac{dV}{dt}&=&-V+W\cdot\Upsilon(V)\\\
V(0)&=&V_{0},\end{array}\right.$ (11)

where the function $s$ is defined by (2).

######  Proof 2.10.

The proof is a an immediate consequence of Proposition 2.5 in [65].

##  3\. Stationary solutions and bifurcations thereof

In this paper, we focus on the stationary (independent of time) solutions to
(11), the steady-states. 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 $\gamma$ and
$\varepsilon$ in (2), $\mu_{s}$, $\nu_{s}$, $\alpha_{s}$, $\beta_{s}$ in (4),
$\xi_{c}$ in (7), and $\mu_{c}$, $\nu_{c}$, $\alpha_{c}$, $\beta_{c}$ in (8)
vary. In effect, we will concentrate on the first parameter, the nonlinear
gain $\gamma$, 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 $\gamma$, 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 $\gamma$ varies of
the solutions to the equation

$F(V,\gamma)=0,$ (12)

where the operator $F$ is defined by

$F\stackrel{{\scriptstyle
def.}}{{=}}\left\\{\begin{array}[]{ccc}\mathbb{W}^{m,2}(\Omega)\times\mathbb{R}&\longrightarrow&\mathbb{L}^{2}(\Omega)\\\
(V,\,\gamma)&\longrightarrow&-V+W\cdot\Upsilon\left(V\right).\end{array}\right.$
(13)

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 $\Omega$. 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 $\mathbb{L}^{2}(\Omega)$ associated to the functions $w_{s}$ of (4)
and $w_{c}$ of (6) as follows

$W_{s}\cdot
U(r,c)=\int_{\Omega_{s}}w_{s}(r-r^{\prime})U(r^{\prime},c)\,dr^{\prime}\quad\text{and}\quad
W_{c}\cdot
U(r,\varrho,\varphi)=\frac{1}{2}\int_{\Omega_{c}}w_{c}(\varrho,\varphi,\varrho^{\prime},\varphi^{\prime})U(r,\varrho^{\prime},\varphi^{\prime})\,d\varrho^{\prime}\,d\varphi^{\prime}.$

The separability of the function $w$ in the space and color variables is
reflected in the notation $W=W_{c}\otimes W_{s}$. Fubini’s Theorem dictates
that $W$ is also equal to $W_{s}\otimes W_{c}$.

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
$W_{c}$ is itself separable, $W_{c}=W_{c,m}\otimes W_{c,a}=W_{c,a}\otimes
W_{c,m}$:

$W_{c,m}\cdot
U(r,\varrho,\varphi)=\frac{1}{2}\int_{0}^{1}w_{c,m}(\varrho,\varrho^{\prime})U(r,\varrho^{\prime},\varphi)\,d\varrho^{\prime}\quad\text{and}\quad
W_{c,a}\cdot
U(r,\varrho,\varphi)=\int_{0}^{1}w_{c,a}(\varphi-\varphi^{\prime})U(r,\varrho,\varphi^{\prime})\,d\varphi^{\prime}$
(14)

where the convolution in the right hand side of the definition of $W_{c,a}$ is
a periodic convolution.

##  4\. The spectrum of $W$ in $\mathbb{L}^{2}(\Omega)$

The reader can verify that, given the symmetry properties of the functions
$w_{c}$ and $w_{s}$, the operator $W$ is symmetric in
$\mathbb{L}^{2}(\Omega)$, i.e. satisfies

$\langle U_{1},\,W\cdot U_{2}\rangle_{\mathbb{L}^{2}(\Omega)}=\langle W\cdot
U_{1},\,U_{2}\rangle_{\mathbb{L}^{2}(\Omega)}.$

Another important property of $W$ is that it is _c_ ompact, i.e. given any
bounded sequence $(U_{n})_{n\geq 0}$ of functions in $\mathbb{L}^{2}(\Omega)$
the sequence $(W\cdot U_{n})_{n\geq 0}$ contains a converging subsequence for
the norm of $\mathcal{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}$.

###  4.1. The spectrum of $W_{s}$ in $\mathbb{L}^{2}(\Omega_{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 $\mathcal{L}$ generated by the two vectors
$k_{1}$ and $k_{2}$ where $(k_{1},k_{2})$ is the canonical basis of
$\mathbb{R}^{2}$ and commands that we quotient $\mathbb{R}^{2}$ by
$\mathcal{L}$ thus obtaining the 2-torus
$\tilde{\Omega}_{s}=\mathbb{R}^{2}/\mathcal{L}$. The dual lattice
$\mathcal{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:

######  Hypothesis 1.

The biologically relevant solutions to (1) are $\mathcal{L}$-periodic i.e.
satisfy

$V(r+\boldsymbol{l},c,t)=V(r,c,t)\quad\forall\boldsymbol{l}\in\mathcal{L}$

We now work on the Hilbert space $\mathbb{L}^{2}(\tilde{\Omega}_{s})$ of
$\mathcal{L}$-periodic functions with the same inner product as before. We
note $\tilde{W}_{s}$ the operator $W_{s}$ restricted to this space i.e.
defined by

$\tilde{W}_{s}\cdot
V(r)=\int_{\Omega_{s}}w_{s}\left(r-r^{\prime}\mod\mathcal{L}\right)V(r^{\prime})dr^{\prime}\quad\forall
r\in\Omega_{s},$

Note that the spatial convolution in (3) is now a _periodic_ convolution.
$\tilde{W}_{s}$ is clearly a symmetric compact operator on
$\mathbb{L}^{2}(\tilde{\Omega}_{s})$.

It is easy to characterize the spectrum of $\tilde{W}_{s}$:

######  Proposition 4.1.

The eigenvalues $\lambda^{s}_{m,n}$, $m,\,n\in\mathbb{Z}$, of
$\tilde{W}_{s}\in\mathbb{L}^{2}(\tilde{\Omega}_{s})$ are the (real) Fourier
coefficients of the $\mathcal{L}$-periodic even function $\tilde{w}_{s}$:

$\lambda^{s}_{m,n}=\int_{\Omega_{s}}w_{s}(r)\cos 2\pi\langle
k_{m,n},\,r\rangle\,dr,$

where

$k_{m,n}\stackrel{{\scriptstyle def.}}{{=}}mk_{1}+nk_{2}.$

Because $w_{s}$ is even, so is the sequence of eigenvalues:
$\lambda^{s}_{m,n}=\lambda^{s}_{-m,-n}$, $m,\,n\in\mathbb{Z}$. For
$(m,n)\neq(0,0)$ the eigenspaces of $\tilde{W}_{s}$ are of even dimension,
larger than or equal to 2. If $(m,n)=(0,0)$, because of (5),
$\lambda^{s}_{0,0}=0$ and the dimension of the kernel of $\tilde{W}_{s}$ is 1.
Given an eigenvalue $\lambda^{s}_{m,n}$ the corresponding eigenspace is
generated by the functions:

$\cos 2\pi\langle k_{p,q},\,r\rangle\quad\text{and}\quad\sin 2\pi\langle
k_{p,q},\,r\rangle,$

for all $(p,q)$ such that $\lambda^{s}_{p,q}=\lambda^{s}_{m,n}$.

######  Proof 4.2.

This follows for example from the fact that $\tilde{W}_{s}$ is diagonalized by
the Fourier Transform.

###  4.2. Spectrum of $W_{c}$ in $\mathbb{L}^{2}(\Omega_{c})$

Define the operator $\tau_{c}$ acting on the $\varphi$-periodic functions of
period $1$ as

$\tau_{c}\cdot u(\varrho,\varphi)=u(\varrho,-\varphi\mod 1)$ (15)

We have the following Proposition.

######  Proposition 4.3.

The linear operator $W_{c}$ on $\mathbb{L}^{2}(\Omega_{c})$ is symmetric and
compact, it commutes with the operator $\tau_{c}$ defined by (15).

######  Proof 4.4.

Obvious from the definitions.

From this we have

######  Corollary 4.5.

The spectrum of $W_{c}$ is real and at most countable. The eigenfunctions are
separable w.r.t. $\varrho$, $\varphi$.

######  Proof 4.6.

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

The eigenspaces of $W_{c,a}$ are generated by the functions $\sin 2\pi
n\varphi$ and $\cos 2\pi n\varphi$, $n\in\mathbb{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 4.8.

This follows immediatly from the fact that the operator $W_{c,a}$ is a
$1$-periodic convolution and the $1$-periodic function
$\mu_{c}e^{-2\pi\alpha_{c}|(\varphi-\frac{1}{2}\mod
1)-\frac{1}{2}|}-\nu_{c}e^{-2\pi\beta_{c}|\varphi\mod 1)-\frac{1}{2}|}$ is
even.

Regarding $W_{c,m}$, we have the following Proposition.

######  Proposition 4.9.

The eigenspaces of $W_{c,m}$ corresponding to the non zero eigenvalues are
one-dimensional. Its kernel is reduced to the null function. The eigenvalues
$\lambda_{n}^{c}$, $n\in\mathbb{N}$, are obtained from the countable solutions
$x_{n}$ to the transcendental equation

$\tan x=\frac{2\xi_{c}x}{x^{2}-\xi_{c}^{2}}\quad x>0$ (16)

from the relation

$\lambda_{n}^{c}=\frac{\xi_{c}}{\xi_{c}^{2}+x_{n}^{2}}$ (17)

The eigenfunction $e_{\lambda_{n}}^{c,m}$ corresponding to the eigenvalue
$\lambda_{n}\neq 0$ is given by

$e_{\lambda_{n}}^{c,m}(\varrho)=\xi_{c}\sin x_{n}\varrho+x_{n}\cos
x_{n}\varrho$ (18)

######  Proof 4.10.

We have, according to (9) and (14)

$W_{c,m}\cdot
u(\varrho)=\frac{1}{2}\int_{0}^{1}e^{-\xi_{c}|\varrho-\varrho^{\prime}|}u(\varrho^{\prime})\,d\varrho^{\prime}$

Deriving twice w.r.t. $\varrho$

$\left(W_{c,m}\cdot u\right)^{(2)}(\varrho)=\xi_{c}^{2}W_{c,m}\cdot
u(\varrho)-\xi_{c}u(\varrho),$

so that if $u$ is an eigenfunction corresponding to the eigenvalue
$\lambda^{c}$ we obtain

$\lambda^{c}u^{(2)}(\varrho)=\xi_{c}\left(\xi_{c}\lambda^{c}-1\right)u(\varrho),$

which shows that if $\lambda^{c}=0$, $u=0$. If $\lambda^{c}\neq 0$ we define
$K:=\xi_{c}\left(\xi_{c}-\frac{1}{\lambda^{c}}\right)$. The eigenfunctions are
solutions to

$u^{(2)}(\varrho)=Ku(\varrho)$ (19)

which are of the form $u(\varrho);=ae^{\sqrt{K}\varrho}+be^{-\sqrt{K}\varrho}$
for some constants $a$ and $b$, real or complex. $\lambda^{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 $\lambda^{c}$.

It can be verified that if $K\geq 0$, i.e. if
$\lambda^{c}\geq\frac{1}{\xi_{c}}$, there are no solutions. Hence we must
assume $\lambda^{c}<\frac{1}{\xi_{c}}$ and hence $K<0$. The solutions to (19)
write

$u(\varrho)=a\sin\sqrt{-K}\varrho+b\cos\sqrt{-K}\varrho$ (20)

A symbolic computation system shows that

$W_{c,m}\cdot u(\varrho)=\lambda^{c}u(\varrho)+\\\
\frac{\lambda^{c}}{2\xi_{c}}\Big{(}e^{\xi_{c}\varrho}e^{-\xi_{c}}\left(-a\left(\sqrt{-K}\cos\sqrt{-K}+\xi_{c}\sin\sqrt{-K}\right)+b\left(\sqrt{-K}\sin\sqrt{-K}-\xi_{c}\cos\sqrt{-K}\right)\right)+\\\
e^{-\xi_{c}\varrho}\left(a\sqrt{-K}-b\xi_{c}\right)\Big{)}$

The two conditions

$-a\left(\sqrt{-K}\cos\sqrt{-K}+\xi_{c}\sin\sqrt{-K}\right)+b\left(\sqrt{-K}\sin\sqrt{-K}-\xi_{c}\cos\sqrt{-K}\right)=0,$

and

$a\sqrt{-K}-b\xi_{c}=0$ (21)

are necessary and sufficient to guarantee that $W_{c,m}\cdot
u(\varrho)=\lambda^{c}u(\varrho)$ for all $0<\varrho<1$. Since $a$ and $b$ are
not both equal to 0 (otherwise $u=0$ and hence $\lambda^{c}=0$) we must have

$2\xi_{c}\sqrt{-K}\cos\sqrt{-K}+(\xi_{c}^{2}+K)\sin\sqrt{-K}=0,$

which is equivalent to

$\tan\sqrt{-K}=-\frac{2\xi_{c}\sqrt{-K}}{\xi_{c}^{2}+K}$

Letting $x=\sqrt{-K}$, the eigenvalues are found by solving for $x>0$ the
equation

$\tan x=\frac{2\xi_{c}x}{x^{2}-\xi_{c}^{2}}.$

For each value of $x$, the corresponding value of $\lambda^{c}$ is found to be

$\lambda^{c}=\frac{\xi_{c}}{\xi_{c}^{2}+x^{2}}$

The function $x\to\frac{2\xi_{c}x}{x^{2}-\xi_{c}^{2}}$ decreases monotonically
from $+\infty$ for $x=\xi_{c}^{+}$ to 0 when $x\to+\infty$. The corresponding
curve therfore intersects the curve representing the $\pi$-periodic function
$x\to\tan x$ at a countable number of points $x_{n}$, $n\geq 0$, yielding the
eigenvalues $\lambda^{c}_{n}$ of $W_{c,m}$.

According to (20) and (21) the corresponding eigenfunction is

$e_{\lambda^{c}_{n}}^{c,m}(\varrho)=\xi_{c}\sin x_{n}\varrho+x_{n}\cos
x_{n}\varrho$

##  5\. Symmetries of the model and equivariant bifurcations of the solutions

It follows from (5) that $\overline{V}=0$ is a solution to (12) and (13) for
all values of $\gamma$. We are interested in the problem of determining how
this solution bifurcates when $\gamma$ increases, while allowing us to change
somewhat the value of the threshold $\varepsilon$ in (2).

Because of (12), (13) and (2) we have

$dF(0,\gamma)=-{\rm Id}+\gamma{\rm Sig}(-\varepsilon)W,$

Let $\lambda_{p}$ be the largest positive eigenvalue of $W$. When one
increases $\gamma$ from 0 to $\gamma_{p}\stackrel{{\scriptstyle
def.}}{{=}}1/(\lambda_{p}{\rm Sig}(-\varepsilon))$ the solution
$\overline{V}=0$ bifurcates to another solution.

We investigate the properties of the operator $F$ defined by (13) under the
action of the Euclidean group $E(2)$ restricted to the lattice $\mathcal{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_{\varphi_{c}}^{c}$ acting on the $\varphi$ angle by
$\varphi\to\varphi+\varphi_{c}\mod 1$ and the reflection
$\tau_{c}:\varphi\to-\varphi\mod 1$. An element $\psi$ of this group is of the
form

$\psi=R_{\varphi_{c}}^{c}\tau_{c}^{m},\,m\in\\{0,\,1\\}.$ (22)

The action of $E(2)$ on the space of $\mathcal{L}$-periodic functions is best
understood by considering separately the translations and the rotations. Since
the translations of $\mathbb{R}^{2}$ leave the set of $\mathcal{L}$-periodic
functions invariant and translations in $\mathcal{L}$ fix all
$\mathcal{L}$-periodic functions, the effective action of the group of
translations of $\mathbb{R}^{2}$ is as the 2-torus
$T^{2}=\mathbb{R}^{2}/\mathcal{L}$ which is compact. For the rotations, recall
that the holohedry $H$ of the lattice $\mathcal{L}$ is the largest subgroup of
$O(2)$ that leaves $\mathcal{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 $\mathcal{L}$ invariant is the
compact semi-direct sum $E(\mathcal{L}):=D_{4}\overset{\cdot}{+}T^{2}$.
Therefore, we are interested in the action on $F$ of the compact group
$\mathcal{G}=E(\mathcal{L})\times O(2)$. We have the following Proposition.

######  Proposition 5.1.

The operator $F$ defined in (13) is equivariant w.r.t. the action of the group
$\mathcal{G}$.

######  Proof 5.2.

The action of the element $g=(\varpi,\psi)\in\mathcal{G}$ on $F$, with $\psi$
given by (22), is defined as follows

$g\cdot
F(V,\gamma)(r,\varrho,\varphi)=F(V,\gamma)(\varpi^{-1}r,\varrho,(-1)^{m}(\varphi-\varphi_{c})),\
m\in\\{0,\,1\\},$

and, because of Proposition 4.3, the reader will verify that $F$ is
equivariant w.r.t. the action of $\mathcal{G}$, i.e. that we have

$g\cdot F(V,\gamma)(r,c)=F(g\cdot V,\gamma)(r,c)\ \forall g\in\mathcal{G}.$
(23)

At this point, branches of planforms are usually obtained by applying the
equivariant branching lemma [60, 17, 30, 14] as follows. The kernel
$\mathcal{V}$ of $dF(0,\gamma_{p})$ is $\mathcal{G}$-invariant. We fix an
isotropy subgroup $\Sigma$ of $\mathcal{G}$ (i.e. for which there exists
$v_{0}\in\mathcal{V}$ such that $g\cdot v_{0}=v_{0}$ for all $g\in\Sigma$) and
compute the dimension of the fixed point subspace ${\rm
Fix}_{\mathcal{V}}(\Sigma)$ of $\mathcal{V}$ associated with $\Sigma$ where

${\rm Fix}_{\mathcal{V}}(\Sigma)\stackrel{{\scriptstyle
def.}}{{=}}\left\\{v\in\mathcal{V},\,\sigma\cdot v=v,\
\forall\sigma\in\Sigma\right\\}.$

The equivariant branching lemma states that if $\Sigma$ is an axial subgroup
of $G$, i.e. such that

${\rm dim}\ {\rm Fix}_{\mathcal{V}}(\Sigma)=1,$ (24)

then generically, there is a branch of steady-state solutions to (1) with
symmetries $\Sigma$. 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 $\mathcal{G}$ and the corresponding
fixed point subspaces in the following Proposition.

######  Proposition 5.3.

The isotropy subgroups $\Sigma$ of $\mathcal{G}=E(\mathcal{L})\times O(2)$
fall in two classes: those which do not contain the color reflection
$\tau_{c}$ are those of $\mathcal{H}:=E(\mathcal{L})\times SO(2)$, and those
that do contain $\tau_{c}$. In the first case the corresponding fixed
subspaces are those of $\Sigma$ as a subgroup of $\mathcal{H}$. In the second
case the fixed subspaces are those of $\Sigma\cap\mathcal{H}$ projected on
$\mathcal{V}_{1}$ along $\mathcal{V}_{-1}$, the eigenspaces of $\tau_{c}$
corresponding to the eigenvalues 1 and -1, respectively.

######  Proof 5.4.

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(\mathcal{L})$ symmetry [55, 23]. We note
that the fixed point subspaces have even dimension since $SO(2)$ commutes with
$E(\mathcal{L})$ and all real finite dimensional non trivial representation of
$SO(2)$ are of even dimension. All fixed point subspaces of $\mathcal{V}$ for
the action of $\mathcal{H}$ are therefore of even dimension.

If we now extend $SO(2)$ to $O(2)$, i.e. $\mathcal{H}$ to $\mathcal{G}$, given
an isotropy subroup $\Sigma$ of $\mathcal{G}$, either it is an isotropy
subgroup of $\mathcal{H}$ and we are done, or it is not and the color
reflection operator $\tau_{c}$ is not reduced to the identity on ${\rm
Fix}_{\mathcal{V}}(\Sigma)$. But $\tau_{c}$ is a projector. It has therefore
two eigenspaces $\mathcal{V}_{\pm 1}$ corresponding to its eigenvalues $\pm 1$
and $\mathcal{V}$ is the direct sum of the corresponding eigenspaces. Let
$v\in{\rm Fix}_{\mathcal{V}}(\Sigma\cap\mathcal{H})$. We write
$v=v_{1}+v_{-1}$ where $v_{\pm 1}\in\mathcal{V}_{\pm 1}$. We have

$\tau_{c}\cdot v=\tau_{c}\cdot v_{1}+\tau_{c}\cdot v_{-1}=v_{1}-v_{-1},$ (25)

so that the projection $p_{1}\left({\rm
Fix}_{\mathcal{V}}(\Sigma\cap\mathcal{H})\right)$ of ${\rm
Fix}_{\mathcal{V}}(\Sigma\cap\mathcal{H})$ on $\mathcal{V}_{1}$ along
$\mathcal{V}_{-1}$ is included in ${\rm Fix}_{\mathcal{V}}(\Sigma)$.

Conversely, let $v\in{\rm Fix}_{\mathcal{V}}(\Sigma)$. The condition $v\in{\rm
Fix}_{\mathcal{V}}(\Sigma)$ and (25) impose $\tau_{c}\cdot v=v$, i.e.
$\tau_{c}\cdot v_{-1}=0$, i.e. $v_{-1}=0$. so that $v=v_{1}$ and ${\rm
Fix}_{\mathcal{V}}(\Sigma)\subset p_{1}\left({\rm
Fix}_{\mathcal{V}}(\Sigma\cap\mathcal{H})\right)$ .

######  Remark 5.5.

In fact, writing $O(2)=SO(2)\oplus\mathbb{Z}_{2}$, with
$\mathbb{Z}_{2}=\\{1,-1\\}$ the (closed) subgroups of $O(2)$ fall in three
classes [36]

I:

Closed subgroups of $SO(2)$.

II:

Closed subgroups containing $-1$: they are of the form
$\Sigma\oplus\mathbb{Z}_{2}$, $\Sigma$ 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 5.3 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 $\Sigma$
are determined by pairs $K\subset H$ of subgroups of $SO(2)$ such that $K$ has
index 2 in $H$, see also [46]. If $p:SO(2)\oplus\mathbb{Z}_{2}\to SO(2)$,
$p(\psi)=det(\psi)\psi$, is the projection, then $H=p(\Sigma)$ is isomorphic
to $\Sigma$ and $\Sigma=K\cup\psi K$ for any $-\psi\in H\backslash 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(\mathcal{L})\times O(2)$ with a nonzero
fixed subspace, basically because the color rotations act on $\mathcal{V}$ by
multiplications with a magnitude 1 complex exponential.

Since $\mathcal{V}$ is $\mathcal{G}$-invariant and $\mathcal{G}$ is compact,
we can write $\mathcal{V}$ as a direct sum of $\mathcal{G}$-irreducible
subspaces (subspaces such that their only $\mathcal{G}$-invariant subspaces
are 0 or the whole subspace), see e.g. [31][Theorem 1.22]:

$\mathcal{V}=\mathcal{V}_{1}\oplus\cdots\oplus\mathcal{V}_{\rho}.$

It follows easily that

${\rm Fix}_{\mathcal{V}}(\Sigma)={\rm
Fix}_{\mathcal{V}_{1}}(\Sigma)\oplus\cdots\oplus{\rm
Fix}_{\mathcal{V}_{\rho}}(\Sigma).$

Hence if ${\rm dim}({\rm Fix}_{\mathcal{V}}(\Sigma))=1$ then ${\rm
Fix}_{\mathcal{V}}(\Sigma)={\rm Fix}_{\mathcal{V}_{j}}(\Sigma)$ for some $j$
so that the first step in classifying the planforms associated with a fixed
lattice $\mathcal{L}$ is to enumerate each $\mathcal{G}$-irreducible subspace
of $\mathcal{V}$. Moreover $\mathcal{V}_{j}$ is of the form
$\mathcal{V}^{s}_{j}\times\mathcal{V}^{c}_{j}$ where $\mathcal{V}^{s}_{j}$ is
an eigenspace of $W_{s}$ irreducible under the action of $E(\mathcal{L})$ and
$\mathcal{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(\mathcal{L})$ by several authors
including [21]. Dropping the upper index $j$ for simplicity, the irreducible
representations of $\mathcal{V}^{s}$ must be of the form

$\mathcal{V}^{s}=\mathcal{V}_{K_{1}}^{s}\oplus\cdots\oplus\mathcal{V}_{K_{\rho}}^{s},$

where $\mathcal{V}_{K_{j}}=Vect\left\\{ze^{2\pi i\langle
K_{j},r\rangle}+c.c.,\ z\in\mathbb{C}\right\\}$ , $j=1,\cdots,\rho$.
$K_{1},\cdots,K_{\rho}$ are elements of $\mathcal{L}^{*}$ and the set $\\{\pm
K_{1},\cdots,\pm K_{\rho}\\}$ is an orbit in $\mathcal{L}^{*}$ of the
holohedry $D_{4}$. As proved in [21], there exists only one such
representation of dimension 4 and a countable number of them of dimension 8 if
we add the further constraint that they have to be translation free. A
representation is translation free if there are no (nontrivial) translations
in $E(\mathcal{L})$ that act trivially on (28). This requirement ensures that
we have found the finest lattice, $\mathcal{L}$, that supports the neutral
modes (28) [21].

The first one corresponds to

$\left\\{\begin{array}[]{ccc}K_{1}&=&k_{1}\\\
K_{2}&=&k_{2},\end{array}\right.$ (26)

the second to

$\left\\{\begin{array}[]{ccc}K_{1}&=&uk_{1}+vk_{2}\\\
K_{2}&=&-vk_{1}+uk_{2}\\\ K_{3}&=&vk_{1}+uk_{2}\\\
K_{4}&=&-uk_{1}+vk_{2},\end{array}\right.$ (27)

where $u$ and $v$ are relatively prime strictly positive integers such that
$u+v$ is odd.

With our notations, we have

$\mathcal{V}=\\\ \left\\{u(r,\varrho,\phi)=\sum_{j=1}^{k}z_{j}e^{2\pi i\langle
K_{j},\,r\rangle}e^{2\pi in\varphi}e^{h,m}(\varrho)+w_{j}e^{-2\pi i\langle
K_{j},\,r\rangle}e^{2\pi
in\varphi}e^{h,m}(\varrho)+c.c,\,z_{j},\,w_{j}\in\mathbb{C}\right\\},\\\
k=2,\,4,$ (28)

which is isomorphic to $\mathbb{C}^{2k}$ by
$u\to(z,w):=(z_{1},\cdots,z_{s},\,w_{1},\cdots,w_{s})$.

If we now consider an isotropy subgroup $\Sigma$ of $\mathcal{H}$, such that
${\rm dim\ Fix}_{\mathcal{V}}(\Sigma)=2$, Proposition 7.2 in Chapter XVI in
[30] asserts that $\Sigma$ is a twisted group, i.e. there is a pair of
subgroups $K\subset G$ of $E(\mathcal{L})$ and a unique homomorphism
$\Theta:G\to SO(2)$ such that $\Sigma=G^{\Theta}:=\\{(g,\Theta(g))\in
E(\mathcal{L})\times SO(2),\ g\in G\\}.$ $K$ is the kernel of $\Theta$:
$K=\ker\Theta$. 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 5.3 we can obtain the corresponding subgroups of $\mathcal{G}$ and
their corresponding fixed point subspaces.

In detail $D_{4}$ is generated by the $\pi/2$ rotation, noted $R_{\pi/2}^{s}$
and the reflection, noted $\tau_{r_{1}}$ through the $r_{1}$-axis. The action
of $E(\mathcal{L})$ on $\mathcal{V}$ induces an action on $\mathbb{C}^{2k}$.
The reader will verify that this action is given by

$R_{\pi/2}^{s}\cdot(z,w)=\left\\{\begin{array}[]{cc}(z_{2},w_{1},w_{2},z_{1})&k=2\\\
(z_{2},w_{1},w_{2},z_{1},z_{4},w_{3},w_{4},z_{3})&k=4\end{array}\right.,$ (29)

and

$\tau_{r_{1}}\cdot(z,w)=\left\\{\begin{array}[]{cc}(z_{1},w_{2},w_{1},z_{2})&k=2\\\
(w_{4},w_{3},w_{2},w_{1},z_{4},z_{3},z_{2},z_{1})&k=4\end{array}\right..$ (30)

Similarly, for the action of $O(2)$, we have

$R_{\varphi_{c}}^{c}\cdot(z,w)=e^{2\pi in\varphi_{c}}(z,w),\ k=2,\,4$ (31)

and

$\tau_{c}\cdot(z,w)=(\bar{w},\bar{z}),\ k=2,\,4.$ (32)

Writing

$(z_{j},w_{j})=\frac{1}{2}(z_{j}+\bar{w}_{j},\bar{z}_{j}+w_{j})+\frac{1}{2}(z_{j}-\bar{w}_{j},-\bar{z}_{j}+w_{j}),\,j=1,\cdots,k,$

shows that the projection of $(z,w)$ on $\mathcal{V}_{1}$ is
$\frac{1}{2}(z+\bar{w},\bar{z}+w)$.

The results in [23] allow us to determine all the axial subgroups of
$\mathcal{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
$\frac{1}{2}(k_{1}+k_{2})$, $e^{s}$ the identity of $D_{4}$ and $e^{c}$ the
identity of $O(2)$.

| Generators of $G^{\theta}$ |  | $p_{1}\left({\rm Fix}\left(G^{\Theta}\right)\right)$
---|---|---|---
|  | $\operatorname{dim}\mathcal{V}=4$ |
$\mathrm{SR}$ | $((-e^{s},0),0),((e^{s},\frac{1}{2}),R^{c}_{\frac{1}{2}})$ |  | $z_{1}=w_{1}{\rm real},\,z_{2}=w_{2}=0$
$\mathrm{\leavevmode\nobreak\ S}2$ | $((R_{\pi/2}^{s},0),e^{c}),\,((\tau_{v_{1}},0),e^{c}),\,((e^{s},v_{d}),R_{\frac{1}{2}}^{c})$ |  |
| $((R_{\pi/2}^{s},0),\tau_{c}),\,((\tau_{v_{1}},0),\tau_{c}),\,((e^{s},v_{d}),R_{\frac{1}{2}}^{c}\tau_{c})$ |  | $z_{1}=z_{2}=w_{1}=w_{2},\,z_{1}{\rm real}$
|  | $\operatorname{dim}\mathcal{V}=8$ |
$\mathrm{S}2_{u,v}$ | $((R_{\pi/2}^{s},0),e^{c}),\,((\tau_{v_{1}},0),e^{c}),\,((e^{s},v_{d}),R_{\frac{1}{2}}^{c})$ |  |
| $((R_{\pi/2}^{s},0),\tau_{c}),\,((\tau_{v_{1}},0),\tau_{c}),\,((e^{s},v_{d}),R_{\frac{1}{2}}^{c}\tau_{c})$ |  | $z_{1}=z_{2}=z_{3}=z_{4},w=z,\,z_{1}{\rm real}$
$\mathrm{S}4_{u,v}$ | $((R_{\pi/2}^{s},0),e^{c}),\,((\tau_{v_{1}},v_{d}),e^{c}),\,((e^{s},v_{d}),R_{\frac{1}{2}}^{c})$ |  |
| $((R_{\pi/2}^{s},0),\tau_{c}),\,((\tau_{v_{1}},v_{d}),\tau_{c}),\,((e^{s},v_{d}),R_{\frac{1}{2}}^{c}\tau_{c})$ |  | $z_{1}=z_{2}=-z_{3}=-z_{4},w=z,\,z_{1}{\rm real}$
Table 2. The axial subgroups of $\mathcal{G}$ in dimension 4 and two of them
in dimension 8, as well as their corresponding fixed spaces.

We can now apply the equivariant branching lemma . In detail we have the
following Proposition.

######  Proposition 5.6.

Let us denote by $(\lambda_{n})$ the eigenvalues of
$W\in\mathcal{L}(\mathbb{L}^{2}(\Omega))$ and $\gamma_{n}:=\frac{1}{{\rm
Sig}^{\prime}(-\varepsilon)\lambda_{n}}$ the bifurcation points. The smooth
$\mathcal{G}$-equivariant mapping
$F:\mathbb{W}^{m,2}(\Omega)\times\mathbb{R}\to\mathbb{L}^{2}(\Omega)$
satisfies the assumptions of the equivariant branching lemma. For each axial
subgroup $\Sigma\subset\mathcal{G}$, there exists a unique branch of solutions
to $F(V,\gamma_{n})=0$ emanating from $(0,\gamma_{n})$ where the symmetry of
the solution is $\Sigma$. Moreover each of these branches stems from a
pitchfork bifurcation.

######  Proof 5.7.

The proof is a direct application of [14][Theorem 2.3.2] around
$\gamma\approx\gamma_{n}$. From Section 3 the mapping $F$ is smooth. Its
Jacobian at $V=0$ is the linear operator

$dF(0,\gamma_{n})=-{\rm Id}+\gamma_{n}{\rm Sig}^{\prime}(-\varepsilon)W,$

where ${\rm Sig}$ is defined in (2). As $\mathbb{W}^{m,2}(\Omega)$ is
compactly included in $\mathbb{L}^{2}(\Omega)$ and
$W\in\mathcal{L}(\mathbb{L}^{2}(\Omega))$ is compact, it follows that $W$ is
compact from $\mathbb{W}^{m,2}(\Omega)$ to $\mathbb{L}^{2}(\Omega)$. This
implies that $dF(0,\gamma_{n})$ is a Fredholm operator with index $0$.

From the previous linear analysis, $0$ is an isolated eigenvalue of
$dF(0,\gamma_{n})\in\mathcal{L}(\mathbb{L}^{2}(\Omega))$ with finite
multiplicity.

Let us denote the reduced equation by $f(v,\gamma)=0,\ v\in\mathcal{V}$: it is
given by the Lyapunov-Schmidt method, see e.g. [29]. It remains to show that
$d^{2}_{v\gamma}f(0,0)\neq 0$ to complete the proof of the generic
existence/uniqueness of a bifurcating branch.

We restrict ourselves to the case $\varepsilon=0$ for which
$s^{\prime\prime}(0)=0$. The methods described in e.g.[29][Chapter VII, page
295], allow us to find, using the fact that $dF$ is hermitian, that
$d^{2}_{v\gamma}f(0,0)=\left\langle\zeta,\partial_{\gamma}dF(0,0)\cdot\zeta\right\rangle$
where $\zeta\in\mathcal{V}$ is normalised. A direct computation gives
$d^{2}_{v\gamma}f(0,0)=\frac{1}{\gamma_{n}}\neq 0$.

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^{\theta}$ of
$\mathcal{G}$ which appear in Table 2 we verify (using the relation
$\tau_{c}R_{\frac{1}{2}}^{c}=R_{\frac{1}{2}}^{c}\tau_{c}$) that the element
$g_{0}:=((e^{s},0),R_{\frac{1}{2}}^{c})$ is in the normalizer $N(G^{\theta})$
of $G^{\theta}$ since it satisfies $g_{0}gg_{0}^{-1}=g$ for all generators of
$G^{\theta}$ glven in Table 2 and acts as $-\mathbf{1}$ on $p_{1}\left({\rm
Fix}(G^{\theta})\right)$.

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 ${\rm dim}({\rm Fix}_{\mathcal{V}}(\Sigma))>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, 16, 15].

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

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

##  6\. 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 ${\rm dim}\mathcal{V}=4$, we
have

$SR(r,\varrho,\varphi)=\cos 2\pi r_{1}e^{h,m}(\varrho)\cos(2\pi n\varphi),$
(33)

and

$S2(r,\varrho,\varphi)=\left(\cos 2\pi r_{1}+\cos 2\pi
r_{2}\right)e^{h,m}(\varrho)\cos(2\pi n\varphi).$ (34)

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 7.2. Similarly
we show in Figure 3 an example of spots described by equation (34).

In the case ${\rm dim}\mathcal{V}=8$ we have the two families of planforms

$S2_{u,v}(r,\varrho,\varphi)=\Big{(}\cos 2\pi(ur_{1}+vr_{2})+\cos
2\pi(-vr_{1}+ur_{2})+\\\ \cos 2\pi(vr_{1}+ur_{2})+\cos
2\pi(-ur_{1}+vr_{2})\Big{)}e^{h,m}(\varrho)\cos(2\pi n\varphi),$ (35)

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

and

$S4_{u,v}(r,\varrho,\varphi)=\Big{(}\cos 2\pi(ur_{1}+vr_{2})+\cos
2\pi(-vr_{1}+ur_{2})-\\\ \cos 2\pi(vr_{1}+ur_{2})-\cos
2\pi(-ur_{1}+vr_{2})\Big{)}e^{h,m}(\varrho)\cos(2\pi n\varphi),$ (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}$.

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

######  Remark 6.1.

To determine the color represented at cortical location $(r_{1},r_{2})$, and
produce the images in Figures 2-5, we compute $c_{\rm max}(r_{1},r_{2})={\rm
argmax}_{\varrho,\varphi}\Upsilon\left(V^{\rm
eq}(r_{1},r_{2},\varrho,\varphi)\right)\in\Omega_{c}$ (corresponding to the
maximum activity of the neural mass with spatial coordinates $(r_{1},r_{2})$),
$V^{\rm eq}\in\\{SR,\,S2,\,S2_{u,v},\,S4_{u,v}\\}$, and define the Luminance
to be the corresponding activity level $\Upsilon\left(V^{\rm
eq}(r_{1},r_{2},c_{\rm max}(r_{1},r_{2}))\right)$. 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]:

$\displaystyle H$ $\displaystyle=2\pi\varphi$ $\displaystyle S$
$\displaystyle=\sqrt{\varrho}$ $\displaystyle L$
$\displaystyle=\Upsilon\left(V^{\rm eq}(r_{1},r_{2},c_{\rm
max}(r_{1},r_{2}))\right)$

######  Remark 6.2.

In Figures 2-5 the value of $n$ in Proposition 4.7 is equal to 6. The value of
$\xi_{c}$ in Proposition 4.9 is equal to 2 and $x_{n}$ is the fourth strictly
positive root of equation (16). This implies that the product
$e^{h,m}(\varrho)\cos(2\pi n\varphi)$ has several minima and maxima so that
there are several values of $(\varrho,\varphi)$ where the maximum of
$\Upsilon(SR(r,\varrho,\varphi))$(respectively of
$\Upsilon(S2(r,\varrho,\varphi))$, $\Upsilon(S2_{u,v}(r,\varrho,\varphi))$ or
$\Upsilon(S4_{u,v}(r,\varrho,\varphi))$) with respect to $(\varrho,\varphi)$
is reached. As in Section 7.1 we have chosen arbitrarily one of them.

##  7\. 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_{\rm 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,\gamma_{n})$ in $\mathcal{F}\times\mathbb{R}$. The size of this
neighborhood quantifies the predictability of the local theory and is bounded
by $\min(|\gamma_{n+1}-\gamma_{n}|,|\gamma_{n-1}-\gamma_{n}|)$. 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_{\rm 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
$\varphi_{0}$ (between 0 and 1 to be consistent with our notations). A point
on this diameter is characterized by its polar coordinates
$(2\pi\varphi_{0},\rho)$ or $(2\pi(\varphi_{0}+\frac{1}{2}\mod 1),\rho)$,
$\rho\in[0,1]$. The relevant color coordinate is therefore
$c:=\pm\rho\in[-1,\,1]$.

Note that the problem associated has now the symmetry group
$(D_{4}\overset{\cdot}{+}T_{2})\times\mathbb{Z}_{2}$ instead of
$(D_{4}\overset{\cdot}{+}T_{2})\times O(2)$. This has been extensively studied
[56, 22]. A key remark is that the additional reflection symmetry has no
effect on the bifurcation problems associated with the square lattice666This
was pointed out to us by Pascal Chossat..

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.

###  7.1. 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^{\rm 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})={\rm
argmax}_{c}\Upsilon\left(V^{\rm eq}(r_{1},r_{2},c)\right)\in[-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}):=\Upsilon\left(V^{\rm
eq}(r_{1},r_{2},c(r_{1},r_{2}))\right)$. 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]:

$\left\\{\begin{array}[]{lcl}H&=&2\pi\varphi_{0}+\frac{\pi}{2}\times(1-{\rm
sign}(c(r_{1},r_{2})))\\\ S&=&|c(r_{1},r_{2})|\\\
L&=&a(r_{1},r_{2})\end{array}\right.$ (37)

where

${\rm sign}(x)=\left\\{\begin{array}[]{ll}+1&\text{if }x\geq 0\\\
-1&\text{otherwise}\end{array}\right.$

and obtain the result shown in Figure 6-Middle. $\varphi_{0}=\frac{1}{8}$ is
defined in Appendix A.

######  Remark 7.1.

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\to\Upsilon(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 7.2.

The visual cortex of several species (like monkeys) has the property of being
retinotopically organized, e.g. [53, 48]. 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.

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.

###  7.2. Connectivity functions

In all subsequent examples, we use the following spatial connectivity function
$w_{s}(r)$ of the form (4) whose parameters $\mu_{s}$, $\nu_{s}$ and
$\beta_{s}$ are chosen so that it satisfies (5):

$w_{s}(r)=A\left(\mu_{s}e^{-\|r\|_{2}^{2}/2\pi^{2}}-\nu_{s}e^{-\|r\|_{2}^{2}/2(\beta_{s}\pi)^{2}}\right)\mathbf{1}_{\|r\|\leq
B}.$ (38)

The specific values of these parameters corresponding to different numerical
experiments are shown in the captions of Figures 8\- 11. A cross-section is
shown in Figure 7-Left: it has the traditional Mexican hat shape.

Similarly $w_{c}(c,c^{\prime})$ is defined by

$w_{c}(c,c^{\prime}):=\mu_{c}e^{-\alpha_{c}|c-c^{\prime}|}-\nu_{c}e^{-\beta_{c}|c+c^{\prime}|}$

with the following values of the parameters:

$\displaystyle\alpha_{c}=0.3\quad$ $\displaystyle\beta_{c}=0.4$
$\displaystyle\mu_{c}=0.6,\quad$ $\displaystyle\nu_{c}=0.69.$

A heatmap of this function is shown in Figure 7-Right.

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^{\prime})$ used in the numerical experiments.

###  7.3. 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 $\gamma$, 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 [8, 7]).

In order to compute the bifurcating branches, the reduced equations (see [30]
or proof of proposition 5.6) 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 $\approx 100$, due
to the symmetries, in order to compute the eigenvalues, and this limited the
size $N_{r_{1}}\times N_{r_{2}}\times N_{c}$ of the discretisation.

###  7.4. First bifurcation diagram

Figure 8. Parameters
$\varepsilon=0,\,\Omega_{s}=(-20\pi,20\pi)^{2},\,\Omega_{c}=(-1,1),\,\beta_{s}=1.9,\,B=3\pi,\,A=1.0$.
Equilibria as functions of the nonlinear gain $\gamma$. The stable parts of
the branches are indicated with thick lines. The dots indicate bifurcation
points. The discretization is $256\times 256\times 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 8 represents the variation of the equilibria (represented by their
$\mathbb{L}^{1}$ norm) as functions of the slope $\gamma$ of the activation
function $s$, see (2). The trivial equilibrium $\bar{V}=0$ loses stability at
a first bifurcation point around $\gamma\approx 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).

###  7.5. Second bifurcation diagram

Figure 9.
$\varepsilon=0.5,\,\Omega_{s}=(-20\pi,20\pi)^{2},\,\Omega_{c}=(-1,1),\,\beta_{s}=20.9,\,B=6\pi,\,A=1.0.$
Equilibria are shown as functions of the nonlinear gain $\gamma$. 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.

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 the previous diagram. Using aBS (see
Section 7.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.

###  7.6. Third bifurcation diagram, changing the criticality

Figure 10.
$\varepsilon=2.3,\,\Omega_{s}=(-20\pi,20\pi)^{2},\,\Omega_{c}=(-1,1),\,\beta_{s}=1.9,\,B=3\pi,\,A=9.$
Equilibria as function of the nonlinear gain $\gamma$. 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.

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]

$\left\\{\begin{array}[]{l}\dot{z}_{1}=z_{1}\left(\alpha+\beta_{1}\left|z_{1}\right|^{2}+\beta_{2}\left|z_{2}\right|^{2}\right)\\\
\dot{z}_{2}=z_{2}\left(\alpha+\beta_{1}\left|z_{2}\right|^{2}+\beta_{2}\left|z_{1}\right|^{2}\right)\end{array}\right.$

where $z_{i}\in\mathbb{C},\alpha,\beta_{i}\in\mathbb{R}$, $i=1,\,2$, and (for
example for $\beta_{1}$, see [62] )

$\beta_{1}/\gamma_{0}^{3}\lambda_{0}=\gamma_{1}s_{2}^{2}A+s_{3}/2$

where
$A\in\mathbb{R},s_{2}=Sig^{(2)}(\varepsilon),s_{3}=Sig^{(3)}(\varepsilon)$,
$\gamma_{0}$ is the value of $\gamma$ at the bifurcation point and
$\lambda_{0}$ is the associated eigenvalue of $W$. We recall that for a
threshold $\varepsilon=0$, $s_{2}=0,s_{3}<0$ and thus the Pitchfork is
supercritical. When varying the threshold $\varepsilon$, the ratio
$\frac{s_{3}}{s_{2}^{2}}$ takes on any value, thus we can change the
criticality by altering the threshold. Compared to the previous diagram, we
use $\varepsilon=2.5$ and scale the connectivity to restrict the value of
$\gamma$ at the first bifurcation point.

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 saddle-node 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.

###  7.7. 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
$\epsilon=\gamma_{1}-\gamma_{1}^{SN}$ where $\gamma_{1}$ is the value of
$\gamma$ at the 4D bifurcation point and $\gamma_{1}^{SN}$ 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 $\varepsilon$. By performing codim 2 continuation of the
saddle-node bifurcation and the first primary bifurcation point in the
$(\gamma,\varepsilon)$ plane, we were able to select a better value of the
threshold $\varepsilon$ 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].

Figure 11.
$\varepsilon=2.9,\,\Omega_{s}=(-20\pi,20\pi)^{2},\,\Omega_{c}=(-1,1),\,\beta_{s}=1.9,\,B=3\pi,\,A=9.$
Equilibria as functions of the nonlinear gain $\gamma$. 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.
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.

##  8\. 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
[65, 57] 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 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.

##  Appendix A Color space

It is not important for this work to precisely define the color space we use
because the phenomena described in this paper are robust to changes of this
definition as long as the color space can be thought of as the product of a
luminance axis with a 2D opponent chromaticity space which is assumed to be a
disc of radius one centered at the origin. There are many possible choices for
the details of such space. For convenience we have adopted the HSL space for
Hue, Saturation and Luminance [37]. The (Hue, Saturation) coordinates are
precisely the polar coordinates in the 2D chromaticity space, $2\pi\varphi$ is
the polar angle, varying between $0$ and $2\pi$ since $\varphi$ varies between
0 and 1, and $\rho$ is the distance to the origin, varying between 0 and 1. In
Section 7 for minimizing computer memory we restrict ourselves to a diameter
of the unit disc, i.e. $\varphi=\varphi_{0}=\frac{1}{8}$. The conversion from
the HSL space to the usual RGB space is nonlinear and can be found in
textbooks such as [69, 37, 40].

## Acknowledgments

The authors thank Pascal Chossat for useful discussions. They are also
grateful to the OPAL infrastructure from UniversitÃ© CÃ´te d’Azur for
providing resources and support.
The third author would like to thank Hannes Uecker for his help.

## References

  * [1] R. Adams, _Sobolev spaces_ , Pure and Applied Mathematics, Series of Monographs and Textbooks, vol. 65, Academic Press, Inc., New York, San Francisco, London, 1975.
  * [2] R. Adams, J. Fournier, “Sobolev Spaces”, 2003.
  * [3] S.-I. Amari, “Dynamics of pattern formation in lateral-inhibition type neural fields”, _Biological Cybernetics_ 27 (1977), no. 2, p. 77-87.
  * [4] D. Avitabile, D. Lloyd, J. Burke, E. Knobloch, B. Sandstede, “To snake or not to snake in the planar Swift–Hohenberg equation”, _SIAM Journal on Applied Dynamical Systems_ 9 (2010), p. 704.
  * [5] R. Barret, M. Berry, T. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine, H. van der Vonst, _Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods_ , SIAM, 1994.
  * [6] M. Berthier, “Geometry of color perception. Part 2: perceived colors from real quantum states and Hering’s rebit.”, _J. Math. Neurosc._ 10 (2020), no. 14.
  * [7] T. Besard, V. Churavy, A. Edelman, B. De Sutter, “Rapid software prototyping for heterogeneous and distributed platforms”, _Advances in Engineering Software_ 132 (2019), p. 29-46.
  * [8] T. Besard, C. Foket, B. De Sutter, “Effective Extensible Programming: Unleashing Julia on GPUs”, _IEEE Transactions on Parallel and Distributed Systems_ (2018).
  * [9] P. Bressloff, “Spatiotemporal dynamics of continuum neural fields”, _Journal of Physics A: Mathematical and Theoretical_ 45 (2012).
  * [10] P. Bressloff, J. Cowan, M. Golubitsky, P. Thomas, M. Wiener, “Geometric visual hallucinations, Euclidean symmetry and the functional architecture of striate cortex”, _Phil. Trans. R. Soc. Lond. B_ 306 (2001), no. 1407, p. 299-330, http://dx.doi.org/doi:10.1098/rstb.2000.0769.
  * [11] F. Busse, “Patterns of convection in spherical shells”, _Journal of Fluid Mechanics_ 72 (1975), no. 1, p. 67-85.
  * [12] S. Chapman, G. Kozyreff, “Exponential asymptotics of localised patterns and snaking bifurcation diagrams”, _Physica D: Nonlinear Phenomena_ 238 (2009), no. 3, p. 319-354.
  * [13] S. Chatterjee, K. Ohki, R. Reid, “Chromatic micromaps in primary visual cortex”, _Nature Communications_ 12 (2021), no. 2315, https://doi.org/10.1038/s41467-021-22488-3.
  * [14] P. Chossat, R. Lauterbach, _Methods in Equivariant Bifurcations and Dynamical Systems_ , World Scientific Publishing Company, 2000.
  * [15] P. Chossat, R. Lauterbach, I. Melbourne, “Steady-state bifurcation with 0 (3)-symmetry”, _Archive for Rational Mechanics and Analysis_ 113 (1990), no. 4, p. 313-376.
  * [16] P. Chossat, “Solutions avec symétrie diédrale dans les problèmes de bifurcation invariants par symétrie sphérique”, _Comptes rendus des séances de l’Académie des sciences. Série 1, Mathématique_ 297 (1983), no. 12, p. 639-642.
  * [17] G. Cicogna, “Symmetry breakdown from bifurcation”, _Lettere al Nuovo Cimento (1971-1985)_ 31 (1981), no. 17, p. 600-602.
  * [18] C. Colby, J. Duhamel, M. Goldberg, “Oculocentric spatial representation in parietal cortex”, _Cereb. Cortex_ 5 (1995), p. 470-481.
  * [19] S. Coombes, P. beim Graben, R. Potthast, “Tutorial on Neural Field Theory”, in _Neural Fields: Theory and Applications_ (S. Coombes, P. beim Graben, R. Potthast, J. Wright, eds.), Springer Berlin Heidelberg, Berlin, Heidelberg, 2014, p. 1-43.
  * [20] S. Denève, C. K. Machens, “Efficient codes and balanced networks”, _Nature neuroscience_ 19 (2016), no. 3, p. 375-382.
  * [21] B. Dionne, M. Golubitsky, “Planforms in two and three dimensions”, _ZAMP_ 43 (1992), p. 36-62.
  * [22] B. Dionne, M. Silber, A. Skeldon, “Stability results for steady, spatially periodic planforms”, _Nonlinearity_ 10 (1997), p. 321.
  * [23] B. Dionne, M. Golubitsky, M. Silber, I. Stewart, “Time-periodic spatially periodic Planforms in Euclidean equivariant partial differential equations”, _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ 352 (1995), p. 125-168.
  * [24] B. Ermentrout, “Neural networks as spatio-temporal pattern-forming systems”, _Reports on Progress in Physics_ 61 (1998), p. 353-430.
  * [25] G. Ermentrout, J. Cowan, “Temporal oscillations in neuronal nets”, _Journal of mathematical biology_ 7 (1979), no. 3, p. 265-280.
  * [26] P. E. Farrell, Á. Birkisson, S. W. Funke, “Deflation Techniques for Finding Distinct Solutions of Nonlinear Partial Differential Equations”, _SIAM Journal on Scientific Computing_ 37 (2015), no. 4, p. A2026-A2045, https://epubs.siam.org/doi/abs/10.1137/140984798, Publisher: Society for Industrial and Applied Mathematics.
  * [27] O. Faugeras, R. Veltz, F. Grimbert, “Persistent neural states: stationary localized activity patterns in nonlinear continuous n-population, q-dimensional neural networks”, _Neural Computation_ 21 (2009), no. 1, p. 147-187.
  * [28] S. Funahashi, C. Bruce, P. Goldman-Rakic, “Mnemonic coding of visual space in the monkey’s dorsolateral prefrontal cortex”, _J. Neurophysiol._ 61 (1989), p. 331-349.
  * [29] M. Golubitsky, D. Schaeffer, _Singularities and Groups in Bifurcation Theory_ , vol. I, Springer, 1984.
  * [30] M. Golubitsky, I. Stewart, D. Schaeffer, _Singularities and Groups in Bifurcation Theory_ , vol. II, Springer, 1988.
  * [31] M. Golubitsky, I. Stewart, _The Symmetry Perspective_ , Birkhauser Verlag, 2000.
  * [32] E. Hering, _Outlines of a theory of the light sense._ , Harvard University Press, 1964.
  * [33] R. Hoyle, _Pattern formation: an introduction to methods_ , Cambridge Univ Pr, 2006.
  * [34] D. H. Hubel, “Blobs and color vision”, in _Bioscience at the Physical Science Frontier_ , Springer, 1986, p. 91-102.
  * [35] D. Hubel, _Eye, Brain, and Vision_ , vol. 22, Scientific American Library, 1995\.
  * [36] E. Ihrig, M. Golubitsky, “Pattern selection with O(3) symmetry”, _Physica D: Nonlinear Phenomena_ 13 (1984), no. 1-2, p. 1-33.
  * [37] G. H. Joblove, D. Greenberg, “Color Spaces for Computer Graphics”, _SIGGRAPH Comput. Graph._ 12 (1978), no. 3, p. 20-25, https://doi.org/10.1145/965139.807362.
  * [38] E. R. Kandel, T. M. Jessell, J. H. Schwartz, S. A. Siegelbaum, A. Hudspeth, _Principles of neuroscience_ , 5th ed., McGraw-Hill, 2013.
  * [39] E. Knobloch, “Spatial Localization in Dissipative Systems”, _Annual Review of Condensed Matter Physics_ 6 (2015), no. 1, p. 325-359, http://www.annualreviews.org/doi/10.1146/annurev-conmatphys-031214-014514.
  * [40] J. J. Koenderink, _Color for the Sciences_ , The MIT Press, 2010.
  * [41] J. S. Lund, A. Angelucci, P. C. Bressloff, “Anatomical Substrates for Functional Columns in Macaque Monkey Primary Visual Cortex”, _Cerebral Cortex_ 12 (2003), p. 15-24.
  * [42] E. Miller, C. Erickson, R. Desimone, “Neural mechanisms of visual working memory in prefrontal cortex of the Macaque”, _J. Neurosci._ 16 (1996), p. 5154-5167.
  * [43] P. Monnier, “Standard definitions of chromatic induction fail to describe induction with S-cone patterned backgrounds”, _Vision Research_ 48 (2008), p. 2708-14.
  * [44] P. Monnier, S. K. Shevell, “Chromatic induction from S-cone patterns”, _Vision Research_ 44 (2004), p. 849-56.
  * [45] J. B. Nganou, “How rare are subgroups of index 2?”, _Mathematics Magazine_ 85 (2012), no. 3, p. 215-220.
  * [46] M. Olive, “Effective computation of SO(3) and O(3) linear representation symmetry classes”, _Mathematics and Mechanics of Complex Systems_ 7 (2019), no. 3, p. 203-237.
  * [47] G. Oster, “Phosphenes.”, _Scientific American_ 222 (1970), no. 2, p. 82.
  * [48] J. Petitot, _Elements of neurogeometry_ , Springer, 2017.
  * [49] R. Potthast, P. Beim Graben, “Existence and properties of solutions for neural field equations”, _Mathematical Methods in the Applied Sciences_ 33 (2010), no. 8, p. 935-949.
  * [50] E. Provenzi, “Geometry of color perception. Part 1: structures and metrics of a homogeneous color space.”, _J. Math. Neurosc._ 10 (2020), no. 7.
  * [51] P. Rabinowitz. “Some Global Results for Nonlinear Eigenvalue Problems.”, emph Journal of Functional Analysis 7, (1971)
  * [52] Y. Saad, _Numerical methods for large eigenvalue problems_ , rev. ed ed., Classics in applied mathematics, no. 66, Society for Industrial and Applied Mathematics, Philadelphia, 2011.
  * [53] E. Schwartz, “Spatial mapping in the primate sensory projection: Analytic structure and relevance to perception”, _Biological cybernetics_ 25 (1977), no. 4, p. 181-194.
  * [54] R. Siegel, “Hallucinations.”, _Scientific American_ (1977).
  * [55] M. Silber, E. Knobloch, “Hopf bifurcation on a square lattice”, _Nonlinearity_ 4 (1991), p. 1063.
  * [56] M. Silber, E. Knobloch, “Pattern selection in steady binary-fluid convection”, _Physical Review A_ 38 (1988), no. 3, p. 1468.
  * [57] A. Song, O. Faugeras, R. Veltz, “A neural field model for color perception unifying assimilation and contrast”, _PLoS Comput Biol_ 15 (2019), no. 6.
  * [58] D. Y. Ts’o, M. D. Zarella, G. Burkitt, “Whither the hypercolumn?”, _The Journal of physiology_ 587 Pt 12 (2009), p. 2791-805.
  * [59] H. Uecker, D. Wetzel, “Numerical Results for Snaking of Patterns over Patterns in Some 2D Selkov–Schnakenberg Reaction-Diffusion Systems”, _SIAM Journal on Applied Dynamical Systems_ 13 (2014), no. 1, p. 94-128, https://epubs.siam.org/doi/abs/10.1137/130918484.
  * [60] A. Vanderbauwhede, “Local Bifurcation and Symmetry”, PhD Thesis, Rijksuniversiteit Gent, 1980.
  * [61] N. Vattuone, T. Wachtler, I. Samengo, “Perceptual spaces and their symmetries: The geometry of color space”, _Mathematical Neuroscience and Applications_ 1 (2021), https://mna.episciences.org/7680.
  * [62] R. Veltz, “Nonlinear analysis methods in neural field models”, PhD Thesis, Université Paris Est, 2011, ftp://ftp-sop.inria.fr/neuromathcomp/publications/phds/veltz-11.pdf.
  * [63] R. Veltz, “BifurcationKit.jl”, 2020, https://hal.archives-ouvertes.fr/hal-02902346.
  * [64] R. Veltz, P. Chossat, O. Faugeras, “On the Effects on Cortical Spontaneous Activity of the Symmetries of the Network of Pinwheels in Visual Area V1”, _J. Math. Neurosc._ 5 (2015), no. 11.
  * [65] R. Veltz, O. Faugeras, “Local/Global Analysis of the Stationary Solutions of Some Neural Field Equations”, _SIAM Journal on Applied Dynamical Systems_ 9 (2010), no. 3, p. 954-998, https://epubs.siam.org/doi/pdf/10.1137/090773611.
  * [66] C. Ware, W. B. Cowan, “Changes in perceived color due to chromatic interactions.”, _Vision Research_ 22 11 (1982), p. 1353-62.
  * [67] H. Wilson, J. Cowan, “Excitatory and inhibitory interactions in localized populations of model neurons”, _Biophys. J._ 12 (1972), p. 1-24.
  * [68] H. Wilson, J. Cowan, “A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue”, _Biological Cybernetics_ 13 (1973), no. 2, p. 55-80.
  * [69] G. Wyszecki, W. S. Stiles, _Color Science: Concepts and Methods, Quantitative Data and Formulas_ , John Wiley and Sons, Inc., New York, London, Sidney, 1967.
  * [70] Y. Xiao, A. Casti, J. Xiao, E. Kaplan, “Hue maps in primate striate cortex”, _NeuroImage_ 35 (2007), no. 2, p. 771 - 786, http://www.sciencedirect.com/science/article/pii/S1053811906011608.