Continuous Attractors with Morphed/Correlated Maps

Continuous attractor networks are used to model the storage and representation of analog quantities, such as position of a visual stimulus. The storage of multiple continuous attractors in the same network has previously been studied in the context of self-position coding. Several uncorrelated maps of environments are stored in the synaptic connections, and a position in a given environment is represented by a localized pattern of neural activity in the corresponding map, driven by a spatially tuned input. Here we analyze networks storing a pair of correlated maps, or a morph sequence between two uncorrelated maps. We find a novel state in which the network activity is simultaneously localized in both maps. In this state, a fixed cue presented to the network does not determine uniquely the location of the bump, i.e. the response is unreliable, with neurons not always responding when their preferred input is present. When the tuned input varies smoothly in time, the neuronal responses become reliable and selective for the environment: the subset of neurons responsive to a moving input in one map changes almost completely in the other map. This form of remapping is a non-trivial transformation between the tuned input to the network and the resulting tuning curves of the neurons. The new state of the network could be related to the formation of direction selectivity in one-dimensional environments and hippocampal remapping. The applicability of the model is not confined to self-position representations; we show an instance of the network solving a simple delayed discrimination task.


Introduction
The ability to keep an internal representation of a continuous variable in the absence of sensory stimuli, is a crucial requirement in order to succeed in what can be considered trivial day to day actions or experimenter designed tasks. For instance one may think about the eye position between successive saccades [1], the angle of stimulus presentation in an oculomotor delayed protocol [2], the spatial position or the head direction in a dark environment [3][4][5], or the phase of the recently discovered grid fields [6,7].
A widely used class of models for this kind of working memory is constituted by attractor neural networks. The temporary maintenance of an item in memory corresponds to a specific network pattern of activity which is stabilized via strengthened recurrent connections between the active neurons in the pattern [8][9][10][11]. These connections are usually imposed, or trained, as the outcome of some form of Hebbian learning. The attractor is called continuous when the stable states form a continuous manifold which can be parametrized by the state variables. This outcome is obtained under certain conditions on the synaptic connection, for example when the connections between neurons are lateralinhibition like (e.g. Mexican hat) [12][13][14]. The underlying idea is that each neuron is assigned a location on an abstract map. The synaptic weights (encoding) depend on the location of the pre-and post-synaptic neurons. By means of Turing instability, the network dynamics creates a localized pattern of activity (or bump) on the map [15]. The external input links the position on the map to the state variable, forming a representation.
Continuous attractors have been used to explain the maintenance of various analog quantities close or far from the primary sensory and motor regions. For instance, the orientation tuning in the visual cortex [16,17], hippocampal place fields in one [18,19] and two dimensions [20,21], eye position [1,22], head direction tuning in the postsubiculum [23,24] and entorhinal grid fields [25,26].
The simple picture of a single continuous attractor can be naturally extended to the case of multiple attractors. The encoded maps can then be assumed to be either uncorrelated or correlated, and in particular to exhibit some structure (e.g. deriving from a morphing procedure). Assuming a complete lack of correlations between maps is not realistic, though useful for obtaining analytic results [27]. In this contribution, we analyze the network representations arising from the storage of two maps, with a varying degree of correlation between them, and from the storage of a morph sequence between two uncorrelated maps. We are interested in finding the conditions under which the network representation can provide some information about the state variables. Surprisingly, even when the correlation between two maps is very high, under conditions which will be clarified later it is possible for the network to maintain separate representations of the state variables.

Multiple maps
Multiple state variables can be encoded in the same network. An example is offered by the place representations of several environments [20,21]. To each environment corresponds a neural map which is encoded in the synaptic efficacies. Sensory inputs would then select the correct representation, i.e. both the environment and the position in the environment. The selected map wins the competition with the other maps stored in the network, and a localized pattern appears. In this case the network only maintains information about one of the several encoded state variables.
A more peculiar property of multiple continuous attractors, is their ability to represent simultaneously the values of several state variables. This property was explored in [28], where two partially overlapping neural populations (representing discrete features), are assigned two uncorrelated maps. Another example is provided in the study of [29], where a single network stores and represents simultaneously a continuous and discrete attractors.
In principle, given the existence of multiple representations in different brain regions (either one per region, or many in one region), a brain area downstream would necessarily encode several state variables. In light of a Hebbian interpretation on how this encoding takes place, it seems natural to distinguish between two cases. When multiple representations provide a simultaneous input to a region, the result is probably encoded multiplicatively [29], or, in general, non-linearly. For inputs happening non concurrently, as for instance when walking through several rooms sequentially, an additive encoding of each room is expected [21]. In the following we will analyze additive encoding.

Correlations
The present contribution addresses the issue of encoding correlated maps. The motivations come from recent experimental results on place cells recording in morphed environments [30][31][32], where place fields remapping along a sequence of morphed arenas is experimentally tested, and from theoretical and experimental studies concerning the morphing of discrete attractors [33][34][35].
In general, we would consider the encoding of p manifolds X f , each of dimension d f , where f~1 . . . p. We will refer to a single manifold as a map, once a coordinate system x f is chosen. The use of uppercase (e.g. X ) or lowercase (e.g. x) will distinguish between the whole map and a single point on it respectively. Given a pre-synaptic neuron indexed by x f , and a post-synaptic x' f , the encoding of a single map is obtained using a synaptic matrix W f x f ,x' f ð Þ, and is such that a continuous attractor representation would arise if it were the only map. We assume, as mentioned above, that the complete encoding arises from a linear superposition of the p matrices, W~1 p X f W f . The statistical properties of the maps, and in particular the correlation between them, can be fully specified by providing the probability density n x f f g ð Þ. The general problem is too difficult to be studied analytically. Some results can be obtained for the case of uncorrelated maps on the same manifold [27], though the system can be explored by simulating the full microscopic networks (see e.g. [21] for the uncorrelated case and [36] for simulation results of the correlated case).
In order to simplify the analysis, while retaining the basic structure of the problem, we focus on the case of p~2 representations, on a 1-dimensional circular manifold (i.e. the ring model [16,37]). The correlation between the maps is constructed by limiting the distance between the single neuron locations on the two maps. We devise a simple method to generate a morph sequence between two uncorrelated maps, by linearly modifying the neurons locations between the original maps. This method also suggests a way to test the network response to the exposure of intermediate maps between the two stored correlated maps.
For concreteness, one could think about maps of two similar circular arenas, and reason in term of spatial coding. In this context, we are interested in clarifying how the information about the position in the current environment is represented by the network, when varying the constitutive parameters of the model; And how the representation changes when the network is exposed to environments along a morph sequence.
In the following we will describe with mean-field (MF) theory the attractor landscape of a network, i.e. the stable solutions in absence of any place specific input. We then consider the behavior of the solutions when a spatially tuned input is present. We will establish the approximate relationship between two strongly correlated maps and the encoding of a morph sequence between two reference rings, and study the behavior of the solutions in presence of a tuned input varying along the sequence. Finally we will verify the results with microscopic simulations of finite networks. The network properties can be tested experimentally to confirm (or falsify) the attractor hypothesis.

Let us consider two circular environments
The factor J 1 is a measure of the amplitude of the map specific interaction, while J 0 v0 ð Þis a uniform inhibitory term. This form of connectivity can be thought as arising from the first two terms

Author Summary
How is your position in an environment represented in the brain, and how does the representation distinguish between multiple environments? One of the proposed answers relies on continuous attractor neural networks. Consider the web page of your campus map as a network of pixels. Every pixel is a neuron, and nearby pixels excite each other, while distant pairs are inhibited. As a result of their interactions, a bunch of close-by pixels will light up, indicating your current position as suggested by your web-cam (the sensory input). When you travel to another campus, the common assumption holds that pixels are completely scrambled and the excitatory/inhibitory pattern of connections is summed to the existing one. Now these connections and the sensory input will activate the pixels corresponding to your location in the new campus. The active pixels will look like noise in the old map. But what if the campuses are similar, i.e. the pixels are not completely scrambled? We show that the network has a novel way of distinguishing between the environments, by lighting up distinct subsets of pixels for each campus. This emergent selectivity for the environment could be a mechanism underlying hippocampal remapping and directional selectivity of place cells in 1D environments.
of Fourier series of a more general coupling. The rate dynamics for the network activitym m h A ,h B ð Þis [38]: where we assumed a threshold-linear transfer function for the neurons, ½x z~x when xw0 and 0 otherwise. The external afferent current is denoted by I, and it is assumed uniform in the current Section. We build the maps with a simple procedure which induces a correlation between them.   Note that a change in the strength of the applied uniform input I produce no changes in the order parameters (see Methods -Reduced dynamics and Eqs. 5). Several examples of network activity (Eq. 3), corresponding to different representative choices of the order parameters, are shown in Fig. 2. The various panels show the network activity in the two two-dimensional maps H A ,H B ð Þ and H,R ð Þ, and the onedimensional projections of the activity to H A and H B . Note that not all the choices of order parameters corresponds to actual solutions of the dynamics (which are determined by the parameters J 1 ,J 0 ,m ð Þand the initial conditions), as will be shown later.
The meaning of the order parameters can be read out from Eq. 3. The variable a represents a scaling factor for the amplitude of the network activity, which in turn is proportional to the uniform input I.
The variable s is a measure of the spatial size of the activity profile, i.e. of the region in either H A ,H B ð Þor H,R ð Þin which the network activity m is strictly positive (Eq. 3). The activity profile is also referred to as a bump. For instance the case s~{1 would correspond to absence of activity (the current in the thresholdlinear transfer function would always be negative), while s~1 would make all the neurons in the network active.
The order parameter c tells us how much the network representation ''favors'' one of the two maps. By its definition (Methods -Reduced dynamics, Eq. 24), the possible range for c is ½{1,1. The two extreme cases c~+1 correspond to a network activity localized in either map H B or H A . For instance, the network activity Eq. 3 for c~{1 reads where in the last equality we used Eq. (1). From here we see that the position of the bump peak is located at y A :y z {y { ; the same derivation, with c~1, would give us y B :y z zy { . The network representation is in this case a bump of activity localized in one map, and does not have any spatial modulation in the other map, as exemplified in Fig. 2A,III. In the case c~0, from the explicit expression of the activity m we get This representation exhibits an equal amount of spatial modulation in both maps H A and H B , i.e. the solution represents equally the two stored maps, Fig The quantities y z and y { m identify respectively the location of the maximum of the network activity in the h,r ð Þ coordinates, which is uniquely mapped to the maximum in h A ,h B ð Þvia Eqs. 1. In the following, we will show that the network activity examples depicted in Fig. 2 are possible solutions of the dynamics described by Eq. 2. We refer to each of these classes as double ring ( Fig. 2A,C,) single ring (Fig. 2B) and cylinder (Fig. 2D), for reasons that will be clarified in the next Section. The cylinder class represents an interesting novel regime (simultaneous localized projections in both environments), and we will devote most of the paper to describe the properties of this particular solution.
Phase diagram of the model In this Section we analyze the fixed point solutions of the system, and heuristically describe the region of stability of these solutions. A more rigorous description of the stability can be found in Methods -Stability.
In Methods -Reduced dynamics we derive the dynamics of the order parameters from Eqs. 2. We report here the result where the function g h,r ð Þ is defined as i.e. the rescaled steady state activity profile Eq. 3. Note that y z can be eliminated from the right hand sides of the Eqs. 5, rotating the integration variable h. This is possible because there is no spatial dependence in the external input to the network. The first four equations in Eqs. 5 can then be solved independently of the fifth one, since the right hand sides do not depend on y z . We show in Methods -Solutions properties that, once we have the solution for the variables (c,a,s,y { ), the last equation reduces to _ y y z~0 . We can thus restrict the analysis to four out of five equations in Eqs. 5. The elimination of one angular degree of freedom is a consequence of the rotation invariant structure of the encoding, and is the hallmark of continuous attractors arising from spontaneous symmetry breaking. On the other hand, the integrals over r in Eqs. 5 are not over the whole circle and we cannot rotate y { away.

Homogeneous solution
Before analyzing the fixed point solutions of the system described by Eqs. 5, we briefly mention an uninteresting region in the parameters space which can be found also in the classical ring model. This region corresponds to the homogeneous solution, i.e. all the neurons in the network are active at a constant level, and can be obtained from Eq. 2. The expression corresponding to the line of separation in the plane J 1 ,m ð Þ between the homogeneous solution and the spatially localized bump (see where sinc x ð Þ~s in x ð Þ x . This result is obtained in Methods -Stability, see also below.

Fixed point equations for the localized activity state
Let us start by imposing y {~0 , a restriction that will be addressed later on. The first tree equations at steady state from Eqs. 5 become then equations for the three order parameters a,c,s ð Þ: The first two equations determine the shape of the bump c,s ð Þ. Given the map specific modulation in the coupling and the distance between the maps J 1 ,m ð Þ, we can derive from the first two equations the size of the bump s and the order parameter c, representing how close the network representations are to the stored environments A and B. The last equation gives us the amplitude of the network activity a, which also depends on the parameter J 0 .
As mentioned in Results -Phase diagram of the model, the order parameter y z can be chosen arbitrarily, due to the rotation invariance of the problem; for simplicity we choose y z~0 .

Amplitude instability
We deal first with the equation concerning the amplitude of the solution. Given that the activity can be rescaled by changing the value of the applied external current I, we are not interested in actually solving the equation. The only requirement is that a §0 in order for the solution to be meaningful, i.e. no negative amplitudes are allowed. This requirement translates to a constraint on the inhibition J 0 : We show with stability analysis (Methods -Stability) that the critical value J C 0 , obtained by choosing the equality in the previous expression, corresponds to the onset of amplitude instability; given a choice for the parameters J 1 ,m ð Þ, which specifies the bump shape c,s ð Þ, for values of the inhibition weaker than J C 0 the solution grows to infinity. This qualitative behavior was present also in the classical ring model. Fig. 3B shows the values of J C 0 as a function of J 1 for various choices of m. In order to stabilize the solutions, the inhibition must grow with increasing J 1 and decreasing m. Note that it is reasonable to consider the previously mentioned homogeneous solution as a bump with maximal size s~1 ð Þ. In this case the critical J 0 can be explicitly computed, and turns out to be J C 0~1 .

Single ring solution
Now we focus on the possible solution c~0. It is easy to see that when c~0, the second of Eqs. 8 is automatically satisfied due to the symmetry of the integrand in h (and r); This means that the solution c~0 exists everywhere in the parameter space.
The steady state activity Eq. 3 with c~0 (and y {~0 , our initial assumption) reads which corresponds to a packet of activity localized in the h coordinate, and modulated in r, see Fig. 2B for a plot of the activity profile. The remaining fixed point equation can be used to obtain s~s J 1 ,m ð Þ. We refer to the case c~0,y {~0 ð Þas a single ring solution; the ring is spanned by the freedom of choice in the angle y z . In this regime of activity the network is not able to represent separately the environments A and B, but only the middle environment described by h. Even though the solution exists everywhere, it is destabilized in some regions of the parameter space, as shown in the phase diagram ( Fig. 3A, Single ring region).
By looking at the maximal bump size s~1, we can expect to reproduce the curve separating the homogeneous solution from the single ring. Inserting s~1 in the first of Eqs. 8, it is possible in this case to compute explicitly the integral, which in fact yields Eq. 7.

Double ring solution
In order to find the region of existence of the solutions with c=0, we can solve numerically Eqs. 8 in the parameters plane J 1 ,m ð Þ. The result is shown in Fig. 4, where the color code represents c for a given choice of the parameters. It can be seen that there is only a narrow region of high m (low correlation) and low J 1 where such a solution exists.
It is important to note that the equations used to find c are invariant under the symmetry c?{c. This means that both solutions (+c ? ) representing map H A or H B are possible. The steady state activity profile in this case looks like: Given the freedom of choice for the phase y z , each of this solutions lives on a ring; we call the solution c=0,y {~0 ð Þ , double ring. An instance of the network activity in this regime is shown in Fig. 2C.
The curve separating representations preferring one of the two maps (c=0), and c~0, can be obtained by expanding the second of Eqs. 8 to first order in c: where H x ð Þ is the Heaviside step function, H xw0 ð Þ~1 and H xv0 ð Þ~0. Dividing by c, we get rid of the c~0 solution. By finding the zeros of the integral, we select the curve in the parameter space corresponding to the onset of existence of the double ring solution. This curve is shown in Fig. 4. We have found that the stability of the double ring solution coincides, empirically, with the region of existence of such solution (compare the phase diagram in Fig. 3A, Double ring region with Fig. 4).

Cylinder solution
Finally, we examine the meaning of the equation for y { , the order parameter linked to the location of the maximum of the bump in r. We have assumed y {~0 for simplicity, given that a rotation in the integrands in Eqs. 5 is in general not viable due to the restricted range of integration in r. Note though, that when the size of the bump s is small enough, it is possible to perform the rotation without affecting the value of the integrals; the only requirement is that the rotation keeps the bump from touching the boundaries r~+ p 2 .
In Methods -Solutions properties we verify that there are no solutions with both c and y { different from 0. We can therefore set c~0 in the steady state activity Eq. 3, and impose the activity itself to be zero on the boundary r~+ p 2 to find cos m p 2 ~{ s: This equation corresponds to the curve of separation in the plane m,J 1 ð Þ (using the relationship J 1~J1 m,s ð Þ, Eq. 8) between the single ring solution and a cylinder solution (Fig. 3A, curve surrounding the C region). In this regime, in addition to the freedom of choice for the location of the bump in H, the solution is also partially marginal in y { . The bump can be freely moved on a segment and a circle, defining a cylinder; the activity profile in this case is described by Eq. 4, see an instance in Fig. 2D. This region extends in the high J 1 limit and covers the whole range of correlations.
Despite the fact that each of the maps H A and H B defines a ring, it shouldn't come as a surprise that the topology of the attractor is a cylinder instead of a torus. The correlation between maps gives rise by definition to a cylinder structure, as can be seen for instance by inspecting Fig. 2B,II. It can be shown that when m~1 the cylinder solution degenerates in a torus; the bump of activity can be in any location of the h,r ð Þ coordinates (hence, also in (h A ,h B )). This regime is linked to the observation of an activity bump simultaneously localized in two environments in network simulations [39], and the study in [28]. Fig. 3 summarizes the results obtained so far. When J 1 is low, the only solutions is a constant level of activity which spreads over the whole network (Homogeneous region). As J 1 is increased, the interplay between the short range excitation and long range inhibition creates a pattern of localized activity in the middle map H (Single ring, see also Fig. 2B) or, if the correlation between maps is small enough, a localized pattern in either H A or H B (Double ring, Fig. 2C). Intuitively, the network ''remembers'' the two maps separately (y {~0 ,c=0, two solutions) if they are weakly correlated (mw*0:9). When the maps are more similar, the network represents just an average between them (y {~0 ,c~0).

Phase diagram
The bump size decreases with increasing J 1 . When J 1 is further increased, instead of having a reduced size of the localized activity in just one of the maps, the presence of two stored maps in the synaptic structure and the inhibition J 0 produce a packet of activity which looks localized in both maps (Cylinder solution, Fig. 2D).
Three particular values of the distance m deserve a special mention. The case m~0, corresponding to the encoding of two identical maps, can be shown to be identical to the ring model [37], as expected. In particular, besides the homogeneous solution and the amplitude instability region, the system can only exhibit the single ring solution.
The case m~1, corresponding to the encoding of two uncorrelated maps, does not have the single ring regime as a possible solution. The double ring solution in this case is depicted in Fig. 2A, where it can be seen that the bump is perfectly localized in either maps H A or H B , lacking any spatial tuning in the other map. This is the desired outcome in the ''multi-chart'' approach of [21].
The third case is m~1 2 . We will see in Results -Morphing maps that this case is closely related to the behavior of a network storing a morph sequence between two uncorrelated maps. As can be seen in the phase diagram, the double ring solution is not possible in this regime. How the environment, and the position in the environment, are represented by the network activity? For the single ring (Eq. 10) and the double ring (Eq. 11) solutions, both characterized by y {~0 , it is evident that the position is coded by the order parameter y z . The identity of the environment can only be represented with the ambiguity in the choice of the sign of c when the network operates in the double ring regime.
In the cylinder regime, it is not clear how the information about the environment is represented in the network, since now the solution is described by y z and y { . The following Section is mainly devoted to explore the link between the state variable (eventually time-dependent) Y t ð Þ in the active environment, and the behavior of the solution in this novel regime, by introducing a spatially tuned external input.

Tuned external input
Until now we considered the condition in which the only external input to the network, I, was steady and uniform. Let us introduce a tuned input, for instance in map H A at position Y t ð Þ: For simplicity we assume the shape of the external input to be I E x ð Þ~cos x ð Þ. The parameter measures the strength of the tuned component of the external input as a fraction of the constant baseline I we adopted so far. In general what we are interested in, and what is experimentally observable, are the tuning curves of the neurons i.e. their profile of activity as a function of the input angle in the active environment. It is easy to see the effect on the dynamics of the order parameters (Eqs. 5) when the location specific external current is inserted in the original dynamics for the network activity, Eq. 2. The dynamics keeps the same form as in Eq. 5, with the exception of the threshold-linear term in g h,r ð Þ, which now reads where m Ã~z m correspond to the choice of map H A in the input, and m Ã~{ m for map H B .

Tuning curves
With the input at a constant location Y t ð Þ~Y, one can see that a solution of Eqs. 5 for the single and double ring regime (y {~0 ), is y z~Y , i.e. the input pinpoints the location of the bump. This implies that, assuming a weak tuned input a vv1, the tuning curve of a neuron h,r ð Þ can be written in the single and double ring regime (from Eqs. 10 respectively. The tuning curve in the single ring regime has a maximum for Y~h (hence h is the preferred angle for a neuron h,r ð Þ), independently of which map m ?~+ m is being used in the external input, as can be seen from Eq. 14. This implies that each neuron has identical tuning curves in both environments, and that the preferred angle of a neuron does not coincide with either the assigned h A or h B but with their average.
For the double ring regime, the preferred angle assumes the form ( Fig. 2D). If several randomly selected external locations Y in one of the maps are presented to the network, once at time and starting from random initial conditions, the tuning curves would be an average over y { : where the allowed range for The cylinder regime extends the region of existence of two tuning curves per neurons to an higher correlation between the stored maps; the difference is that the coding becomes unreliable: during a single exposure to a given value of the input angle Y, a neuron could remain silent even if its average tuning curve would predict a response.
When the representation refers to the location in an environment, it is natural to think about a smoothly varying location Y. With a moving input like Y t ð Þ~v : t,vw0, the tuning curve depends as before on which map is stimulated, but in a novel way. Assume for simplicity to start from a y {~0 initial condition, corresponding to (y A~yB ). A moving input in the map H A would tend to move the bump along that map (i.e. increase the y A of the solution), while keeping y B constant (hence the bump will move to rv0). This movement is possible only until the bump reaches the part of configuration space not occupied by neurons due to the distance between maps m, see Fig. 2D. At that point, the bump will start to move equally along h A and h B , maintaining y { v0, which is proportional to y B {y A , and increasing y z (proportional to y B zy A ). A similar scenario, but with y { w0, is obtained when stimulating the map H B .
If the size of the bump is sufficiently small, this effect has dramatic consequences. The small bump will move along neurons with rv0 when a moving stimulus is presented in environment A, and viceversa neurons with rw0 will be active only when the moving stimulus is presented in environment B. As a consequence, neurons will essentially just have a tuning curve (or field), only in one map, and will be silent in the other one. We refer to this phenomenon as dynamical pattern separation (see Fig. 5 for an example). The separation of the activity patterns is essentially a dynamical phenomenon, dependent on the history of the inputs. The figure shows also the robustness of the dynamical pattern separation behavior to the addition of Gaussian d-correlated noise in the external current (see Methods -Numerical Methods). Note that neurons characterized by r*0 (i.e. h A~hB ), will have tuning curves in the same location. The number of neurons with tuning curves in both environments grows with the size of the bump. Note though that by changing the sign of the velocity in the moving input, the behavior would reverse; neurons with positive (negative) r would be active during a stimulation in map H A (H B ). In order to maintain the dynamical pattern separation and the analogy with place coding, one could think about two circular environments, as we did so far, with the additional constraint that the environments can only be traveled, for instance, in the counter-clockwise direction (CCW). As an alternative, the two environments may be thought as the same circular arena, but traveled clockwise (CW, environment A) and CCW (B); this interpretation would give rise to place fields with directional selectivity (see Discussion).
The dynamical pattern separation is basically dependent on the history of the input (positive or negative velocity), in addition to the identity of the map used in the stimulation. This history dependence is present also for non smooth time-dependent stimuli, as for instance the sequential presentation of stimuli with an intervening delay period. In this case the history dependence gives rise to a memory effect: the current location of the bump following a stimulation depends on the location attained after the previous stimulus presentation. Let us consider a basic example of this phenomenon, where the tuned external input is always presented in map H A . Consider for simplicity the state of the network being characterized by y A~yB~y0 , as a result of the presentation of stimulus y 0 sometime in the past. If we now present a stimulus y 1 , the bump will move, through the shortest arc on the map, to the new location y A~y1 . Depending on the stimuli, this movement can happen in two ways. If the shortest arc from y 0 to y 1 is directed CCW, the bump will move with a positive velocity vw0 and will end up being located in the rv0 region (as we previously saw in the case of moving tuned input). If the shortest arc is directed CW, then the movement will happen with a negative velocity, and the final location of the bump will be in the rw0 region. Hence, by looking at the activity resulting from the presentation of y 1 , we know whether the shortest way on the ring to it from y 0 is CW or CCW. A similar result can be obtained if the stimulus presentation alternates between map H A and H B .
If we vary the manifolds on which the maps live, for example to segments instead of circles, the history dependence changes accordingly. For instance, on segments the activity would give us information about the second stimulus being greater/smaller than the first one (see Discussion). In the next section we present a simple (albeit artificial) delayed discrimination task which the network can perform by exploiting the memory effect.

An application of the memory effect
Let us suppose to have a screen with a circle on it. A first stimulus (a dot) appears on the circle at some random location (described by an angle, y 0 ), for the duration of 1s. This first stimulus is then removed for a delay period of 2s. Then a second stimulus appears at another random angle y 1 ; the subject's task is to determine whether the shortest path on the circle from angle y 0 to y 1 is CW or CCW. The basic idea is that it is enough to look at the network activity (location of the bump in the r axis), to determine the relationship between the first and the second stimulus (see Results -Tuned external input for a description of the idea).
To test the ability of the network to solve this task, we numerically solve the dynamics for the order parameter (Results -Phase diagram of the model) with an external input (Results -Tuned external input) mimicking the presentation of the stimuli, for a sequence of 50 trials. We used no inter-trial interval, i.e. the presentation of the second stimulus in the k-th trial is immediately followed by the presentation of the first stimulus in trial kz1. The time courses of the bump location on the r axis (y { ) in two example trials for which y 0~0 , are shown in Fig. 6A. When looking at the location of the bump in the r axis at the end of a trial, there is a clear difference between the two cases of shortest CW, corresponding to positive y { (in the specific example y 1~3 p 2 ), or CCW arcs (y { v0, where y 1~p 2 ). Fig. 6B shows that the bump location at the end of trial, can be used to easily discriminate between the two possible answers (except for the cases in which the first and second stimuli are relatively close to each other). Note that this result has been obtained without any activity reset to new initial conditions during the inter-trial intervals.

Morphing maps
How do the results described so far change when, instead of storing just two correlated maps, the network encodes a sequence of maps gradually morphed between two uncorrelated ones? Let us start by constructing two random uncorrelated maps, H {1 and H 1 . We would like to define the intermediate maps as gradual rotations between the two extreme ones; since we are dealing with circles, the rotation should be performed along the shortest arc between h {1 and h 1 (see Eq. 21, Methods -Inverse transformation). We assume here to have already transformed the variables in such a way that we can write directly where f [ ½{1,1 indexes the maps along the morph sequence. Hence a neuron with label h {1 in the first map, will rotate along the sequence to its location h 1 on the last map, following the shortest path on the circle. With this choice of the morphing procedure, each neuron is still characterized by just two quantities, its labels in the extreme maps. We store the whole morph sequence by a superposition of the synaptic structures generated in each map separately, as for the case of two correlated maps previously described. For the sake of analytical tractability, we study the resulting coupling in the limit p??
Introducing the definition of two uncorrelated maps (Eq. (1) with m~1) into Eq. (16) The first term of the infinite product in the Euler formula, or the first term in the limit sum, gives us cos r{r' 2 . Comparing the coupling in Eq. (18), and the one derived for two maps, Eq. (2), we see that to first order, the synaptic coupling induced by the storage of the whole morph sequence, is equivalent to the storage of two correlated maps with m~1 2 .
In Fig. 7, we compare the network activity generated by the approximated coupling and the full result of Eq. 18, when the external input is constant. The results are qualitatively similar but the full morph case reaches the cylinder regime for lower J 1 compared to the m~1 2 case. Note that the network storing the morph sequence shows the same dynamical pattern separation observed in the two maps case (Fig. 8) In the experiment of [32], the rat is trained until it develops two separate place coding for a single arena with different light configurations (representing two distinct environments). The advantage of this setup is that it allows, for instance, to slowly morph the light configuration between the two environments familiar to the rat. The experimental results shows a sharp transition around the middle of the light morphing (lasting  Fig. 8 for the approximated whole morph sequence storage, for two slightly correlated maps in the cylinder region of the parameter range and for the double ring regime. For each run we show the dynamics of the relevant order parameter for the regime under consideration, c for the double ring case and y { for the cylinder solution. In addition, we numerically solve the dynamics for a moving stimulus in either environment A or B. We use this as a reference for computing, at each time step, the correlation coefficient between the network activity during the morphing protocol and the activity in the fixed environment. The transition is sharpest for the storage of two slightly correlated maps. Note that similar results would be obtained by testing the network separately in each environment of the sequence (see e.g. [31]). The sharp transition is maintained when increasing the amplitude of the external tuned input, because a small tilt in the tuned input towards either map A or B is sufficient to generate the dynamical pattern separation described in the previous Section. The transition in the cylinder regime occurs few seconds later than the one occurring in the double ring regime, which in turn happens in the middle of the morphing ( This delay is due to the time required for the bump to move from the region of rv0 to rw0, or viceversa (see also Fig. 5B). This result could be compared with the experimental results of [32]. The delay does not occur when testing the network in separate environments along the morph sequence. There are two additional observations to be made (data not shown). The first one is related to the sharpness of the transition in the double ring regime; by further reducing the amplitude of the external input, the mean-field dynamics can produce a sharp transition between the environments representations, which is also delayed compared to the middle of the morphing period. The delay gets longer as the external input gets weaker, in extreme cases it happens just before the end of the morphing procedure. This sharp and delayed transition is not observed in microscopic simulations with up to 2000 neurons, since the weak input is not able to overcome the local inhomogeneities in which the bump is trapped (see e.g. [18]). It is possible that in larger networks the transition can be observed. The fine-tuning of the external input strength required to have the transition around the middle of the sequence, makes the double ring regime a weaker candidate explanation for the experimental results of [32] compared to the cylinder regime.
The second observation concerns the dependency of the transition parameters on the velocity of the moving external input. We have noticed that the transition becomes smoother and closer to the middle as the velocity of the simulated animal is reduced. The details of the transition in a realistic setting would depend on the velocity history of the animal.

Comparison with simulations
In order to verify that the results obtained in the previous Sections are not artifacts coming from our assumptions of having an infinite number of neurons (and maps, referred to the morphing procedure) we compare some of the MF predictions to simulations of networks with a finite number N of neurons. Each neuron is assigned a random pair of labels (h (i) ,r (i) , for the i-th neuron), from which we create either two maps with distance m, or a finite number p of maps (h (i) f k , for the k-th map) along the morph sequence between two uncorrelated references (see Methods -Numerical Methods).
In Fig. 9 we compare the order parameters from MF and estimated from simulations, at a fixed value of the distance between the maps and inhibition. Varying J 1 , the solution goes through the double ring, single ring and cylinder regime. The order parameter c is particularly sensitive to the finite size of the network (and the randomized maps, see [18]). Fig. 10 shows the time evolution of a network storing few maps from a morph sequence. This is the best example to show dynamical pattern separation at finite size, since it is less intuitive than the case of two correlated maps. From an arbitrary initial position, the bump of activity starts moving first towards negative r (increasing angles in map A), then along increasing h without changing its location in r. Note that, despite the presence of neurons everywhere in the (h A , h B ) plane, the bump moves along an invisible barrier resulting from the storage of the morph sequence.
We have also verified that all the qualitative behaviors, number and type of solutions, unreliable coding, dynamical pattern separation and memory effect, are maintained when moving from maps on rings, to segments (either two correlated maps or morphed), as studied e.g. in [18,37] for the single map (data not  shown). Instead of having neurons arranged on a cylinder in the Þcoordinates, as for the ring case (see e.g. Fig. 2B,II), the geometry resulting from two correlated linear maps would be an infinite strip. A strong enough map-specific interaction would produce a bump localized in both maps. An external moving input in one of the maps would move the bump on the strip up to the boundary, and then the bump would crawl along such boundary. Depending on the direction of the moving input or the identity of the stimulated environment, the bump can settle either in ''upper'' of ''lower'' part of the strip as in the cylinder regime.

Discussion
We have studied a continuous attractor network model storing a pair of correlated maps. The storage of a morph sequence between two uncorrelated maps falls in this class of model, since it is approximately equivalent to the storage of two strongly correlated maps. The other relevant parameter for describing the possible network behaviors, beside the correlation between the maps, is the strength of map-specific interaction between neurons.
The analysis of the solutions of the system with a weak tuned external input, reveals several interesting behaviors. When the correlation between the maps is weak, neurons have two different tuning curves corresponding to the stimulus presentation in different maps. The representation is reliable, in that the single neuron response is consistent between presentations. This is the operating regime which is usually considered useful in place coding applications.
For higher correlations between the maps and weak mapspecific interactions, each neuron possesses only one tuning curve, irrespectively of the stimulated map. In contrast to the previous regime, this one is rendered useless by the inability to represent fully the state of the external world, i.e. the identity of the environment in the context of place coding analogy.
We find another, novel regime for strong interactions and for any amount of correlation between maps. The surprising aspect of this regime is that the state of the world does not uniquely determine the state of the network; there is an additional degree of freedom in the network representation.
To a closer look, this additional freedom found in the novel regime is rich of consequences. When the external input location is randomly varied between presentations in one map, we can define the response of a neuron to a particular location as an average of the neuron activity over external input presentations in that location. In this context each neuron has different tuning curves relative to the different maps used in the stimulation, but the price to pay is unreliable coding; a neuron which should be active during a particular state of the world, could remain silent.
When the location of the external input changes smoothly in time on one map, some neurons develop a selectivity to the direction of change. When the increase happens on the other map, another subset of neurons fires. The overlap between the two subsets may be arbitrarily small, depending on the parameters choice. Neurons active in both maps would have tuning curves around similar values of the external input location. We refer to this phenomenon as dynamical pattern separation. There is an ambiguity in the network representation, due to the fact that the subset of neurons activating with the increase of the external location in map A, will also activate with a decrease of the location in map B. There are three possible experimental contexts in which this ambiguity does not arise.
A simple experimental context would arise if the input is tuned in only one of the two maps and the only parameter changing is the location of the external input. Given some state variable, like size and orientation of objects, or frequency of sound waves for instance, our model would produce respectively tuning for expansion/contraction, CW/CCW rotation and upward/downward frequency sweeps (all experimentally observed, see e.g. [40,41]).
Our model provides a unique way for producing selectivity for the direction of change of a state variable, given a selectivity for the variable itself. Both kind of responses give rise to another interesting phenomenon: The current representation of the state of the world is influenced by the preceding one, even with an intervening delay. It is possible to read out from the network the direction of change of the state variable. This property may be exploited when solving delayed discrimination tasks (see [42] for data analysis and modeling in terms of remapping for a somatosensory discrimination task).
A second experimental context is related to place coding; the two environments should be considered as two distinct circular arenas which can be traveled only in one direction. Experimental observations show that when an animal is exposed to two environments, the majority of place cells have a place field in only one of the two environments (see e.g. [43,44]). A possible experiment to test the model would consist in training the animals in two well differentiated environments. After measuring the distance between preferred locations for neurons having tuning curves (place fields) in both environments, one could train the animals in intermediate environments, which would correspond to the storage of the morph sequence in the model. For the novel regime of the model, the disappearance of the place fields in one of the environment would be predicted for neurons with very different preferred locations, and the remaining fields will converge to a common representation. Alternatively the training could be performed by using the initial two environments, and then slowly changing them across several training days to increase their similarity. This would correspond to the storage of two correlated environments.
A third experimental context is related to direction selectivity in place cells. Animals trained to shuttle back and forth in a onedimensional track (a segment or a circle), have place cells showing selectivity to the direction of motion. For instance a cell could be active in a certain region of the circular environment when the animal is moving clockwise, while being completely silent when the animal moves counterclockwise. The link with our model is provided by the simple observation that the same 1D track, but walked in opposite directions, correspond to two different environments. Dynamical pattern separation would produce directional selective neurons, while a neuron having place fields in both environments would have similar preferred locations. In [45], place cells recorded from rats trained in a circular environment indeed showed bi-directional place fields in similar locations. There was however a systematic bias in the difference between the preferred locations in the CW and CCW directions of the majority of the bi-directional cells: place fields were displaced backward with respect to the direction of motion of the animal. We believe that this result, termed by the authors ''prospective misalignment'', could be obtained in the context of our model in more than one way. One possibility is the introduction of an asymmetry in the synaptic connections (following [46]), with the asymmetry determined by the emerging direction selectivity of the neurons. The spread of activity due to the asymmetry would activate neurons earlier compared to the symmetric case, reproducing the prospective misalignment. A similar result could be obtained with short-term synaptic plasticity, which is known to produce a moving bump of activity ( [47]). A third option could be the introduction of a systematic shift between the maps, possibly resulting from Hebbian learning of the configurations generated by the suggested asymmetry mechanisms.
In the experiments of [32], two environments correspond to two different light configurations in the same arena. A slow linear morph between light configurations results in a sharp transition from the population representation for one environment to the other. This is a promising experimental technique which is able to probe with unprecedented flexibility the dynamics of remapping between two environments or along a morph sequence [32], and could serve as a fertile ground for our model's predictions, hence for testing the attractor hypothesis. We show that, in agreement with the experiment, the slow morph protocol produces sharp transitions due to dynamical pattern separation. This result is even more significant considering the acknowledged difficulties in reproducing sharp transitions between correlated maps in a ''traditional'' setting [36]. The model predicts a transition between representations slightly delayed compared to half of the morphing period; it remains to be seen whether this occurs also in the experiment.
Our results can be related to experimental observations about changes in place representation between distinct environments. Two major classes of remapping have been observed when an animal is tested in two distinct environments: rate remapping, in which cells maintain the positions of their firing fields while differentially changing their amplitudes, and global remapping, where changes in firing location are observed in addition to firing rate modifications (see e.g. [43]). Based on these properties, we could associate the double ring regime to the global remapping and the cylinder regime to the rate remapping.
The model results can also be compared to experiments with sequences of continuously morphed environments. When animals explored intermediate environments, both sharp and smooth transitions in representations were observed in different experiments (see [31] and [30] correspondingly). Our model exhibits both sharp transitions between the place representations corresponding to intermediate environments (cylinder regime) and smooth transitions (double ring regime).
The linkage of cylinder and double ring regimes to sharp and smooth transitions respectively, taken together with the above mentioned association between these two model regimes with global and rate remapping, would be against the hypothesis made in [30] that related global remapping and sharp transitions on one hand, and rate remapping with smooth transitions on the other. In the present form, our model cannot be made compatible with this hypothesis. Since both the recordings of [31] and [30] contained populations of neurons exhibiting different transition behaviors, we speculate that the introduction of an additional selectivity for the environments (see below) could help in resolving the contradiction. Rate remapping would then correspond to a mixed single ring-cylinder regime (different subsets of the network would exhibit the different regimes), while global remapping would resemble a mix of the double ring and cylinder regimes.
A future extension of the model would include neurons with some form of selectivity for the context; each neuron would then be characterized not only by its location on the two maps, but also by selectivity indexes measuring its ''preference'' for the maps (e.g. [17]). This more realistic setting including selectivity would produce silent neurons and place fields with variable peak rates/ widths even when storing a single map.
A second issue to be addressed is how the network can learn the synaptic structure from its inputs. The long-term plasticity (e.g. [33,34]), could bring the network through various operating regimes depending on the training protocol. This could impose additional constraints on the model and provide additional predictions.
Finally, with the introduction of short-term plasticity [48][49][50][51][52][53], the network could exhibit an even richer repertoire of dynamics. This extension of the model would be an important step towards the experimental results of [32]. In this study, it was observed that when there is a fast switch between the two light configurations, the population vector sometime oscillates between the place representation of the environments, before settling on the current one. Preliminary results coming from the introduction of short term facilitation and depression in a network exhibiting a double ring solution, show that is indeed possible to observe oscillations between place representations. A detailed analysis of this behavior will be matter for a future report.

Numerical methods
To solve numerically the MF dynamics described by Eq. 5, we The integrals in the rhs of the equations were estimated using a trapezoidal method. The system of ODEs were integrated with an adaptive 4-th order Runge Kutta scheme.
The simulation of the microscopic networks, whose results are reported in Figs. (9,10), were performed by solving numerically the system of ODEs where j~1 Á Á Á N indexes the neurons. The matrix W ij is built by summing the single map encoding W ij~1 p To obtain the p labels h f i characterizing each neuron, we first randomly generated a h i ,r i and used Eq. 1 for p~2 or Eq. 16 for pw2. For the comparison of the simulation with the MF results in Fig. 9 from which we constructed the estimates for the order parameters, using Eq. 24.
For the noisy simulations shown in Fig. 5B, we used a currentbased version of the dynamics described by Eqs. 19: We then estimated the order parameters via Eqs. 20, using the firing rates m j~hj Â Ã z . The noise was introduced as an additional term in the current where B t ð Þ is a zero average, unit variance Gaussian d-correlated noise. We used N~0 :05 for the results in Fig. 5. The numerical solution was obtained using the Euler-Maruyama integration scheme.
The simulations performed in Fig. 7, for a network storing the whole morph sequence, were carried out as follows. Substituting the synaptic coupling obtained in Eq. 18 with the one in Eq. 2, it is possible to derive a dynamics for the ''order function'' The 2p rotation in the first equation is just needed to select the shortest distance between two maps on a ring, and it is transparent for the connectivity given its periodicity. This rotation was implicitly assumed when defining the neurons locations along the morph sequence, Eq. 16.

Reduced dynamics
A first reduction of the dynamics described by Eq. 2 is done using the first two Fourier components of the activity m with respect to the two correlated maps h A and h B , rewritten in terms of center map h and the distance r using Eq. 1. In line with [37] we define the following variables From Eqs. 23, after some algebra, it is finally possible to obtain the dynamics of the new order parameters, Eq. 5.

Solutions properties
In Methods -Phase diagram of the model, we mentioned that the equation _ y y z~0 from Eqs. 5, i.e.
where C~ð ð DhDr cos 2 h cos 2 mr ð Þ: Therefore, two conditions must be satisfied for the solution to be stable: J 0 v1 (amplitude instability) and J 1 v 1 C . Evaluating the integral in C explicitly, we get the line of separation between the homogeneous solution and the localized bump (Turing instability), expressed in Eq. 7. For the single, double ring and cylinder solution we can linearize directly Eq. 5, posing y z~y{~0 . The matrix associated with the dynamics of the vector dc,da,ds ð Þ , after using the fixed points equations (Eq. 8), is We define Let us examine the single ring and cylinder solution, c ?~0 . Given the symmetry in the integrand, A C~ASC~0 in this case. The stability matrix from Eq. 25 becomes then It is immediately seen that the eigenvalue corresponding to a destabilization of c ?~0 changes sign when J 1 A S 2~1. Substituting for J 1 the expression in Eq. 27, it is easy to verify that this reproduces the curve of separation between the single and double ring regime described by Eq. 12. We analyze the remaining two eigenvalues by looking at the trace T and the determinant D of the 2|2 sub-matrix in the a,s subspace (from 28): Given that A C §0, we see immediately that the eigenvalues have the same sign for a ? w0, and one of them changes sign when a ? v0. Recall that a ?~0 is the onset of amplitude instability we introduced without proof in Eq. 9. If we find that when Dw0 the trace is negative, then we know that a ?~0 correspond to a destabilization of the solution. Using Eq. 27, we see that imposing Dw0 is equivalent to , and that the trace T satisfies The numerator in the first term is non-negative (A C §0). The denominator is simply ÐÐ DhDrg, non-negative by definition. The denominator in the second term is J {1 1 w0, and we can write the numerator as ðð DhDrh h,r ð Þ cos h cos mr ð Þzs ? ð Þ 2 w0: Finally, we numerically verified that the region of existence of the double ring solution coincides with its stability region.