Equilibrium States and Their Stability in the Head-Direction Ring Network

Head-direction cells have been found in several areas in the mammalian brains. The firing rate of an ideal head-direction cell reaches its peak value only when the animal's head points in a specific direction, and this preferred direction stays the same regardless of spatial location. In this paper we combine mathematical analytical techniques and numerical simulations to fully analyze the equilibrium states of a generic ring attractor network, which is a widely used modeling framework for the head-direction system. Under specific conditions, all solutions of the ring network are bounded, and there exists a Lyapunov function that guarantees the stability of the network for any given inputs, which may come from multiple sources in the biological system, including self-motion information for inertially based updating and landmark information for calibration. We focus on the first few terms of the Fourier series of the ring network to explicitly solve for all possible equilibrium states, followed by a stability analysis based on small perturbations. In particular, these equilibrium states include the standard single-peaked activity pattern as well as double-peaked activity pattern, whose existence is unknown but has testable experimental implications. To our surprise, we have also found an asymmetric equilibrium activity profile even when the network connectivity is strictly symmetric. Finally we examine how these different equilibrium solutions depend on the network parameters and obtain the phase diagrams in the parameter space of the ring network.


INTRODUCTION
Head-direction cells were first reported in several brain areas related to the limbic system in the rodents (Taube, 2007) and later in other mammalian species such as monkeys (Robertson et al., 1999) and bats (Finkelstein et al., 2015). A stereotypical head-direction cell increases its firing rate when the animal's head is facing in a specific direction in a world-centered coordinate system regardless of the animal's spatial location, and the firing rate decreases to its baseline level as the animal's head turns away from the preferred direction (Taube et al., 1990). It has been proposed that the head-direction cells may form a ring network that allows an activity bump to be selfsustained by attractor dynamics, and the peak position of the activity bump is updated by selfmotion information and calibrated by learned landmarks (Skaggs et al., 1995;Redish et al., 1996;Zhang, 1996). Multiple versions of the ring network have been studied for the head-direction cells (Goodridge and Touretzky, 2000;Arleo and Gerstner, 2001;Sharp et al., 2001;Stringer et al., 2002;Xie et al., 2002;Song and Wang, 2005) as well as for a variety of applications beyond the original head-direction system (Ben-Yishai et al., 1995;Pouget et al., 1998;Hahnloser et al., 2000;Kakaria and de Bivort, 2017;Zhang et al., 2019). Besides head-direction cells, attractor networks have been used as a general theoretical framework for modeling other types of spatial cells in the hippocampus and related systems (Knierim and Zhang, 2012).
The equilibrium state of the head-direction ring network is often visualized as a single bump of activity whose peak position corresponds to the animal's current heading direction (Figure 1, top and middle rows). While this picture is compelling and highly intuitive, it is not the only theoretical possibility for explaining the experimental data. For instance, imagine that the ring network can sustain two activity bumps instead of one (Figure 1, bottom row), then if one records from an individual cell in the ring, one would still find a head-direction cell with a perfectly normal, single-peaked tuning curve, assuming that the activity bumps now rotates at half of the speed as the single activity bump. Indeed, if we focus on a single cell corresponding to north, we see that in both situations, the cell fires at maximal rate only when the animal is facing north (N).
Despite the functional equivalence, the structures of the two ring networks are different. More specifically, unlike the standard single-peaked network, the double-peaked network has strong connections between cells in opposing directions although they are not as strong as the connections between neighboring cells. In the consideration above, we assume that the rotation speed is halved in the double-peaked network. If the rotation speed is kept the same, then the network should generate tuning curves with two peaks that are 180 • apart. In fact, double-peaked head-direction tuning curves have been reported in the retrosplenial cortex (Jacob et al., 2017), although the phenomenon could be attributed to a single preferred FIGURE 1 | Head-direction ring network. In the classic view, the head direction of a rat (Top) is represented by the peak location of the activity bump (Middle: red shades) in a ring of head-direction cells. An alternative possibility is a ring network that allows two stable activity bumps that rotate at half of the speed (Bottom). direction switching back and forth in time rather than implying a truly double-peaked activity pattern (Page and Jeffery, 2018). The consideration above can be generalized readily to activity patterns with three or more peaks. The possible existence of multi-peaked activities in the head-direction system together with their potential functional significance has motivated us to examine the equilibrium states in the ring network model in greater depth. This paper is aimed at a thorough analysis of the equilibrium states in the ring network, with a focus on the exact conditions for the existence of activity patterns with multiple peaks. We will use the simple continuous ring network to simplify the mathematical analysis. The rotational symmetry of the system allows Fourier analysis to be used effectively. We strive to derive the exact analytical conditions whenever possible, and the analytical treatments are complemented by systematic numerical simulations. Once the exact expressions of all different kinds of equilibrium states are obtained, we will employ small perturbations and eigenanalysis of the linearized system to determine the stabilities of these equilibrium states. We will examine the dependence of various equilibrium states on the network parameters and summarize the results by the phase diagrams. Our analysis may provide a necessary step for extending the application and analysis of the ring network beyond the classic singlepeaked condition.

MATERIALS AND METHODS
We consider a continuous formulation of the head-direction system which has a continuous ring structure (Zhang, 1996). Such continuous formulation has a long history in neural modeling (Wilson and Cowan, 1972;Amari, 1977;Bressloff, 2012). The standard simplified time evolution continuous dynamics is governed by the equation (1) where the convolution is defined by (2) In this system, u(θ , t) represents the state of voltage of a unit with θ as its preferred direction, w(θ − ϕ, t) represents the synaptic weight between units with θ − ϕ being the difference of their preferred directions, g(u) is a monotone increasing and sigmoid gain function, I(θ , t) represents external inputs, τ is a time constant, θ is head-direction, and variable t is time. So the whole head-direction system is given by 3. RESULTS

Boundedness of Solutions
According to the form of Equation (3), we multiply integrating factor e t τ of both sides of this equation. It is easy to get the solution of Equation (3) as follows: Obviously, as t → +∞ the first term u(θ , 0)/e t τ tends to zero. If w(θ , t) and g(u) are bounded function, then the second term is bounded. This is because: where M is the maximum value of the integrand. At the same time we have: So when t → +∞ the second term of solutions (4) is bounded. In general the external inputs I(θ , t) is also bounded, thus we know that the third term of solutions (4) is bounded, too. Therefore if synaptic weight w(θ , t), gain function g(u) and external inputs I(θ , t) are bounded, then all solutions of system (3) are tending to a bounded domain rather than wandering in the whole space. Namely as variable t tends to infinity, the average net state u(θ , t) for each head-direction cell is bounded and changes in a bounded domain. Meanwhile we find that if there is no external input (i.e., I(θ , t) = 0), then |u(θ , t)| is less than the maximum of the product of |w(θ , t)| and |g(u(θ , t))| as t → +∞. So the output of system (3) is under control by synaptic weight and gain function.
(10) Notice that when external input I(θ , t) is a constant and the Fourier series coefficients of synaptic weight w(θ , t) just have finite term, i.e., then when |n| > m the Fourier series coefficients of solution u(θ , t) is: i.e., Obviously the Fourier series coefficients of solution u(θ , t) will reduce to zero (t → +∞) when |n| > m. When the synaptic weight w(θ , t) and sigmoid function g(u) are chosen in some special forms, we can derive the form and properties of the solutions, especially the equilibrium solutions.

The Equilibrium Solutions
Until now, as far as we know the general solutions of the ring attractor network (3) can not be solved, but we can get special solutions, such as equilibrium solutions. Once the equilibrium states are determined, we can further obtain the properties of other solutions near the equilibrium solutions by local stability analysis. Sometimes we even determine the tendency of all solutions in the solution space. Base on the definition of an equilibrium solution, we let This equation shows that the equilibrium solutions is independent of time t. In other words, the activity u(θ , t) does not change with time. In the head-direction neural network, different equilibrium solutions are related to different equilibrium states. We write u(θ , t) = u(θ ) which depends only on variable θ . An equilibrium solution satisfies: Here we assume I(θ , t) represents a fixed external input current (i.e., I(θ , t) = I), so the system has equilibrium solutions if and only if the synaptic weight w(θ , t) is also independent of time t. That means w(θ , t) = w(θ ) and the equilibrium solution satisfies As mentioned in section 3.1.2 we know that w(θ , t), u(θ , t) and I(θ , t) are periodic with period 2π on θ . According to Equation (16) we conclude that the equilibrium solution u(θ ) is always rotation-invariant. It means that if u(θ ) is an equilibrium solution of system (3), then u(θ − θ 0 ) is also an equilibrium solution. Therefore from the viewpoint of symmetry, we show that every head-direction cell in the ring network has similar equilibrium states. A mathematical verification of this conclusion is in the next paragraph. Because u(θ ) is an equilibrium solution which satisfies (Equation 16) for any θ 0 , we have: Set ϕ = φ − θ 0 and plug this relation into the integration, then we get: That is, Thus u(θ − θ 0 ) also is an equilibrium solution. According to this property there is no need to consider any shift of head direction, and we only need to focus on the mathematical forms all the equilibrium solutions.
Therefore we find that the form of the equilibrium states of the ring network depends heavily on the form of the synaptic weights.

Lyapunov Function
In this section we consider the stability of system (14) by constructing a continuous version of a Lyapunov function for symmetric networks (Cohen and Grossberg, 1983;Hopfield, 1984). The Lyapunov function or energy function is as follows: Obviously, the energy function is bounded and the time derivative of Equation (25) is: We find if w(θ ) is even, i.e., w(θ − ϕ) = w(ϕ − θ ), then Thus its time derivative becomes: Since V = g(u) is a monotonically increasing gain function, its inverse function u = g −1 (V) is also a monotonically increasing function, i.e., Therefore when the synaptic weight w(θ ) is even, the gain function g(u) is monotonically increasing, and the synaptic weight and the gain function are all bounded, all solutions of system (3) are convergent to the corresponding equilibrium states as t → +∞. If the system has one, and only one, equilibrium solution, then this equilibrium solution must have global stability. That means all flows of system (3) converge to this stable state no matter what the initial state is.

An Example of Head-Direction Ring Network
Based on the above analysis, in order to get all equilibrium solutions we just need to pay attention to the Fourier form of the synaptic weight. Here we choose the synaptic weight as w(θ ) = a + b cos θ + c cos 2θ which only has three terms. Of course, the conclusions and methods can be extended to more general cases as long as the synaptic weight has finite terms.
According to the analysis in section 3.1.3, we know that the equilibrium solutions of system have the following form: (30) Once synaptic weight is chosen, we can determine the form of equilibrium solutions. According to the solution (29), the equilibrium solution u(θ ) has a similar form as synaptic weight w(θ ). Therefore if the synaptic weight has no more than two peaks, then all equilibrium solutions have no more than two peaks. That is, an equilibrium of the system may be a flat solution, a single-peaked solution, or a double-peaked solution, and on the basis of relationship (30) all equilibrium solutions are dependent on control parameters a, b and c. In fact we can obtain all the equilibrium solutions when the sigmoid function g(u) is also chosen.
The gain function is often described by the sigmoid: where k > 0 is the gain and u 0 is the threshold. For convenience in this paper we set the threshold u 0 = 0. When the gain k is larger enough or as k → +∞, the form of the gain function converges to the Heaviside step function, and the derivative of the gain function reduces to the Dirac δ function.
Here we choose I(θ , t) = 0, g(u) =H(u) =Heaviside(u), and w(θ ) = a + b cos θ + c cos 2θ . Now the ring network becomes: The equilibrium solution equation becomes: (33) Parameter a only moves the equilibrium state u(θ ) up and down, without changing its shape. So we set parameter a = 0 in this paper. Figure 2A shows synaptic weight function w(θ ) and Figure 2B shows the weight matrix of the ring network, with parameters a = 0, b = 3 and c = 2.

All Possible Equilibrium States
The general form of the equilibrium solutions in system (32) is: Using basic trigonometric formula, we rewrite this form as: and 2θ 2 = arctan c 2 c 1 . As mentioned in section 3.1.3 for any θ 0 if u(θ ) is a solution of system (32), then u(θ − θ 0 ) is also a solution. So we should pay attention only to the form of the equilibrium solutions and ignore the influence of the phaseshift θ 0 . Meanwhile according to Equation (35) we know that any equilibrium solution is a linear combination of functions cos(θ − θ 1 ) and cos 2(θ − θ 2 ). Next we will discuss all possible linear combinations of these three terms and solve for the exact equilibrium states respectively.
The first situation is that the equilibrium solution is a constant. If u(θ ) = A is a solution of system (32), then Since a = 0, the constant A = a 2π 0 g(u(ϕ))dϕ is always equal to zero. In other words, we have u(θ ) = 0 as the only constant solution of the system. The general form of non-zero equilibrium solution for head-direction neural network (32) is: where B and C are unknown constants and θ 1 and θ 2 are unknown angles.
Due to rotation-invariance, we just need to find one solution of this form u(θ ) = B cos θ (B > 0). Plugging u(θ ) = B cos θ (B > 0) into (Equation 33) we have the following equality: Figure 3A. Figure 3B shows numerical simulation with initial state u 0 is 0.1 b π cos(θ − π) with additional small perturbations, and as t → +∞, the state approaches this equilibrium solution, suggesting that the equilibrium solution u(θ ) = b π cos(θ − θ 0 ) is likely to be stable. We will return this stability problem in section 3.2.2.
perturbation. The result indicates that this equilibrium solution is probably stable.
When the equilibrium solutions have form B cos(θ − θ 1 ) + C cos 2(θ − θ 2 ), we change the form of equilibrium solution as follows: where φ = θ − θ 1 and φ 0 = θ 2 − θ 1 . Due to the property of equilibrium solutions we just need to find the basic equilibrium solutions By numerical simulation we find that the phase diagram of u(φ) has two situations as parameters B, C and φ 0 are changing.
One situation is there is only one positive domain above φ-axis in period 2π, another situation is there are two positive domains above φ-axis in period 2π. Because positive domain decides the value of Heaviside function, we have to discuss these two cases respectively. In this section we mainly consider the first situation in which u(φ) just has one positive domain above φ-axis in period 2π. that is shown in Figure 5A.
As a equilibrium solution, we plug u i.e., For any value φ (Equation 43) must hold, which means that the corresponding coefficients must be equal. So we get sin( α ′ +α 2 ) = 0, i.e., sin(α ′ + α) = 0 and sin 2φ 0 = 0. This result shows that in this case if the system has an equilibrium solution, then its form must be u(φ) = B cos φ ± C cos 2φ. And the parameters B and C must satisfy the following relationship: 4c 2 > 0, and the parameter sb and c must satisfy 2b ≥ c > 0. Comprehensively, in the first situation if the system has the equilibrium solution u(φ) = B cos φ ± C cos 2φ, then parameters b and c must satisfy the condition 0 < b < c ≤ 2b. So we have the following conclusion: If 0 < b < c ≤ 2b, then for any θ 0 the system has an equilibrium solution of the form: For parameters a = 0, b = 1 and c = 1.5, the profile of the synaptic weight is shown in Figure 5B. By numerical simulation we find two of this kind solutions as shown in Figures 5C,E, while Figures 5D,F are time evolution with initial state u 0 = 0.1u(θ ) with additional small perturbations. The results suggest that this type of solutions may be stable too. Now we consider the second situation two domains above φaxis in period 2π as shown in Figure 6A. As above analysis we know the basic form of this kind equilibrium solution is u(φ) = B cos φ + C cos 2(φ − φ 0 ). Plug u(φ) into the equilibrium solution Equation (33), we get Compare the corresponding relationship we have cos α −cos α ′ + cos β − cos β ′ = 0, where α, α ′ , β and β ′ are four roots of equation u(φ) = 0. That means we have to solve equation u(φ) = B cos φ + C cos 2(φ − φ 0 ) = 0. Using trigonometric formulas to expand u(φ) = 0 the equation u(φ) = 0 becomes: B cos φ + C(2 cos 2 φ − 1) cos 2φ 0 + 2C sin φ cos φ sin 2φ 0 = 0.
If 0 < b < 2c and b c = 1 2 , then for any θ 0 the system has equilibrium solutions which are where When parameters set a = 0, b = 3 and c = 2 and 0 < b < 2c, Figures  Obviously, these two kinds of equilibrium solutions keep immovability when t tends to infinity.
Solve above equation and get the roots cos α = 0, sin α = 0 and cos 2 α = b 2c . For B > 0 and C > 0 we choose cos 2 α = b 2c , other roots are given up. Of course, the parameter must meet 0 < b 2c < 1, i.e., 0 < b < 2c. Therefore we have sin α = 2c−b 2c . Now we get another form equilibrium solution which is In a words, if 0 < b < 2c, then for any θ 0 the system has equilibrium solutions which are: When parameters are chosen a = 0, b = 3 and c = 2, which satisfy 0 < b < 2c, you must get this equilibrium solution. Figure 8A is time evolution of equilibrium solutions (63).
(4) 2φ 0 = 3π 2 . Actually when 2φ 0 = 3π 2 the analysis and calculating process is in the same way as 2φ 0 = π 2 . Here we ignore the complex calculating process and give the follow conclusion: If 0 < b < 2c, then for any θ 0 the system has equilibrium solution which are We find that this kind solution is mirror-symmetry with equilibrium solution (63). Actually in equilibrium solution (63) if we replace θ − θ 0 with −(θ − θ 0 ), then we get equilibrium solution (64). So we can image the graph of equilibrium  solution (64) according to equilibrium solution (63). To our surprise, we find that equilibrium solution (63) and equilibrium solution (64) are two asymmetric equilibrium activity pattern even when the network connectivity pattern is strictly symmetric.
(66) So as to further consider the perturbation equation of equilibrium solutions. We plug Taylor's expansion of g(u(θ , t)) into Equation (65), then we have Since u(θ ) is the equilibrium solution, we have U = 1 n WG 1 . Thus head-direction neural network (67) can be simplified as follows: and the linear part is where J = −I+ 1 n WG is the Jacobian matrix and I is n×n identity matrix. The general solution to this linear ordinary differential equations is: E = c 1 e λ 1 t V 1 + c 2 e λ 2 t V 2 + · · · + c n e λ n t V n , where λ 1 , λ 2 , · · · , λ n are eigenvalues of the Jacobian matrix, V 1 , V 2 , · · · , V n are corresponding eigenvectors, and c 1 , c 2 , · · · , c n are any constant. It is easy to obtain that the stability of perturbation (Equation 69) depends on the eigenvectors of the Jacobian matrix. Since synaptic weight function w(θ ) is even and periodic, the corresponding matrix W is a circulatory and symmetric matrix. According to properties of a circulatory and symmetric matrix (Gray, 2005;Tee, 2007) we can get the normalized and orthogonal eigenvectors of 1 n W. They are where ρ j = exp( 2π i n j) (j = 0, 1, · · · , n − 1) are the n-th roots of unity and i is the imaginary unit. The corresponding eigenvalues6 are given by Since thus ρ n−k j + ρ k j = ρ k j + ρ k j = 2 cos( 2kπ n j), j = 0, 1, · · · , n − 1. (75) So we change the form of eigenvalues as follows: i.e., λ j = 1 n n−1 k=0 w(k 2π n ) cos(k 2π n j), j = 0, 1, · · · , n − 1.
(88) So we can further determine the stability of other equilibrium solutions.
When the equilibrium is u Therefore characteristic equation is For convenient we set k = c b , the Equation (97) becomes Thus the eigenvalues are λ 1 = · · · = λ n−2 = −1, λ n−1,n = −1 Obviously when parameters satisfy 0 < b < c ≤ 2b, i.e., 1 < k ≤ 2, all eigenvalues are negative. So the equilibrium cos 2θ is stable. By the same analysis method we also get that solution u cos 2θ is stable too. According to the Equation (86), we know that at least n − 2 eigenvalues of the four remaining equilibrium solutions are equals to −1 and at most two eigenvalues are different. So we just need to consider these different eigenvalues in each case. And the two different eigenvalues satisfy equation thus all the eigenvalues are: λ 1 = · · · = λ n−2 = −1, λ n−1,n = −1 + Obviously if λ max < 0, all eigenvalues are negative, so that the corresponding equilibrium solution is stable. If λ max > 0, at least one eigenvalue is positive, so that the corresponding equilibrium solution is unstable. Therefore putting the equilibrium solution into Equation (100) allows allows us to determine the stability of the equilibrium solution.
For example, in the equilibrium solution u(θ ) = b π 2c−b 2c cos(θ − θ 0 ) + b 2π sin 2(θ − θ 0 ) (63) we set the parameters as a = 0, b = 3, c = 2, and θ 0 = 0 to obtain u(θ ) = 3 2π (1 + 2 sin θ ) cos θ . Plugging this equilibrium solution into (Equation 100), we find that the maximum eigenvalues of the equilibrium solutions are positive. This means that the equilibrium solution (63) are unstable. That is, adding small perturbations to the equilibrium solutions (63) will make the state moving away from the equilibrium solutions (63), as shown in Figure 8B. In fact the state approaches other stable equilibrium solution. As illustrated in Figure 8B we find that the state converges to the one-peaked stable equilibrium solution given by u θ = 3 π cos(θ − θ 0 ). This result is confirmed in Figure 8C where the blue line is the equilibrium solution (63), the red line represents the state of the system after 100, 000 time steps starting from the equilibrium solution (63) as the initial state, and the green star line is the one-peaked stable equilibrium solution u(θ ) = 3 π cos(θ − 6.5π 42 ). We can see that the system starts from an unstable equilibrium state but converges gradually to this stable equilibrium.

Phase Diagram in the Parameter Space of the Ring Network
So far we have obtained all kinds of equilibrium solutions and determined their stability when the synaptic weight is w(θ ) = a+b cos θ +c cos 2θ and sigmoid function is g(u) = Heaviside(u). We find that the form of these solutions are strongly dependent on weight coefficients a, b and c. The values of these parameters determine not only the form of the equilibrium solutions, but also their stability. Setting parameter a = 0, we summarize the results in section 3.2.1 with the following conclusion. Proposition: (i) For any b and c, system (32) has one, and only one, constant equilibrium solution u(θ ) = 0; (ii) When b > 0, system (32) has symmetric and one-peaked equilibrium solutions u(θ ) = b π cos(θ − θ 0 ); (iii) When c > 0, system (32) has symmetric and two-peaked equilibrium solutions u(θ ) = c π cos 2(θ − θ 0 ); (iv) When 0 < b < c < 2b, the system (32) has the equilibrium solutions (v) When 0 < b < 2c, the system (32) has four kinds of special equilibrium solutions: and The distribution of all possible equilibrium solutions in system (32) is shown in Figure 9.

Shifting Mechanism
In head-direction ring network, one important biological characteristics is that the attractor bumps can shift in time in response to a head turn. In Figure 10 we select N = 500 headdirection cells, and time step t = 0.01, and use numerical simulation to demonstrate the shifting mechanism. Following the derivative rule (Zhang, 1996), now the synaptic weight becomes: where the original synaptic weight is w(θ ) = a + b cos θ + c cos 2θ , α is the speed of shifting. Here we choose α = 0.2 and initial values are equilibrium solutions plus small disturbances. Figures 10A,B show the shifting mechanism for the one-peaked attractor state u(θ ) = b π cos θ and the two-peaked attractor state u(θ ) = b π cos 2θ , with parameters a = 0, b = 3, and c = 2. Figure 10C shows the shifting mechanism for the stable special two peaked attractor state u cos 2θ , with parameters a = 0, b = 1 and c = 1.5.

DISCUSSION
We have analyzed the equilibrium states of the head-direction ring network and found multiple solutions. Even for the simplest network with step gain function and only two Fourier terms in the weight distribution profile, there are a rich variety of equilibrium solutions. Some of the equilibrium solutions are well  known, such as the flat solution and the single-peaked solution, while other solutions are unexpected, such as the asymmetric solutions. In particular, the equilibrium states with two peaks can be generated under many parameter combinations. A necessary condition for the two peaked solution is that the weight profile must have a Fourier component with two peaks. Our analysis reconfirms that it is possible for the ring network to have two stable activity bumps as illustrated in Figure 1. To determine how these different types of equilibrium states depend on the parameters, we have calculated the phase diagram in the parameter space of the ring network with step gain function. Our method can be extended other gain functions such as the standard sigmoidal gain function, as discussed below.
The simple head-direction ring network has some essential dynamical features such as boundedness of state, convergence to stable equilibrium states, and strong dependence of the equilibrium states on the synaptic weight. In this paper our analytical treatment relies on several simplifying assumptions such as truncated Fourier series and step gain function. To generate the equilibrium states, we assume an even weight function w(−θ ) = w(θ ) which is equivalent to symmetric reciprocal connection weights because the existence of a Lyapunov function guarantees the stability in this situation. We have only briefly considered the case of asymmetric weights by adding a derivative of the weight profile in order to shift the activity bumps. The parameters in synaptic weight w(θ ) are the main control parameters that determine the form of an equilibrium state and its stability. As mentioned in section 3, parameter a determines the position of the constant solution, parameter b determines the existence of one-peaked equilibrium solutions, parameter c determines the existence of two-peaked equilibrium solutions, and the values of parameters a, b and c together determine the exact form of an equilibrium state as well as its stability. We find that if the head-direction ring network has one equilibrium solution, then there must exist at least one stable equilibrium solution, i.e., one attractor solution.
The step gain function or Heaviside function may be regarded as the limit of the standard sigmoidal gain function as the slope approaches infinity. When the slope of the sigmoidal gain function is not too small, the behaviors of the ring network are qualitatively quite similar to the network with Heaviside function. For example, consider a ring network with the gain function g(u) = 1 1+e −ku with k = 2 and the weight function w(θ ) = a + b cos θ + c cos 2θ . The parameters a, b and c also determine which kind of equilibrium states the ring network has as well as their stabilities in a manner similar to the network with Heaviside function. The equilibrium state of the system may allow a flat solution, a one-peaked solution and a two-peaked solution, as shown in Figure 11. For k = 2, a = 0, b = 3.5, and c = 3.5, the equilibrium state is constant 0, as shown Figure 11A. Actually because parameter a = 0, we know that u(θ ) = 0 is the one and the only constant solution of the system. Figure 11B shows the final state reaching a one-peaked equilibrium state for parameters a = 0, b = 4.5, and c = 3.5. Figure 11C shows a two-peaked equilibrium solution for parameters a = 0, b = 3.5 and c = 4.5. In addition, according to section 3.2.2 we obtain that when g ′ (0)b 2 < 1 and g ′ (0)c 2 < 1, i.e., b < 8 k and c < 8 k , the constant solution u(θ ) = 0 is stable. Furthermore b = 8 k and c = 8 k demarcate the boundary of different parameter domains which contain different types of equilibrium solutions. Figure 12A shows the parameter space for a = 0 and k = 2. Such phase diagram is obtained by repeated numerical simulation with random initial states. The one-peaked state u 0 = cos θ , the two-peaked state u 0 = cos 2θ , and the flat state u(θ ) = 0 are all possible equilibrium state, depending on the parameters. In the red domain, u(θ ) = 0 is the FIGURE 11 | Equilibrium states in the ring network with sigmoidal gain function g(u) = 1 1+e −2u . The number of head-direction cells N = 50, τ = 1, and time step t = 0.001. The synaptic weight function is also w(θ) = b cos θ + c cos 2θ and parameter a = 0. (A) When b = 3.5 and c = 3.5, starting from a random initial state, the final state converges to a flat solution u(θ) = 0. (B) When b = 4.5 and c = 3.5, starting from a random initial state, the final state converges to a one-peaked equilibrium solution. (C) When b = 3.5 and c = 4.5, starting from a random initial state, the final state converges to a two-peaked equilibrium solution.
FIGURE 12 | The phase diagram for the distribution of equilibrium states in the ring network with the sigmoidal gain function g(u) = 1 1+e −ku . The number of head-direction cells N = 50, τ = 1 and time step t = 0.001. The synaptic weight function is still w(θ) = b cos θ + c cos 2θ and parameter a is 0. The equilibrium state of the system is strongly dependent on parameters b, c and k. (A) The parameter space of head-direction ring network when k is set to 2 and the initial values are random, u 0 = cos θ and u 0 = cos 2θ. The red domain represents flat equilibrium states, the blue domain represents one-peaked equilibrium states, the green domain represents two-peaked equilibrium states, and black domain allows both one-peaked equilibrium states and two peaked equilibrium states. (B) The parameter space with different gain k, where k = 2, k = 4, k = 8, and k = +∞ and the initial state is same as (A). stable equilibrium state. In the blue domain, the state of system converges to the one-peaked equilibrium state. While in the green domain, the state of system converges to the two-peaked equilibrium state, and in black domain the state may converge to either a one-peaked equilibrium state or a two-peaked state, depending on the initial state. In other words, gain k determines the boundary of different domains. The phase diagrams for different gain k are shown in Figure 12B, where are four cases with k = 2, k = 4, k = 8, and k = +∞. We can see that phase diagram changes gradually as the gain slope k changes.
Our analysis reveals a diverse set of equilibrium states of the ring network. Although an unstable equilibrium state is not as robust as a stable equilibrium state, it might be useful for generating slow dynamics that slows down around these special states. The existence of stable activity pattern with multiple peaks provides a theoretical foundation for future study of the headdirection system so that new data analysis methods and new experimental designs could be developed to distinguish different computational mechanisms.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/supplementary material.

AUTHOR CONTRIBUTIONS
This paper was completed and written by CW under the guidance of KZ.

FUNDING
KZ was supported partially by NIH grant U01NS111695.