Tracking Changing Stimuli in Continuous Attractor Neural Networks
C. C. Alan Fung, K. Y. Michael Wong Department of Physics The Hong Kong University of Science and Technology Clear Water Bay, Hong Kong, China alanfung@ust.hk, phkywong@ust.hk Si Wu  Department of Informatics, University of Sussex, Brighton, United Kingdom siwu@sussex.ac.uk

Abstract
Continuous attractor neural networks (CANNs) are emerging as promising models for describing the encoding of continuous stimuli in neural systems. Due to the translational invariance of their neuronal interactions, CANNs can hold a continuous family of neutrally stable states. In this study, we systematically explore how neutral stability of a CANN facilitates its tracking performance, a capacity believed to have wide applications in brain functions. We develop a perturbative approach that utilizes the dominant movement of the network stationary states in the state space. We quantify the distortions of the bump shape during tracking, and study their effects on the tracking performance. Results are obtained on the maximum speed for a moving stimulus to be trackable, and the reaction time to catch up an abrupt change in stimulus.

1 Introduction
Understanding how the dynamics of a neural network is shaped by the network structure, and consequently facilitates the functions implemented by the neural system, is at the core of using mathematical models to elucidate brain functions [1]. The impact of the network structure on its dynamics is twofold: on one hand, it decides stationary states of the network which leads to associative memory; and on the other hand, it carves the landscape of the state space of the network as a whole which may contribute to other cognitive functions, such as movement control, spatial navigation, population decoding and object categorization. Recently, a type of attractor networks, called continuous attractor neural networks (CANNs), has received considerable attention (see, e.g., [2, 3, 4, 6, 7, 8, 9, 10, 11, 12, 13, 5]). These networks possess a translational invariance of the neuronal interactions. As a result, they can hold a family of stationary states which can be translated into each other without the need to overcome any barriers. Thus, in the continuum limit, they form a continuous manifold in which the system is neutrally stable, and the network state can translate easily when the external stimulus changes continuously. Beyond pure memory retrieval, this large-scale stucture of the state space endows the neural system with a tracking capability. This is different from conventional models of associative memory, such as the Hopfield model [14], in which the basin of each attactor is well separated from the others.
 Current address: Institute of Neuroscience, Shanghai Institutes for Biological Sciences, Chinese Academy of Sciences - China.

1


The tracking dynamics of a CANN has been investigated by several authors in the literature (see, e.g., [3, 4, 5, 8, 11]). These studies have shown that a CANN has the capacity of tracking a moving stimulus continuously and that this tracking property can well justify many brain functions. Despite these successes, however, a detailed analysis of the tracking behaviors of a CANN is still lacking. These include, for instance, 1) the conditions under which a CANN can successfully track a moving stimulus, 2) the distortion of the shape of the network state during the tracking, and 3) the effects of these distortions on the tracking speed. In this paper we will report, as far as we know, the first systematic study on these issues. We hope this study will help to establish a complete picture about the potential applications of CANNs in neural systems. We will use a simple, analytically-solvable, CANN model as the working example. We display clearly how the dynamics of a CANN is decomposed into different distortion modes, corresponding to, respectively, changes in the height, position, width and skewness of the network state. We then demonstrate which of them dominates the tracking behaviors of the network. In order to solve the dynamics which is otherwise extremely complicated for a large recurrent network, we develop a time-dependent perturbation method to approximate the tracking performance of the network. The solution is expressed in a simple closed-form, and we can approximate the network dynamics up to an arbitory accuracy depending on the order of perturbation used. We expect that our method will provide a useful tool for the theoretical studies of CANNs. Our work generates new predictions on the tracking behaviors of CANNs, namely, the maximum tracking speed to moving stimuli, and the reaction time to sudden changes in external stimuli, both are testable by experiments.

2 The Intrinsic Dynamics of CANNs
We consider a one-dimensional continuous stimulus being encoded by an ensemble of neurons. The stimulus may represent, for example, the moving direction, the orientation, or a general continuous feature of an external object. Let U (x, t) be the synaptic input at time t to the neurons with preferred stimulus of real-valued x. We will consider stimuli and responses with correlation length a much less than the range of x, so that the range can be effectively taken to be (-, ). The firing rate r(x, t) of these neurons increases with the synaptic input, but saturates in the presence of a global activity-dependent inhibition. A solvable model that captures these features is given by r(x, t) = U (x, t)2 d , 1 + k x U (x , t)2 (1 )

where  is the neural density, and k is a small positive constant controlling the strength of global inhibition. The dynamics of the synaptic input U (x, t) is determined by the external input Iext (x, t), the network input from other neurons, and its own relaxation. It is given by d dU (x, t)  = Iext (x, t) +  x J (x, x )r(x , t) - U (x, t), (2 ) dt where  is the time constant, which is typically of the order 1 ms, and J (x, x ) is the neural interaction from x to x. The key characteristic of CANNs is the translational invariance of their neural interactions. In our solvable model, we choose Gaussian interactions with a range a, namely,  J (x, x ) = exp[-(x - x )2 /(2a2 )]J / 2 a2 . (3 ) CANN models with other neural interactions and inhibition mechanisms have been studied [2, 3, 4, 7, 9]. However, our model has the advantage of permitting a systematic perturbative improvement. Nevertheless, the final conclusions of our model are qualitatively applicable to general cases (to be further discussed at the end of the paper). We first consider the intrinsic dynamics of the CANN model in the absence of external stimuli. For  0 < k < kc  J 2 /(8 2 a), the network holds a continuous family of stationary states, which are - , (x - z )2 ~ U (x|z ) = U0 exp (4 ) 4 a2  where U0 = [1 + (1 - k /kc )1/2 ]J /(4  ak ). These stationary states are translationally invariant among themselves and have the Gaussian bumped shape peaked at arbitrary positions z . 2


Height v0
-2 1 0.5 0 -1 0 -0.5 -1 1 2 -2

v1

0.5 0.25

Position

0 -1 0 -0.25 -0.5

1

2

Width v2
-2 1 0.5 0 -1 0 -0.5 -1 1 2 -2

Skew v3
1 0.5 0 -1 0 -0.5 -1 1 2

Figure 1: The first four basis functions of the quantum harmonic oscillators, which represent four distortion modes of the network dynamics, namely, changes in the height, position, width and skewness of a bump state. The stability of the Gaussian bumps can be studied by considering the dynamics of fluctuations. ~ Consider the network state U (x, t) = U (x|z ) +  U (x, t). Then we obtain d d   U (x, t) = x F (x, x ) U (x , t) -  U (x, t), (5 ) dt dJ where the interaction kernel is given by F (x, x ) =  x (x, x ) r(x )/ U (x ).

2.1 The motion modes

To compute the eigenfunctions and eigenvalues of the kernel F (x, x ), we choose the wave functions of the quantum harmonic oscillators as the basis, namely, exp(- 2 /2)Hn ( ) ( , (6 ) vn (x|z ) = 2 )1/2 an!2n  where   (x - a)/( 2a) and Hn ( ) is the nth order Hermite polynomial function. Indeed, the ground state of the quantum harmonic oscillator corresponds to the Gaussian bump, and the first, second, and third excited states correspond to fluctuations in the peak position, width, and skewness of the bump respectively (see Fig. 1). The eigenvalues of the kernel F are calculated to be 0 = 1 - (1 - k /kc )1/2 ; n = 1/2n-1 , for n  1. (7 ) The eigenfunctions of F can also be analytically calculated, which turn out to be either the basis functions vn (x|z ) or a linear combination of them. Here we only listhe first four of t them, wh1ch are u0 (x|z ) = v0 (x|z ), u1 (x|z ) = v1 (x|z ), u2 (x|z ) = 1/( 2D0 )v0 (x|z ) + i 1 (1 1 2 - - k /kc )6 D0 v2 (x|z ), with D0 = [(1 - 2 / - k /kc )2 + 1/2]1/2 and u3 (x|z ) = /7v1 (x, z ) + /7v3 (x, z ).

The eigenfunctions of F correspond to the various distortion modes of the bump. Since 1 = 1 and all other eigenvalues are less than 1, the stationary state is neutrally stable in one component, and stable in all other components. The first two eigenfunctions are particularly important. (1) The eigenfunction for the eigenvalue 0 is u0 (x|z ), and represents a distortion of the amplitude of the bump. As we shall see, amplitude changes of the bump affect its tracking performance. (2) Central to the tracking capability of CANNs, the eigenfunction for the eigenvalue 1 is u1 (x|z ) and is neutrally stable. We note that u1 (x|z )   v0 (x|z )/ z , corresponding to the shift of the bump position among the stationary states. This neutral stability is the consequence of the translational invariance of the network. It implies that when there are external inputs, however small, the bump will move continuously. This is a unique property associated with the special structure of a CANN, not shared by other attractor models. Other eigenfunctions correspond to distortions of the shape of the bump, for example, the eigenfunction u3 (x|z ) corresponds to a skewed distortion of the bump. 3


2 1.5

U(x)

1 0.5 0 -2 0 x 2

Figure 2: The canyon formed by the stationary states of a CANN projected onto the subspace formed by b1 |0 , the position shift, and b0 |0 , the height distortion. Motion along the canyon corresponds to the displacement of the bump (inset). 2.2 The energy landscape It is instructive to consider the energy landscape in the state space of a CANN. Since F (x, x ) is not symmetric, a Lyapunov function cannot be derived for En . (5). Nevertheless, for each peak position q (1 - n )bn |2 /2, where bn |z is the overlap z , one can define an effective energy function E |z = z ~ between U (x) - U (x|z ) and the nth eigenfunction of F centered at z . Then the dynamics in Eq. (5) can be locally described by the gradient descent of E |z in the space of bn |z . Since the set of points bn |z = 0 for n = 1 traces out a line with E |z = 0 in the state space when z varies, one can envisage a canyon surrounding the line and facilitating the local gradient descent dynamics, as shown in Fig. 2. A small force along the tangent of the canyon can move the network state easily. This illustrates how the landscape of the state space of a CANN is shaped by the network structure, leading to the neutral stability of the system, and how this neutral stability shapes the network dynamics.

3 The Tracking Behaviors
We now consider the network dynamics in the presence of a weak external stimulus. Suppose the neural response at time t is peaked at z (t). Since the dynamics is primarily dominated by the translational motion of the bump, with secondary distortions in shape, we may develop a time-dependent perturbation analysis using {vn (x|z (t))} as the basis, and consider perturbations in increasing orders of n. This is done by considering solutions of the form n ~ an (t)vn (x|z (t)). (8 ) U (x, t) = U (x|z (t)) +
=0

Furthermore, since the Gaussian bump is the steady-state solution of the dynamical equation in the absence of external stimuli, the neuronal interaction term in Eq. (2) can be linearized for weak stimuli. Making use of the orthonormality and completeness of {vn (x|z (t))}, we obtain from Eq. (2) expressions for dan /dt at each order n of perturbation, which are d 1 a U(   In dz 1 - n 1/2 a nan-1 - n + 1an+1 + - 0 2 ) n1 + n= dt   2a dt (  1r n + 2r)! (-1)r an+2r , (9 ) +  =1 n! 2n+3r-1 r! where In (t) is the projection of the external input Iext (x, t) on the nth eigenfunction. Determining z (t) by the center of mass of U (x, t), we obtain the self-consistent condition . I n  !!/(n - 1)!!In + a1 1+ dz 2a n=3,o dd ( ( =  dt  U0 2 )1/2 a + n=0,even n - 1)!!/n!!an

(1 0 )

Eqs.(9) and (10) are the master equations of the perturbation method. We can approximate the network dynamics up to an arbitary accuracy depending on the choice of the order of perturbation. In practice, low order perturbations already yield very accurate results. 4


4 3 2 1 0 0 s

(a)

1 0.8 0.6 s 0.4 0.2
50 100 150 200 250

(b)

t

0 0

vmax 0.01 0.02 v 0.03 0.04

Figure 3: (a) The time dependence of the separation s starting from different initial values. Symbols: simulations with N = 200 and v = 0.025. Lines: n = 5 perturbation. Dashed lines: s1 (bottom) and s2 (top). (b) The dependence of the terminal separation s on the stimulus speed v . Symbols: simulations with N = 200. Dashed line: n = 1 perturbation. Parameters:  = 0.05, a = 0.5,   = 1, k = 0.5,  = N/(2 ), J = 2 a2 . 3.1 Tracking a moving stimulus Consider the external stimulus consisting of a Gaussian bump, namely, Iext (x, t) = U0 exp[-(x - z0 )2 /4a2 ]. Perturbation up to the order n = 1 yields a1 (t) = 0, [d/dt + (1 - 0 )/ ]a0 = ( U0 2 )1/2 a exp[-(z0 - z )2 /8a2 ]/ , and R - dz  (z0 - z )2 (t)-1 , (1 1 ) = (z0 - z ) exp dt  8 a2 t where R(t) = 1 +  - (dt / ) exp[-(1 - 0 )(t - t )/ - (z0 - z (t ))2 /8a2 ], representing the ratio of the bump height relative to that in the absence of the external stimulus ( = 0). Hence, the dynamics is driven by a pull of the bump position towards the stimulus position z0 . The factor R(t) > 1 implies that the increase in amplitude of the bump slows down its response. The tracking performance of a CANN is a key property that is believed to have wide applications in neural systems. Suppose the stimulus is moving at a constant velocity v . The dynamical equation becomes identical to Eq. (11), with z0 = v t. Denoting the lag of the bump behind the stimulus by s = z0 - z we have, after the transients, 1 -1 2 2 2 2 e-s /8a ds se-s /8a + = v - g (s); g (s)  . (1 2 ) dt  1 - 0

The value of s is determined by two competing factors: the first term represents the movement of the stimulus, which tends to enlarge the separation, and the second term represents the collective effects of the neuronal recurrent interactions, which tends to reduce the lag. Tracking is maintained when these two factors match each other, i.e., v = g (s); otherwise, s diverges.  The function g (s) is concave, and has the maximum value of gmax = 2a/( e) at s = 2a. This means that if v > gmax , the network is unable to track the stimulus. Thus, gmax defines the maximum trackable speed of a moving stimulus. Notably, gmax increases with the strength of the external signal and the range of neuronal recurrent interactions. This is reasonable since it is the neuronal interactions that induce the movement of the bump. gmax decreases with the time constant of the network, as this reflects the responsiveness of the network to external inputs. On the other hand, for v < gmax , there is a stable and unstable fixed point of Eq. (12), respectively denoted by s1 and s2 . When the initial distance is less than s2 , it will converge to s1 . Otherwise, the tracking of the stimulus will be lost. Figs. 3(a) and (b) show that the analytical results of Eq. (12) well agree with the simulation results. 3.2 Tracking an abrupt change of the stimulus Suppose the network has reached a steady state with an external stimulus stationary at t < 0, and the stimulus position jumps from 0 to z0 suddenly at t = 0. This is a typical scenario in experi5


ments studying mental rotation behaviors. We first consider the case that the jump size z0 is small compared with the range a of neuronal interactions. In the limit of weak stimulus, the dynamics is described by Eq. (11) with R(t) = 1. We are interested in estimating the reaction time T , which is the time taken by the bump to move to a small distance  from the stimulus position. The reaction time increases logarithmically with the jump size, namely, T  ( /) ln(|z0 |/).
400 (a) 300 200 100 0 0 T
Simulation "n=1" perturbation "n=2" perturbation "n=3" perturbation "n=4" perturbation "n=5" perturbation

2 (b) 1.5 U(x)
2 2.5 3

1 0.5

0.5

1

1.5 z0

0

-2

0 x

2

Figure 4: (a) The dependence of the reaction time T on the new stimulus position z0 . Parameters: as in Fig.3. (b) Profiles of the bump between the old and new positions at z0 =  /2 in the simulation. When the strength  of the external stimulus is larger, improvement using a perturbation analysis up to n = 1 is required when the jump size z0 is large. This amounts to taking into account the change of the bump height during its movement from the old to new position. The result is identical to Eq. (11), with R(t) replaced by . - + - t  (1 - 0 ) (1 - 0 ) (z0 - z (t ))2 dt R(t) = 1 + ex p t  ex p (t - t ) - 1 - 0   8 a2 0 (1 3 ) Indeed, R(t) represents the change in height during the movement of the bump. Contributions from the second and third terms show that it is highest at the initial and final positions respectively, and lowest at some point in between, agreeing with simulation results shown in Fig. 4(b). Fig. 4(a) shows that the n = 1 perturbation overcomes the insufficiency of the logarithmic estimate, and has an excellent agreement with simulation results for z0 up to the order of 2a. We also compute the reaction time up to the n = 5 perturbation, and the agreement with simulations remains excellent even when z0 goes beyond 2a. This implies that beyond the range of neuronal interaction, tracking is influenced by the distortion of the width and the skewed shape of the bump.

4 The Two-Dimensional Case
We can straightforwardly extend the above analysis to two-dimensional (2D) CANNs. Consider a neural ensemble encoding a 2D continuous stimulus x = (x1 , x2 ), and the network dynamics satisfies Eqs. (1-3) with x and x being replaced by x and x , respectively. We can check that the network holds a continuous family of stationary states given by - , (x - z)2 ~ U (x|z) = U0 exp (1 4 ) 4 a2 where z is a free parameter indicating the position of the network state in a 2D manifold, and (x - z)2 = (x1 - z1 )2 + (x2 - z2 )2 the Euclidean distance between x and z. By applying the stability analysis as in Sec. 2, we obtain the distortion modes of the bump dynamics, which are expressed as the product of the motion modes in the 1D case, i.e., um,n (x|z) = um (x1 |z1 )un (x2 |z2 ), f o r m, n = 0 , 1 , 2 , . . . (1 5 )

The eigenvalues for these motion modes are calculated to be 0,0 = 0 , m,0 = m , for m = 0, 0,n = n , for n = 0, and m,n = m n , for m = 0 and n = 0. 6


The mode u1,0 (x|z) corresponds to the position shift of the bump in the direction x1 and u0,1 (x|z) the position shift in the direction x2 . A linear combination of them, c1 u1,0 (x|z) + c2 u0,1 (x|z), corresponds to the position shift of the bump in the direction (c1 , c2 ). We see that the eigenvalues for these motion modes are 1, implying that the network is neutrally stable in the 2D manifold. The eigenvalues for all other motion modes are less than 1. Figure 5 illustrates the tracking of a 2D stimulus, and the comparison of simulation results on the reaction time with the perturbative approach. The n = 1 perturbation already has an excellent agreement over a wide range of stimulus positions.

400 (b)
(a) 1.2 1 0.8 0.6 0.4 0.2 0 -3 -2 -1 x

300 200 100 00 T

Simulation Theory

U(x,y)

0

1

2

3 -3

-2

1 0y -1

2

3

0.5

1

1.5 2 |z0 - z(0)|

2.5

3

Figure 5: (a) The tracking process of the network; (b) The reaction time vs. the jump size. The simulation result is compared with the theoretical prediction. Parameters: N = 40 × 40, k = 0.5,  a = 0.5,  = 1, J = 2 a2 ,  = N/(2 )2 and  = 0.05.

5 Conclusions and Discussions
To conclude, we have systematically investigated how the neutral stability of a CANN facilitates the tracking performance of the network, a capability which is believed to have wide applications in brain functions. Two interesting behaviors are observed, namely, the maximum trackable speed for a moving stimulus and the reaction time for catching up an abrupt change of a stimulus, logarithmic for small changes and increasing rapidly beyond the neuronal range. These two properties are associated with the unique dynamics of a CANN. They are testable in practice and can serve as general clues for checking the existence of a CANN in neural systems. In order to solve the dynamics which is otherwise extremely complicated for a large recurrent network, we have developed a perturbative analysis to simplify the dynamics of a CANN. Geometrically, it is equivalent to projecting the network state on its dominant directions of the state space. This method works efficiently and may be widely used in the study of CANNs. The special structure of a CANN may have other applications in brain functions, for instance, the highly structured state space of a CANN may provide a neural basis for encoding the topological relationship of objects in a feature space, as suggested by recent psychophysical experiments [15, 16]. It is likely that the distance between two memory states in a CANN defines the perceptual similarity between the two objects. Interestingly to note that the perceptual similarity measured by the psychometric functions of human subjects in a categorization task has a similar logarithimic nature as that of reaction times in a CANN [17]. To study these issues theoretically and justify the experimental findings, it is important for us to have analytic solutions of the state space and the dynamical behaviors of CANNs. We expect the analytical solution developed here will serve as a valuable mathematical tool. The tracking dynamics of a CANN has also been studied by other authors. In particular, Zhang proposed a mechanism of using asymmetrical recurrent interactions to drive the bump, so that the shape distortion is minimized [4]. Xie et al. further proposed a double ring network model to achieve these asymmetrical interactions in the head-direction system [8]. It is not clear how this mechanism can be generated in other neural systems. For instance, in the visual and hippocampal systems, it is often assumed that the bump movement is directly driven by external inputs (see, e.g., [5, 19, 20]), 7


and the distortion of the bump is inevitable (indeed the bump distortions in [19, 20] are associated with visual perception). The contribution of this study is on that we quantify how the distortion of the bump shape affects the network tracking performance, and obtain a new finding on the maximum trackable speed of the network. Finally, we would like to remark on the generality of the results in this work and their relationships to other studies in the literature. To pursue an analytical solution, we have used a divisive normalization to represent the inhibition effect. This is different from the Mexican-hat type of recurrent interactions used by many authors. For the latter, it is often difficult to get a closed-form of the network stationary state. Amari used a Heaviside function to simplify the neural response, and obtained the boxshaped network stationary state [2]. However, since the Heaviside function is not differentiable, it is difficult to describe the tracking dynamics in the Amari model. Truncated sinusoidal functions have been used, but it is difficult to use them to describe general distortions of the bumps [3]. Here, by using divisive normalization and the Gaussian-shaped recurrent interactions, we solve the network stationary states and the tracking dynamics analytically. One may be concerned about the feasibility of the divisive normalization. First, we argue that neural systems can have resources to implement this mechanism [7, 18]. Let us consider, for instance, a neural network, in which all excitatory neurons are connected to a pool of inhibitory neurons. Those inhibitory neurons have a time constant much shorter than that of excitatory neurons, and they inhibit the activities of excitatory neurons in a uniform shunting way, thus achieving the effect of divisive normalization. Second, and more importantly, the main conclusions of our work are qualitatively indpendent of the choice of the model. This is because our calculation is based on the fact that the dynamics of a CANN is dominated by the motion mode of position shift of the network state, and this property is due to the translational invariance of the neuronal recurrent interactions, rather than the inhibition mechanism. We have formally proved that for a CANN model, once the recurrent interactions are translationally invariant, the interaction kernel has a unit eigenvalue with respect to the position shift mode irrespective of the inhibition mechanism (to be reported elsewhere).

References
[1] P. Dayan and L. Abbott, Theoretical Neuroscience: Computational and Mathematical Modelling of Neural Systems, (MIT Press, Cambridge MA, 2001). [2] S. Amari, Biological Cybernetics 27, 77 (1977). [3] R. Ben-Yishai, R. Lev Bar-Or and H. Sompolinsky, Proc. Natl. Acad. Sci. USA, 92 3844 (1 9 9 5 ). [4] K.-C. Zhang, J. Neurosicence 16, 2112 (1996). [5] A. Samsonovich and B. L. McNaughton, J. Neurosci. 17, 5900 (1997). [6] B. Ermentrout, Reports on Progress in Physics 61, 353 (1998). [7] S. Deneve, P. Latham and A. Pouget, Nature Neuroscience, 2, 740 (1999). [8] X. Xie, R. H. R. Hahnloser and S. Seung, Phys. Rev. E 66, 041902 (2002). [9] A. Renart, P. Song and X. Wang, Neuron 38, 473 (2003). [10] C. Brody, R. Romo and A. Kepecs, Current Opinion in Neurobiology, 13, 204-211 (2003) [11] S. Wu and S. Amari, Neural Computation 17, 2215 (2005) [12] B. Blumenfeld, S. Preminger, D. Sagi and M. Tsodyks, Neuron 52, 383 (2006). [13] C. Chow and S. Coombes, SIAM J. Appl. Dyn. Sys. 5, 552-574, 2006. [14] J. Hopfield, Proc. Natl. Acad. Sci. USA, 79 2554 (1982). [15] J. Jastorff, Z. Kourtzi and M. Giese, J. Vision 6, 791 (2006). [16] A. B. A. Graf, F. A. Wichmann, H. H. Bulthoff, and B. Scholkopf, Neural Computation 18, ¨ ¨ 1 4 3 (2 0 0 6 ). [17] J. Zhang, J. Mathematical Psychology 48, 409 (2004) [18] D. Heeger, J. Neurophysiology 70, 1885 (1993). [19] M. Berry II, I. Brivanlou, T. Jordon and M. Meister, Nature 398, 334 (1999). [20] Y. Fu, Y. Shen and Y. Dan, J. Neuroscience 21, 1 (2001).

8