James Rankin

Emilien Tlapale

Romain Veltz

Olivier Faugeras

Pierre Kornprobst

Bifurcation analysis applied to a model of motion integration with a multistable stimulus

Keywords: motion, perception, multistability, visual cortex, barber pole, bifurcation

A computational study into the motion perception dynamics of a multistable psychophysics stimulus is presented. A diagonally drifting grating viewed through a square aperture is perceived as moving in the actual grating direction or in line with the aperture edges (horizontally or vertically). The different percepts are the product of interplay between ambiguous contour cues and specific terminator cues. We present a dynamical model of motion integration that performs direction selection for such a stimulus and link the different percepts to coexisting steady-states of the underlying equations. We apply the powerful tools of bifurcation analysis and numerical continuation to study changes to the model's solution structure under the variation of parameters. Indeed, we apply these tools in a systematic way, taking into account biological and mathematical constraints, in order to fix model parameters. A region of parameter space is identified for which the model reproduces the qualitative behaviour observed in experiments. The temporal dynamics of motion integration are studied within this region; specifically, the effect of varying the stimulus gain is studied, which allows for qualitative predictions to be made.

Analyse de bifurcation appliquée à un modèle d'intégration de mouvement en présence d'un stimulus multistable

Résumé : Nous présentons une analyse d'un modèle de la dynamique de la perception du mouvement, en présence d'un stimulus multistable. Un réseau en mouvement dans la diagonale présenté derrière une ouverture carrée est perçu comme se déplaçant dans la direction du réseau ou dans la direction des cotés de l'ouverture (horizontalement ou verticalement). Ces différents percepts sont le produit de l'interaction entre les indices ambigus des contours du réseau et les indices de terminaisons de lignes. Nous présentons un modèle dynamique de l'intégration du mouvement qui effectue la sélection de la direction pour un tel stimulus, en reliant les différents percepts à des états d'équilibre coexistants des équations sous-jacentes. Nous appliquons les puissants outils d'analyse de bifurcations et de continuation pour étudier les changements de la structure des solutions du modèle en fonction des paramètres. Nous appliquons ces outils de façon systématique afin de fixer les paramètres du modèle, en prenant en compte les contraintes biologiques et mathématiques. Une région de l'espace des paramètres est identifiée, région dans laquelle le modèle reproduit le comportement qualitatif observé dans les expériences. La dynamique temporelle de l'intégration du mouvement est étudiée au sein de cette région; plus précisément, nous étudions l'effet produit par la modification du gain du stimulus ce qui permet d'effectuer des prédictions qualitatives.

Mots-clés : mouvement, perception, multistabilité, cortex visuel, barber pole, bifurcation

Introduction

The interesting and long-studied aperture problem constitutes an important ambiguity that must be resolved by the visual system in order to attribute an accurate direction of motion to moving objects [START_REF] Wallach | Über visuell wahrgenommene Bewegungsrichtung[END_REF][START_REF] Wuerger | on the visually perceved direction of motion" by hans wallach: 60 years later[END_REF]. The motion of a uniform contour is consistent with many possible directions in the absence of terminator information provided by line endings.

The ambiguous contour information is referred to as a 1D cue and the specific terminator information, which can be intrinsic to the object or produced by an occluding aperture, is referred to as a 2D cue.

In this article, we focus on a classical psychophysics stimulus used to probe the interactions between 1D and 2D motion cues, the so-called "barber pole" illusion (Hildreth, 1983, Chapter 4) . A a diagonally drifting grating viewed through an elongated rectangular aperture is perceived as drifting in the direction of the long edge of the aperture. The illusion is generated by the larger number of unambiguous (but misleading) 2D cues parallel to the long edges of the aperture. Mitigating the effect of the 2D cues can break the illusion as demonstrated in experiments by introducing a depth separation between grating and aperture [START_REF] Shimojo | Occlusion and the solution to the aperture problem for motion[END_REF], or, by introducing indentations on the aperture edges [START_REF] Kooi | Local direction of edge motion causes and abolishes the barberpole illusion[END_REF]. More interesting is the fact that for certain stimulus parameters, the barber pole can elicit a multistable perceptual response. The effect of the ratio of the lengths of the aperture edges, or terminator ratio, was investigated in (Castet et al, 1999;[START_REF] Fisher | The directional tuning of the barber-pole illusion[END_REF]. For large ratios the direction of the long aperture edge dominates perception. However, for ratios close to 1 : 1 the stimulus is multistable with a tri-modal response for the perceived direction; the main percepts either agree with the direction perpendicular to the grating's orientation or with one of the two aperture edge directions . Another interesting aspect in resolving the aperture problem is the temporal dynamics of the integration of 1D and 2D cues. The barber pole illusion was studied in ocular following experiments [START_REF] Masson | Temporal dynamics of motion integration for the initiation of tracking eye movements at ultra-short latencies[END_REF]; early tracking responses were shown to be initiated in the contour direction with a later response in the direction of the long axis of the aperture. The existence of an early response dominated by 1D cues, which is later refined when 2D cues are processed is supported by further studies in ocular following [START_REF] Barthélemy | Temporal dynamics of 2d motion integration for ocular following in macaque monkeys[END_REF], psychophysics [START_REF] Lorenceau | Different motion sensitive units are involved in recovering the direction of moving lines[END_REF] and physiology [START_REF] Pack | Temporal dynamics of a neural solution to the aperture problem in visual area MT of macaque brain[END_REF].

The middle temporal area (MT) of the visual pathway, that receives the majority of its synaptic inputs from the primary visual cortex (V1), plays a key role in the perception of moving objects and, more specifically, the solution of the aperture problem. MT is characterised by direction-selective neurons that are organised in a columnar fashion, similar to the organisation of orientationselective neurons in V1 [START_REF] Diogo | Electrophysiological imaging of functional architecture in the cortical middle temporal visual area of cebus apella monkey[END_REF]. For an extensive discussion of the function of MT, see the following review articles [START_REF] Britten | The middle temporal area: Motion processing and the link to perception[END_REF][START_REF] Born | Structure and function of visual area MT[END_REF]. Cortical responses of MT have been linked specifically to perception of motion; see again the review article [START_REF] Britten | The middle temporal area: Motion processing and the link to perception[END_REF] and, more recently, the paper [START_REF] Serences | The representation of behavioral choice for motion in human visual cortex[END_REF].

Several models of motion integration have been proposed in the literature to solve the aperture problem, providing some insights into the underlying neural mechanisms. Building on the first linear/non-linear models [START_REF] Chey | Neural dynamics of motion processing and speed discrimination[END_REF][START_REF] Simoncelli | A model of neuronal responses in visual area MT[END_REF], several approaches added extensions to modulate Inria the motion integration stages: feedback between hierarchical layers [START_REF] Grossberg | Neural dynamics of motion integration and segmentation within and across apertures[END_REF][START_REF] Bayerl | Disambiguating visual motion through contextual feedback modulation[END_REF], inclusion of input form cues [START_REF] Ben-Yishai | Laminar cortical dynamics of visual form and motion interactions during coherent object motion perception[END_REF][START_REF] Bayerl | Disambiguating visual motion by form-motion interaction -a computational model[END_REF], luminance diffusion gating (Tlapale et al, 2010b), or depth cues [START_REF] Beck | Interactions of motion and form in visual cortexa neural model[END_REF]. Although these models reproduce the predominant percepts in a wide range of stimuli, in none of the articles describing them are multistable results depicted. Furthermore, although they show some limited individual case of dynamical behaviour just on the level of simulations, there is no rigorous analysis of the dynamical behaviour and no comprehensive parameter studies that fully explore all possible dynamical behaviour. In summary, the questions of what mechanisms are behind multistable motion perception and what dynamical processes are involved have largely been overlooked from the modelling point of view.

The focus of this paper will be the analysis of a mathematical model of motion integration with an input that generates a multistable perceptual response.

To do so, we study cortical behaviour at the population level; see [START_REF] Pouget | Information processing with population codes[END_REF] for a discussion of how information can be encoded at a population level. We work within the neural fields formalism, a mathematical framework that was originally studied in [START_REF] Wilson | Excitatory and inhibitory interactions in localized populations of model neurons[END_REF][START_REF] Amari | Characteristics of randomly connected threshold element networks and neural systems[END_REF][START_REF] Amari | Characteristics of random nets of analog neuron-like elements[END_REF]; see (Ermentrout and Terman, 2010, Chapter 11) for various derivations of the equations. Instead of looking at the spiking behaviour of individual, interconnected neurons, the neural field approximation represents the mean firing rate of a neural population at the continuum limit and activity levels are represented in a spatially continuous way. Neural fields equations have been successfully applied to the study of motion in, e.g., [START_REF] Giese | Dynamic Neural Field Theory for Motion Perception[END_REF][START_REF] Deco | The role of multi-area interactions for the computation of apparent motion[END_REF] and (Tlapale et al, 2010a). In the later, the complex model presented, describing behaviour of multiple cortical layers and their feedforward and feedback connections, was capable of performing motion integration on both natural image sequences and classical psychophysics presentations. In terms of multistable stimuli, the 1 : 1 barber pole discussed above lead to coexisting steady states in the model, but the temporal dynamics of multistable perception were not investigated. We aim to develop a tractable model of manageable complexity that allows for a detailed study of the temporal dynamics of multistable motion perception using powerful tools from dynamical systems theory.

A natural tool for the study of dynamical systems for which multiple steadystate solutions co-exist is bifurcation analysis. Throughout the manuscript when the term solution is used this refers to a steady-state solution. In dynamical systems theory, a bifurcation is a critical point encountered under the variation of one or more parameters at which there is a change in the stability and number of solutions. Indeed, under the variation of a parameter a solution of a dynamical system will vary in state space; when the solution is plotted in terms of some norm against the parameter, it lies on a solution branch. A dynamical system may have multiple solution branches that are characterised by, for example, different spatial and stability properties. At special points when varying a parameter, solution branches can meet and bifurcate from each other. At these so-called bifurcation points, the number of solutions and their stability changes. In order to know what type of behaviour a model can produce it is necessary to gain full understanding of the type of bifurcations that occur, the types of solutions that are involved and at which parameter values. For a general introduction to bifurcation analysis of finite dimensional systems see [START_REF] Strogatz | Nonlinear dynamics and chaos[END_REF][START_REF] Kuznetsov | Elements of Applied Bifurcation Theory[END_REF], and of infinite dimensional systems with symmetries see [START_REF] Chossat | Methods in Equivariant Bifurcations and Dynamical Systems[END_REF][START_REF] Haragus | Local bifurcations, center manifolds, and normal forms in infinite dimensional systems[END_REF]. Bifurcation analysis has been used to study pattern formation in a number of different settings [START_REF] Ermentrout | Large scale spatially organized activity in neural nets[END_REF][START_REF] Bressloff | Nonlocal ginzburg-landau equation for cortical pattern formation[END_REF][START_REF] Coombes | Bumps, breathers, and waves in a neural network with spike frequency adaptation[END_REF]. More specifically, a spatialised model of V1 has been used to investigate hallucinatory visual patterns [START_REF] Bressloff | Geometric visual hallucinations, euclidean symmetry and the functional architecture of striate cortex[END_REF][START_REF] Golubitsky | Bifurcation on the visual cortex with weakly anisotropic lateral coupling[END_REF][START_REF] Bressloff | Nonlocal ginzburg-landau equation for cortical pattern formation[END_REF], localised patterns have been studied in models of working memory [START_REF] Laing | Multiple bumps in a neuronal model of working memory[END_REF][START_REF] Guo | Existence and stability of standing pulses in neural networks: I. existence[END_REF][START_REF] Faye | Localized states in an unbounded neural field equation with smooth firing rate function: a multi-parameter analysis[END_REF] and in a model of texture perception [START_REF] Faye | Analysis of a hyperbolic geometric model for visual texture perception[END_REF]. In all of these studies only spontaneous activity is studied, that is, in the absence of any cortical input.

In this article, we propose a spatialised ring model of direction selection, where the connectivity in the direction space and physical space is closely related to the Mexican-hat type connectivity typically used in the ring model of V1 [START_REF] Ben-Yishai | Laminar cortical dynamics of visual form and motion interactions during coherent object motion perception[END_REF][START_REF] Somers | An emergent model of orientation selectivity in cat visual cortical simple cells[END_REF][START_REF] Hansel | Modeling feature selectivity in local cortical circuits[END_REF]. We apply analytical tools such as stability analysis and normal form computations in order to identify and categorise bifurcations in our model. These tools have been used successfully to study the neural field equations, see, for example [START_REF] Curtu | Pattern formation in a network of excitatory and inhibitory cells with adaptation[END_REF][START_REF] Coombes | Modeling electrocortical activity through local approximations of integral neural field equations[END_REF][START_REF] Roxina | How effective delays shape oscillatory dynamics in neuronal networks[END_REF]. However, in the presence of an input to the model, these analytical tools are no longer applicable; although certain perturbation-type methods can be applied if the input is considered to have a specific, simple spatial structure and to be small [START_REF] Veltz | Local/global analysis of the stationary solutions of some neural field equations[END_REF]Ermentrout et al, 2010;[START_REF] Kilpatrick | Response of traveling waves to transient inputs in neural fields[END_REF]. Note that in the study [START_REF] Veltz | Local/global analysis of the stationary solutions of some neural field equations[END_REF], large inputs with a simple spatial structure were also studied with numerical continuation. In this paper, we first investigate the model's behaviour in the absence of an input using analytical techniques. Then building on the knowledge gained we apply the tool of numerical continuation to track solution branches under parameter variation and detect bifurcation points; effectively continuation provides a computational tool for performing bifurcation analysis. For an introduction to continuation algorithms see [START_REF] Krauskopf | Numerical Continuation Methods for Dynamical Systems[END_REF]. For the problem we study here the continuation module LOCA of the numerical tools package Trilinos is well suited [START_REF] Heroux | An overview of the Trilinos project[END_REF]. Numerical bifurcation tools have not previously been used to study a neural fields model in the presence of a large, spatially complex input. Furthermore, the application of bifurcation analysis and numerical continuation to the study of a model of motion integration is new. These methods allow us to build a complete picture of the model's possible dynamics in terms of parameter regions exhibiting qualitatively different behaviour and to identify the boundaries between these regions. In this way, we are able to ensure that parameter regions in which a desired behaviour is present are not isolated and that the behaviour is robust with respect to small changes in the model set up. The typical approach of a simple numerical search in parameter space often misses important parameter regions and does not provide information about robustness of behaviour with respect to parameter variation. The technical parts of this paper form the basis of the strategy employed to fix model parameters based on biological and mathematical constraints. We perform a systematic study of the model's solution structure depending on key parameters including the stimulus strength and identify a parameter region for which the model has steady-state solutions corresponding to the three predominant percepts (diagonal, horizontal, vertical) observed in experiment [START_REF] Fisher | The directional tuning of the barber-pole illusion[END_REF]) . An extended study of Inria the system's temporal dynamics allows for experimental predictions to be made regarding the distributions of the different percepts seen for different length presentations. In particular we predict that for long enough presentations only the horizontal and vertical percepts should be seen.

Model of direction selection

The model described here uses a neural fields description of the firing-rate activity of a population of neurons in middle temporal area (MT) over a physical (cortical) space and a feature space of motion direction. Two essential mechanisms are represented by the model: direction selectivity in the feature space and spatial diffusion of activity across the physical space. Stimulus input to the model is represented as preprocessed motion direction signals from V1 complex cells. Such a representation is comparable to classic motion detectors such as the output of elaborated Reichardt detectors [START_REF] Van Santen | Temporal covariance model of human motion perception[END_REF][START_REF] Bayerl | Disambiguating visual motion through contextual feedback modulation[END_REF] or of motion filters. The output of the model is the time evolution of activity levels across the physical-direction space.

The functionality of the model is encoded by the connectivity across directionspace and physical-space, which is processed by a nonlinearity. In the direction space, the connectivity is based on a narrowly tuned excitation with broadly tuned inhibition, as described for MT in [START_REF] Grunewald | Orthogonal motion after-effect illusion predicted by a model of cortical motion processing[END_REF]. Such a Mexican-hat-type connectivity in the orientation-space has been used previously for the ring model of orientation selection [START_REF] Ben-Yishai | Laminar cortical dynamics of visual form and motion interactions during coherent object motion perception[END_REF][START_REF] Bressloff | Geometric visual hallucinations, euclidean symmetry and the functional architecture of striate cortex[END_REF]. Here, we assume that the inhibition tuning in MT is broad enough to be approximated by uniform lateral inhibition; with uniform lateral inhibition when cells selective for a particular direction are active, cells selective for the opposite direction are inhibited. This is consistent with known direction selectivity properties of MT [START_REF] Albright | Direction and orientation selectivity of neurons in visual area MT of the macaque[END_REF][START_REF] Diogo | Electrophysiological imaging of functional architecture in the cortical middle temporal visual area of cebus apella monkey[END_REF] and with recordings from MT with transparent drifting dot patterns moving in opposite directions [START_REF] Snowden | The response of area MT and V1 neurons to transparent motion[END_REF]. Furthermore, the uniform lateral inhibition connectivity has the convenient property of fixing the first non-trivial Fourier mode of the connectivity to be the largest, which is necessary for the model to produce tuning-curve solutions. In the physical space, diffusion is captured by an inverted Mexican hat connectivity, which has been used in a number of neural fields models with delays; see, for example, [START_REF] Hildreth | Pattern formation in intracortical neuronal fields[END_REF][START_REF] Venkov | Dynamic instabilities in scalar neural field equations with space-dependent delays[END_REF] and (Veltz, 2011, Chapter 6). As motivated in [START_REF] Venkov | Dynamic instabilities in scalar neural field equations with space-dependent delays[END_REF], for cortical tissue the principal pyramidal cells are often surrounded by inhibitory interneurons and their long range connections are typically excitatory [START_REF] Gilbert | Spatial integration and cortical dynamics[END_REF][START_REF] Salin | Corticocortical connections in the visual system: structure and function[END_REF][START_REF] Mcguire | Targets of horizontal connections in macaque primary visual cortex[END_REF]. The inverted Mexican hat connectivity propagates activity outwards from stimulated regions and is consistent with a model output that describes a coherent motion across the physical space.

Representation of the stimulus

The stimulus that we consider here is a single drifting grating viewed through a square aperture, which has been shown to exhibit multistability of the viewer's percept during the first 2 seconds of presentation (Castet et al, 1999;[START_REF] Fisher | The directional tuning of the barber-pole illusion[END_REF]. For longer presentations, the stimulus is known to exhibit perceptual switching of the kind studied in, for example, binocular rivalry experiments [START_REF] Blake | A primer on binocular rivalry, including current controversies[END_REF]. However, the focus in this article will be the early dynamics after onset of the stimulus. The stimulus is shown in Fig. 1(a) and the three dominant percepts are labelled as diagonal D (the direction perpendicular to the grating's orientation at 0 • ) , horizontal H (in line with the horizontal aperture edge at -45 • ) or vertical V (in line with the vertical aperture edge at 45 • ). Note that here in the description of the stimulus and later in the numerical results of Sec. 4 we give angles in degrees; however, for ease of presentation in the technical parts that follow in this section and in Sec. 3 angles are given in radians. One of the most significant simplifications in the model is that we consider a 1D approximation of the cortex for physical space, which restricts the way in which the stimulus can be represented as input to the model, as discussed below. We consider the different cues at different points on a 1D cut across the stimulus, as shown in Fig. 1(b). On the interior of the aperture the input is directionally ambiguous because 1D motion cues centred around the direction of the drifting grating are received. At the aperture edges the input is specific because 2D motion cues parallel to the aperture edges are received. It is the interplay between these two types of competing motion signals that produces the multistable percept. In this paper, we use the terms stable and unstable in the sense of dynamical systems theory to reflect whether a solution is attracting or repelling in state-space; the term multistable refers to the coexistence of different possible percepts for a visual stimulus.

The inclusion of 2D physical space would allow for many variations to the stimulus to be considered. However, there are a number modifications that are possible with the present 1D representation that we now discuss. For example, rectangular apertures are often considered with a specific ratio between the different edges [START_REF] Fisher | The directional tuning of the barber-pole illusion[END_REF] and by taking a 1D cut across the stimulus this information would be lost. However, a similar effect could be achieved by giving a stronger weighting to the longer aperture edge. Another feature often investigated in psychophysics experiments is the angle of the grating with respect to the aperture edges (Castet et al, 1999). Again, although certain information is lost in the 1D cut, such a change to the stimulus could be considered by placing the 2D elements as shown in Fig. 1(c) asymmetrically about v = 0 (but still orthogonal to each other). Furthermore, the choice of the cut made Fig. 1(b) is quite arbitrary. For example, moving the cut closer to the line between the top left and bottom right corners would widen the stimulated region in x; in terms of the models response, we would expect the further separation of the 2D cue elements to delay the time it takes for the 2D cues to affect on the dynamics. The specific case where a cut is taken between two corners of the aperture would result in there being competition between 2D different cues at the aperture edge points. The stimulus shown in Fig. 1(c) could be modified to incorporate this by adding points at both v = ±45 • at each aperture edge x = ±0.75 as opposed to v = 45 • at x = -0.75 and v = -45 • at x = 0.75 as is currently the case. Overall the 1D approximation used here is well suited to the specific stimulus studied with its particular symmetry properties, that is, with a square aperture and a grating that forms an equal angle with each aperture edge.

Inria

D ( v = 0 • ) H (v = -45 • ) V (v = 45 • ) (a) (b) (c) x ❅ ❅ | x = 0 v x

Neural field equation

We now introduce the neural field equation that describes the neuronal population's activity, detailing how the different elements discussed so far enter into the model. The firing rate activity of a single population of neurons in MT is denoted

p : (x, v, t) ∈ Ω × V × R + → p(x, v, t) ∈ R + , ( 1
)
where Ω is the spatial domain (a bounded subset of R) and V = R/2πZ is the direction space. Note that we impose the restriction that p ≥ 0, because it is not physically relevant to have a negative firing rate. We assume that x and v are independent variables: every direction v is represented at each physical position x. The neural field equation is given by

   ∂p ∂t = -µp + S(λ [Jp + T + kI ext ]), p(x, v, 0) = p 0 (x, v), (2)
where µ is the decay rate, J is the connectivity operator, k is the input gain (analogous to contrast) for the stimulus I ext (Fig. 1(c)) and T is the constant threshold parameter. The parameter λ is the stiffness (slope at x = 0) of the sigmoid nonlinearity S(x) = 1 1+exp(-x) . The choice of a sigmoid function, which is smooth and infinitely differentiable, facilitates the study of steady-state behaviour and allows for the application of numerical continuation [START_REF] Veltz | Local/global analysis of the stationary solutions of some neural field equations[END_REF]. The initial condition at t = 0 is p 0 (x, v). The connectivity J has the following form

Jp = ν 1 (G E p) -ν 2 (G I p) -ν 3 p.
(3)

The first two terms represent a difference of 2D Gaussian functions G E and G I , where ν 1 and ν 2 represent the relative strength of excitation and inhibition, respectively. The convolutions G E and G I are defined by:

G E p = Ω V g E x (x -y)g E v (v -w)p(y, w, t)dwdy, (4)
G I p = Ω V g I x (x -y)b I v (v -w)p(y, w, t)dwdy, (5)
where

g * x (x) = 1 √ 2πσ * x exp(-|x| 2 2σ * x ) for σ * x = σ I x , σ E
x and similarly for g E v (v). In physical space, we have an inverted Mexican hat type connectivity so σ E

x > σ I x . The width of excitation in direction space is set by σ E v and the uniform lateral inhibition is represented by the box-function b I v (v) = 1 2π . In order to illustrate the shape of the connectivity function across (x, v)-space we first plot in Fig. 2

(a) the function M x (x) = ν 1 g E x -ν 2 g I x with σ E x > σ I x and ν 2 > ν 1 along with its Fourier transform Mx (j). Similarly in panel (b) we plot M v (v) = ν 1 g E v -ν 2 b I v
and its Fourier transform Mv (k). Furthermore, the full connectivity

M (x, v) = ν 1 g E x (x)g E v (v) -ν 2 g I x (x)b I v (v)
is plotted in panel (c); functions with the shape of M x and M v are recovered by taking 1D slices across this surface. The parameter ν 3 is ignored here as this serves only to shift the surface up or down. Restrictions on the choice of parameters describing the connectivity function are discussed An additional linear threshold term ν 3 p is used to tune the maximum firing rate of direction-selected solutions, whereas the constant threshold T is used to tune the firing rate of homogeneous, non-direction-selected solutions. It is convenient to include the linear threshold term in the connectivity operator as this simplifies the stability computations in Sec. 3.2. The choice of all parameters is discussed in Sec. 4.

Inria in Sec. 3.2. . . -1 0 1 -2 -1.5 -1 -0.5 0 0.5 -2 0 2 -1.5 -1 -0.5 0 0.5 1 x, j v, k M x Mx M v Mv ✻ ✏ ✏ ✶ ❅ ❅ | M x v
3 Analytic results in the absence of stimulus, k = 0

In this section we discuss inherent properties of Eq. ( 2) and its solutions without a stimulus input. The results we obtain analytically provide a foundation of knowledge about the different types of solution the model can produce, the role of key parameters and a means to set appropriate values of parameters based on mathematical and biological constraints. We begin our study by looking at the symmetry properties satisfied by the connectivity and the governing equation in Sec. 3.1. We show that Eq. ( 2) with the connectivity J as described above is equivariant with respect to a certain symmetry group. This important property dictates the types of solution that can be produced by the model. Furthermore, it determines the type of bifurcations that occur.

In the single population model that we consider here, the only types of solutions that we encounter are steady states (or, persistent states). Given an initial condition p 0 , the time evolution of the equations can be computed numerically; the particular initial condition chosen will determine which steadystate solution the system converges to. It is important to note that the transient dynamics encountered before the system converges to a steady state can also be greatly affected by the initial condition. In Sec. 3.2 we calculate analytically the steady states that have the additional property of being independent of both the physical and direction space. For these spatially independent solutions a constant level of activity persists across (x, v)-space; this type of solution can be thought of as the baseline activity that we would see in the absence of a stimulus (or below the contrast threshold). We give an expression that allows us to compute these solutions depending on the system parameters. Also in Sec. 3.2, we compute the eigenvalues and eigenvectors of the connectivity operator J (a spatial-mode decomposition), which allows us in turn to compute the stability of the steady-state solutions dependent on the nonlinearity stiffness λ. We show that for small enough λ the steady-state solutions are stable. As λ is increased, the most destabilising mode of J, as determined by its largest eigenvalue will become unstable at a critical value of λ. This critical value is the system's principal bifurcation point.

We determine the type of the principal bifurcation in Sec. 3.3. Furthermore, given the mode of J that loses stability in this bifurcation, and given the symmetry properties of the governing equation, we are able to characterise the spatially dependent solutions produced by the model. A normal form computation determines the way in which the transition from spatially homogeneous solutions to spatially dependent solutions occurs in the model.

Symmetry group

Here we discuss the symmetry properties of Eq. ( 2), which will play an important role in determining the type of bifurcation that the model produces. The general concept is to specify the group of translations and reflections for which the governing equation is equivariant. The same group of translations and reflections, when applied to a solution of the equations, will produce coexisting solutions; for example, we will see in Sec. 3.3 that translational invariance in v means that a direction-selected solution associated with one specific direction can be translated by any angle to give direction-selected solutions associated with all other possible directions. Note that when a stimulus is introduced, the symmetry group of the equations will be in some way reduced and it is, therefore, important to first identify the full symmetry group before its introduction.

In order to simplify subsequent calculations we impose periodicity on the spatial domain so that Ω = R/cZ for some c ∈ R; this simplification does not affect the stability properties of the system [START_REF] Faugeras | Abolute stability and complete synchronization in a class of neural fields models[END_REF]. Let us consider the group, denoted Γ, of translations of R/cZ: a one parameter group parametrised by α ∈ R/cZ. An element Γ α of this group acts in the following way on the variables (x, v, t):

Γ α • (x, v, t) = (x + α, v, t). ( 6
)

Inria

In general the action of the group Γ on the function p is

Γ • p(x, v, t) def = p Γ -1 (x, v, t) , (7)
and more specifically for a translation in x:

Γ α • p(x, v, t) = p(x -α, v, t). (8)
Let

F (p(x, v, t)) = ∂p(x, v, t) ∂t + µp(x, v, t) (9) -S(λ[Jp(x, v, t) + T ]);
it can be shown that F is equivariant with respect to the group Γ or, equivalently

Γ α F (p) = F (Γ α p).
Furthermore, the function F is equivariant with respect to the reflection group generated by R which has the following action on the variables and activity:

R • (x, v, t) = (-x, v, t), (10)
R • p(x, v, t) = p(-x, v, t). (11)
If we denote H x the group generated by Γ and R and since we have

         Γ α1 Γ α2 = Γ α1+α2 ∀α 1 , α 2 , ∈ R/cZ, RΓ α = Γ -α R ∀α ∈ R/cZ, Γ 0 = Id, R 2 = Id, (12)
the group H x is isomorphic to O(2), the group of two-dimensional orthogonal transformations. Furthermore, the equation F is equivariant with respect to the similarly defined group H v generated by translation and reflection in v. Therefore, F is equivariant under the action of the symmetry group

H = H x × H v which is isomorphic to O(2) × O(2).

Spatially homogeneous solutions and their stability

Steady-state solutions are those for which ∂p ∂t = 0. We first consider solutions that are independent in both the physical space x and velocity space v and the level of activity across the population p is equal to a constant value p ∈ R + . To find these solutions we set the right hand side of Eq. ( 2) equal to 0 and, thus, search for solutions p to the following equation

µp = S (λ [J p + T ]) .
Given that for the normalised Gaussian functions we have G E p = p and G I p = p, we can further write

µp = S (λ [(ν 1 -ν 2 -ν 3 )p + T ]) , (13)
an implicit expression for p.

Next, we wish to determine the linear stability of Eq. ( 2) at the solution p depending on system parameters. The first task is to compute the spectrum of the operator J. In order to do this we decompose J into Fourier modes by obtaining its eigenvalues and the associated eigenvectors; this information will determine exactly which (Fourier) modes of J have the greatest destabilising effect. The eigenvalues ζ (j,k) of J are given by the following relation:

ζ (j,k) = ν 1 ĝE x (j) ĝE v (k) -ν 2 ĝI x (j) bI v (k) -ν 3 , (14)
where ĝE x (j), ĝE v (k) and ĝI x (j) are Fourier coefficients of the respective periodically extended Gaussian functions. The coefficients bI v (k) are defined as follows:

bI v (k) = 1 for k = 0, 0 otherwise. ( 15
)
Due to the functions g E x , g I x , g E v and b I v being even, their Fourier coefficients are real positive and even. Hence we have

ζ (±j, ±k) = ζ (j,k) .
The dimension of the eigenspace E (j,k) associated with each eigenvalue ζ (j,k) depends on the indices j and k. Here we let the indices (j, k) be positive numbers. The eigenvectors χ (j,k) are given by:

χ (j,k) =                {1} j = 0, k = 0, {e ikv , e -ikv } j = 0, k > 0, {e 2πijx c , e -2πijx c } j > 0, k = 0, {e i( 2πjx c +kv) , e i( 2πjx c -kv) , e i(kv-2πjx c ) , e -i( 2πjx c +kv) } j > 0, k > 0. ( 16
)
The eigenvectors could equivalently be represented as combinations of sin and cos functions.

Using the modal decomposition of J, we can obtain an expression for the eigenvalues associated with the solution p, for each Fourier mode of J. The sign of the eigenvalue for each mode will tell us whether it is stable (-) or unstable (+). We define S 1 to be the linear coefficient in the Taylor expansion of S at the fixed point p; note that we Taylor expand about λ[(ν 1ν 2ν 3 )p + T ] and S 1 depends on the values of several other system parameters. By linearising about the solution p of Eq. ( 2), we obtain the following expression for the eigenvalues

(j,k) = -µ + λS 1 ζ (j,k) . ( 17
)
By identifying the mode of J with the largest eigenvalue ζ (j,k) , we can find the smallest value λ for which p is unstable. Indeed for small enough λ the solutions are stable as (j,k) ≈ 0. For the values of σ E x , σ I x and σ E v used in this paper (see Table 1) and imposing certain restrictions on the values of ν 1 , ν 2 and ν 3 , we can identify exactly which modes (j, k) are the most destabilising. If we impose ν 1 > 0, 0 ≤ ν 3 < ν 1 and ν 2 > ν 1 then the following properties hold:

• The mode (0, 0) is stable because ζ (0,0) = ν 1 -ν 2 -ν 3 < 0.
• The ν 2 ĝI x (j) bI v (k) term ensures that all modes for which k = 0 are stable.

Inria

• The positive ν 1 ĝE x (j) ĝE v (k) term produces the destabilising contribution, which is greatest for j = 0.

• Further, this destabilising contribution is greatest for k closest to 0 and then diminishes for increasing k.

Therefore, the largest eigenvalue ζ (j,k) of J corresponds to the mode (0, 1) followed by the subsequent modes with increasing k. Accordingly, in the analysis that follows it is convenient to drop the subscript j and to assume that it is zero, such that

ζ k = ζ (0,k)
. The largest eigenvalue is ζ 1 ; therefore, the smallest value of λ for which (j,k) = 0 is given by

λ c = µ S 1 ζ 1 . ( 18
)
This value λ c is the system's principal bifurcation point, which we study in the next section. The term S 1 depends on λ and p, but values of λ c can be found by solving the following system for the pair (p c , λ c ):

pc = S(λc[(ν1-ν2-ν3) pc+T ]) µ , λ c = µ S1(λc[(ν1-ν2-ν3) pc+T ])ζ1 . (19)
By taking advantage of the equality S = S(1 -S), it was proved in [START_REF] Veltz | Local/global analysis of the stationary solutions of some neural field equations[END_REF]) that given ζ 0 < 0 and ζ 1 > 0, the pair (p c , λ c ) is unique. These two inequalities for the eigenvalues hold given the restrictions on ν 1 , ν 2 and ν 3 discussed above.

Bifurcation points associated with other modes that occur as λ is increased beyond λ c can be found in a similar fashion, however, it is the branch of solutions that are born from the principal bifurcation that will determine the types of spatially dependent solutions that the model will produce.

Normal form of the principal bifurcation point

In this section we classify the principal bifurcation point by first, applying the center manifold theorem and secondly, giving the appropriate change of variables to reduce the system's dynamics into a normal form; we now introduce these concepts. In the previous section we performed a modal decomposition of the connectivity and computed the linear stability of Eq. ( 2) with respect to perturbations in the different modal components. For the different modal components, or eigenvectors χ (j,k) (Eq. ( 16)), the sign of the associated eigenvalue (j,k) (Eq. ( 17)) gives the linear stability. In the case when (j,k) = 0 the stability is neutral; the parameter value for which this occurs is the bifurcation point. At this bifurcation point it is necessary to also consider nonlinear terms in some parameter neighbourhood in order to capture the local dynamics. A center manifold reduction allows us to compute these nonlinear terms by means of a leading order Taylor approximation; the center manifold theorem allows us to prove rigorously that the computed reduced system accurately captures the local dynamics. A normal form computation is a change of variables that classifies the type of bifurcation present in our system and allows for the dynamics local to the bifurcation point to be seen clearly. The coefficients found in the normal form computation provide important information about the direction of bifurcating branches in terms of the bifurcation parameter and the stability of these branches.

We prove in Appendix A that the relevant hypotheses for the centre manifold hold in our case. This computation depends both on the symmetry properties discussed in Sec. 3.1 and the fixed point stability analysis from Sec. 3.2. Indeed, in the previous section we identified the system's principal bifurcation point as given by the pair (p c , λ c ), solutions to the system (19). We now define respectively the first, second and third order coefficients in the Taylor expansion of S at pc to be S 1 , S 2 and S 3 . We drop the subscript notation for the eigenvectors χ = e iv and χ = e -iv , which span the two-dimensional eigenspace E 1 associated with the eigenvalue ζ 1 . The eigenvalues ζ 0 (for the homogeneous mode) and ζ 2 (for the j = 0, k = 2 mode) will also appear in the analysis that follows.

Here we define a centre manifold on the two dimensional eigenspace of ζ 1 . This centre manifold will be independent of physical space x, therefore, the manifold must be equivariant with respect to the reduced symmetry group H v , which is isomorphic to O(2). Here we introduce P the real valued solutions on the center manifold, which it is convenient to express in terms of a complex variable w. We decompose P into linear components on the eigenspace E 1 and nonlinear components orthogonal to E 1 . We set

P = pc + w • χ + w • χ + Ψ(w, w, λ -λ c ), w(t) ∈ C, (20)
where Ψ is a grouping of nonlinear terms called the center manifold correction.

A simple change of coordinates can be used to eliminate the constant term pc .

From [START_REF] Haragus | Local bifurcations, center manifolds, and normal forms in infinite dimensional systems[END_REF], Chapter 2), we have in this case a pitchfork bifurcation with O(2) symmetry and the reduced equation for the dynamics of w has the following normal form equation

dw dt = aw(λ -λ c ) + bw|w| 2 + O(|w|(|λ -λ c | 2 + |w| 4 )), (21)
together with the complex conjugated equation for w. Indeed, due to the fact that the reduced equation for the dynamics of w must also be O(2)-equivariant, even powered terms in the reduced equation are prohibited. We define the parameter dependent linear operator L λ = -µId + λS 1 J, which represents the linearisation of the right hand side of Eq. ( 2) at pc . The linear coefficient a can be determined by considering the action of L λ on the eigenvector χ (equivalently, on the linear terms of P ) at the bifurcation point given by Eq. ( 18):

L λ • χ = (-µ + λS 1 J) • χ, = (-µ + λS 1 ζ 1 ) • χ, = (-µ + λ µ λ c ) • χ, = µ λ c (λ -λ c ) • χ,
which, by comparing with Eq. ( 21), gives a = µ λc > 0. The fact that a is positive means that the principal solution branch is stable before the bifurcation point. Note that the linear term disappears at the bifurcation point λ = λ c and it is necessary to consider higher order terms in order to quantify the dynamics. The sign of the cubic coefficient b determines the direction of the bifurcating branch; Inria see Fig. 3 for the two possibilities when a > 0. We now use the expression for solutions P on the center manifold to determine the coefficient b in terms of our model parameters. Substituting the expression (20) for P into Eq. ( 2) we obtain

dP dt = (-µ + λ c S 1 J) L λc (P -pc ) + λ 2 c S 2 2 (J(P -pc )) 2 + λ 3 c S 3 6 (J(P -pc )) 3 + O(P 4 ), (22)
where

J(P -pc ) = ζ 1 w • χ + ζ 1 w • χ + JΨ
and the linear terms are collected in the operator L λc . Higher order terms in the Taylor expansion may be neglected given the form of Eq. ( 21) assuming that a and b are non-zero. It remains to determine the coefficient and b by matching terms between Eqs. ( 21), ( 22).

In order to compute b we Taylor expand Ψ at λ = λ c :

Ψ(w, w, 0) = ∞ p,q=2 ψ pq w p w q ,
where the coefficients ψ pq are orthogonal to the eigenvectors χ and χ. We identify terms with common powers in the Taylor expansion in order to obtain the following equation to be solved for b:

L λc Ψ 21 + λ 2 c S 2 2 (2JΨ 11 ζ 1 χ + 2JΨ 20 ζ 1 χ) + λ 3 c S 3 6 3ζ 3 1 χ 2 χ = bχ.
After some calculations given in Appendix B, we obtain the following expression for b:

b = λ 3 c ζ 3 1 S 3 2 + λ c S 2 2 µ ζ 0 (1 -ζ0 ζ1 ) + ζ 2 2(1 -ζ2 ζ1 ) . ( 23
)
As we can see from Eq. ( 23), the criticality of the pitchfork bifurcation, as determined by the sign of b, depends on all system parameters in a complex way. We briefly discuss the implications of this criticality in our model.

• b < 0: The bifurcation is supercritical [START_REF] Haragus | Local bifurcations, center manifolds, and normal forms in infinite dimensional systems[END_REF]. The new branch of solutions exist after the bifurcation, i.e. for λ > λ c . Furthermore the branch of solutions will be stable (attracting) local to the bifurcation. In our model, for increasing λ passing the bifurcation point, there would be smooth transition (smooth change in activity levels) from spatially homogeneous solutions to a direction-selected solution.

• b > 0 The bifurcation is subcritical [START_REF] Haragus | Local bifurcations, center manifolds, and normal forms in infinite dimensional systems[END_REF]. The new branch of solutions exist before the bifurcation, i.e. for λ < λ c . Furthermore the branch of solutions will be unstable (repelling) local to the bifurcation. In our model, for increasing λ passing the bifurcation point, there would be non-continuous transition (jump in activity levels) from spatially homogeneous solutions to a direction-selected solution (see below for explanation).

Sketches of the bifurcation diagrams for the two cases are shown in Fig. 3. In

RR n°7822

the second case, and from the analysis in [START_REF] Veltz | Local/global analysis of the stationary solutions of some neural field equations[END_REF], we know that the unstable branch existing for λ < λ c must also have a fold bifurcation λ f at some point 0 < λ < λ c . Thus the unstable branch of solutions existing for λ < λ c will change direction and stability at λ f . For λ f < λ < λ c there are coexisting stable solutions and for λ > λ c the direction-selected solution is the only stable one. Therefore, passing the bifurcation at λ c results in a jump from spatially homogeneous solutions to a direction-selected solution.

. Ultimately we need to choose a value of λ such that 1. In the absence of a stimulus there are no direction selected solutions;

. b < 0 b > 0 λ λ (a) (b) λ c λ c λ f |.| |.| ✻ ✻
2. When a stimulus is introduced a solution is selected that is intrinsically present in the model.

For λ too large the model will produce direction-selected solutions in the absence of a stimulus. For λ too small the solutions will be purely driven by the stimulus once it is introduced, the connectivity that dictates the solutions intrinsically present in the model will not play a role. Therefore, for the case b < 0 we must choose a value close to but still less than λ c . For the case b > 0 we must choose a value close to but still less than λ f . The bifurcating branch of solutions is characterised by the mode involved ((j, k) = (0, 1)), that is, the solutions will be uniform in physical space and will have a single maximum in v-space; we can think of the maximal point as being centred at a selected direction. Hence, we refer to solutions on the bifurcated branch as direction-selected solutions. Secondly, due to the O(2) symmetry in the absence of a stimulus, taking a direction-selected solution, we know that it will still be a solution under any angular translation in v; i.e. there continuum (or ring) of solutions representing all possible selected directions.

Numerical results

In order to progress in the study of our model, we make use of computational tools. In Sec. 4.1 we take advantage of the analytical results from Sec. 3 in order to explore properties of the model dependent on parameters and to determine relevant ranges of the parameters. In Sec. 4.2 we build on the normal form computation from Sec. 3.3 by computing the bifurcated branch of solutions using numerical continuation. We compute the relationship between the activity levels along the bifurcated branch and certain parameters in order to fix their values. An in depth study of the solution structure of the model varying two parameters is given in Sec. 4.3. We identify a region of interest in parameters space, which we study in more detail in Sec. 4.4. In particular we investigate the temporal dynamics of the model which provides insight into the relation between the model's solution structure and different stimulus driven responses.

The default parameter values that we use in Secs. 4.3 and 4.4 are given in Table 1. We can arbitrarily set the decay parameter to µ = 2 and the diffusion parameter to ν 1 = 3. Equally, these could be set to equal 1, however, the chosen values lead to all system parameters being of roughly the same order of magnitude, which facilitates the computations with Trilinos.

The activity levels for steady states of Eq. ( 2) are bounded by [0, 1 µ ]. In order to simplify the presentation of the results it is convenient to give activity levels in terms of either the max-norm p max or L 2 -norm |p| as a percentage of the maximum value 1 µ . For the spatial connectivity, the surround width is set by σ 1

x which is three times that of the centre width region σ 2 x [START_REF] Tadin | Perceptual consequences of centresurround antagonism in visual motion processing[END_REF]. The spatial extent of the stimulus is roughly 1.5 times larger than the surround width. The width of the excitation in direction space is set by σ 1

v to be about 20 • . This relatively tight tuning allows for the model to produce direction selected solutions (tuning curves) of widths appropriate for distinguishing between the main percepts H, V and D; tuning widths are discussed further in 4.2. Finally, it remains to set values for T , ν 2 and ν 3 , which we do in Secs. 4.1 and 4.2.

Parameter tuning

Ultimately, we will operate the model at a value of λ close to the either λ c or λ f , depending on the coefficient b, as discussed in Sec. 3.3. We wish to determine an appropriate firing rate for the spontaneous activity. We do this by investigating the value of pc , which will be approximately the same as the activity before the bifurcation. In the experimental study [START_REF] Sclar | Coding of image contrast in central visual pathways of the macaque monkey[END_REF] the contrast response curves of neurons at various stages of visual processing, including MT, were found. The data were fitted with the standard Naka-Rushton function characterising low responsiveness at low contrast, followed by a region of high sensitivity where the firing increases rapidly, which eventually plateaus out at some maximal firing rate. From this data it appears that at very low contrast (in the absence of a stimulus), the firing rate is at approximately 5%-20% of the maximum firing rate. It is therefore appropriate to set model parameters such that 5% < pc < 20%. Figure 4(a) and Fig. 4(b) show the dependence of pc and b on the threshold parameter T , used to set the spontaneous activity level, and ν 2 , the inhibition strength; all other parameters are set to the values in Table 1. Note that for each value of pc there is an associated unique value of λ c given by Eq. ( 19); it is possible vary T and ν 2 linearly such that the argument of S in the first equation is fixed, which explains the affine relationship in Fig. 4(a). The same property carries through to the Taylor coefficients S 2 and S 3 and to the values of b given by Eq. ( 23).

For a given value of ν 2 , Fig. 4(a) provides the necessary range of T for which we have 5% < pc < 20%. Furthermore, for values of (ν 2 , T ) in this range, we always have the case b > 0 corresponding to a subcritical pitchfork bifurcation; this implies that it will be necessary to operate the model close to λ f . In order to set the final values of T , ν 2 and ν 3 , it is necessary to consider the level of activity and tuning width of the direction selected branch of solutions. In the next section we use the information presented here, along with the dependence of these solution properties on ν 2 and ν 3 to determine the final values of the parameters used.

Solution structure in the absence of stimulus

We now look at numerically computed solutions in the absence of stimulus (k = 0). In Sec. 3.3 it was shown that the system's principal bifurcation is a pitchfork with O(2)-symmetry; furthermore, it was shown in the previous section that this pitchfork bifurcation is subcritical for the parameter values specified in Table 1 . We now employ numerical continuation to compute directly the bifurcated solution branch under the variation of λ. We use the continuation package LOCA, part of the Trilinos package of numerical tools. We use a discretisation of 37 points for both the physical space x and direction space v, verifying that the integrals involved in the solution of the equations satisfy suitable error bounds; such a discretisation gives rise to a system of 1369 ordinary differential equations (ODEs), which resulted in manageable computations with Trilinos. Increasing the discretisation further results in a drastic increase in computation times.

Figure 5(a) shows a bifurcation diagram in the parameter λ in the absence of a stimulus. As predicted from the analytical results herein and [START_REF] Veltz | Local/global analysis of the stationary solutions of some neural field equations[END_REF], the pitchfork bifurcation P 1 from the homogeneous (lower) branch of solutions is subcritical and the bifurcated (upper) branch undergoes a fold bifurcation at F 1 . Here we have λ c = 22.1 and λ f = 15.4. Therefore, it will be necessary to operate the model with λ < 15.4 such that direction selected solutions do not exist in the absence of a stimulus. Due to the O(2)symmetry properties of the solutions on the bifurcated branch, the direction selected solutions are invariant to translations in v (there is a continuum, or ring, of solutions corresponding to all possible directions). Panels (b) and (c) show a direction-selected solution from the upper branch centred at an arbitrary direction; the same solution exists for any translation in v. These directionselected solutions are modulated in v, but are still uniform in physical space x; it is only later, with the introduction of an input, that the model produces responses that are also modulated in x (see Fig. 8).

In order to set values of the parameters ν 2 and ν 3 we investigate how certain properties of the direction-selected solutions change with respect to these parameters. Firstly, the solutions should have a suitable tuning width in v-space and, secondly, the activity should not be close to saturation. Given that we will operate the model at a value of λ close to λ f , it is convenient to study these properties at the fold point F 1 , as shown in Fig. 5(a), and how they change under variation of the ν 2 and ν 3 . Note that as we vary either ν 2 or ν 3 , the value of λ for which the fold bifurcation occurs also changes. Therefore, we perform two-parameter continuation in λ and k, varying one parameter to satisfy the condition that the system be at a steady state and a second parameter to satisfy the condition that the system also be at a fold bifurcation. We define w 1 2 to be the tuning width in v of the direction selected solutions at half-height. Figure 6(a) shows the variation of w 1 2 at the fold point F 1 with respect to ν 2 , for convenience of presentation we do not show the variation of λ. As ν 2 is increased (greater lateral inhibition) the tuning width decreases. Introducing the linear threshold term (ν 3 = 1.5) does not change this characteristic, but shifts the curves up. It is necessary to set a value of ν 2 such that the tuning width is less than 45 • in order to make the distinction between solutions moving diagonally or vertically/horizontally as shown in Fig. 1. At ν 2 = 66, which gives w 1 2 ∈ (30 • , 40 • ), the maximum level of activity for the direction-selected solution is shown in Fig. 6(b). Imposing the condition that p max < 50% at F 1 , the system should still exhibit suitable sensitivity to changes in the input gain when the stimulus is introduced. We arrive at a choice of (ν 2 , ν 3 ) = (66, 1.5), which gives w 1 2 ≈ 40 • and p max ≈ 48%. Overall, in the last two sections, we have explored the relationship between the three parameters T , ν 2 and ν 3 and properties of the homogeneous and direction-selected solutions. This has allowed us to set suitable values of these parameters such that the solutions produced by the model satisfy important biological and mathematical constraints; refer back to Table 1 for the parameter values used in the subsequent sections.

Introduction of stimulus and two-parameter analysis

From the results presented in the previous section we know that the sigmoidal slope λ should be set at a value for which the model cannot spontaneously produce direction-selected solutions. This gives rise to the requirement that λ < 15.4 so as to be at a λ-value less than the first fold point labelled F 1 in Fig. 5. We now investigate the solutions when a stimulus is introduced to the model by increasing the stimulus gain parameter k from 0; this is analogous to increasing the stimulus contrast. Recall that the stimulus described in Sec. 1 has the form shown in Fig. 1(c) when represented in the (x, v)-plane.

Figure 7 shows one-parameter bifurcations diagrams in k initialised at the values of λ indicated in the panels. In order to see the solution structure clearly the diagrams are shown with k on a logarithmic scale. In panels (a), (b) and (c) solutions are plotted in terms of the the norm |p|; note that in panel (c) this allows the main solution branch at a lower activity level to be distinguished from the solution branch associated with the percept D. In panel (d) we plot solutions in terms of the average direction v in order to distinguish between the solutions associated the percepts H and V .

F H,V 1 P 1 P 2 F H,V 1 F 2 F 3 F 4 F 5 P 1 P 2 F V 1 F H 1 P 1 P 2 |p| |p| |p| v k k k k λ = 12.5 λ = 13 λ = 14 λ = 14 (a) (b) (c) (d)

Inria

At λ = 12.5, as shown in panel (a), there is a single branch of stable solutions on which no bifurcations are encountered. The solutions effectively mimic the input stimulus along the solution branch and the level of activity increases with k.

At λ = 13, for increasing k bifurcations are encountered at the fold points labelled F H 1 and F V 1 , the subcritical pitchfork P 1 and the supercritical pitchfork P 2 . With the introduction of the stimulus, the O(2)-symmetry of the pitchfork for k = 0, as discussed in the previous section, is broken. We now find that the pitchfork bifurcation at P 1 gives rise to a pair of direction-selected-solution branches, as opposed to a continuum of solutions at all possible directions when k = 0 as discussed in Sec. 4.2 . The two branches are associated with the directions H and V and the branches undergo a fold bifurcations F H 1 and F V 1 , respectively. These two fold points coincide when plotted in terms of the solution norm |p|. The pair of direction-selected branches reconnect at the supercritical (for k decreasing) pitchfork P 2 . The main solution branch (that mimics the input stimulus and has a lower level of activity than the bifurcated branch) is unstable between P 1 and P 2 ; it coexists with the direction-selected solutions.

At λ = 14, as shown in Fig. 7(c), the same bifurcations are encountered as for λ = 13, but with F H 1 , F V 1 and P 1 at lower values of k, and P 2 at a larger value of k. The significant difference is the introduction of two pairs of fold bifurcations F 2 , F 3 and F 4 , F 5 on the main unstable branch of solutions. This series of bifurcations gives rise to an unstable direction-selected solution between F 2 and F 5 . For values of k between F 2 and F 5 the inset shows the largest eigenvalue on the unstable direction-selected branch (grey) and the two symmetrical direction-selected branches (black). Close to F 2 , the largest eigenvalue is positive but very close to zero implying that the solution is weakly unstable in this region. Figure 7(d) shows the pair of symmetrical direction-selected branches with average direction v as the solution measure. The unstable solution corresponding to D (v = 0 • ) is shown in Fig. 8(a) and the stable solutions corresponding to V (v ≈ 45 • ) and H (v ≈ -45 • ) are shown in Fig. 8(b) and Fig. 8(c), respectively. We note that a slightly elevated level activity persists outside of the aperture (stimulated region), which is a consequence of the underlying solutions encoded by the connectivity as shown in Fig. 5 that have a single direction selected across the entire physical space. It is suggested that this (subthresold) activity could facilitate further recruitment of direction selection beyond the aperture edges. See, for example, [START_REF] Shimojo | Occlusion and the solution to the aperture problem for motion[END_REF] where experiments were carried out using gratings masked by multiple apertures.

The pitchfork and fold bifurcations encountered in the one-parameter bifurcation diagrams presented thus far are codimension-one bifurcations that lie on curves in the parameter plane and, can therefore, be tracked under variation of two parameters. Indeed, the necessary routines to achieve this are implemented in the package Trilinos. Figure 9 shows the locus of pitchfork bifurcation P and the loci of fold bifurcation F H , F V , F D and F T over the (λ, k)-plane; note that F H and F V coincide and that F D is associated with the unstable D solution. In order to illustrate how solutions are organised we consider three one-parameter slices of the parameter plane. Indicated by vertical dashed lines in Fig. 9, the slices correspond to the three one-parameter cases shown in Fig. 7. Firstly, for the trivial case at λ = 12.5 the corresponding slice in Fig. 9 does not intersect the pitchfork or fold curves. At λ = 13, for increasing k the first (simultaneous) intersection is with F H,V (corresponding to F H 1 and F V 1 ), followed by P RR n°7822

.

. (corresponding to P 1 ), and P again (corresponding to P 2 ); compare Fig. 9 and Fig. 7(b). Finally, at λ = 14, for increasing k the intersections occur in the following order F H,V , F D , P , F T , F T , F D and P ; these intersections correspond to the bifurcation points F H,V

1

, P 1 , F 2 , F 3 , F 4 , F 5 and P 2 in Fig. 7(c). We are now able to summarise the type of solutions that exist in different regions of the (λ, k)-plane. Outside of the region bounded by P and below F H,V , there are no direction selected solutions. One can think of F H,V representing a contrast threshold in k that changes with respect to λ. For λ too small, the solutions are purely stimulus driven. In the region between F H,V and P the homogeneous solution is stable and coexists with the two direction selected solutions. Inside the area bounded by P , the system's homogeneous solutions are unstable and there are always two direction selected solutions corresponding to H and V. For points in the region bounded by P and also to the right of F D , there exists stable solutions associated with H and V, and an unstable solution associated with D. We identify this as a region of interest for which we expect the model to qualitatively reproduce the behaviour observed in experiment.

Temporal dynamics

By means of an in-depth, two-parameter bifurcation analysis we have identified a region of interest in parameter space for which the model produces two sta-Inria . . 12.5 13 13.5 14 0.01

0.1 1 k λ F H,V P F D F T F 1 P 1 P 2 F 1 F 3 F 4 F 5 P 2 P 1 F 2 ❇ ❇ ❇ ▼ ❅ ❅ ❅ ■ Figure 9: A two-parameter bifurcation diagram is shown in the (λ, k)-plane
where the locus of pitchfork bifurcation P is a black curve and the loci of fold bifurcations F H,V , F D and F T are a grey curves. The point where the coinciding fold curves F H,V terminate at an intersection with the pitchfork curve P is indicated by black point. One-parameter slices indicated by vertical dashed lines correspond to the bifurcation diagrams for fixed λ and varying k shown in Fig. 7; the labels along these slices correspond with the bifurcation points in Fig. 7. ble direction-selected solutions corresponding to the horizontal (H) and vertical (V) percepts, and one unstable direction-selected solution corresponding to the diagonal (D) percept. We know that in this parameter region the system will converge to one of the stable percepts; in this section we study how the unstable solution D plays a role in the temporal dynamics. At the chosen parameter values in Table 1 with λ = 14, the model produces a low level of homogeneous activity in the absence of stimulus k = 0. We take this homogeneous state as an initial condition for simulations with a random perturbation drawn from a standard uniform distribution in order to introduce a stochastic element. The equations remain deterministic, but we introduce variability in the initial conditions for each simulation. First, we present an example simulation and define some quantities that characterise the temporal dynamics; next, we study how the temporal dynamics change with respect to the strength of the input gain k.

The reader may find it helpful to refer back to the bifurcation diagrams Fig. 7(c) and (d), which show the solution branches corresponding to H, D and V.

Figure 10 shows an example of the temporal dynamics produced by the model in the region of parameter space with solutions corresponding to H, D and V, (see Sec. 4.3). The behaviour can be broken into three phases: 1) initially no specific direction is selected but there is a slightly higher level of activity for stimulated directions, 2) the first direction-selected solution to appear corresponds to D but the system diverges from this unstable solution, 3) . . the final direction-selected solution corresponds to either H or V and the system will remain at this stable solution after convergence. The point in time of the transition from phase 1 to phase 2 is denoted t 1 , phase 2 to phase 3 is denoted t 2 (see panels (b), (c) and accompanying caption for definition of t 1 and t 2 . The time-mark t 1 coincides with the first appearance of a direction-selected solution, in this case D. The time-mark t 2 coincides with the convergence of the solution to either H or V. These two quantities characterise the temporal dynamics of the system as transitions are made between the solutions associated with the different percepts.

Inria

We now investigate the way in which the characteristic quantities t 1 and t 2 vary depending on initial conditions; as such we study their distribution in time based on many simulations initiated with random initial conditions. Furthermore, we study the dependence of these distributions on the parameter k. For a range of discrete values of k ∈ {0.2, 0.7}, we run at each k-value N = 500 simulations. Figure 11 shows a summary of the statistics of the distributions of t 1 and t 2 depending on k. The results show that as k increases, the mean value of t 1 decreases (D is perceived sooner) as does the mean value of t 2 -t 1 (the switch from D to either H or V occurs sooner). Another important property is that for smaller values of k, there is a large spread for t 2 with a long tail extending off to large values of t. This can be explained from the bifurcation results presented earlier, specifically, the unstable branch of solutions corresponding to D and the pair of stable solution branches corresponding to H and V as shown in Fig. 7(c),(d). The inset of Fig. 7(c) shows the largest eigenvalues along the solution branches for the parameter values considered here. For the unstable solution corresponding to D the largest eigenvalue is indeed positive, but still very close to 0. Therefore, very long transients can be observed resulting in a large time being spent at the D solution. As k increases the eigenvalue increases and the transient times decrease resulting in less time being spent near D. Finally, for large values of k the D solution is barely seen at all, there is convergence directly to H or V.

We now look more closely at how the unstable solution D affects the dynamics at k = 0.2 and k = 0.5, in particular in terms of the distributions of the average direction v taken at different time snapshots. At k = 0.2 the average t1 = 280 and at k = 0.5 the average t1 = 34; these two average times represent when the rate of change of the system's activity is greatest. Figure 12 shows histograms of v for time snapshots before t1 (first row) at t1 (second row) and after t1 (third row). In the first case k = 0.2 there is a tri-modal short-term response dominated by D (panel (a)), there is a tri-modal medium-term response spread between H, D and V (panel (c)) and there is a late-term tri-modal response dominated by H and V (panel (e)). In the second case k = 0.5, we see a shift from a uni-modal, short-term response (panel (b)), that becomes a bimodal response for the medium-and late-term with no peak at v = 0 (panels (d), (f)). In the second case the unstable solution D does not have an distinguishable effect on the dynamics, there is a rapid convergence to one of the two stable solutions.

Discussion

In this paper we presented a spatialised model of direction selection that was used to study the dynamical behaviour of multistable responses to the 1 : 1 ratio barber pole (a diagonally drifting grating viewed through a square aperture). The aim was to reproduce, at the cortical level, firing rate activity that can be related to qualitative behaviour observed in psychophysics experiments: the stimulus exhibits multistability of perception for short presentations where the dominant percepts are either 1) the diagonal direction of the grating (D) 2) in agreement with horizontal (H) aperture edge, or 3) in agreement with the vertical (V) aperture edge (Castet et al, 1999;[START_REF] Fisher | The directional tuning of the barber-pole illusion[END_REF]. The model reproduced this multistable behaviour, where the percepts H and V correspond to stable solutions and the percept D to an unstable solution. The temporal dy-namics were investigated and it was shown that early responses were dominated by the diagonal percept, midterm responses were tri-modal between the three percepts and later responses were dominated by the two stable percepts H and V. This behaviour is consistent with experimental findings that show an early response dominated by 1D cues, which is later refined by 2D cues [START_REF] Barthélemy | Temporal dynamics of 2d motion integration for ocular following in macaque monkeys[END_REF][START_REF] Lorenceau | Different motion sensitive units are involved in recovering the direction of moving lines[END_REF][START_REF] Pack | Temporal dynamics of a neural solution to the aperture problem in visual area MT of macaque brain[END_REF]. One of the main predictions to be made from these results is that the percept D is only seen as a transient behaviour; for long term presentations on the order of seconds either H or V will be seen or perceived.

One of the main advantages of the model used here is its simplicity; the philosophy was to reproduce interesting behaviour observed in experiments with a minimalistic set of features to perform motion integration. Given the particular stimulus studied in this paper and its inherent symmetry properties, it was important to utilise a framework, such as neural fields, where these symmetry properties can be preserved. Another positive aspect of working with a relatively simple model is the small set of parameters that must be determined, in contrast to, for example, the study [START_REF] Chey | Neural dynamics of motion processing and speed discrimination[END_REF] where a huge number of parameters must be determined heuristically. The strategy employed for setting model parameters took into account a number of important biological and mathematical constraints. One of the main ideas was to ensure that the model is operating close to its principal bifurcation where it will be most sensitive to subtleties of the stimulus input [START_REF] Veltz | Nonlinear analysis methods in neural field models[END_REF]; the bifurcation analysis and other analytical results helped to ensure the model was operated in the right parameter regime. For the remaining parameters, we used numerical continuation to study the relationships between a given parameter and biologically relevant properties of the model's solutions. We looked at these relationships over a wide range of parameter values and ensured that appropriate values were set. The general principles applied to tune parameters are applicable to a broad class of models that covers the various possible extensions proposed below. In existing studies of motion integration, behaviour was studied at fixed parameters [START_REF] Chey | Neural dynamics of motion processing and speed discrimination[END_REF][START_REF] Simoncelli | A model of neuronal responses in visual area MT[END_REF][START_REF] Bayerl | Disambiguating visual motion through contextual feedback modulation[END_REF]Tlapale et al, 2010b); in certain cases the influence of a single term is tested by setting its weight to zero. Here the two-parameter continuation analysis allow us to determine an entire region of interest in parameter space. Further, the extended investigation in this region looked at the dynamical behaviour taking into account the effect of changing initial conditions. The two-parameter investigation also tells us that the behaviour produced is robust over entire regions of parameter space.

An important question in studies of motion perception is identifying exactly what happens in the first few hundred milliseconds of presentation. The most insightful work so far has come from ocular following experiments [START_REF] Masson | Temporal dynamics of motion integration for the initiation of tracking eye movements at ultra-short latencies[END_REF] and physiological recordings from individual neurons [START_REF] Pack | Temporal dynamics of a neural solution to the aperture problem in visual area MT of macaque brain[END_REF]. In order to support the results in this paper we propose a psychophysics experiment that investigates the same stimulus described in this paper with incrementally increasing short presentation time in order to identify the qualitative change in percept distributions that were found in our model. We found a uni-modal distribution in the short-term dynamics and a tri-modal distribution for the medium-term dynamics. Furthermore, we predict that the D percept can only be seen as a transient and for long enough presentation times only the percepts H or V would be seen.

The model and results presented here not only captures a number of important aspects seen in experiment, but also forms a solid basis for further study of motion integration. One natural extension would be to lift the one dimensional approximation of the physical space. This would first of all allow for the validity of the original approximation to be tested, further, it would allow for other stimulus parameters such as the terminator ratio to be tested for further comparison with the results of [START_REF] Fisher | The directional tuning of the barber-pole illusion[END_REF]. More complex aperture arrangements could also be considered such as the cross-shaped stimulus studied in [START_REF] Castet | Long-range interactions in the spatial integration of motion signals[END_REF]. A directional bias that varied from subject to subject but generally towards the horizontal [START_REF] Fisher | The directional tuning of the barber-pole illusion[END_REF] is a feature of experimental results that was not captured by our model. The most straightforward way to investigate this would be to consider an asymmetry in the model's input that assumes stronger (or more numerous) inputs at the preferred direction. The model in (Tlapale et al, 2010b), implemented in a similar framework to the one studied here, considers a much more detailed description taking into account filtering stages applied to the input along with feedforward and feedback interactions between MT and V1. The methods used in this paper would be applicable to models of much greater complexity though at the cost of increasing the intricacy of analytical computations and scale of numerical computations. The focus of this article was multistable perception in the first few seconds after stimulus onset, a similar model and stimulus could be used to study perceptual switches that are known to occur for extended presentations by considering an adaptation dynamic on a slow timescale such as the one studied in [START_REF] Curtu | Pattern formation in a network of excitatory and inhibitory cells with adaptation[END_REF]. defined for all times t ≥ 0 and continuously differentiable with respect to t for all t ≥ 0.

Regarding the applicability of the parameter dependent centre manifold theorem (Haragus and Iooss, 2010, Theorem 3.3, Chapter 2), we consider the linear operator L λ defined in Eq. ( 22):

L λ = -µId + λS 1 J.
We note that J can be written Kν 3 Id, where K is a convolution operator (in effect a linear combination of two convolutions with bounded convolution kernels). From these remarks it follows that L λ is in L(G, F), the set of linear continuous operators from G to F, and that it is a compact operator, hence that its spectrum is discrete and the eigenvalues have finite multiplicity.

It is also not difficult to prove that (iωId -L λ ) -1 L(F ) ≤ 2 |ω| and (iωId -L λ ) -1 L(G) ≤ 2 |ω| for ω large enough. Therefore, all hypotheses necessary for applying the parameter dependent centre manifold theorem [START_REF] Haragus | Local bifurcations, center manifolds, and normal forms in infinite dimensional systems[END_REF] are satisfied and the computation of the normal form Eq. ( 21) is justified.

B Normal form coefficient

In order to compute b we need to solve:

L λc Ψ 21 + λ 2 c S 2 2 (2JΨ 11 ζ 1 χ + 2JΨ 20 ζ 1 χ) + λ 3 c S 3 6 3ζ 3 1 χ 2 χ = bχ. ( 24
)
It is first necessary to obtain expressions for Ψ 11 and Ψ 20 by solving the following equations:

L λc Ψ 20 + λ 2 c S 2 2 ζ 2 1 χ 2 = 0, (25)
L λc Ψ 11 + λ 2 c S 2 ζ 2 1 χχ = 0. (26)
In order to solve Eq. ( 25) we decompose Ψ 20 into components on the eigenspace v k e ikv , for coefficients v k to be determined. From

L λc Ψ 20 = - λ 2 c S 2 2 ζ 2 1 χ 2 ,
we find that v k = 0 for k = ±2 and that

v 2 = v -2 = - λ 2 c S 2 ζ 2 1 2(-µ + λ c S 1 ζ 2 ) .
The ζ 2 coefficient comes from the action of J on Ψ. Finally we have a general expression for Ψ 20 :

Ψ 20 = - λ 2 c S 2 ζ 2 1 2(-µ + λ c S 1 ζ 2 ) e i2v + α 20 χ + α 20 χ.
We obtain a similar expression for Ψ 11 :

Ψ 11 = - λ 2 c S 2 ζ 2 1 -µ + λ c S 1 ζ 0 + α 11 χ + α 11 χ.
The ζ 0 coefficient comes from the action of J on Ψ. These expressions are further simplified by the fact that at the bifurcation point ζ 1 = µ λcS1 , so we have: .

Ψ 20 = λ 2 c S 2 ζ 2 1 2µ(1 -ζ2 ζ1 ) e i2v + α 20 χ + α 20 χ, (27)
Finally, we obtain the expression for b given in Eq. ( 23).

Figure 1 :

 1 Figure 1: Description of the model stimulus. (a) A sinusoidal luminance grating viewed through a square aperture (grey); a large arrow indicates the diagonal grating direction D (v = 0 • ). Indicated by small arrows are the horizontal H (v = 45 • ) and vertical V (v = -45 • ) directions corresponding to the aperture edges. (b) A cut across the aperture is used to represent the stimulus in the 1D physical space x (black line). At points on the interior of the aperture (black line is dashed) many directions are stimulated, on the edges of the aperture a single direction is stimulated parallel to the edge (black points) and outside the aperture no direction is stimulated (black line is solid). (c) Representation of the stimulus I ext in physical-direction space (x, v); at each point x each of the possible directions v ∈ (-180 • , 180 • ) is either stimulated (white) or not stimulated (black). Multiple directions are stimulated on the interior of the aperture x ∈ (-0.75, 0.75), unique directions are stimulated at the aperture edges x = -0.75 and x = 0.75. Note that the 1D cut in panel (b) is illustrative and the actual interior and exterior regions used in the model are those shown in panel (c).

Figure 2 :

 2 Figure 2: Profile of the connectivity function. (a) Form of the connectivity M x in x-space and the corresponding Fourier transform Mx . (b) Form of the connectivity M v in v-space and the corresponding Fourier transform Mv . (c) Full, coupled connectivity M . The functions M x and M v correspond to the 1D slices at v = 0 (black curve) and x = 0 (white curve), respectively.

Figure 3 :

 3 Figure 3: Sketch of bifurcation diagrams in λ for some norm |.| of a (a) supercritical or (b) subcritical pitchfork. Lower solution branch corresponds to spatially homogeneous solutions, upper solution branch corresponds to direction-selected solutions; stable parts are black, unstable parts are grey. Pitchfork occurs at λ c and in (b) there is a fold bifurcation on the direction-selected branch λ f ; arrows indicate temporal dynamics when passing λ c on the lower branch.

Figure 4 :

 4 Figure 4: Tuning the parameters T and ν 2 . (a) Grey scale map of activity level pc at bifurcation point λ c ; plotted in both panels are the contours 5% and 20% (grey) along with the contour 99% (dashed black). (b) Grey scale map of the normal form coefficient b as given by equation Eq. (23); plotted in both panels are the contours b = 0 (black); positive values correspond to a subcritical bifurcation and negative values to a supercritical bifurcation. For |T | large (white region in bottom left of (a), grey region in bottom left of (b)) the activity level is saturated at pc = 100% and b ≈ 0 because the Taylor coefficients S * ≈ 0.

Figure 5 :

 5 Figure 5: Model solutions in absence of stimulus. (a) A bifurcation diagram varying λ for k = 0 and the parameter values given in Table 1, where the solution norm is |p|; branches of steady-state solutions are plotted with stable sections black and with unstable sections grey. The lower branch of solutions are spatially uniform in both x and v ; for increasing λ the branch loses stability in a subcritical pitchfork bifurcation at P 1 . The bifurcated branch represents direction-selected solutions (modulated in v but still uniform in x); a change in stability occurs at the fold bifurcation F 1 . (b) Activity p plotted over the (x, v)-plane for the direction selected solution indicated by a diamond marker on the upper branch in panel (a). (c) Cross section in v at an arbitrary value of x from panel (b).

Figure 6 :

 6 Figure 6: Tuning of parameters ν 2 and ν 3 . (a) Dependence of velocity tuning half-width w 12 on ν 2 at fold point F 1 ; black curve computed for ν 3 = 0, grey curve computed for ν 3 = 1.5. (b) Dependence of maximum activity level for direction-selected solutions p max on ν 3 computed for ν 2 = 66.

Figure 7 :

 7 Figure 7: Branches of steady-state solutions are plotted for varying bifurcation parameter k (shown on a logarithmic scale) as computed for different values of λ; stable sections of solution branches are black and unstable sections are grey. Changes in stability occur at fold bifurcations (black points) and pitchfork bifurcations (black stars). In panels (a), (b) and (c) solutions are plotted in terms of the L 2 -norm of the solution vector |p| for λ = {12.5, 13, 14}, respectively. The inset in panel (c) shows the largest eigenvalue ζ max along the two upper solution branches from the main panel over the range of k as indicated by vertical dashed lines. Panel (d) shows the same solution branches as in panel (c), but in terms of the average direction v.

Figure 8 :

 8 Figure 8: Three possible direction-selected solutions. Panels (a), (b) and (c) show, for λ = 14, k = 0.3, three possible direction-selected solutions corresponding to the unstable diagonal (D), stable vertical (V) and stable horizontal (H) percepts , respectively.

Figure 10 :

 10 Figure10: Example of temporal dynamics for (λ, k) = (14, 0.3) initialised from random initial conditions. (a) Grey scale map showing the evolution activity levels p(x, v, t) with time on the vertical axis. The state vector p is indexed locally by position in physical space x and globally by position in direction space v and it is, therefore, convenient to show v-position on the horizontal axis. At the time-point labelled t 1 (horizontal white dashed line, star in panels (b) and (c)) there is a transition to the diagonal direction-selected solution. At the time-point labelled t 2 (horizontal white line, point in panels (b) and (c)) there is a transition to the final steady-state, direction-selected solution. (b) Time-evolution of |p|; the time-mark t 1 is the point of steepest increase for |p|. (c) Time evolution of the average direction v with horizontal dashed black lines indicating the two possible final steady-state solutions. The time-mark t 2 is taken as the first point where the v is within 3 • of the final selected direction.

  The next technical step is to do the change of variables p = pc +P , λ = λ c +η and to write an equation for P as a function of the parameter η. In detail we have∂P ∂t = ∂p ∂t = -µ(p c + P ) + S((λ c + η)[J(p c + P ) + T ]).We subtract from this equation0 = -µp c + S(λ c [J pc + T ]),to obtain∂P ∂t = -µP + S((λ c + η)[J(p c + P ) + T ]) -S(λ c [J pc + T ]).We then write∂P ∂t = L λc P + S((λ c + η)[J(p c + P ) + T ]) -S(λ c [J pc + T ])λ c S 1 JP. If we define G(p, λ) = S(λ[Jp+T ]), it should be clear that λ c S 1 JP = D p G(p c , λ c )P .From this follows the fact thatR(P, η) = S((λ c + η)[J(p c + P ) + T ]) -S(λ c [J pc + T ])λ c S 1 JP = G(p c + P, λ c + η) -G(p c , λ c ) -D p G(p c , λ c )P,is such that (use first order Taylor's expansion with integral remainder) R(0, 0) = 0, D p R(0, 0) = 0, and R is infinitely differentiable, because the sigmoid S is.

E 1 Ψ

 1 || 20 = α 20 χ + α 20 χ, for some arbitrary coefficients α 20 and α 20 , and components orthogonal to the eigenspace Ψ 20 = |k| =1

  α 11 χ + α 11 χ. (28) Inria Now to solve for b we take the inner product of Eq. (24) with χ which gives b = LΨ 21 , χ =0 +λ 2 c S 2 ζ 1 ( JΨ 11 χ, χ + JΨ 20 χ, χ ) +

Table 1 :

 1 Default parameter values used in the numerical studies in Secs. 4.3 and 4.4.

	Inria

RR n°7822

route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgements

The authors would like to thank Guillaume Masson, Andrew Meso and Gregory Faye for their insight and for their helpful criticism of the work presented here. This research work received funding from CNRS, the European Community (through FACETS, IST-FET, Sixth Framework, No 025213, the ERC grant No 227747 (NERVI) and the Agence Nationale de la Recherche (ANR, NATSTATS)).

Inria

A Centre manifold

Here, we show that the necessary hypotheses are satisfied in order to apply the centre manifold theorem as given in (Haragus and Iooss, 2010, Chapter 2). We introduce some functional spaces in which the problem defined by Eq. ( 2) is well-posed. We note F the Hilbert space L 2 (Ω × V, R), where Ω = R/cZ and V = R/2πZ , endowed with the usual inner product

We also note G the space L ∞ (Ω × V, R) ⊂ F. It is then easy to show, using techniques similar to those in [START_REF] Faugeras | Abolute stability and complete synchronization in a class of neural fields models[END_REF] and (Veltz, 2011, Chapter 4), that for any initial condition p 0 (x, v) ∈ F there is a unique solution p(x, v, t), Inria