# Cortical origins of MacKay-type visual illusions: A case for the non-linearity

Cyprien Tamekue, Dario Prandi, Yacine Chitour

► To cite this version:

Cyprien Tamekue, Dario Prandi, Yacine Chitour. Cortical origins of MacKay-type visual illusions: A case for the non-linearity. 22nd World Congress of the International Federation of Automatic Control (IFAC 2023), Jul 2023, Yokohama, Japan. pp.476-481, 10.1016/j.ifacol.2023.10.1613. hal-03894606

HAL Id: hal-03894606

https://u-paris.hal.science/hal-03894606v1

Submitted on 29 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Cortical origins of MacKay-type visual illusions: A case for the non-linearity*

Cyprien Tamekue*. Dario Prandi*. Yacine Chitour*.

* Université Paris-Saclay, CNRS, CentraleSupélec, Laboratoire des Signaux et Systèmes, 91190, Gif-sur-Yvette, France (e-mail: {cyprien.tamekue, dario.prandi, yacine.chitour}@centralesupelec.fr).

**Abstract:** To study the interaction between retinal stimulation by redundant geometrical patterns and the cortical response in the primary visual cortex (V1), we focus on the MacKay effect (Nature, 1957) and Billock and Tsou's experiments (PNAS, 2007). We use a controllability approach to describe these phenomena starting from a classical biological model of neuronal field equations with a non-linear response function. The external input containing a localised control function is interpreted as a cortical representation of the static visual stimuli used in these experiments. We prove that while the MacKay effect is essentially a linear phenomenon (i.e., the nonlinear nature of the activation does not play any role in its reproduction), the phenomena reported by Billock and Tsou are wholly nonlinear and depend strongly on the shape of the nonlinearity used to model the response function.

**Keywords:** Control in neuroscience, modelling of biological systems, Neuronal field equations, non-linear systems, human visual system, perception and psychophysics, MacKay effect.

## 1. INTRODUCTION

In many situations, humans perceive an illusory component that is not physically present in a visual stimulus. Helmholtz (1867) is probably the first to be interested in the visual effect induced by the presentation of a pattern consisting of black-and-white zones. In particular, he related the perception of rotating darker and brighter radial zones after viewing a pattern consisting of black and white concentric rings to the fluctuation of eye accommodation. In this direction, MacKay (1957) reported striking after-effect of visual stimulation by regular geometrical patterns with highly redundant information (see Fig. 1 for the so-called “MacKay rays”) and attributed the phenomena to some part of the visual cortex, which might profit from such redundancy. In these experiments, an illusory contour consisting of a pattern of white and black concentric rings (tunnel pattern) is evoked by all observers as the after-image induced by a pattern consisting of white and black fan shape (funnel pattern) with high redundant information in the fovea (the centre of the visual field). Due to the retino-cortical map<fnref n="1" /> between the retina and the primary visual cortex (V1, henceforth), the after-images in V1 are superimposed patterns consisting of orthogonal horizontal and vertical stripes. This indicates that neuronal response

<figure><img src="image_2.png" /><figcaption>Fig. 1. MacKay effect: the presentation of the stimulus to the left (“MacKay rays”) induces an illusory perception of the image on the right. Adapted from MacKay (1957) and Zeki et al. (1993).</figcaption></figure>

in V1 tends to favour directions at a right angle to the visual stimulus.

Even more striking visual effects have been obtained in the psychophysical experiments reported by Billock and Tsou (2007), see Fig. 2. As in the case of the MacKay effect, they found that biasing stimuli could induce orthogonal responses in the visual cortex. More precisely, a funnel pattern localised at the fovea (resp. in the periphery) with a background flicker induces the perception of a tunnel pattern in the periphery (resp. at the fovea). The spatial interaction is localized, meaning the hallucination does not extend through the physical stimulus nor into empty non-flickering regions.

This work is concerned with a theoretical description of MacKay and Billock and Tsou's illusory phenomena in V1. This is achieved by studying the properties of the Amari-type equation (Amari, 1977, Eq. (3)) describing the dynamics of the activity $a: \mathbb{R}_+ \times \mathbb{R}^2 \to \mathbb{R}$ on V1:

$$ \frac{\partial a}{\partial t} = -a + \mu\omega * f(a) + I. \qquad (\text{NF}) $$

$$ re^{i\theta} \mapsto (x_1, x_2) := (\log r, \theta). $$

<footnote>* This work has been supported by the ANR-20-CE48-0003. The first author was supported by a grant from the “Fondation CFM pour la Recherche”.</footnote>

<footnote n="1">¹ Let $(r, \theta) \in [0, \infty) \times [0, 2\pi)$ denote polar coordinates in the visual field (or in the retina) and $(x_1, x_2) \in \mathbb{R}^2$ Cartesian coordinates in V1. The retino-cortical map (see, e.g., Tamekue et al. (2022) and references within) is given by</footnote>
Here, * denotes the spatial convolution operation, $\omega: \mathbb{R}^2 \to \mathbb{R}$ is an interaction kernel modelling cortical connections in V1, $f$ is a sigmoid non-linearity, and $\mathrm{I}: \mathbb{R}^2 \to \mathbb{R}$ is the cortical representation of the presented static visual stimulus, that is assumed to be time-independent. Finally, $\mu > 0$ is a parameter measuring the strength of intra-neuron connectivity. In this work, following Tamekue et al. (2022), we assume that the parameter $\mu$ is smaller than the threshold parameter $\mu_c$ where cortical patterns (e.g., funnels, tunnels, spirals, checkerboards, cobwebs, etc...) spontaneously emerge in V1 (see, e.g. Ermentrout and Cowan (1979); Bressloff et al. (2001)). Neurophysiologically seeing, this corresponds to considering an unaltered state where no spontaneous hallucinations emerge.

By an asymptotic analysis of the properties of (NF) we describe why the after-image in the MacKay effect consists of illusory contours in the background of the physical visual stimulus. In particular, the result we provide here implies that a motion in the after-image moves at a right angle to the stimulus pattern. This is because the static physical stimulus and the after-image in V1 are superimposed horizontal and vertical stripes combined with the fact that the inverse retino-cortical map conserves this opponency in the retina.

Our main finding is that while the MacKay effect is essentially a linear phenomenon, Billock and Tsou's experiments are completely non-linear phenomena that strongly depend on the shape of the non-linear function used to model the neuronal response after an activation. Moreover, due to the equivariance of equation (NF) with respect to the plane Euclidean group $\mathbf{E}(2)$, we find that the MacKay effect results from the highly redundant information in the visual stimulus aiming to break its plane Euclidean symmetry. The same is true for Billock and Tsou's phenomenon, where symmetry-breaking arises due to the localization of the visual stimulus in the visual field.

We conclude this section by mentioning that, up to our knowledge, the only other attempt to describe these phenomena theoretically is due to Nicks et al. (2021). There, the authors study a different model of neuronal fields equation containing an adaptation variable and a feedback (state-dependent) external input. Their theoretical result relies on bifurcation and multi-scale analysis, which can be applied only for values of $\mu$ near the threshold parameter $\mu_c$ and in the presence of fully distributed external inputs. In particular, their analysis does not apply to the range of $\mu$ that we consider, nor to localized inputs, such as those used by MacKay and Billock and Tsou. Nevertheless, they provide numerical results showing the capability of their model to reproduce Billock and Tsou's experiments.

**Notation.** In the following, $d \in \{1, 2\}$ is the dimension of $\mathbb{R}^d$ and $|x|$ denote the Euclidean norm of $x \in \mathbb{R}^d$. For $p \in \{1, \infty\}$, $L^p(\mathbb{R}^d)$ is the Lebesgue space of class of real-valued measurable functions $u$ on $\mathbb{R}^d$ such that $|u|$ is integrable over $\mathbb{R}^d$ if $p = 1$, and $|u|$ is essentially bounded over $\mathbb{R}^d$ when $p = \infty$. We endow these spaces with their standard norms $\|u\|_1 = \int_{\mathbb{R}^d} |u(x)|dx$ and $\|u\|_{\infty} = \sup_{x \in \mathbb{R}^d} |u(x)|$. We let $S(\mathbb{R}^d)$ be the Schwartz space of rapidly-decreasing $C^{\infty}(\mathbb{R}^d)$ functions, and $S'(\mathbb{R}^d)$ be its dual space, i.e., the space of tempered distributions. Then, $S(\mathbb{R}^d) \subset L^p(\mathbb{R}^d)$ and $L^p(\mathbb{R}^d) \subset S'(\mathbb{R}^d)$ continuously. The Fourier transform

<figure><img src="image_3.png" /><figcaption>Fig. 2. Billock and Tsou's phenomena: the presentation of funnel pattern in the centre induces an illusory perception of tunnel pattern in surround after a flickering (image on the left). We have a reverse effect on the right. Taken from Billock and Tsou (2007).</figcaption></figure>

of $u \in L^1(\mathbb{R}^2)$ is defined by

$$ \hat{u}(\xi) := \mathcal{F}\{u\}(\xi) = \int_{\mathbb{R}^d} u(x)e^{-2\pi i \langle x, \xi \rangle} dx, \quad \forall \xi \in \mathbb{R}^d. $$

Since $S(\mathbb{R}^d) \subset L^1(\mathbb{R}^2)$, one can extend the above by duality to $S'(\mathbb{R}^d)$, and in particular to $L^\infty(\mathbb{R}^d)$. Finally the convolution of $u \in L^1(\mathbb{R}^d)$ and $v \in L^p(\mathbb{R}^d)$, $p \in \{1, \infty\}$, is

$$ (u * v)(x) = \int_{\mathbb{R}^d} u(x - y)v(y)dy, \quad \forall x \in \mathbb{R}^d. $$

## 2. NEURONAL FIELDS EQUATIONS

In their seminal paper, Ermentrout and Cowan (1979) develop a theory describing (spontaneous) geometric visual hallucinations perceived in the retina. More precisely, using bifurcation techniques near a static Turing-like instability, they found that a simplified biological model of neuronal fields equation suffices to describe the (spontaneous) formation of cortical patterns (horizontal, vertical and oblique stripes, square, hexagonal and rectangular patterns etc.) in V1. Then, applying the inverse retino-cortical map to these patterns, they obtained some of the geometric visual hallucinations or form constants that Klüver (1966) had meticulously classified. We refer to Fig. 3 for a visual illustration concerning funnel patterns. In their considerations, V1 is treated as a sheet of isotropically interconnected excitatory and inhibitory neurons. A more biologically realistic model of neuronal fields, including the anisotropic properties of cortical connections in V1 (orientation preference of “simple” cells, see Hubel and Wiesel (1959)), was done in Bressloff et al. (2001). The authors were then able to describe all of Klüver's form constants.

Due to the success of the Ermentrout and Cowan model in describing simple patterns, such as funnel patterns, we expect that a similar model (i.e. without orientations preference) should be sufficient to describe sensory hallucinations (visual illusions) induced by these patterns. We, therefore, consider in this work that neuronal activity in V1 evolves according to equation (NF). It models the average membrane potential $a(x, t)$ of a neuron located at $x \in \mathbb{R}^2$ at time $t \ge 0$. In the next section, we present the assumptions on the parameters involved in this equation.

### 2.1 Assumptions on model parameters

Throughout the following we assume the response function $f$ to be an odd non-decreasing function of class $C^2(\mathbb{R})$ such
that $f'(0) = \max_{s \in \mathbb{R}} f'(s) = 1$. Unless explicitly stated otherwise, $f$ is a nonlinear sigmoidal function.

The kernel $\omega$ is taken to be a DoG distribution (difference of Gaussians, also called “Mexican hat” distribution). Namely, we let for all $x \in \mathbb{R}^2$

$$ \omega(x) = [2\pi\sigma^2]^{-1}e^{-\frac{|x|^2}{2\sigma^2}} - [2\pi\kappa^2\sigma^2]^{-1}e^{-\frac{|x|^2}{2\kappa^2\sigma^2}}, \quad (1) $$

where $\kappa > 1$ and $0 < \sigma < 1$. Clearly, $\omega$ is radial and $\omega$ belongs to the Schwartz space $S(\mathbb{R}^2)$. Moreover, its Fourier transform is explicitly given by

$$ \hat{\omega}(\xi) = e^{-2\pi^2\sigma^2|\xi|^2} - e^{-2\pi^2\sigma^2\kappa^2|\xi|^2}, \quad \forall \xi \in \mathbb{R}^2, $$

and $\hat{\omega}$ reaches its maximum at every vector $\xi_c \in \mathbb{R}^2$ such that $|\xi_c| = \sqrt{\log \kappa/\pi^2\sigma^2(\kappa^2 - 1)} =: q_c$.

Observe that with this choice, we fall into the framework of Bressloff et al. (2001), i.e., there exists a critical interaction parameter $\mu_c := \hat{\omega}(\xi_c)^{-1}$ around which spontaneous cortical patterns in V1 emerge.

**Remark 1.** The kernel $\omega$ satisfies the balance<fnref n="2" /> condition $\hat{\omega}(0) = 0$ between excitation and inhibition. Moreover, $\hat{\omega}(\xi) \ge 0$ for all $\xi \in \mathbb{R}^2$ and therefore $\|\hat{\omega}\|_{\infty} = \hat{\omega}(\xi_c)$. Nevertheless, this condition is just for mathematical convenience since it is not explicitly required in our study. Indeed, the following kernel $\omega$ works as well

$$ \omega(x) = [2\pi\sigma_1^2]^{-1}e^{-\frac{|x|^2}{2\sigma_1^2}} - \kappa[2\pi\sigma_2^2]^{-1}e^{-\frac{|x|^2}{2\sigma_2^2}}, \quad x \in \mathbb{R}^2, \quad (2) $$
where $\kappa \ge 1$, $0 < \sigma_1 < \sigma_2$ and $\sigma_1\sqrt{\kappa} < \sigma_2$.

## 2.2 Mathematical preliminaries

We briefly recall some useful results related to equation (NF). It is straightforward to show that the r.h.s. of equation (NF) is a Lipschitz continuous map on $L_t^\infty(\mathbb{R}) \times L_x^\infty(\mathbb{R}^2)$. Thus, it is standard to obtain that, for every external input $I \in L^\infty(\mathbb{R}^2)$ and any initial datum $a_0 \in L^\infty(\mathbb{R}^2)$, equation (NF) admits a unique solution $a \in C([0, +\infty); L^\infty(\mathbb{R}^2))$. Let us recall the following.

**Definition 1.** (Stationary state). Let $a_0$, $I \in L^\infty(\mathbb{R}^2)$. A stationary state $a_I \in L^\infty(\mathbb{R}^2)$ to equation (NF) is a time-invariant solution, viz.

$$ a_I = \mu\omega * f(a_I) + I. \quad (3) $$

Via the contraction mapping principle, one obtains the existence of a unique stationary solution whenever $\mu < \mu_0$, see (Tamekue et al., 2022, Proposition 1). Here, we let

$$ \mu_0 := \|\omega\|_1^{-1} \le \mu_c. $$

This implies in particular that if $\mu < \mu_0$, the map $\Psi: L^\infty(\mathbb{R}^2) \to L^\infty(\mathbb{R}^2)$ associating to each external input $I$ its corresponding stationary state is well-defined and bi-Lipschitz continuous. Observe that $\Psi$ is defined by

$$ \Psi(I) = I + \mu\omega * f(\Psi(I)), \quad \forall I \in L^\infty(\mathbb{R}^2). $$

It is then immediate that $\Psi$ and $\Psi^{-1}$ are $\mathbf{E}(2)$-equivariant, see e.g. (Ermentrout and Cowan, 1979, Appendix A).

**Remark 2.** As a consequence of the $\mathbf{E}(2)$-equivariance of $\Psi$, a subgroup $\Gamma \subset \mathbf{E}(2)$ is a symmetry group of the external input $I \in L^\infty(\mathbb{R}^2)$ if and only if it is a symmetry

<figure><img src="image_4.png" /><figcaption>Fig. 3. Funnel pattern with $\lambda = 1$ in (4): Pattern in the retina (left), corresponding one in V1 (right) after applying the retino-cortical map.</figcaption></figure>

group of the output stationary state $\Psi(I)$. E.g., if $I(x) = I(x_1)$ then $\Psi(I)(x) = \Psi(I)(x_1)$, for all $x = (x_1, x_2) \in \mathbb{R}^2$.

## 2.3 Binary representation of patterns

Due to the retino-cortical map, funnel, and tunnel patterns are respectively given in Cartesian coordinates $x := (x_1, x_2) \in \mathbb{R}^2$ of V1 by

$$ P_F(x) = \cos(2\pi\lambda x_2), \quad P_T(x) = \cos(2\pi\lambda x_1), \quad \lambda > 0. \quad (4) $$

This choice is motivated by analogy with the (spontaneous) geometric hallucinatory patterns described in Ermentrout and Cowan (1979) and Bressloff et al. (2001).

Given the above representation of funnel and tunnel patterns in cortical coordinates, to see how they look in terms of images, we represent them as contrasting white and black regions, see Fig 3. More precisely, define the binary pattern $B_h$ of a function $h: \mathbb{R}^2 \to \mathbb{R}$ by

$$ B_h(x) = \begin{cases} 0, & \text{if } h(x) > 0 \quad (\text{black}) \\ 1, & \text{if } h(x) \le 0 \quad (\text{white}). \end{cases} $$

It follows that $B_h$ is essentially determined by the zero level-set of $h$. Since stimuli involved in the MacKay effect and Billock and Tsou experiments are binary patterns, our strategy in describing these phenomena consists in characterising the zero level-set of output patterns. That is, we are mainly devoted to studying the qualitative properties of patterns by viewing them as binary patterns.

## 3. MACKAY EFFECT

In (Tamekue et al., 2022, Theorem 1), we proved that highly redundant information is needed in the funnel and tunnel patterns for equation (NF) to reproduce the MacKay effect if $\mu < \mu_0/2$. More precisely, if the external input $I = P_F$ or $I = P_T$ in equation (NF), $I$ and $\Psi(I)$ have the same binary pattern, then the same geometric shape in terms of images. Moreover, in the case where the response function $f$ is linear<fnref n="3" />, we obtain that $\Psi(I)$ is even proportional to $I = P_F$ and $I = P_T$. These correspond respectively to $\xi_0 = (0, \lambda)$ and $\xi_0 = (\lambda, 0)$ in the following.

**Proposition 1.** Assume that $I \in L^\infty(\mathbb{R}^2)$ is given by $I(\cdot) = \cos(2\pi\langle\xi_0, \cdot\rangle)$, for some $\xi_0 \in \mathbb{R}^2$ and that the response function $f$ is linear. Then, if $\mu < \mu_0$, it holds

$$ a(\cdot, t) \xrightarrow{t \to \infty} \frac{I(\cdot)}{1 - \mu\hat{\omega}(\xi_0)}, \quad \text{exponentially in } L^\infty(\mathbb{R}^2). $$

<footnote n="2">² For a homogeneous NF equation (i.e., if $I = 0$), this condition enforces the existence of a unique stationary state $a_0 = 0$ even if $f(0) \ne 0$. It was assumed, for instance, in Nicks et al. (2021) for the derivation of the amplitude equation.</footnote>

<footnote n="3">³ Observe that if $f$ is linear, the fact that $f'(0) = 1$ implies that $f(s) = s$. This is non-restrictive since the slope of $f$ can always be factored in the parameter $\mu$.</footnote>
**Proof.** The statement is a direct consequence of the fact that $I(\cdot) = \cos(2\pi\langle\xi_0, \cdot\rangle)$ satisfies $\omega * I = \hat{\omega}(\xi_0)I$, and of the exponential convergence of $a$ to the stationary state.

We deduce that the system's resulting symmetry (the underlying Euclidean symmetry of the interaction kernel) restricts the geometrical shape of visual stimuli that can induce an illusory perception in the after-image. So in describing the MacKay effect with the considered model of neuronal activity (NF) (provided we are far from the bifurcation point, i.e. $\mu < \mu_0 \le \mu_c$), it is necessary to break the Euclidean symmetry of the funnel and tunnel patterns by localising redundant information.

Concerning, e.g. the description of the MacKay effect related to funnel pattern, the “MacKay rays” (see the image on the left of Fig 1), we use the following cortical representation:

$$I(x) = P_F(x) + \varepsilon H(-x_1) = \cos(2\pi\lambda x_2) + \varepsilon H(-x_1). \quad (5)$$

Here $\varepsilon > 0$, and $H$ is the Heaviside step function which would model redundant information in the funnel pattern.

The following result shows that equation (NF) with a linear response function $f$ suffices to describe illusory contours perceived in the after-image induced by the “MacKay rays” having (5) as a V1 analytical representation.

For simplicity, we assume in the rest of this section that parameters in the kernel $\omega$ defined in (1) are such that $2\pi^2\sigma^2 = 1$ and $\kappa^2 = 2$. We also choose $\mu := 1 < 2 = \mu_0$.

**Theorem 2.** Assume the response function $f$ is linear and the input $I$ is given by (5). Then, the unique stationary state to equation (NF) is given for all $(x_1, x_2) \in \mathbb{R}^2$ by

$$a_I(x_1, x_2) = \frac{\cos(2\pi\lambda x_2)}{1 - \hat{\omega}(\xi_0)} + \varepsilon g(x_1), \quad \xi_0 := (0, \lambda). \quad (6)$$

Here $g: \mathbb{R} \to \mathbb{R}$ has a discrete set of zeroes on $(0, +\infty)$.

The hypothesis on $f$ implies that the stationary equation (3) is linear. It follows that the first term in the r.h.s. of (6) is the stationary state associated with the input $P_F$ provided by Proposition 1 and $g$ is the stationary state associated to the second term in the r.h.s. of (5). Therefore, $g$ is the solution of the 1-D stationary equation

$$b(x) = H(-x) + (\omega_1 * b)(x), \quad x \in \mathbb{R}. \quad (7)$$

where the 1-D kernel $\omega_1$ is given for all $x \in \mathbb{R}$ by

$$\omega_1(x) = [\sigma\sqrt{2\pi}]^{-1}e^{-\frac{x^2}{2\sigma^2}} - [\kappa\sigma\sqrt{2\pi}]^{-1}e^{-\frac{x^2}{2\kappa^2\sigma^2}}.$$

Consequently, Theorem 2 follows from the following.

**Proposition 3.** The solution $b \in L^\infty(\mathbb{R})$ of (7) is given, for $x > 0$, by

$$e^{\pi x} \sqrt{\frac{2\pi}{3}} b(x) = \frac{\sqrt{3}}{\pi} \cos \left( \frac{\pi}{3} + \pi x \sqrt{\frac{2\pi}{3}} \right) + \frac{R(x)}{\pi x}, \quad (8)$$

$$\text{where } |R(x)| \le \frac{2}{\pi\sqrt{10\pi}}.$$

Moreover, letting $(\theta_k)_k$ and $(\tau_k)_k$ be respectively zeroes and extrema of $x \mapsto \cos(\pi/3 + \pi x\sqrt{2\pi/3})$ for $x > 0$, the zeroes of $b$ in $(0, +\infty)$ are a countable sequence $(\rho_k)_k$ such that $\rho_k$ is unique in the interval $J_k := ]\tau_k, \tau_{k+1}[$ for all $k \in \mathbb{N}^*$ and

$$|\theta_{k+1} - \rho_k| \le \frac{\sqrt{3}}{\pi\sqrt{2\pi}} \arcsin \left( \frac{2}{\pi(3k-1)\sqrt{5}} \right).$$

The proof of Proposition 3 relies on harmonic and complex analysis techniques. Due to space constraints, we only present a sketch of the proof in the following. We refer to a work in preparation for the complete details.

**Proof of Proposition 3.** Set $I(x) := H(-x)$ in (7). Applying the Fourier transform in $S'(\mathbb{R})$ to both sides of (7), we obtain

$$\hat{b}(\xi) = (1 + \hat{K}(\xi))\hat{I}(\xi), \quad \xi \in \mathbb{R}, \quad (9)$$

where we have set

$$\hat{K}(\xi) := \frac{\hat{\omega}_1(\xi)}{1 - \hat{\omega}_1(\xi)}.$$

Obviously one has $\hat{K} \in S(\mathbb{R})$ so that its inverse Fourier transform $K \in S(\mathbb{R})$ can be computed for all $x \in \mathbb{R}$ by

$$K(x) = \int_{-\infty}^{+\infty} e^{2i\pi\xi x} \hat{K}(\xi) d\xi = \int_{-\infty}^{+\infty} \frac{e^{2i\pi\xi x} \hat{\omega}_1(\xi)}{1 - \hat{\omega}_1(\xi)} d\xi.$$

By a precise analysis on $\hat{K}$, applying the residue Theorem, we can show the existence of a function $S \in L^\infty(\mathbb{R})$ such that for all $x \in \mathbb{R} \setminus \{0\}$, it holds

$$\frac{e^{\pi|x|\sqrt{\frac{2\pi}{3}}}}{2\sqrt{\pi}} K(x) = \cos \left( \frac{\pi}{12} + \pi|x|\sqrt{\frac{2\pi}{3}} \right) + \frac{S(x)}{|x|}, \quad (10)$$

$$|S(x)| \le \frac{2}{\pi\sqrt{6\pi}}, \quad \forall x \in \mathbb{R} \setminus \{0\}.$$

Taking now the inverse Fourier transform of equation (9) in the space $S'(\mathbb{R})$, we find

$$b(x) = I(x) + (K * I)(x), \quad x \in \mathbb{R} \setminus \{0\}. \quad (11)$$

Letting $I(x) = H(-x)$, we obtain (8). Standard arguments based on the intermediate value Theorem, allow us to prove the second part of the statement.

**Remark 3.** Theorem 2 implies that if the external input is the V1 representation of the “MacKay rays” defined by (5), then the associated stationary state corresponds to the V1 representation of the after-image reported by MacKay (1957). Moreover, we have the exponential convergence of $a(\cdot, t)$ on the stationary state when $t \to \infty$. It follows that equation (NF) theoretically describes the MacKay effect associated with the “MacKay rays” at the cortical level. Due to the retino-cortical map, we deduce the theoretical description of the MacKay effect for the “MacKay rays” in the retina.

We refer to Section 5 for numerical results.

**Remark 4.** One can prove that no illusory contours are present in the after-image using a Gaussian kernel to model synaptic interactions (i.e., discarding the excita-tory/inhibitory nature of interactions). Indeed, the corre-sponding kernel $K$ in (10) is positive on $\mathbb{R} \setminus \{0\}$ and thus letting $I(\cdot) = H(-\cdot)$ in (11) shows that $b$ is non-negative.

## 4. BILLOCK AND TSOU'S EXPERIMENTS

In Tamekue et al. (2022), we exhibited numerical results showing the capability of equation (NF) to reproduce Billock and Tsou's experiments. The stimuli used in these experiments are funnel or tunnel patterns localised at the fovea or periphery. Due to the retino-cortical map, this corresponds to taking as external inputs in equation (NF), $I = \varepsilon P_F v$ or $I = \varepsilon P_T v$, where $\varepsilon > 0$ and $v$ is a localised function either in the left or in the right area of the cortex.
In particular, the external inputs in these experiments do not fill all of the visual field. The function $v$ can then be thought of as a localised control aiming to break the global plane Euclidean symmetry of stimuli patterns.

In this section, we prove that the equation (NF) with a linear response function $f$ cannot describe these phenomena: In contrast to the MacKay effect, the phenomena reported by Billock and Tsou are completely nonlinear. We mention that the numerical experiments of Section 5 will also highlight that the shape of the nonlinearity is crucial and, in particular, they suggest that an asymmetric nonlinearity is essential for the description.

We will focus on the funnel pattern localised at the fovea. In V1, it corresponds to the following external input.

$$I(x) = \cos(2\pi\lambda x_2)H(-x_1), \quad \lambda > 0, \quad x \in \mathbb{R}^2, \quad (12)$$

where $H$ is the Heaviside step function. For ease of notation, we assume the kernel $\omega$ in (1) is such that $\kappa^2 = 2$. If $\delta$ is the Dirac distribution at zero and $h \in S'(\mathbb{R}^2)$,

$$h_*^0 := \delta, \quad h_*^j := h * h * \dots * h, \quad j \in \mathbb{N}^*.$$

By Newton binomial formula, for all $h, g \in S'(\mathbb{R}^2)$, $n \in \mathbb{N}^*$

$$(h + g)^n_* = \sum_{j=0}^{n} \binom{n}{j} h_*^{n-j} * g_*^j. \quad (13)$$

Since the convolution of Gaussians with zero mean remains a Gaussian with zero mean, the following is a direct consequence of (13).

**Lemma 4.** Let $I \in S'(\mathbb{R}^2)$. Assume that the response function $f$ is linear. If $\mu < \mu_0$, the stationary state $a_I \in S'(\mathbb{R}^2)$ of (NF) is given by

$$a_I = I + \sum_{n=1}^{\infty} \mu^n \sum_{j=0}^{n} \binom{n}{j} (-1)^j g_{n,j} * I.$$

Here $g_{n,j}$ is the Gaussian defined for $x \in \mathbb{R}^2$, $n \in \mathbb{N}^*$ by

$$g_{n,j}(x) = \frac{1}{2\pi(n+j)\sigma^2} e^{-\frac{|x|^2}{2(n+j)\sigma^2}}, \quad j \in [|0, n|].$$

**Proposition 5.** Assume that the response function $f$ is linear and $I \in L^\infty(\mathbb{R}^2)$ is defined by (12). If $\mu < \mu_0$, there exists a function $R \in L^\infty(\mathbb{R}^2)$ such that the stationary state of (NF) is given for all $x \in \mathbb{R}^2$, by

$$a_I(x) = [H(-x_1) + R(x_1)] \cos(2\pi\lambda x_2).$$

**Proof.** By Lemma 4, one computes for all $x \in \mathbb{R}^2$,

$$(g_{n,j} * I)(x) = \frac{\cos(2\pi\lambda x_2)}{\sigma\sqrt{2\pi(n+j)}} e^{-\frac{\lambda^2(n+j)\sigma^2}{2}} \int_{-\infty}^{0} e^{-\frac{(x_1-y_1)^2}{2(n+j)\sigma^2}} dy_1 \\ = \frac{e^{-\frac{\lambda^2(n+j)\sigma^2}{2}}}{2} \operatorname{erfc}\left(\frac{x_1}{\sigma\sqrt{2(n+j)}}\right) \cos(2\pi\lambda x_2),$$

where erfc is the complementary error function.

**Remark 5.** Proposition 5 shows that the output $a_I$ associated with (12) has a contribution in the right area of the cortex given by

$$a_{I,r}(x) = R(x_1) \cos(2\pi\lambda x_2), \quad x_1 > 0.$$

Since $a_{I,r}$ depends on the factor $\cos(2\pi\lambda x_2)$, the retinocortical map tells us that a visual stimulus consisting of fan shapes in the centre induces an after-image containing fan shapes in the periphery instead of concentric rings only,

as Billock and Tsou reported. Equation (NF) (provided we are far from the bifurcation point, i.e. $\mu < \mu_0 \le \mu_c$) with a linear response function cannot describe Billock and Tsou's experiments. Therefore, these phenomena depend fundamentally on the presence of the nonlinearity $f$.

## 5. NUMERICAL RESULTS

The numerical implementation is performed with Julia and is available at https://github.com/dprn/MacKay-Billock_Tsou-2022/. Given an input $I$, the stationary state $a_I$ is numerically implemented via an iterative fixed-point method. The cortical data is defined on a square $(x_1, x_2) \in [-L, L]^2$, $L = 10$ with steps $\Delta x_1 = \Delta x_2 = 0.01$.

For the reproduction of the MacKay effect, parameters in the kernel $\omega$ given by (1) are $2\pi^2\sigma^2 = 1$ and $\kappa^2 = 2$. We exhibit in Fig. 4 the MacKay effect for “MacKay rays” $I(x) = \cos(5\pi x_2) + \varepsilon H(2 - x_1)$, $\varepsilon = 0.025$. Here, we use a linear response function ($f(s) = s$). We stress that the phenomenon can be reproduced with any odd sigmoidal function, see e.g. (Tamekue et al., 2022, Fig. 3).

In (Tamekue et al., 2022, Figs. 5 and 6), we illustrated the capability of equation (NF) to reproduce Billock and Tsou experiments with the nonlinear response function $f(s) = (1+\exp(-s+0.25))^{-1}-(1+\exp(0.25))^{-1}$. In Section 4, we proved that a linear response function does not reproduce these phenomena. We exhibit in Figs. 5 and 6 Billock and Tsou's experiments for a funnel-like stimulus localised at the periphery and fovea, respectively. As images, we have a fan shape pattern at the periphery (resp. fovea) and white in the fovea (resp. periphery). We use the kernel $\omega$ defined in (2) with $\sigma_1 = 0.1$, $\sigma_2 = 0.5$ and $\kappa = 4.56$. In Fig. 5, the stimulus is $I(x) = \cos(4\pi x_2)H(x_1 - 6)$, and the nonlinearity is $f(s) = \max(-0.2, \min(1, 1.7s))$. In Fig. 6, the stimulus is $I(x) = \cos(4\pi x_2)H(6 - x_1)$ and the nonlinearity is $f(s) = \max(-0.2, \min(1, 1.2s))$. In Fig. 7, the stimulus is $I(x) = \cos(4\pi x_2)H(6 - x_1)$: On the left, the nonlinearity used is $f(s) = \max(-0.2, \min(1, 1.2s))$, and on the right, we use $f(s) = \max(-1.2, \min(1, s))$. In the after-image on the left of Fig. 7, the fan shape does not extend to the periphery, whereas on the right, the fan shape extends through the periphery.

To see how the shape of the nonlinearity is involved in the reproducibility, we consider the family of nonlinear functions $f_{m\alpha}(s) = \max(-m, \min(1, \alpha s))$. We exhibit in Fig. 8 the range of parameters $(m, \alpha) \in \{(k/10, \ell/10) | k = 0, \dots, 20, \ell = 1, \dots, 20\}$ where equation (NF) with the nonlinear response function $f_{m\alpha}$ reproduces Billock and Tsou's experiments or not, for a funnel-like stimulus $I(x) = \cos(4\pi x_2)H(6 - x_1)$. The magenta (resp. black) region corresponds to the value of $(m, \alpha)$ where $f_{m\alpha}$ reproduces (resp. do not reproduce) the phenomenon and the yellow region corresponds to the value of $(m, \alpha)$ where $f_{m\alpha}$ reproduces the phenomenon, but the stimulus extends through the periphery. According to Billock and Tsou's experimental description, the good range of parameters $(m, \alpha)$ is that of the magenta region.

## 6. CONCLUSIONS

We showed that breaking the Euclidean symmetry of stimuli is necessary for theoretically describing the MacKay
<figure><img src="image_5.png" /><figcaption>Fig. 4. Initial stimulus (“MacKay rays”, *left*) inducing the MacKay effect (*right*).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Billock and Tsou's experiments: funnel stimulus localised at the periphery (*left*) and after-image (*right*).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Billock and Tsou's experiments: funnel stimulus localised at the fovea (*left*) and after-image (*right*).</figcaption></figure>

effect and Billock and Tsou's experiments using equation (NF) whenever the parameter $\mu$ is smaller than $\mu_0$. To do this, we used a localised function. The issue of characterising all possible localized perturbations (i.e., “controls”) of the funnel and tunnel inputs that give rise to MacKay-type observations is still an open problem.

Although the Gaussian kernel is usually used in image processing and computer vision tasks due to its proximity to the visual system, it is unable to reproduce the phenomena described here. A physiological reason for this is that we used a one-layer model of the NF equation. It is not then biologically realistic to model synaptic interactions with a Gaussian, which would model only excitatory-type interactions between neurons.

Numerical results indicate that the anisotropic nature of cortical connections of “simple” cells in V1 need not be integrated into equation (NF) to reproduce Billock and Tsou’s experiments. Moreover, the model reproduces the phenomena without a temporal flicker of the complementary region where the stimulus is not localized.

Further work is ongoing to theoretically recover the numerical results in Fig. 8 for general nonlinear response functions. In particular, the numerical results suggest that to reproduce Billock and Tsou’s phenomena, an asymmetric (i.e., not odd nor even) nonlinearity is required.

<figure><img src="image_8.png" /><figcaption>Fig. 7. After-images in Billock and Tsou's experiments: funnel stimulus localised at the fovea as in Fig. 6 (right). On the *right* (resp. *left*) stimulus extends (resp. does not) to the periphery.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Range of ($m, \alpha$) where $f_{m\alpha}$ reproduce Billock and Tsou's experiments or not, see text for details.</figcaption></figure>

## REFERENCES

* Amari, S. (1977). Dynamics of pattern formation in lateral-inhibition type neural fields. *Biol. Cybern.*, 27(2), 77-87.

* Billock, V.A. and Tsou, B.H. (2007). Neural interactions between flicker-induced self-organized visual hallucinations and physical stimuli. *PNAS*, 104(20), 8490-8495.

* Bressloff, P.C., Cowan, J.D., Golubitsky, M., Thomas, P.J., and Wiener, M.C. (2001). Geometric visual hallucinations, euclidean symmetry and the functional architecture of striate cortex. *Philos. Trans. R. Soc. Lond., B, Biol. Sci.*, 356(1407), 299-330.

* Ermentrout, G.B. and Cowan, J.D. (1979). A mathematical theory of visual hallucination patterns. *Biol. Cybern.*, 34(3), 137-150.

* Helmholtz, H.L.F. (1867). *Optic physiologique*. Masson.

* Hubel, D.H. and Wiesel, T.N. (1959). Receptive fields of single neurones in the cat's striate cortex. *J. Physiol.*, 148(3), 574.

* Klüver, H. (1966). Mescal and mechanisms of hallucinations. *Chicago: University of Chicago*.

* MacKay, D.M. (1957). Moving visual images produced by regular stationary patterns. *Nature*, 180, 849-850.

* Nicks, R., Cocks, A., Avitabile, D., Johnston, A., and Coombes, S. (2021). Understanding sensory induced hallucinations: From neural fields to amplitude equations. *SIAM J. Appl. Dyn. Syst.*, 20(4), 1683-1714.

* Tamekue, C., Prandi, D., and Chitour, Y. (2022). Reproducing sensory induced hallucinations via neural fields. In *2022 IEEE-ICIP*, 3326-3330.

* Zeki, S., Watson, J.D., and Frackowiak, R.S. (1993). Going beyond the information given: the relation of illusory visual motion to brain activity. *Proc. Royal Soc. B: Biol. Sci.*, 252(1335), 215-222.