J Rankin

A I Meso

G S Masson

O F A U G E R A S *

P Kornprobst

Bifurcation Study of a Neural Fields Competition Model with an Application to Perceptual Switching in Motion Integration

Keywords: multistability, competition, perception, neural fields, bifurcation, motion

The phenomenon of perceptual multistability in which alternate interpretations of a fixed stimulus are perceived intermittently is an active area of research and the underlying mechanisms that gate perception are little understood. Numerical tools from bifurcation analysis are applied to the study of a competition model posed as a feature-only neural field equation. In the absence of input the model with a spike-frequency adaptation mechanism has been shown to produce an array of complex spatio-temporal dynamics local to a Bogdanov-Takens point. Here, we demonstrate how, with the introduction of an input, the organisation of solutions in parameter space changes via symmetry breaking. The model is then used to investigate a more complex stimulus in the context of motion integration that is multistable in terms of its perceived direction of motion, the so-called multistable barber pole, which has been the subject of concurrent psychophysics experiments. We bring the model to an operating regime where known physiological response properties are reproduced whilst also working close to bifurcation. We find that in this regime the model is able to account for characteristic behaviour from experiment in terms of the type of switching observed and changes in the rate of switching with respect to contrast. In this way, the modelling study sheds light on the underlying mechanisms that drive perceptual switching in different contrast regimes. The general approach presented is applicable to a broad range of perceptual competition problems in which spatial interactions play a role.

Etude des bifurcations d'un modèle de compétition en champs neuronaux avec une application aux changements perceptifs dans l'intégration du mouvement

Résumé : Le phénomène de multistabilité perceptive au cours duquel plusieurs interprétations d'un même stimulus alternent est un domaine de recherche actif et les mécanismes sous-jacent qui guident cette perception restent mal compris. Dans ce papier, les outils numériques d'analyse de bifurcations sont appliqués à l'étude d'un modèle de compétition décrit dans le cadre des champs neuronaux. En l'absence d'entrée et avec un mécanisme d'adaptation de type spike-frequency,l e modèle produit un ensemble de dynamiques spatio-temporelles complexes au voisinage d'un point de Bogdanov-Takens. Dans ce papier, nous montrons comment, avec l'introduction d'une entrée, l'organisation des solutions dans l'espace des paramètres change via la cassure de symétrie. Le modèle est ensuite utilisé pour étudier un stimulus plus complexe dans le cadre de l'intégration de mouvement qui est multistable en terme de direction perçue du mouvement. Ce stimulus est le barber pole et une étude psychophysique a été menée simultanément. Nous plaçons le modèle dans un régime de fonctionnement où les propriétés physiologiques des réponses sont reproduites, tout en travaillant également à proximité de la bifurcation. Nous trouvons que dans ce régime, le modèle est capable de rendre compte des comportements caractéristiques observés expérimentalement, en terme de nature des changements perceptifs et du taux de changements en fonction du contraste. Cette étude met en lumière les mécanismes sous-jacents qui entrainent les changements perceptifs en fonction du contraste. L'approche générale présentée est applicable àunlargeév en taildeproblèmesdecompétitionauniv eauperceptifpourlesquelslesin teractions spatiales jouent un rôle.

Mots-clés : multistabilité, compétition, perception, champs neuronaux, bifurcation, mouvement

Introduction

The link between neural activity and conscious perception is possibly among the least understood and more interesting aspects of neuroscience. Perception can evolve dynamically for fixed sensory inputs and so-called multistable stimuli have been the attention of much recent experimental and computational investigation. Working within a mathematical framework where neural activity is described at the population level, and using powerful computational tools from the field of dynamical systems, we aim to gain a deeper understanding of the neural mechanisms that underpin the complex computations driving perception. The neural field equations provide an established framework for studying the dynamics of cortical activity, represented as an average membrane potential or mean firing rate, over a spatially continuous domain. Since the seminal work by [START_REF] Amari | Characteristics of randomly connected threshold element networks and neural systems[END_REF]; Wilson andCowan (1972, 1973) a broad range of mathematical tools have been developed for their study; see reviews by [START_REF] Ermentrout | Neural networks as spatio-temporal pattern-forming systems[END_REF]; [START_REF] Coombes | Waves, bumps, and patterns in neural fields theories[END_REF]; [START_REF] Bressloff | Spatiotemporal dynamics of continuum neural fields[END_REF] along with Ermentrout and Terman (2010, Chapter 11) for a derivation of the equations. The equations describe the dynamical evolution of activity of one or more connected populations of neurons, each defined in terms of a spatial domain that can represent either physical space (on the cortex), an abstracted feature space (orientation, direction of motion, texture preference, etc.), or, some combination of the two. This framework has proved especially useful in the study of neuro-biological phenomena characterised by complex spatio-temporal patterns of neuronal activity, such as, orientation tuning in the primary visual cortex V1 [START_REF] Ben-Yishai | Theory of orientation tuning in visual cortex[END_REF][START_REF] Somers | An emergent model of orientation selectivity in cat visual cortical simple cells[END_REF][START_REF] Hansel | Methods in neuronal modeling: From ions to networks[END_REF][START_REF] Veltz | Lo cal/global analysis of the stationary solutions of some neural field equations[END_REF], cortical waves [START_REF] Coombes | Waves and bumps in neuronal networks with axo-dendritic synaptic interactions[END_REF][START_REF] Folias | Breathing pulses in an excitatory neural network[END_REF][START_REF] Laing | Spiral waves in nonlocal equations[END_REF], binocular rivalry [START_REF] Kilpatrick | Binocular rivalry in a competitive neural network with synaptic depression[END_REF][START_REF] Bressloff | Neural field model of binocular rivalry waves[END_REF], and motion integration [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]. Ac ruc ialtoolinthes tudyoftheneur alfie lde quations,anddynam ic alsy ste m singe ne ral,is bifurcation analysis. Bifurcations are special points at which there are qualitative changes to the underlying solution structure of a dynamical system, under the variation of model parameters. The model's dynamics is governed by this solution structure in terms of its coexisting solutions, their spatial properties and stability. The typical approach taken in many parametric studies of as i m p l en u m e r i c a ls e a r c hc a no f t e nm i s si m p o r t a n tp a r a m e t e rr e g i o n sa n dd o e sn o tp r o v i d ei nformation about robustness of model behaviour with respect to parameter variation. Bifurcation analysis provides a solution to this problem and has shown to be a powerful tool in a very broad set of applications including, but in no way limited to physics, engineering and the biological sciences; for a general introduction see [START_REF] Strogatz | Nonlinear dynamics and chaos. Addison-Wesley Reading[END_REF]; [START_REF] Kuznetsov | Elements of Applied Bifurcation Theory[END_REF]. The methods have also been widely applied to a range single neuron models [START_REF] Ermentrout | Foundations of Mathematical Neuroscience[END_REF]. In the study of the neural field model presented here, an understanding of specific methods applicable to infinite dimensional (spatially continuous) dynamical systems in the presence of symmetries is required; 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]. Within the framework of neural fields, 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 with interesting symmetry properties 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 (Laing et al, 2002;[START_REF] Guo | Existence and stability of standing pulses in neural networks: I. existence[END_REF][START_REF] Faye | Lo calized radial bumps of a neural field equation on the euclidean plane and the poincaré disk[END_REF] and in a model of texture perception [START_REF] Faye | Analysis of a hyperbolic geometric model for visual texture perception[END_REF]. The computational counterpart to the analytical tools of bifurcation analysis is numerical continuation [START_REF] Krauskopf | Numerical Continuation Methods for Dynamical Systems[END_REF]. These numerical schemes allow one to locate, classify and track bifurcation points, which provides a means to map out regions in parameter space with qualitatively different behaviour and accurately compute the boundaries between these regions; e.g. in a model of working memory, numerical continuation was used to map out parameter regions for which the application-relevant patterns of spatially localised activity are observed as opposed to non-localised spatially periodic patterns [START_REF] Faye | Lo calized radial bumps of a neural field equation on the euclidean plane and the poincaré disk[END_REF]. This kind of information forms a basis for tuning a model's parameters; indeed, it is possible to ensure that parameter regions in which a desired behaviour is present are not isolated and ensure robustness with respect to small changes in the model set up. Numerical continuation has been used in general studies of the neural field equations (Veltz and Faugeras, 2010;[START_REF] Veltz | Lo cal/global analysis of the stationary solutions of some neural field equations[END_REF], to investigate localised states [START_REF] Laing | PDE methods for nonlocal models[END_REF][START_REF] Faye | Lo calized radial bumps of a neural field equation on the euclidean plane and the poincaré disk[END_REF] and to study rotating and spiral waves [START_REF] Laing | Spiral waves in nonlocal equations[END_REF][START_REF] Owen | Bumps and rings in a two-dimensional neural field: splitting and rotational instabilities[END_REF]. One key advantage of using bifurcation and continuation techniques is that they allow for a model to be brought into an operating regime, close to bifurcation, where the model is most sensitive to changes in its input and where the combination of mechanisms involved in performing complex computations can be revealed. This general philosophy has been used to great effect in studies of orientation tuning in V1 (Veltz, 2011, Chapter 9), simplified rate models of neural competition [START_REF] Shpiro | Dynamical characteristics common to neuronal competition models[END_REF]Theodoni et al, 2011b) and studies of decision making (Theodoni et al, 2011a). Furthermore, it is possible to incorporate known biological and experimental observations as done with the orientation tuning width of V1 responses (Veltz, 2011, Chapter 9).

In this paper we are interested in certain ambiguous visual stimuli for which two or more distinct interpretations are present, but where only one of these interpretations, or percepts, can be held at a time. Not only can the initial percept be different from one short presentation to the next, but for extended presentations, the percept can change, or switch, dynamically. This phenomenon of multistability has been observed and investigated with a number of different experimental paradigms including both static and motion displays, e.g. ambiguous figures [START_REF] Necker | Observations on some remarkable optical phaenomena seen in Switzerland; and on an optical phaenomenon which occurs on viewing a figure of a crystal or geometrical solid[END_REF][START_REF] Rubin | Visuell wahrgenommene figuren: studien in psychologischer analyse[END_REF], binocular rivalry experiments [START_REF] Levelt | On binocular rivalry[END_REF][START_REF] Blake | A neural theory of binocular rivalry[END_REF][START_REF] Blake | A primer on binocular rivalry, including current controversies[END_REF], random-dot rotating spheres (kinetic depth effect) [START_REF] Wallach | The kinetic depth effect[END_REF][START_REF] Sperling | Early vision and beyond, Bradford Books, chap Depth from motion[END_REF][START_REF] Stonkute | The role of attention in ambiguous reversals of structure-from-motion[END_REF], apparent motion [START_REF] Ramachandran | Bifurcation analysis applied to a model of motion integration with a multistable stimulus[END_REF], motion plaids that are bistable [START_REF] Hupé | The dynamics of bi-stable alternation in ambiguous motion displays: a fresh look at plaids[END_REF] or tristable [START_REF] Hupé | The initial phase of auditory and visual scene analysis[END_REF] and the multistable barberpole illusion [START_REF] Castet | The extrinsic/intrinsic classification of two-dimensional motion signals with barber-pole stimuli[END_REF][START_REF] Fisher | The directional tuning of the barber-pole illusion[END_REF]Meso et al, 2012b). During extended presentations of these stimuli, the dominant percept switches randomly and the dominance durations between switches have been shown to fit certain distributions dependent on the experimental paradigm [START_REF] Levelt | On binocular rivalry[END_REF][START_REF] Leopold | Activity changes in early visual cortex reflect monkeys' percepts during binocular rivalry[END_REF][START_REF] Logothetis | What is rivalling during binocular rivalry?[END_REF][START_REF] Lehky | Binocular rivalry is not chaotic[END_REF][START_REF] Zhou | Perceptual dominance time distributions in multistable visual perception[END_REF][START_REF] Rubin | Binocular rivalry, The MIT Press, chap Dynamics of perceptual bistability: Plaids and binocular rivalry compared[END_REF]. The focus of many perceptual competition modelling studies has been to reproduce the switching behaviour observed in experiment and provide insight into the underlying mechanisms [START_REF] Laing | A spiking neuron model for binocular rivalry[END_REF][START_REF] Freeman | Multistage mo del for bino cular rivalry[END_REF][START_REF] Kim | Stochastic resonance in binocular rivalry[END_REF][START_REF] Shpiro | Dynamical characteristics common to neuronal competition models[END_REF][START_REF] Moreno-Bote | Noise-induced alternations in an attractor network model of perceptual bistability[END_REF][START_REF] Borisyuk | Visual perception of ambiguous figures: synchronization based neural models[END_REF][START_REF] Ashwin | A low-dimensional model of binocular rivalry using winnerless competition[END_REF]. Bifurcation analysis and numerical continuation have already proved effective in the study of reduced rate models where the competing states are represented by individual space-clamped neuronal populations [START_REF] Shpiro | Dynamical characteristics common to neuronal competition models[END_REF]; Theodoni et al (2011b). Two commonly proposed mechanisms that drive the switching behaviour in these models are adaptation and noise. In [START_REF] Shpiro | Balance between noise and adaptation in competition models of perceptual bistability[END_REF] a strong argument is made that a balance of adaptation and noise accounts best for experimental findings across different model architectures and different adaptation mechanisms.

Here we will take advantage of the neural fields formalism and the strengths of bifurcation and continuation methods in order to study neural competition in a model with a continuous feature space where adaptation and noise are implemented as mechanisms that can drive activity switches. The model describes the mean firing rate of a populations of feature selective neurons. Indeed, this feature-only model with spike frequency adaptation has been studied previously [START_REF] Curtu | Pattern formation in a network of excitatory and inhibitory cells with adaptation[END_REF], but in the absence of an input and noise. A key difference with existing rivalry models is that, instead of representing competing percepts by discrete populations, the competing percepts form tuned responses in a continuous feature space. The more general model presented here allows for perceptual transitions to occur in a smooth way as opposed to discrete switches between two isolated percepts. Starting from the results presented Inria in [START_REF] Curtu | Pattern formation in a network of excitatory and inhibitory cells with adaptation[END_REF], we will introduce first a simple input and investigate how the various types of solutions previously found in their work are modified. We find that although the boundaries between parameter regions featuring different types of responses are gradually distorted with increasing input strength, much of the global structure is preserved. This allows for all possible types of behaviour, and parameter regions for which it can occur, to be comprehensively described across a wide range of model parameters controlling input gain, adaptation gain and the shape of the firing rate function. For a simple input we are able to match the models output to known response properties from the literature before considering the introduction of more complex inputs that give rise to multistable behaviour.

As an example of application of the proposed model, we will investigate the temporal dynamics of perception for a particular multistable motion stimulus. The so-called multistable barberpole illusion was investigated in complementary psychophysics experiments (Meso et al, 2012b); some of these results will be presented alongside the modelling work. We will demonstrate how the general neural fields model can reproduce the main dynamical characteristics of the perceptual switches observed in the experiments. The model also captures the relationship between contrast and the rate of switching. Importantly, the two contrast regimes identified experimentally, one in which the rate increases with contrast, the other in which the rate decreases with contrast, are captured by different features in the model. This provides insight into the mechanisms that drive switching at different contrast regimes. Although a combination of noise and adaptation drive the switching, the dominant mechanism changes with contrast. Furthermore, we are able to quantify this in an experimentally testable way: the distribution of dominance durations fit different statistical distributions in each contrast regime.

In Sec. 2 section we give a mathematical description of the model before presenting general results that map out the model's possible behaviours across parameter space in Sec. 3 and then applying the model to the study multistable perception in Secs. 4 and 5.

Competition model with continuous feature space

In this section we describe a general neural competition model that implements adaptation and noise as mechanisms that, for an ambiguous input, drive switches between selected states. The model represents the time-evolution of a neuronal population defined across a continuous feature space in which a selected state corresponds to a tuning curve. A spatial connectivity is chosen such that a winner-takes-all mechanism leads to a definitive tuned response at any given time instant; the connectivity produces mutual inhibition between competing tuned responses. Over time, shifts between tuned responses are driven by a combination of adaptation and noise.

Model equations

We will consider a single neuronal population p(v, t),d e fi n e da c r o s st h ec o n t i n u o u s ,p e r i o d i c feature space v 2 [-π, π), whose evolution depends on time t. The variable p(v, t) takes values in [0, 1] representing activity as a proportion of a maximal firing rate normalised to 1.W ea l s o define secondary adaptation α(v, t) and stochastic X(v, t) variables. The time evolution of p(v, t) is described by the following coupled system of integro-differential equations:

d dt p(v, t)=-p(v, t)+S λ[J(v) ⇤ p(v, t) -k α α(v, t) + k X X(v, t)+k I I(v) -T ] , (1)
τ α d dt α(v, t)=-α(v, t)+p(v, t). (2)
The principal equation Eq. ( 1) evolves with cortical timescale 1ms and has a standard decay term -p.As m o o t h ,n o n l i n e a rs i g m o i d a lfi r i n gr a t ef u n c t i o n

S(x)= 1 1+exp(-x) (3)
is used as plotted in Fig. 1(a). The slope and threshold of the firing rate function are controlled by the parameters λ and T , respectively. The firing rate function processes lateral connections described by J and inputs from adaptation α,a d d i t i v en o i s eX and a time independent input I; the respective input gain parameters are k α , k X and k I . The connectivity in the feature space v is represented by a convolutional operator J that approximates a Mexican hat connectivity (local excitation, lateral inhibition). As in [START_REF] Curtu | Pattern formation in a network of excitatory and inhibitory cells with adaptation[END_REF], we use a 3-mode Fourier expansion and J takes the form

J(v)=J 0 + J 1 cos(v)+J 2 cos(2v), (4)
see Fig. 1(b). The adaptation dynamic in Eq. ( 2) describes linear spike frequency adaptation that evolves on a slow time scale τ α . The additive noise X(v, t) variables is an Ornstein-Uhlenbeck process that evolves on the same slow timescale τ α as the adaptation, has mean hX(v, t)i =0 , variance Var(X)=1and no feature correlation; see Appendix A for further details. The input I depends only on the feature v and the so-called simple input studied in Sec. 3 is shown in Fig. 1(c).

. .

Parameter values, initial conditions and numerical computations

The parameter values used in the numerical computations in each section of the paper are given in Table 2 we solve the system described by Eqs. ( 1)-(2) using the ODE23t solver in Matlab with default settings except the relative tolerance, which was set to 10 -6 .A200-point discretisation of v was used that satisfies error bounds for the computation of the integral term J ⇤ p with a standard trapezoid method. For the simulations with noise (k X =0.0025), as presented in Sec. 5, the same discretisation is used with a standard Euler-Maruyama method and a fixed timestep of 0.5ms. In all simulations the initial conditions are set to a low level of activity p 0 (v)=0.1 (10% of the maximal firing rate) with a small randomised perturbation. Initial conditions are given by p(v, 0) = p 0 (v),α (v, 0) = 0,X (v, 0) = 0.

In order to carry out a bifurcation analysis of the system (1)-(2) we use a numerical continuation package AUTO [START_REF] Doedel | Auto 97: Continuation and bifurcation software for ordinary differential equations (with homcont[END_REF]) that allows us to compute branches of steady state and oscillatory solutions, and to detect and track bifurcations of these solutions. These computations are forcibly done in the absence of noise (k X =0)a n di nt h i sc a s ew ec a nt a k ea d v a n t a g eo ft h e 3-mode approximation of J (4) by expressing p(v, t) and α(v, t) in terms of the same modes plus some orthogonal components b p ? and b α ? :

p(v, t)=b p 0 (t)+b p 1 (t) cos(v)+b p 2 (t)sin(v) + b p 3 (t) cos(2v)+b p 4 (t)sin(2v)+b p ? , (5)
α(v, t)=b α 0 (t)+b α 1 (t) cos(v)+b α 2 (t)sin(v) + b α 3 (t) cos(2v)+b α 4 (t)sin(2v)+b α ? . (6)
In Veltz and Faugeras (2010) it was proved that as t !1the orthogonal components decay to 0. Therefore, we can study steady-state and oscillatory solutions to (1)-(2) by solving a set of 10 ordinary differential equations in b p i and b α i , i =0...4. The integral term J ⇤ p was computed with the same 200-point trapezoid integration scheme used for in the ODE solver. Periodic orbits were typically computed with 150 mesh points (constants NTST=50 and NCOL=3) in AUTO. The reduced description is used only for the computation of the bifurcation diagrams.

General study of competition model

Recall that a bifurcation is a special point at which, under the variation of a model parameter, there is a qualitative change to the types of response produced by the model. The types of response can differ in terms of 1) spatial properties, such as being tuned or untuned, 2) temporal properties, such as being steady or oscillatory and 3) stability, where stable implies responses that persist in time and unstable implies responses that are transient. Each type of response is associated with a solution of the underlying equations and the organisation of these solutions in state/feature space governs the dynamical behaviour. The main types of bifurcation that we encounter in this study are the pitchfork bifurcation and the Hopf bifurcation; these are so-called codimension-one bifurcations, that occur at a point as one parameter is varied. In this model the pitchfork is associated with a transition from a homogeneous state to a tuned state. Hopf bifurcations are associated with transitions from steady to oscillatory responses that can either travel in the feature space (travelling waves), or that remain static but oscillate in amplitude (standing waves). In the parameter plane these codimension-one bifurcations lie on curves. At points where these curves meet or intersect we encounter codimension-two bifurcations that act as organising centres close to which several solution types can be encountered. Due to the presence of translational and reflectional symmetry properties in the underlying equations when there is no input, the bifurcating solutions encountered have the same symmetry properties [START_REF] Haragus | Local bifurcations, center manifolds, and normal forms in infinite dimensional systems[END_REF]. We will give an account of how these symmetry properties break-down with the introduction of an input.

In [START_REF] Curtu | Pattern formation in a network of excitatory and inhibitory cells with adaptation[END_REF] it was shown that in the absence of input (k I =0 )t h e model possesses O(2)-symmetry and as a consequence several types of steady solutions and oscillatory patterns exist in different parameter regimes local to the codimension-two Bogdanov-Takens (BT) point [START_REF] Dangelmayr | The Takens-Bogdanov bifurcation with O(2)-Symmetry[END_REF]. As we will see, this BT point acts as an important organising center in parameter space. First, in Sec. 3.1, we will give a summary of what is already known from [START_REF] Curtu | Pattern formation in a network of excitatory and inhibitory cells with adaptation[END_REF] in terms of solution branches and bifurcation curves that are relevant to this study (the account will not be exhaustive, omitting several bifurcation curves involving exclusively unstable solutions). Next, in Sec. 3.2, we cover new results describing, with the introduction of an input, how the solutions existing in different parameter regimes change and how the boundaries of these parameter regimes shift in parameter space. We investigate the model's behaviour in terms of two parameters, the adaptation strength k α and the sigmoidal slope λ.

No input (k I =0)

Figure 2 shows the different types of dynamical behaviour produced in different regions of the (λ, k α )-parameter plane as demarcated by bifurcation curves. The model simulations shown in panels (b)-(e) were performed at the corresponding points in panel (a). In each case the simulation time was chosen such that the model reaches its stable behaviour during the simulation; the stable behaviour is either a steady state or an oscillatory state. In the white region the model produces an homogeneous (untuned) steady-state response at a low level of activity as shown in panel (d). In the dark-grey region the model produces a steady-state response tuned to an arbitrary direction as shown in panel (e). The boundary between the white and dark-grey regions is a pitchfork curve P ;asλ is increased and the pitchfork bifurcation is encountered the homogeneous steady-state becomes unstable and a ring of tuning curves forms the stable behaviour. In the light-grey region the stable behaviour is a travelling-wave solution with an arbitrary direction in v;t h et r a n s i e n tb e h a v i o u ro b s e r v e db e f o r er e a c h i n gt h i ss t a b l es t a t ec h a n g e sd e p e n d e n to nt h e chosen parameter values; see panels (b) and (c). The boundary between the white region and the light-grey region is the coinciding Hopf-type curves H SW an H TW1 .A sλ is increased and the two coinciding bifurcation points are encountered the homogeneous steady states loses stability and two new branches bifurcate simultaneously: an unstable branch of standing wave solutions and a stable branch of travelling wave solutions; this will be shown explicitly in Sec. 3.3. In panel (b), close to these curves, the unstable standing wave solution is seen as a transient behaviour before eventual convergence to the stable travelling wave solution. The boundary between the dark-grey and light-grey region is H TW2 and as k α is increased and the bifurcation is encountered the stable tuned response becomes spatially unstable and starts to travel in an arbtirary direction. In panel (c) the unstable tuned response is seen as a transient behaviour before starting to travel. 3.2 Simple input with k I =0.001 The most notable result is that much of the structure from the no-input case has been preserved, albeit with subtle changes that are now discussed. In the white region (to the left of H SW1 and b P ) there is now a low-activity response that is weakly tuned to the input centered at v =0;seepanel (d). In the dark-grey region there is still a steady-state, tuned response, but now centered on the stimulus at v =0. In the light-grey region the stable behaviour is still predominantly a travelling wave solution resembling those shown in Fig. 2(b) and (c) but with a slight modulation as the wave passes over the stimulus; the modulated solution will be shown later. Here we highlight a qualitatively different type of travelling wave solution that can be found close to the Hopf curve H TW2 , whereby the wave has been pinned to the stimulated direction, as a so-called slosher state [START_REF] Folias | Nonlinear analysis of breathing pulses in a synaptically coupled neural network[END_REF][START_REF] Ermentrout | Spatiotemporal pattern formation in neural fields with linear adaptation[END_REF]; see panel (c). Furthermore, a new elongated region in parameter space has opened between H SW1 and the coinciding curves H SW2 and H TW1 , in which the stable behaviour is a standing wave, with one of its peaks aligned to the stimulus at v =0 . In order to describe in detail the changes to bifurcation structure that occur when the stimulus is introduced we now consider several one-parameter slices in λ at fixed values of k α taken from the diagrams already shown in Figs. 2(a) and 3(a); these slices are indicated by horizontal lines.

Figure 4 shows one-parameter bifurcation diagrams with zero adaptation gain k α =0 ,fi r s t with no input in panel (a) and with a small simple input (k I =0 .001) in panel (c). In order to best represent the solution branches we plot them in terms of the even, first-order mode of the solutions b p 1 (the cos(v)-component). Figure 4(a) shows that for small λ there is a single, stable solution branch with b p 1 =0; this corresponds to the flat (untuned) response shown in Fig. 2(d). When λ is increased beyond the pitchfork P 0 ,t h i sfl a ts t a t el o s e si t ss t a b i l i t ya n dar i n go f tuned responses are created. Figure 4(b) shows the profiles of the different solutions that exist at λ = 23. The dashed curve b2 is the unstable flat state, the tuned state b1 centered at v =0 corresponds to when b p 1 is largest and the tuned state b3 centered at v = ±180 corresponds to when b p 1 is smallest. Due to the presence of translational symmetry, intermediate states centered at any value of v also exist; discrete examples of these are shown as grey curves, but note these exist on a continuous ring filling in the direction space v. The easiest way to see the effect of introducing the stimulus and the fact that this breaks the translational symmetry is by studying the states that exist in the small input case also at λ = 23 shown in Fig. 4(d). Now the only stable solution is the tuned response d1 centered at v =0 ,t h e r ei sac o u n t e r p a r tu n s t a b l e solution centered at v = ±180 and all of the intermediate states have been destroyed. In the bifurcation diagram Fig. 4(c) the pitchfork bifurcation has been destroyed and there remain two disconnected solution branches. On the unstable branch the unstable solutions d2 and d3 are connected at a fold point F 0 ;t h i sb i f u r c a t i o ni st r a c e do u ta sF in Fig. 3(a). On the (upper) stable branch there is a smooth transition with increasing λ from a weakly-to a highly-tuned response centered at the stimulated direction. It is useful to detect where the increase in b p 1 is Inria .

. steepest as this signifies the transition to a tuned response. We denote this point b P 0 and as this is not strictly a bifurcation point we call it a pseudo pitchfork; it is still possible to trace out where this transition occurs in the (λ, k α )-plane and this is plotted as b P in Fig. 3(a). In order to best represent the solution branches we plot them in terms of the maximum of the sum of the even and odd first-order mode of the solutions max{ b p 1 + b p 2 } (the cos(v) and sin(v)-components). Figure 5(a) shows that two solution branches bifurcate simultaneously off the trivial branch at the twice-labelled point H tw1 H sw . In panel (b) we show one period of a stable travelling wave solution from the branch corresponding to H tw1 ;t h ew a v ec a nt a k ee i t h e r positive or negative (shown) direction in v. In panel (c) we show one period of an unstable standing wave solution from the branch corresponding to H sw1 ;t h es t a n d i n gw a v eo s c i l l a t e ss u c h that 180 • -out-of-phase (in v) maxima form alternatively. The phase in v of the entire waveform is arbtirary and, as for the pitchfork bifurcation, any translation of the whole waveform in v is also a solution. With the introduction of a stimulus, as for the pitchfrok bifurcation discussed earlier, the translational symmetry of the solutions is broken resulting in changes to the solution structure. The bifurcation point H sw shown in Fig. 5(a) splits into two bifurcation points H sw1 and H sw2 in panel (d) and solution profiles on the respective bifurcating branches are shown in panels (e) and (f). On the first branch, which is initially stable close to H sw1 ,o n eo ft h em a x i m a is centered on the stimulated direction v =0 • .F o rt h es e c o n du n s t a b l eb r a n c hb i f u r c a t i n gf r o m H sw2 the maxima are out of phase with the stimulated direction. The travelling wave branch is now a secondary bifurcation from the branch originating at H sw1 and the solutions now have a small modulation when the travelling wave passes over v =0 • (similar to the solution shown in panel (j)). In panel (h), at a lower value of k α ,t h es t e a d y -s t a t es o l u t i o nb r a n c hi si n i t i a l l y weakly tuned and forms a highly-tuned solution after b P 0 . The tuned response becomes unstable at H tw2 ;c l o s et ot h i sb i f u r c a t i o npo i n tt h es o l u t i o ni sp i n n e db yt h es t i m u l a t e dd i r e c t i o n ,a sw a s shown in Fig. 3(c), and as λ is increased further the amplitude in v of these oscillations about the stimulated direction increases, see panel Fig. 5(i). When the point PS is reached in panel (h) the oscillations become large enough such that there is a phase slip and beyond this point we obtain as t a n d a r dt r a v e l l i n gw a v es o l u t i o no n c em o r e ,s e ep a n e l( j ) . N o t et h a tt h et r a v e l l i n gw a v ei s still modulated as it passes over the stimulated direction. We reiterate the qualitative difference between the branches forming at H tw1 and H tw2 :f o rb r a n c h e so ft r a v e l l i n gw a v ef o r m i n gd i r e c t l y from an untuned or weakly tuned steady state as at H tw1 the solutions cannot be pinned to as t i m u l a t e dd i r e c t i o n ,h o w e v e r ,f o rb r a n c h e so ft r a v e l l i n gw a v ef o r m i n gd i r e c t l yf r o mat u n e d Inria steady state as at H tw2 the solutions are pinned to a stimulated direction close to the bifurcation. Furthermore, we note that, with the introduction of a stimulus a region of stable standing wave solutions, in phase with the stimulated direction, are introduced between H sw1 and H tw1 ;s e e branch segment through the point e in Fig. 5(d).

Explanatory one-parameter bifurcation diagrams

k I =0.001,k α =0.01 (a) (b) (c) (d) (e) (f) (h) (i) (j) b c ef i j H tw1 H sw H sw1 X X z H sw2 ⇠ ⇠ 9 H tw1 PS b P 0 H tw2 v v v v v v
The bifurcation analysis with the same input studied throughout this section, as shown in Fig. 1(c), continues in the next section. We go on to present the case k I =0.01 but in the context of a motion stimulus.

Competition model applied to the study of multistable motion

We now extend the general study presented in Sec. 3 and demonstrate how the same mo del can be used to study a specific neural biological phenomenon for which perceptual shifts are observed. We now associate the model's periodic feature space v 2 [-π, π) with motion direction. We assume that the model's activity in terms of time-evolving of firing rates p(v, t) are responses of direction-selective neurons in the middle temporal (MT) visual area. Indeed, MT is characterised by direction-selective neurons that are organised in a columnar fashion [START_REF] Diogo | Electrophysiological imaging of functional architecture in the cortical middle temporal visual area of cebus apella monkey[END_REF]. Here we only consider a feature space of motion direction and, thus, we assume the model responses to be averaged across physical (cortical) space. The chosen connectivity function Eq. ( 4) shown in Fig. 1(b) represents mutual inhibition between sub-populations of neurons associated with competing directions; this type of connectivity naturally gives rise to winnertakes-all responses tuned to one specific direction; there is evidence that competing percepts have mutually inhibitory representations in MT [START_REF] Logothetis | Neuronal correlates of subjective visual perception[END_REF][START_REF] Leopold | Activity changes in early visual cortex reflect monkeys' percepts during binocular rivalry[END_REF]). We use the models tuned response to dynamically simulate the mechanisms driving perception; cortical responses of MT have been linked specifically to perception of motion [START_REF] Britten | The middle temporal area: Motion processing and the link to perception[END_REF][START_REF] Serences | The representation of behavioral choice for motion in human visual cortex[END_REF]. We assume that over time any particular tuned response will slowly be inhibited as represented by the linear spike-frequency adaptation mechanism in the model. Furthermore, we assume there to be a fixed-amplitude stochastic fluctuation in the membrane potential that is modelled by additive noise (note that the noise is only introduced for the simulations presented in Sec. 5). We use as a model input pre-processed direction signals in the form expected from V1 [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]. In Sec. 4.2 the model's response properties in terms of its contrast dependence and direction tuning properties will be matched to what is known about the direction selective behaviour of MT neurons from physiological studies [START_REF] Albright | Direction and orientation selectivity of neurons in visual area MT of the macaque[END_REF][START_REF] Sclar | Coding of image contrast in central visual pathways of the macaque monkey[END_REF][START_REF] Diogo | Electrophysiological imaging of functional architecture in the cortical middle temporal visual area of cebus apella monkey[END_REF].

Definition of motion stimuli

We introduce two classical psychophysical stimuli where a luminance grating drifting diagonally (up and to the right in the example shown) is viewed through an aperture see Figs. 6(a) and (b). In the first case, with a circular aperture, the grating is consistently perceived as moving in the diagonal direction D (v =0 • ). The classical barberpole illusion [START_REF] Hildreth | The detection of intensity changes by computer and biological vision systems[END_REF], Chapter 4) comes about as a result of the aperture problem [START_REF] Wallach | Üb er visuell wahrgenommene Bewegungsrichtung[END_REF][START_REF] Wuerger | On the visually p erceived direction of motion" by Hans Wallach: 60 years later[END_REF], a diagonally drifting grating viewed through an elongated rectangular aperture is perceived as drifting in the direction of the long edge of the aperture. In the second case, with a square aperture, the stimulus has been shown to be multistable for short presentations on the order of 2-3s, where the dominant percepts are vertical [START_REF] Castet | The extrinsic/intrinsic classification of two-dimensional motion signals with barber-pole stimuli[END_REF][START_REF] Fisher | The directional tuning of the barber-pole illusion[END_REF]. We denote this stimulus the multistable barberpole and it has been the subject of complementary psychophysical experiments (Meso et al, 2012b) from which some results will be presented in Sec. 5.2. Complex model input: It has been shown in [START_REF] Barthélemy | Dynamics of distributed 1D and 2D motion representations for short-latency ocular following[END_REF] the motion signals from 1D cues stimulate a broad range of directions when compared with 2D cues that stimulate a more localised range of directions; see arrows in Figs. 6(a) and (b). Based on these properties, it is proposed that the multistability for the square-aperture stimulus is primarily generated by competition between ambiguous 1D motion direction signals along grating contours on the interior of the aperture and more directionally specific 2D signals at the terminators along the aperture edges. We represent the 1D cues by a Gaussian bump I 1D (v)=exp(-v 2 /2σ 2 1D ) with σ 1D = 18 • centred at v =0 • as shown in Fig. 6(c) (which is the same as Fig. 1(c)); on its own we call this a simple input that represents a drifting grating either filling the visual field (without aperture) or with an aperture that has no net effect on perceived direction such as the circular one shown in Fig. 6(a). We represent the 2D cues by two Gaussian bumps 6(d). Note that the functions I 1D and I 2D are normalised such that their maxima are 1 (not their areas). Figure 6(e) shows the complex input I ext represented as a summation of 1D and 2D motion signals with maximum normalised to 1 and a smaller weighting w 1D 2 [0, 1] given to 1D cues:

V (v = 45 • ), horizontal H (v = -45 • )andD(v =0 • ) (Castet
D ( v = 0• ) D ( v = 0• ) H(v = -45 • ) V(v =45 • ) ? 6 w 1D (a) (b) (c) (d) (e) v v v I 1D (v) I 2D (v) I ext (v)
I 2D (v)=exp(-v 2 /2σ 2 2D ) centred at v = 45 • and v = -45 • with width σ 2D =6 • as shown in Fig.
I ext (v)=w 1D I 1D (v)+I 2D (v -45) + I 2D (v + 45). (7)
The weighting w 1D translates the fact that in motion integration experiments 2D cues play a more significant role that 1D cues in driving perceived direction of motion [START_REF] Masson | Temporal dynamics of motion integration for the initiation of tracking eye movements at ultra-short latencies[END_REF];

Inria [START_REF] Barthélemy | Temporal dynamics of 2D motion integration for ocular following in macaque monkeys[END_REF]. Here we represent this weighting in a simple linear relationship, but in future studies it may be relevant to consider the contrast response functions for 1D and 2D cues separately [START_REF] Barthélemy | Dynamics of distributed 1D and 2D motion representations for short-latency ocular following[END_REF]. In the dark grey region containing the point d2 there is a steady, high-activity, tuned response with tuning width δ;seeuppercurv einpanel(d). W edefineδ as the width at half-height of the tuned response. Again, the boundary between these regions is demarcated by b P in Fig. 3(a). In the region containing the point b there is still a standing-wave-type solution, but it has been modified by the input such that there are oscillations between a tuned and an untuned state over time; see Fig. 7(b). In this context there is not a physiological interpretation for this solution, and so we will ensure that the model is operated in a parameter region where it cannot be observed. In the region containing the point c there is still a periodic response with small-amplitude oscillations in v of a tuned response about v =0; see Fig. 7(c). This is a travelling-wave-type slosher solution as described in Sec. 3.2 that is pinned by the input at v =0 • ;n o t et h et r a n s i e n tt h a tm a k e so n e full excursion before being pinned. Closer to the bifurcation curve the solutions are immediately pinned and further away a phase slip can be encountered as shown in Fig. 5(h)-(j). For a simple (unambiguous) input we again operate the model away from this oscillatory region of parameter space but find that when the complex input is introduced it is the slosher-type solutions that produce the desired switching behaviour.

In order to produce results and predictions that can be related directly to the experiments, where contrast was one of the main parameter investigated, contrast should also be represented as a parameter in the model. We argue that motion signals arriving in MT, primarily from V1, are normalised by shifts in the sigmoidal nonlinearity [START_REF] Carandini | Normalization as a canonical neural computation[END_REF] and, therefore we should not vary the input gain k I in Eq. ( 1) with respect to contrast. Accordingly, we found that when the input gain k I is increased, the model ceased to produce switching behaviour. However, by making the slope parameter λ depend on the contrast c 2 [0, 1],w ea r ea b l et or e p r od u c et h e observed switching behaviour.

We now fix k α =0 .01 such that we operating the model away from the oscillatory regions shown in Fig. 7 and describe how the model can be reparametrised in terms of contrast c.F o r some steady state p,w ed e fi n efi r i n gr a t er e s p o n s eR = max{p}-M as the peak firing rate response above some baseline value M ; max{p} is shown as a dashed line for solutions d1 and d2 in Fig. 7(d) and we set M = max{p d1 }.A sd i s c u s s e di nm o r ed e t a i li nt h eA p p e n d i xB ,t h e solution d1 at λ = 13 is consistent with an MT response to a very low contrast input (c<0.01), whereas the solutions d2 at λ = 25 is consistent with a high contrast input (c>0.2). By making λ aspecificfunctionofc we are able to match the model's contrast response to known behaviour for MT neurons. As shown in Fig. 7(e), we match the model's response to an appropriately parametrised Naka-Rushton function, which was used to fit contrast response data across several stages of the visual pathway including MT in [START_REF] Sclar | Coding of image contrast in central visual pathways of the macaque monkey[END_REF]; again, refer to Appendix B for further details. The operating range for the model is indicated by a horizontal line at k α =0.01 for λ 2 [13, 25] in Fig. 7(a). In Appendix B we also show that the tuning widths δ of the model responses are in agreement with the literature [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]. In the following sections an adaptation timescale of τ α = 16.5s is used, which is consistent with reported values from physiological experiments [START_REF] Descalzo | Slow adaptation in fast-spiking neurons of visual cortex[END_REF]. Up until now the results presented have been carried out with the adaptation timescale arbitrarily fixed at τ α = 100ms. The change of timescale does not qualitatively change the bifurcation diagrams shown previously or in the present section. However, it does change the oscillations produced such that they reproduce the type of switching observed in the experiments that will be discussed in Sec. 5.2; indeed, the specific value of τ α was chosen to match experimental data. In order for the bifurcation results described in earlier sections to be relevant it is desirable that we work with as m a l l -a m p l i t u d ea d d i t i v en o i s ea sg o v e r n e db yk X . The single source of noise in the model evolves with its timescale τ X set equal to τ α . This choice was found to have a pronounced affect on the switching dynamics without the need for large values of k X .N o t et h a tt h et i m eu n i t s displayed in figures up to this point are ms, but will be s in the remainder of the paper.

Figure 8 shows the bifurcation diagram for a complex input and the different types of behaviour that are observed in the operating range of the model as defined in the previous section. In the presence of the complex input, we see the same four regions found for the simple input; compare Fig. 7(a) and Fig. 8(a). The top-right-most region of the (λ, k α )-plane in which oscillations about v =0 • are observed has grown significantly. As for the simple input, there is a weakly tuned steady-state response in the region containing the point b and there is a tuned steady-state response in the region containing the point c; see panels (b) and (c). However, the region contaning the point d now shows an altered oscillatory behaviour. We see a model response that is initially centred at the direction D but after 2-3sshiftstoHandproceedstomak e regular switches between H and V, see panel (d). The model's separation of timescales is now seen more clearly; the model spends prolonged periods at H or V during which the adaptation builds up and eventually induces a switch to the opposite state; with τ α = 16.5s switches occur every ⇡ 3s, but the transition itself takes only ⇡ 50ms. Due to the dynamics being deterministic in the absence of noise (k X =0), the switches occur at regular intervals.

As described in the previous section we fix the operating region of the model with k α =0.01 and for λ 2 [13, 25] as indicated by the horizontal line in Fig. 8(a). As λ is increased from λ = 13 at b to λ = 25 at d (equivalently contrast increases from c =0to c>0.2)t h e r ea r et r a n s i t i o n s from a weakly tuned response (b) to a tuned response (c) to an oscillatory response (d). For c>0.2 the model response saturates as shown in Fig. 7(e). We now incorporate another known aspect of contrast dependence in motion processing by varying the relative weighting between 1D and 2D cues in the input. The psychophysics experiments presented in [START_REF] Lorenceau | The influence of terminators on motion integration across space[END_REF]; [START_REF] Lorenceau | Different motion sensitive units are involved in recovering the direction of moving lines[END_REF] show that 1D cues (contour signals) play an important role in motion perception at low contrast that diminishes with increasing contrast. As contrast increases the 2D cues (terminator signals) play a more significant role. Based on these studies we propose that for the complex model input (7) the relative weighting of 1D cues should decrease linearly with contrast

w 1D = W 0 -W 1 c, (8)
where W 0 =0 .5 and W 1 =1 .1. These specific values were chosen in order to match the experiments; see further comments in Sec. 5.3.

Figure 9 shows a one-parameter bifurcation diagram for the model working in the operating regime shown in Figs. 7(a) and 8(a) but now reparameterised in terms of contrast c as described above. At low contrast there is a stable steady-state response tuned to the direction D. The peak response max{p} increases with contrast. This steady-state response loses stability at at r a v e l l i n g -w a v eH o p fi n s t a b i l i t yH tw beyond which there is a stable oscillatory branch. The dependence of w 1D on the contrast affects the solutions in the two following ways. Firstly, the unstable branch associated with the D direction decreases in max{p} at large contrasts; see dashed curve in panel (a). Secondly, the period and amplitude of the oscillations in v does not saturate but continues to increase with contrast as shown in the inset (a1) and panel (b). We also note that close to the bifurcation point H tw there are long transients before the onset of oscillations, see panel (c), and that further from the bifurcation point the onset of oscillations is faster, see panel (d).

5 Comparison of model with experimental results

Experimental results

Figure 10 shows a summary of experimental data obtained in psychophysics experiments using the complex stimulus shown in Fig. 6(b) Meso et al (2012b). Data was recorded from four subjects for 15sp r e s e n t a t i o n so ft h es t i m u l u s . E y em o v e m e n t sw e r er e c o r d e da n dc o n t i n u o u s smooth trajectories estimated after removing blinks, saccades (very fast abrupt movements) and applying a temporal low pass filter. An SR Eyelink 1000 video recorder was used for the eye movement recordings and psychophysics stimuli were presented on a CRT monitor through a Mac computer using Psychtoolbox version 3.0 running stimulus generation routines written in Matlab. For the specific data shown, four healthy volunteers who provided their informed consent were participants, of whom two were naive to the hypothesis being tested. All experiments were carried out with, and following CNRS ethical approval. The presented stimuli covered 10 degrees of visual angle (the size of the side of the square in Fig. 6(b)) and were presented at a distance of 57cm from the monitor. Each task was done over 8 blocks of up 15 minutes over which 36 trials spanning a range of six contrasts were randomly presented each time. In this paradigm, recorded forced choice decisions indicating shifts in perceived direction through the three directions H, D and V and the estimated eye directions were found to be coupled and both indicative of perceived direction. Further details of these experiments can be found in our previous presentation (Meso et al 2012b) and a full description will appear in the experimental counterpart of this manuscript.

The temporal resolution of the eye traces is much higher than that of the reported transitions and allows for a relatively continuous representation of eye movement direction that can be compared with model simulations. Figure 10(a) and (b) show, for two different subjects, time traces of the time-integrated directional average of eye-movements from a single experimental trial at c =0 .08. Switches in perception can be computed from these trajectories by imposing thresholds for the different percepts. Both trials show that the directions H and V are held for extended durations and regular switches occur between these two states. The switches involve sharp transitions through the diagonal direction D. The diagonal direction can be held for extended durations immediately after presentation onset. However, we note that the eye-movement direction during the first 1so fp r e s e n t a t i o nh a sam o r el i m i t e dh i s t o r yi ni t st e m p o r a lfi l t e r i n g . Short presentations of the same stimulus were investigated in a related set of experiments (Meso et al, 2012a) and modelling work (Rankin et al, 2012).

Figure 10 also shows the relationship between the averaged rate of switches between H and Vo v e rar a n g eo fc o n t r a s tv a l u e sc 2{ 0.03, 0.05, 0.08, 0.1, 0.15, 0.2};i np a n e l( c )t h ed a t ai s averaged across the four subjects and in panel (d) it is separated out by subject. The lowest contrast shown c =0 .03 corresponds to the smallest contrast value for which subjects were able to reliably report a direction of motion for the stimulus. For the grouped data, at low contrast (c<0.1) the rate of switching increases with contrast with the rate being maximal at approximately c ⇡ 0.1.B e y o n dt h ep e a k ,f o rc o n t r a s t sc>0.1, the rate of switching decreases with contrast. For the data separated by subject shown in panel (d), the subjects MK and AB have a peak rate around 3.5 switches per 15sp r e s e n t a t i o na n dt h ep e a ko c c u r sa tc ⇡ 0.1.F o r Inria subjects JR and AM the peak rate is lower at around 2.5 switches per 15spresen tationandthere is a less prominent peak occurring at a higher contrast value c>0.1.H o w e v e r ,t h ec o m m o n pattern reveals two qualitatively different regimes with respect to changing contrast. A low contrast regime for which the switching rate increases with contrast and a high contrast regime for which the switching rate decreases with contrast.

Model simulations with noise (k X =0.0025)

We now study the dynamics of the mo del in the presence of additive noise in the main neural field equation. Recall that the stochastic process in the model is operating on the same slow timescale τ α as the adaptation and that the strength of the noise is k X =0 .0025.T w oc a s e s will be studied, first the low contrast case at c =0 .04,c l o s et ot h ec o n t r a s tt h r e s h o l do nt h e steep part of the model's contrast response; see Fig. 7(e). Second, the high contrast case at c =0 .08, which is above the contrast threshold on the saturated part of the contrast response function. In the first case, noise is introduced in a parameter regime where the model is close to bifurcation and oscillations only occur after a long transient, see Fig. 9(c). When operating in an e a r b yp a r a m e t e rr e g i m ec l o s et ob i f u r c a t i o nt h en o i s ec a u s e sr a n d o md e v i a t i o n sa w a yf r o m the direction D and can drive the model into an oscillatory state more quickly. In the second case, noise is introduced in a parameter regime where the model produces an oscillatory response with a short transient behaviour, see Fig. 9(d). In this regime the noise perturbs the regular oscillations either shortening or prolonging the time spent close to H and V. In processing this simulated data we assume that, initially the activity is centred around the direction D, and after some transient period switching will occur primarily between H and V. In order to detect switches between the directions H and V a so-called perceptual threshold (PT)hasbeensetatv = ±10 • . The first switch from D to either H or V is detected the first time the corresponding threshold is crossed. Subsequent switches are only detected the next time the opposite threshold is crossed. Note that although other algorithms could be employed to detect these switches, we found that these do not have had a great effect on the presented results. Across all the examples shown in Fig. 11, the average direction v oscillates in a random fashion and as time progresses the amplitude of these oscillations grows in v.F o rt h ec a s ec =0.04 there is a long transient and the first switch occurs for approximately t 2 [5s, 10s].F o rt h ec a s ec =0.08 the overall amplitude of the oscillations is larger and the first switch occurs for t<3s. Note also that the level of activity shown as an intensity in Fig. 11 is higher in the c =0 .08 case. An important difference between the two contrast cases is that in the low contrast case, the transitions between H and V occur gradually when compared with the abrupt transitions in the high contrast case. This suggests that at low contrast the direction D could be seen during the transitions, where as in the high contrast case the switches occur directly from H to V.

With respect to the experimental data, the model consistently reproduces the characteristic behaviour of regular switches between the H and V. Furthermore, the sharp transitions through the diagonal direction D are also captured well by the model. Compare the second row of Fig. 11 with the two examples shown in Fig. 10

Dependence of switching rate on contrast

Figure 12 shows the relationship between contrast and switching rate as computed with the model where the rate is expressed as the mean number of switches per 15ss i m u l a t i o n . P a n e l (a) shows the relationship without noise (k X =0 ) and with noise (k X =0 .0025). We show the average switching rate at discrete contrasts c 2 [0.02, 0.25] and at each contrast value we plot the switching rate averaged across 500 model simulations.

The deterministic case can be explained in terms of the bifurcation diagram shown in Fig. 9. At low contrast, no switching behaviour is observed as the model can only produce a steady-state response weakly tuned to the direction D. With increasing contrast, the onset of switching is abrupt, occurring just above c =0 .04 after the bifurcation H tw at c ⇡ 0.03. Switching does not begin immediately at the bifurcation point, due to long transients for values of c nearby, see Fig. 9(c). The switching rate remains constant at around 3 switches per 15si n t e r v a l ,a n ds t a r t s to drop off for contrasts c>0.12. The reduction in switching rate for larger contrasts is due to the increasing period of the oscillations as shown in Fig. 9(a1).

With the introduction of noise, there is an overall increase in switching rate, however, for larger contrasts the increase is minimal. At low contrasts, when the model is operating close to bifurcation, the noise has a more pronounced effect on the dynamics. As the contrast increases from c =0.02 there is a smooth increase in the rate of switching, which peaks after the bifurcation point before starting to decrease as in the no-noise case. Mechanistically these results can be interpreted as follows: at larger contrasts the switching rate is governed by an underlying adaptation driven oscillation. When noise is introduced, the individual switching times are randomly distributed as will be discussed in Sec. 5.4, but the average rate is not affected as can be seen by comparing the two cases shown in Fig. 12(a). At lower contrasts the presence of noise alone can drive a deviation from the diagonal direction leading to a switch. However, because the model is operating close to bifurcation the noise can also serve to shorten the transient period before the onset of adaptation-driven switching. Note that the intensity of noise is fixed across all contrasts, it is only at low contrasts that this has a large effect on the underlying dynamics.

The model results with noise are able to accurately capture the two contrast regimes from the experimental data. That is to say, an increase in switching rate at low contrasts with peak and subsequent decrease in switching rate at higher contrasts, compare Fig. 12(a) black curve with Fig. 10(a). The values of W 0 , W 1 and τ α were chosen in order to fit the experimental data, however, the two contrast regimes are robustly produced by the model independent of the specific values chosen. In Fig. 12(b) we show how, in the model, the relationship between switching rate and contrast changes with respect to PT.W h e nPT is low the peak switching rate is highest and occurs at a low contrast value. As PT is increased, the peak rate decreases and also occurs at a higher contrast value; the relationship also appears to flatten out for larger PT. Figure 10(b) shows the reported switching rate curves from the experiments, separated out by individual subject. The data shows a range of peak switching rate between the subjects. For the two subjects with the highest switching rate (MK,AB), the prominent peak occurs at c ⇡ 0.1. For the other two sub jects (JR,AM), the p eak rate is lower, the resp onse is flatter and the p eak rate occurs at a larger value of c.W ec o n c l u d et h a td i ff e r e n c e si np e r c e p t u a lt h r e s h o l db e t w e e n subjects can account for inter-subject differences.

Distribution of switching times

In the previous section we showed example model outputs for which switches between the directions H and V are detected. We found that the times between these switches vary and that, particularly in the low contrast case, the early transient behaviour can be very different from one simulation to the next. In order to investigate the distribution of the switching times we ran 1, 500 model simulations each of 15sa n df o r m e dad a t as e tb ye x t r a c t i n gt h et i m e sbe t w e e n consecutive switches from each simulation.

Figure 13 shows histograms of the computed switching times t sw . In the low contrast case approximately 1, 483 switches were recorded with mean time tsw =3.73sa n dS D =2.89 (Coefficient of Variance COV=0 .56)a n di nt h eh i g hc o n t r a s tc a s e3, 154 switches were reported with mean time tsw =4 .07sa n dS D =2 .08 (COV=0 .51). Although the mean of t sw is smaller in the high contrast case, more switches are detected because there is a shorter transient period before switching begins; the average time to the first switch in the low contrast case is 7.47s (SD=2 .87) compared with an average time of 3.02s( S D =1 .56) in the high contrast case. The aim now is to determine from which distribution the model data could have arisen. We follow the method presented in [START_REF] Shpiro | Balance between noise and adaptation in competition models of perceptual bistability[END_REF] and compare the model data with a Weibull probability distribution function (pdf), a gamma pdf and a log-normal pdf each with parameters Inria chosen using a standard maximum likelihood estimate. By inspection, it appears that the data in the low contrast case are well fitted by a log-normal distribution and that the data in the high contrast case are well fitted by a gamma distribution. In order to confirm this we perform a Kolmogorov-Smirnov goodness-of-fit test. In the low contrast case the log-normal distribution provides best fit (P =0.087)butthegammaandW eibulldistributionscanberejectedatthe5% significance level. In the high contrast case the gamma distribution provides best fit (P =0.635) and both the log-normal and Weibull distributions can be rejected.

Clearly, studying only the mean and standard deviation for the two different contrast cases does not reveal a significant difference. However, we do find a change in the underlying distributions governing the switching times, which is indicative of a change in the dominant mechanism driving the switching. Typically switching behaviour that is driven by adaptation over noise will have a lower peak that occurs later and a smaller spread with shorter tail as characterised by the gamma distribution [START_REF] Shpiro | Balance between noise and adaptation in competition models of perceptual bistability[END_REF], whereas, the peak is higher, earlier and the tail longer as characterised by log-normal when noise plays a more significant role. We also highlight the fact that the first switch occurs much earlier in the high-contrast case, this prediction could easily be tested experimentally.

Conclusions and discussion

The neural fields model at the center of this study has a periodic, continuous feature space and a linear spike frequency adaptation mechanism. In the presence of noise and a complex input for which multiple points in the feature space are stimulated this can be considered as an extension of classical rivalry models where competing states are modelled as discrete populations [START_REF] Laing | A spiking neuron model for binocular rivalry[END_REF][START_REF] Shpiro | Dynamical characteristics common to neuronal competition models[END_REF][START_REF] Moreno-Bote | Noise-induced alternations in an attractor network model of perceptual bistability[END_REF][START_REF] Shpiro | Balance between noise and adaptation in competition models of perceptual bistability[END_REF]. The initial theoretical exploration of the model without noise, and with a simple input provided a basis for the subsequent application to multistable motion perception. Spatially extended neural fields models with a linear implementation of spike-frequency adaptation have been studied both in ring models [START_REF] Hansel | Methods in neuronal modeling: From ions to networks[END_REF][START_REF] Curtu | Pattern formation in a network of excitatory and inhibitory cells with adaptation[END_REF][START_REF] Kilpatrick | Hallucinogen persisting perception disorder in neuronal networks with adaptation[END_REF][START_REF] Ermentrout | Spatiotemporal pattern formation in neural fields with linear adaptation[END_REF] and infinite spatial domains [START_REF] Pinto | Spatially structured activity in synaptically coupled neuronal networks: 1. traveling fronts and pulses[END_REF][START_REF] Ermentrout | Spatiotemporal pattern formation in neural fields with linear adaptation[END_REF]. A version of the model presented in this article, but without an input or noise was studied in [START_REF] Curtu | Pattern formation in a network of excitatory and inhibitory cells with adaptation[END_REF] where the existence of parameter regions with homogeneous, stationary tuned, travelling-wave and standing-wave responses was shown. In the absence of an input, these various different solution types are known to exist close to a so-called Bogdanov-Takens (BT) point in parameter space, which acts as a parametric organising center for different types of dynamics. In the presence of simple inputs it has been shown that new solution types can be produced such as breathers and pinned travelling-wave solutions [START_REF] Hansel | Methods in neuronal modeling: From ions to networks[END_REF][START_REF] Ermentrout | Spatiotemporal pattern formation in neural fields with linear adaptation[END_REF]. In this article we presented an in-depth numerical study of the complex organisation of these various solution types in three-dimensional parameter space. We were able to show that much of the structure local to the BT point is preserved with the introduction of a small, simple input, albeit in a subtly modified form. We give an account of the changes that occur in terms of the complicated series of bifurcations that delineate regions of parameter space exhibiting qualitatively different dynamics. It was found that close to a travelling-wave-type Hopf bifurcation solutions are pinned to the input. Furthermore, it was shown that for a standing-wave-type Hopf bifurcation giving rise to unstable solutions with no input, a region in parameter space with stable standing waves solutions was opened up when an input is introduced. Although great progress has been made analytically in the study of this class of model for simple inputs where a single location in feature space is stimulated [START_REF] Hansel | Methods in neuronal modeling: From ions to networks[END_REF][START_REF] Ermentrout | Spatiotemporal pattern formation in neural fields with linear adaptation[END_REF], the question of more complex inputs with stimulation of multiple locations provides a challenge. The advantage of the numerical approach used here is that we can directly extend earlier results when a complex input is introduced. In a recent study, perceptual multistability has been investigated in a model with synaptic depression and a two-location stimulus in a continuous feature space [START_REF] Kilpatrick | Short term synaptic depression improves information processing in perceptual multistability[END_REF].

We subsequently investigated p erceptual multistability for a stimulus that is multistable in terms of its perceived direction of motion and that has been the subject of recent psychophysics experiments, of which a summary was presented. This specific application allows for the model's continuous feature space to be exploited; it allows for a truly dynamic consideration of perceived direction, which unlike binocular rivalry or ambiguous shapes, is known to be neurally represented on a continuous scale. We study the multistable barber pole, which consists of a diagonally drifting grating viewed through a square aperture. The stimulus is known to be multistable between the diagonal grating direction D and the horizontal H and vertical V aperture-edge directions. The characteristic perceptual response is dominated by D immediately after onset followed by regular switches between horizontal H and vertical V directions. In the model the complex multistable barberpole can be represented by three-bumps in the feature space of motion direction based on experimental insights about the different ways in which 1D and 2D motion cues are processed by the visual system. The simple input is used to tune model parameters and introduce a contrast parameter such that its behaviour matches the known contrast response properties from [START_REF] Sclar | Coding of image contrast in central visual pathways of the macaque monkey[END_REF]. It is found that for a fixed adaptation strength, we are able to select a range of the nonlinear slope parameter such that the model's activity response can be matched to the qualitative and quantitative behaviour close to contrast threshold observed in physiological experiments. Once appropriately parametrised for the simple input, we find that for a complex input the model produces behaviour that is consistent with the characteristic perceptual responses described above.

We investigate the relationship between contrast and the switching behaviour; in the experiments two different regimes can be identified, the first at low contrast where the switching rate is increasing and a second at higher contrasts where the rate is decreasing. In the model we study a low contrast regime where the model is operating close to bifurcation and a high contrast regime above the contrast threshold. For both regimes we find common features in the switching behaviour produced by the model. Initially the percept D is dominant, but after some delay there is a shift to either H or V, typically within the first 1-8s( t h i sbe h a v i o u ri sc o n s i s t e n t with existing studies [START_REF] Castet | The extrinsic/intrinsic classification of two-dimensional motion signals with barber-pole stimuli[END_REF]; [START_REF] Fisher | The directional tuning of the barber-pole illusion[END_REF]; Rankin et al (2012)), after Inria which regular switching occurs between H and V. In the high contrast case this regular switching starts earlier, which we would expect as the 2D cues associated with the aperture edges and H/V directions should be stronger with increased contrast. We also find that typically the transitions between H and V are relatively smooth, passing gradually through the direction D in the low contrast case when compared with the sharper transitions in the high contrast case. Furthermore, by studying the dynamics either with or without noise, we find that at high contrasts the mean rate of switching is governed by the adaptation-driven oscillations. Although the noise produces random deviations in these switching times, the mean rate is unaffected. However, at low contrast, where the model is operating close to bifurcation the noise has a larger effect on the dynamics. At low contrast, what would be a simple threshold effect where switching turns on at some critical value becomes a smooth increase in the rate of oscillations. We conclude that at low contrast the increasing regime is associated with oscillations resulting from noise driving the model into an oscillatory regime. We further quantify these differences by showing that the distribution of switching times is best fitted to a log-normal distribution in the low contrast case and to a gamma distribution in the high contrast case. Typically switching behaviour that is driven by adaptation over noise will have a lower peak that occurs later and a smaller spread with shorter tail as characterised by the gamma distribution [START_REF] Shpiro | Balance between noise and adaptation in competition models of perceptual bistability[END_REF], whereas, the peak is higher, earlier and the tail longer as characterised by log-normal when noise plays a more significant role.

The general approach applied to a neural field model in this paper -making use of bifurcation methods for tuning parameters such that the model operates close to bifurcation whilst simultaneously matching known response properties from physiological studies -will allow for much broader studies of multistable perception. In particular, extensions to models that take into account physical space in conjunction with an abstracted feature space as used here would allow for the particular spatial properties of multistable visual stimuli to be taken into account. For example, more detailed neural fields mo dels taking into account the spatial integration of motion stimuli such as [START_REF] Tlapale | A neural field model for motion estimation[END_REF]) could be used to, within a single model architecture, investigate the underlying motion integration mechanisms that yield multistable perception across a broad range of stimuli, e.g., barberpoles, plaids, moving diamonds. Multistability has not been investigated in this kind of model to which the methods of numerical continuation and bifurcation analysis would be most applicable.

A The Ornstein-Uhlenbeck process

The noise X(v, t) that appears in Eq. ( 1) is a classical Ornstein-Uhlenbeck process, see e.g. [START_REF] Ermentrout | Foundations of Mathematical Neuroscience[END_REF], that is described by the following stochastic differential equation

τ α dX(v, t)=-X(v, t)dt + σdW (v, t), (9)
where W (v, t) is a feature uncorrelated Wiener or Brownian process. Note that the timescale τ α is the same as the one for the adaptation α in Eq. ( 2). The solution to Eq. ( 9) is readily found to be X(v, t)=e -t/τα X 0 (v)+σ Z t 0 e -t-s τα dW (v, s),t ≥ 0 where X 0 (v) is the initial noise distribution, assumed to be independent of the Brownian. Its mean is therefore given by hX(v, t)i = hX 0 (v)ie -t/τα , and its variance by Var(X(v, t)) = Var(X 0 (v))e -2t/τα + τ α σ 2 2 (1 -e -2t/τα ) It is seen that as soon as t becomes larger than the timescale τ α the mean becomes very close to 0 (it is even exactly equal to 0 for all times if the mean of the initial value is equal to 0), and the variance becomes very close to τασ 2 2 .B yc h o o s i n gσ 2 =2/τ α we ensure that Var(X(v, t)) is close to 1 as soon as t becomes larger than the timescale τ α .

B Tuning the model's contrast response

Here we describe the way in which the model is reparameterised such that its response matches known physiological behaviour. In [START_REF] Sclar | Coding of image contrast in central visual pathways of the macaque monkey[END_REF] contrast response functions were computed across several stages of the visual pathway for a drifting sinusoidal grating. At cortical layers, individual cell recordings were made with the stimulus moving in the cell's preferred direction. The Naka-Rushton function was used to fit contrast response data and the response R as a function of contrast c is given by

R(c)=R max c n c n + c n 50 + M, (10)
where R max is the maximal response c n 50 is the contrast at which the response reaches half its maximum value and n is the exponent that determines the steepness of the curve. For visual area MT in the macaque these parameters were estimated to take the following values: R max = 36, c n 50 =0.07 and n =3;theNak a-Rush tonfunctionisplottedforthesev aluesinFig.14(f)(dashed curve). Note that it is also necessary to adjust for the spontaneous firing rate M , which typically takes a value M 2 [0.05, 0.3]. We wish to reproduce the known contrast response properties for a drifting grating stimulus represented by the input in Fig. 14(a). For very low contrasts c<0.01,themodelshouldrespond Inria to a stimulus with a low level of activity that is weakly tuned or not tuned. Increasing c should result in a change in the model's response that is consistent with crossing the contrast threshold:

1. an increasing level of activty consistent with the Naka-Rusthon function Eq. ( 10), 2. a response that is tuned to the model's input with an appropriate tuning width (consistent with values reported in [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]).

We find that, for k α =0.01 and λ = 13,t h em od e l ' sw e a k l y -t u n e d ,s t e a d ys t a t er e s po n s eh a sa n appropriate level activity that we take as the value for M =0.18;seethelowercurveinFig.14(b). Furthermore, changing to λ = 25 we find that the model produces a tuned response that is consistent with a response above contrast threshold c>0.2 with a maximum level of activity that agrees with the value of R max =0 .52 -M =0 .36;s e eu p p e rc u r v ei nF i g .

Figure 1 :

 1 Figure 1: Model features and stimulus. (a): The smooth (infinitely differentiable) sigmoidal firing rate function S(x). (b): The convolutional kernel J is a three-mode approximation of a Mexican hat connectivity. (c): The simple stimulus is a Gaussian bump centered at v =0 • .

.

.

Figure 2 :

 2 Figure 2: Bifurcation diagram for the no-input case; summary of results from Curtu and Ermentrout (2004) in terms of stable behaviour. (a) Bifurcation curves plotted in the (λ, k α )-parameter plane demarcate regions with qualitatively different dynamics. The Hopf-type curves are H SW , H TW1 (coinciding) and H TW2 ,ap i t c h f o r kc u r v ei sP and these curves meet at the Bogdanov-Takens p oint BT .P a n e l s ( b ) -( e ) s h o w t h e a c t i v i t y p(v, t) indicated by intensity for model simulations at parameter values from the corresponding points b-e in panel (a).

Figure 3

 3 Figure3shows a new bifurcation diagram after the introduction of the simple input I 1D shown in Fig.1(c) with input gain k I =0.001.W ea r ei n t e r e s t e dt os e eh o wt h es o l u t i o n si d e n t i fi e di n the previous section change and how their organisation in parameter space has been modified. The most notable result is that much of the structure from the no-input case has been preserved, albeit with subtle changes that are now discussed. In the white region (to the left of H SW1 and b P ) there is now a low-activity response that is weakly tuned to the input centered at v =0;seepanel (d). In the dark-grey region there is still a steady-state, tuned response, but now centered on the stimulus at v =0. In the light-grey region the stable behaviour is still predominantly a travelling wave solution resembling those shown in Fig.2(b) and (c) but with a slight modulation as the wave passes over the stimulus; the modulated solution will be shown later. Here we highlight a qualitatively different type of travelling wave solution that can be found close to the Hopf curve H TW2 , whereby the wave has been pinned to the stimulated direction, as a so-called slosher state[START_REF] Folias | Nonlinear analysis of breathing pulses in a synaptically coupled neural network[END_REF][START_REF] Ermentrout | Spatiotemporal pattern formation in neural fields with linear adaptation[END_REF]; see panel (c). Furthermore, a new elongated region in parameter space has opened between H SW1 and the coinciding curves H SW2 and H TW1 , in which the stable behaviour is a standing wave, with one of its peaks aligned to the stimulus at v =0 .

Figure 3 :

 3 Figure 3: Bifurcation diagram for the simple input case with k I =0.001.( a )B i f u r c a t i o nc u r v e s plotted in the (λ, k α )-parameter plane demarcate regions with qualitatively different dynamics. The Hopf curves are H SW1 , H TW2 , along with the coinciding H SW2 and H TW1 ;af u r t h e rH o p fi n the light-grey region involves only unstable solutions and is not labelled; all the Hopf curves meet at a the double Hopf point DH .T w ob i f u r c a t i o nc u r v e sr e s u l t i n gf r o mt h es y m m e t r yb r e a k i n g of a pitchfork bifurcation are b P and F ;s e eF i g .4a n da c c o m p a n y i n gt e x t . P a n e l s( b ) -( e )s h o w the activity p(v, t) indicated by intensity for model simulations at parameter values from the corresponding points b-e in panel (a).

Figure 4 :

 4 Figure 4: Symmetry breaking of the pitchfork with introduction of a stimulus. (a) and (c) show bifurcation diagrams in λ for the no-input and small input cases, respectively; stable states are solid curves and unstable states are dashed curves. (b) and (d) show the solution profiles in v-space at the labelled points for λ = 23;s t a b l es t a t e sa r es o l i dc u r v e sa n du n s t a b l es t a t e sa r e dashed curves. In panel (b) the solid black curves correspond to the solution b1, for which b p 1 takes its largest value and the solution b3, for which b p 1 takes its smallest value; several intermediate solutions are plotted as grey curves, see text.

Figure 5

 5 Figure 5 shows one-parameter bifurcation diagrams in λ for three different cases: • No input cases with k α =0.03;seeFig.5(a)-(c);correspondstothehorizon tallinethrough the point labelled b in Fig. 2(a). • First small, simple input case with k α =0 .03;s e eF i g .5 ( d ) -( f ) ;c o r r e s p o n d st ot h eh o r izontal line through the point labelled b in Fig. 3(a). • Second small, simple input case with k α =0 .01;s e eF i g .5 ( h ) -( j ) ;c o r r e s p o n d st ot h e horizontal line through the point labelled c slice in Fig. 3(a).

  ,k α =0.03 k I =0.001,k α =0.03

Figure 5 :

 5 Figure 5: Changes to standing-and travelling-wave branches born in Hopf bifurcations with introduction of the stimulus. The first column shows one-parameter bifurcation diagrams in λ where black curves are steady-state branches and grey periodic branches; stable solution branches are solid and unstable branches are dashed. Hopf bifurcations to travelling waves are H tw1 and H tw2 ,a n dt os t a n d i n gw a v e sa r eH sw , H sw1 ,a n dH sw2 .S e c o n da n dt h i r dc o l u m n ss h o wo n ep e r i o d T of the solutions at corresponding points on solution branches from one-parameter diagrams.

Figure 6 :

 6 Figure 6: Simple and complex motion stimuli. (a): A drifting luminance grating viewed through ac i r c u l a rs o f ta p e r t u r e ;t h ed i a g o n a ld i r e c t i o no fm o t i o nDi sc o n s i s t e n t l yp e r c e i v e d . ( b ) : A drifting luminance grating viewed through a square aperture; the dominant percepts are vertical V, horizontal H and diagonal D. (c): Representation of the 1D motion signals in direction space; the simple motion stimulus (a) is equated with I 1D .( d ) :R e p r e s e n t a t i o no ft h e2 Dm o t i o ns i g n a l s in direction space. (e): Summation of the 1D and 2D motion signals with a weighting w 1D =0.5; the complex motion stimulus (b) is equated with I ext .

4. 2

 2 Figure7(a) shows the two-parameter bifurcation diagram in λ and k α for the same simple 1D input as in Fig.3(a), but with input gain increased by a factor of 10 to k I =0 .01. The diagram shows the same organisation of bifurcation curves, but the oscillatory regions have shifted significantly towards the top-right corner. Once again, in the white region containing the point d1 there is a steady, low-activity, weakly tuned response; see lower curve in panel (d). In the dark grey region containing the point d2 there is a steady, high-activity, tuned response with tuning width δ;seeuppercurv einpanel(d). W edefineδ as the width at half-height of the tuned response. Again, the boundary between these regions is demarcated by b P in Fig.3(a). In the region containing the point b there is still a standing-wave-type solution, but it has been modified by the input such that there are oscillations between a tuned and an untuned state over time; seeFig. 7(b). In this context there is not a physiological interpretation for this solution, and so we will ensure that the model is operated in a parameter region where it cannot be observed. In the region containing the point c there is still a periodic response with small-amplitude oscillations in v of a tuned response about v =0; see Fig.7(c). This is a travelling-wave-type slosher solution as described in Sec. 3.2 that is pinned by the input at v =0 • ;n o t et h et r a n s i e n tt h a tm a k e so n e full excursion before being pinned. Closer to the bifurcation curve the solutions are immediately pinned and further away a phase slip can be encountered as shown in Fig.5(h)-(j). For a simple (unambiguous) input we again operate the model away from this oscillatory region of parameter space but find that when the complex input is introduced it is the slosher-type solutions that produce the desired switching behaviour.In order to produce results and predictions that can be related directly to the experiments, where contrast was one of the main parameter investigated, contrast should also be represented as a parameter in the model. We argue that motion signals arriving in MT, primarily from V1, are normalised by shifts in the sigmoidal nonlinearity[START_REF] Carandini | Normalization as a canonical neural computation[END_REF] and, therefore we should not vary the input gain k I in Eq. (1) with respect to contrast. Accordingly, we found that when the input gain k I is increased, the model ceased to produce switching behaviour. However, by making the slope parameter λ depend on the contrast c 2 [0, 1],w ea r ea b l et or e p r od u c et h e observed switching behaviour.We now fix k α =0 .01 such that we operating the model away from the oscillatory regions shown in Fig.7and describe how the model can be reparametrised in terms of contrast c.F o r some steady state p,w ed e fi n efi r i n gr a t er e s p o n s eR = max{p}-M as the peak firing rate response above some baseline value M ; max{p} is shown as a dashed line for solutions d1 and d2 in Fig.7(d) and we set M = max{p d1 }.A sd i s c u s s e di nm o r ed e t a i li nt h eA p p e n d i xB ,t h e solution d1 at λ = 13 is consistent with an MT response to a very low contrast input (c<0.01), whereas the solutions d2 at λ = 25 is consistent with a high contrast input (c>0.2). By making λ aspecificfunctionofc we are able to match the model's contrast response to known behaviour for MT neurons. As shown in Fig.7(e), we match the model's response to an appropriately parametrised Naka-Rushton function, which was used to fit contrast response data across several stages of the visual pathway including MT in[START_REF] Sclar | Coding of image contrast in central visual pathways of the macaque monkey[END_REF]; again, refer to Appendix B for further details. The operating range for the model is indicated by a horizontal line at k α =0.01 for λ 2[13, 25] in Fig.7(a). In Appendix B we also show that the tuning widths δ of the model responses are in agreement with the literature[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].

Figure 7 :

 7 Figure 7: Bifurcation study and contrast response for simple input I 1D (see inset of (a)) with input gain k I =0 .01.( a ) :T w o -p a r a m e t e r b i f u r c a t i o n d i a g r a m i n t e r m s o f s i g m o i d s l o p e λ and adaptation strength k α shows qualitatively the same organisation of bifurcation curves as Fig. 3(a). (b),(c): Time-traces of the activity p(v, t) indicated by intensity as computed at the corresponding points b and c labelled in (a). (d): Steady-state responses in terms of the activity p(v) at the corresponding points d1 and d2 labelled in (a). (e): Contrast response in terms of normalised peak activity R for simple input (solid curve) fitted to a Naka-Rushton function (dashed curve) as described in Appendix B. The line between d1 and d2 in (a) is the operating range of the model. The response d1 shown in panel (d) corresponds to c =0and the response d2 in panel (d) corresponds to c =1.

Figure 8 :Figure 9 :

 89 Figure 8: Bifurcation study for complex input I ext (see inset of (a)) with input gain k I =0 .01. (a): Two-parameter bifurcation diagram in terms of sigmoid slope λ and adaptation strength k α shows qualitatively the same organisation of bifurcation curves as Fig. 3(a) and Fig. 7(a). Line between points labelled b and d in (a) show the operating region of the model as defined in Fig. 7(a). (b)-(d): Time-traces of the activity p(v, t) indicated by intensity as computed at the corresponding points b, c and d labelled in (a).

Figure 11

 11 Figure 11(a)-(d) shows 15stimetracesofthepopulationactivityp for the cases c =0.04 (first row) and c =0.08 (second row). Note that each individual model simulation is quite different due to the noise, but we have selected representative examples that allow us to highlight key features in the model responses and compare the different contrast cases.In processing this simulated data we assume that, initially the activity is centred around the direction D, and after some transient period switching will occur primarily between H and V. In order to detect switches between the directions H and V a so-called perceptual threshold (PT)hasbeensetatv = ±10 • . The first switch from D to either H or V is detected the first time the corresponding threshold is crossed. Subsequent switches are only detected the next time the opposite threshold is crossed. Note that although other algorithms could be employed to detect these switches, we found that these do not have had a great effect on the presented results.

  (a) and (b).

Figure 10 :Figure 11 :

 1011 Figure 10: Summary of results from psychophysics experiments for the complex stimulus shown in Fig. 6(b). (a), (b): Time traces of average direction from eye-movements during two individual stimulus presentations at c =0 .08.E r r o rb a r ss h o wt h es t a n d a r dd e v i a t i o no ft h ec o m p u t e d direction of smooth components over 200 samples; the re-sampled value at 5Hz is the mean. (c): Relation between contrast and mean switching rate in terms of perception (reported by subjects) and as computed from eye-movement traces; grouped data is averaged across the four subjects with standard error shown. (d): Switching-rate data (from perception) separated out by subject with standard error for each subject shown.

Figure 12 :

 12 Figure 12: Mean switching rates computed with the model and recorded from psychophysics experiments. (a): Switching rates computed with the model without noise and with noise. (b): Switching rate curves computed with the model for a range of PT values.

Figure 13 :

 13 Figure 13: Distribution of perceptual switching times and test distributions. Histograms show the distribution of switching times as computed from model simulations with PT = 15 • (see text for details). Candidate distributions are overlayed, where the shape and scale parameters have been chosen to best fit the model data. (a) Low contrast case in which the log-normal distribution provided a better fit. (b) High contrast case in which the gamma distribution provided a better fit.

Figure 14 :

 14 Figure 14: Model response for drifting grating stimulus and tuning of λ-dependence on c.( a ) Stimulus profile I 1D (v) for a visual-field-covering drifting grating. (b) Steady state responses at λ = 13 (grey curve) and at λ = 25 (black curve); dashed horizontal lines indicate the maximum activity for these steady states. The tuning bandwidth, or tuning width at half-height, is δ and the maximum level of activity above the spontaneous level (dashed grey line) is R.( c )T u n i n g width δ of response for λ 2 [λ min ,λ max ].( d )N o r m a l i s e dr e s p o n s eR for λ 2 [λ min ,λ max ].( e ) Variation of λ with respect to c. (f) Contrast response curve with λ varying as in (e) (black) and the Naka-Rushton function with parameters as described in the text (dashed grey).

  1 4 ( b ) . W ea l s o find that varying λ in the range λ 2[13, 25] that the tuning widths at half-height δ fall in the range δ 2 [80 • , 115 • ] (see Fig.14(c)), which is consistent with the tuning widths obtained by individually probing a population neurons in the studies[START_REF] Albright | Direction and orientation selectivity of neurons in visual area MT of the macaque[END_REF] (mean 90 • ,s t a n d a r d deviation SD= 29 • )a n d D i o g oe ta l( 2 0 0 3 )( m e a n104 • ,S D = 35 • ). Figure14(d)shows the model's response R plotted against λ 2 [13, 25] and we see that initially, close to λ min = 13,t h e response has an appropriate shape, but for larger λ the response does not saturate as required. It is therefore necessary to introduce a function λ(c) that saturates smoothly at the value λ max = 25. We reparameterise the mo del so that λ depends on the contrast using a function that increases linearly from λ = λ min and saturates for large c at λ = λ max .W ee ff e c t i v e l yu s et h et o ph a l fo f as i g m o i df u n c t i o nS(x)=1 /(1 + exp(-x)).A f t e rfi x i n gλ min = 13 and λ max = 25,i ts u ffi c e s to vary just the slope parameter µ in the following equation in order to find a good fit. The equation λ(c)=λ min +(λ max -λ min ) Fig. 14(e) with µ = 60.B yc o n s t r u c t i o n ,w eh a v es e td e p e n d e n c eo fλ on c in order to obtain a function of the model response R in terms of c that matches a Naka-Rushton function appropriately parameterised for MT; see Fig. 14(f) (solid curve).

Table 1 :

 1 Parameter values used in the numerical studies.

	.2. For the model simulations without noise (k X =0 ), as presented in Secs. 3 and 4,
	Inria

RR n°8220

Inria

Rankin & Meso & others

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

JR, OF and PK acknowledge funding from ERC grant No 227747 (NERVI) and from the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement no. 269921 (BrainScaleS) and by the EC IP project FP7- ICT-2011-8 no. 318723, (MatheMACS). AM was supported by the Aix-Marseille Université Foreign Post-doctoral Fellowship 2011 from 05/11-04/12.