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

1. 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 self-sustained by attractor dynamics, and the peak position of the activity bump is updated by self-motion 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 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 eigen-analysis 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 single-peaked condition.

2. 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

where the convolution is defined by

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

3.1. General Properties of Head-Direction Ring Network

3.1.1. Boundedness of Solutions

According to the form of Equation (3), we multiply integrating factor $e^{\frac{t}{\tau }}$ 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(\theta ,0)/e^{\frac{t}{\tau }}$ 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.

3.1.2. The Form of Solutions

The structure of head-direction model is a homogeneous and continuous ring network, so the synaptic weight w(θ, t) and external input I(θ, t) are periodic with period 2π. We have the Fourier series expansions of w(θ, t) and I(θ, t) as follows:

According to the form of solutions (4), we know u(θ, t) and g(u(θ, t)) are also periodic with period 2π on θ. So the Fourier series expansions of u(θ, t) and g(u(θ, t)) are given as follows:

According to the convolution theorem, we know that if g(u(θ, t)) and w(θ, t) are in L¹([−π, π]), then the Fourier series coefficients of 2π-periodic convolution of w(θ) and g(u(θ)) are given by:

Therefore the Fourier series coefficients of system (3) have the following relationships:

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.

3.1.3. 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(θ − θ₀) 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 θ₀, we have:

Set φ = ϕ − θ₀ and plug this relation into the integration, then we get:

That is,

Thus u(θ − θ₀) 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.

Since constant I can be absorbed into the constant term of w(θ), we will set I = 0 in the following analysis. Now we choose w(θ) as a Fourier series with finite terms:

and then an equilibrium solution satisfies

Since

where

the equilibrium solution u(θ) has a similar form as w(θ), namely,

Therefore we find that the form of the equilibrium states of the ring network depends heavily on the form of the synaptic weights.

3.1.4. 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⁻¹(V) is also a monotonically increasing function, i.e., $\frac{d{g}^{-1}(V)}{dV}>0$. Since the time parameter τ > 0, we obtain $\frac{dE}{dt}\le 0$, where $\frac{dE}{dt}=0$ if and only if u(θ, t) = u(θ) which is an equilibrium state.

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.

3.2. 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:

where,

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₀ is the threshold. For convenience in this paper we set the threshold u₀ = 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:

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.

3.2.1. All Possible Equilibrium States

The general form of the equilibrium solutions in system (32) is:

Using basic trigonometric formula, we rewrite this form as:

where A = a₀, $B=\sqrt{b_1^2+b_2^2}$, $C=\sqrt{c_1^2+c_2^2}$, ${\theta }_1=arctan\frac{b_2}{b_1}$ and $2{\theta }_2=arctan\frac{c_2}{c_1}$. As mentioned in section 3.1.3 for any θ₀ if u(θ) is a solution of system (32), then u(θ − θ₀) is also a solution. So we should pay attention only to the form of the equilibrium solutions and ignore the influence of the phase-shift θ₀. Meanwhile according to Equation (35) we know that any equilibrium solution is a linear combination of functions cos(θ − θ₁) and cos 2(θ − θ₂). Next we will discuss all possible linear combinations of these three terms and solve for the exact equilibrium states respectively.

Case 1: u(θ) = A.

The first situation is that the equilibrium solution is a constant. If u(θ) = A is a solution of system (32), then

i.e.,

Since a = 0, the constant $A=a{\int }_0^{2\pi }g(u(\phi ))d\phi$ 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 θ₁ and θ₂ are unknown angles.

Case 2: u(θ) = B cos(θ − θ₀), B > 0.

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:

i.e.,

So if b > 0, then $u(\theta )=\frac{b}{\pi }\text{cos}(\theta -{\theta }_0)$ is the equilibrium solution for any θ₀; If b < 0, then there is no equilibrium solution like this form. For parameters a = 0, b = 3 and c = 2, we show the solution $u(\theta )=\frac{b}{\pi }\text{cos}(\theta -\pi )$ in Figure 3A. Figure 3B shows numerical simulation with initial state u₀ is $0.1\frac{b}{\pi }\text{cos}(\theta -\pi )$ with additional small perturbations, and as t → +∞, the state approaches this equilibrium solution, suggesting that the equilibrium solution $u(\theta )=\frac{b}{\pi }\text{cos}(\theta -{\theta }_0)$ is likely to be stable. We will return this stability problem in section 3.2.2.

Case 3: u(θ) = C cos 2(θ − θ₀), C > 0.

Similar to case 2, we just need to find one special equilibrium solution u(θ) = C cos 2θ and its coefficient C > 0. After that we can obtain all the equilibrium solutions with the form u(θ) = C cos 2(θ − θ₀)(C > 0). Plugging u(θ) = C cos 2θ into (Equation 33) we find the following equality:

So we have the following conclusion. If c > 0, then $u(\theta )=\frac{c}{\pi }\text{cos }2(\theta -{\theta }_0)$ is the equilibrium solution for system (32) for any θ₀; If c < 0, then there is no equilibrium solution like this form. Figure 4A shows solution $u(\theta )=\frac{c}{\pi }\text{cos}(\theta -\frac{\pi }{2})$ for parameters a = 0, b = 3 and c = 2. Figure 4B shows the time evolution with initial state $u_0=0.1\frac{c}{\pi }\text{cos }2(\theta -\frac{\pi }{2}))$ with additional small perturbation. The result indicates that this equilibrium solution is probably stable.

Case 4: u(θ) = B cos(θ − θ₁) + C cos 2(θ − θ₂), B > 0 and C > 0.

When the equilibrium solutions have form B cos θ − θ₁) + C cos 2(θ − θ₂), we change the form of equilibrium solution as follows:

where ϕ = θ − θ₁ and ϕ₀ = θ₂ − θ₁. 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 ϕ₀ 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(ϕ) = B cos ϕ + C cos 2(ϕ − ϕ₀) into Equation (33) to obtain

i.e.,

For any value ϕ (Equation 43) must hold, which means that the corresponding coefficients must be equal. So we get $sin(\frac{{\alpha }^{\prime }+\alpha }{2})=0$, i.e., sin(α′ + α) = 0 and sin 2ϕ₀ = 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:

Comparing the corresponding coefficients we have $B=\frac{b}{\pi }sin\alpha$ and $\pm C=\frac{c}{2\pi }sin2\alpha$. Since α is a root of equation u(ϕ) = 0, next we should solve equation u(ϕ) = 0. We find sin α = 0, cos α = 0 and $\text{cos }2\alpha =-\frac{b}{c}$. Because B > 0 and C > 0, we keep $\text{cos }2\alpha =-\frac{b}{c}$ and reject the others. If $\text{cos }2\alpha =-\frac{b}{c}$ is a valid root, we have $\left|\frac{b}{c}\right|<1$, i.e., |b| < |c|. So we obtain $sin\alpha =\sqrt{\frac{c+b}{2c}}$ and $\text{cos }\alpha =\pm \sqrt{\frac{c-b}{2c}}$. At the same time we have $B=\frac{b}{\pi }\sqrt{\frac{c+b}{2c}}$ and $C=\frac{c}{\pi }\sqrt{\frac{c^2-b^2}{4c^2}}$. For the equations to hold there is only one positive domain in period 2π, we have B ≥ C > 0, i.e., $\frac{b}{\pi }\sqrt{\frac{c+b}{2c}}\ge \frac{c}{\pi }\sqrt{\frac{c^2-b^2}{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 θ₀ 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.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(ϕ − ϕ₀). 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. Using trigonometric formulas to expand u(ϕ) = 0 the equation u(ϕ) = 0 becomes:

For convenience let x = cos ϕ and $k=\frac{B}{C}$ we square both side of the Equation (47), then the Equation (47) becomes

i.e.,

According to the four roots of Equation (49), we can find that if and only if 2ϕ₀ is equals to $0,\frac{\pi }{2},\pi$ and $\frac{3\pi }{2}$, then the roots of Equation (47) hold the corresponding relation cos α − cos α′ + cos β − cos β′ = 0. This result can be seen in Figure 6B. Four solid lines represent four roots of Equation (49), and two dash lines are all possible combination of cos α − cos α′ + cos β − cos β′. According to Figure 6B we find that if and only if 2ϕ₀ equals 0, $\frac{\pi }{2}$, π and $\frac{3\pi }{2}$, then four roots of Equation (47) satisfy cos α − cos α′ + cos β − cos β′ = 0, where k = 0.8. By using numerical simulation we find when k is selected other positive number, although the amount of wing flexing for all lines are changing, the positions of all points of intersection on 2ϕ₀-axis are immovability. And the change of picture is successive. That means if and only if 2ϕ₀ is equals to $0,\frac{\pi }{2},\pi$ and $\frac{3\pi }{2}$, then the system may be exist equilibrium solutions which have two domains above ϕ-axis. Let us consider these four situations one by one.

(1) 2ϕ₀ = 0.

If 2ϕ₀ = 0, then system (32)may be have one kind of equilibrium solutions which like u(ϕ) = B cos ϕ + C cos 2ϕ which has two positive domain in period 2π. In this case the roots of u(θ) = 0 are −α, α, β and 2π − β. So plug u(ϕ) = B cos ϕ + C cos 2ϕ into (Equation 33), and then we obtain the relationship as follow

Since −α, α, β and 2π − β are four roots of equation u(θ) = 0, thus we should get the following the following equalities:

Let the first equation of (51) minus the second equation of (51), we have

Since cos 2β − cos 2β = 2(sin β − sin α)(sin β + sin α) and sin α ≠ sin β, thus the above equality becomes:

For −α, α, β and 2π − β are four roots of equation u(ϕ) = B cos ϕ + C cos 2ϕ = 0, we get:

So we have $\text{cos }\alpha ·\text{cos }\beta =-\frac{1}{2}$. Plug $\text{cos }\beta =-\frac{1}{2\text{cos }\alpha }$ into (Equation 53) and compute cos α. In order to keep meaningful during computational process we set 0 < b < 2c and $\frac{b}{c}\ne \frac{1}{2}$, then we have $\text{cos }\alpha =\sqrt{\frac{s-\sqrt{s^2-1}}{2}}$, where $s=\frac{2-\frac{b}{c}+\sqrt{\frac{2b}{c}}}{2}$. Therefore we have the following conclusion:

If 0 < b < 2c and $\frac{b}{c}\ne \frac{1}{2}$, then for any θ₀ the system has equilibrium solutions which are

where

(2) 2ϕ₀ = π.

When 2ϕ₀ = π the analyzing and calculating process is the same as 2ϕ₀ = 0. Because of this we omit the complex calculation and present the conclusion directly.

If 0 < b < 2c and $\frac{b}{c}\ne \frac{1}{2}$, then for any θ₀ the system has equilibrium solutions which are

where

When parameters set a = 0, b = 3 and c = 2 and 0 < b < 2c, Figures 7A,C show the equilibrium solution (55) and equilibrium solution (57). Meanwhile we can see the time evolution of solution (55) and solution (57) from Figures 7B,D. Obviously, these two kinds of equilibrium solutions keep immovability when t tends to infinity.

(3) $2{\varphi }_0=\frac{\pi }{2}$.

If $2{\varphi }_0=\frac{\pi }{2}$, then the equilibrium solutions of system (32) like u(ϕ) = B cos ϕ + C sin 2ϕ. In this case the four roots of u(ϕ) = 0 are −α, $\frac{\pi }{2}$, π + α and $\frac{3\pi }{2}$. As one kind of equilibrium solution we Plug u(ϕ) = B cos ϕ + C sin 2ϕ into equation and then we obtain the relation as follow:

Since −α is root of equation u(θ) = 0, then we have:

i.e.,

Solve above equation and get the roots cos α = 0, sin α = 0 and ${\text{cos}}^2\alpha =\frac{b}{2c}$. For B > 0 and C > 0 we choose ${\text{cos}}^2\alpha =\frac{b}{2c}$, other roots are given up. Of course, the parameter must meet $0<\frac{b}{2c}<1$, i.e., 0 < b < 2c. Therefore we have $\text{sin}\alpha =\sqrt{\frac{2c-b}{2c}}$. Now we get another form equilibrium solution which is

In a words, if 0 < b < 2c, then for any θ₀ 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{\varphi }_0=\frac{3\pi }{2}$.

Actually when $2{\varphi }_0=\frac{3\pi }{2}$ the analysis and calculating process is in the same way as $2{\varphi }_0=\frac{\pi }{2}$. Here we ignore the complex calculating process and give the follow conclusion:

If 0 < b < 2c, then for any θ₀ 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 θ − θ₀ with −(θ − θ₀), 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.

3.2.2. The Local Stability of Equilibrium States

In order to analyze the stability of equilibrium solutions we select n head-direction neural cells which have the preferred direction from 0 to 2π, and they are $\frac{1·2\pi }{n}$, $\frac{2·2\pi }{n}$, …, $\frac{n·2\pi }{n}$. Therefore the head-direction neural network which contains n units is as follows

We set solution u(θ, t) = u(θ) + ε(θ, t), where u(θ) is a equilibrium solution and ε(θ, t) is a small perturbation. Use Taylor's expansion to expand g(u(θ, t)) at u(θ), we have:

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

where

Since u(θ) is the equilibrium solution, we have $U=\frac{1}{n}W{G}_1$. Thus head-direction neural network (67) can be simplified as follows:

and the linear part is

where $J=-I+\frac{1}{n}WG$ is the Jacobian matrix and I is n × n identity matrix. The general solution to this linear ordinary differential equations is:

where λ₁, λ₂, ⋯, λ_n are eigenvalues of the Jacobian matrix, V₁, V₂, ⋯, V_n are corresponding eigenvectors, and c₁, c₂, ⋯, 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 $\frac{1}{n}W$. They are

where ${\rho }_j=exp(\frac{2\pi 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

So we change the form of eigenvalues as follows:

i.e.,

As n tends to infinite the eigenvalues become:

Because of the synaptic weight w(θ) = b cos θ + c cos 2θ, the eigenvalues of matrix $\frac{1}{n}W$ are:

Set V = (V₀, V₁, V₂, ⋯, V_n−1), obviously V is an normal orthogonal matrix and the characteristic equation becomes

where

We have the characteristic equation as follows

Now we can obtain all the eigenvalues from Equation (82) when the equilibrium solution is chosen, and then we can further determine the stability of each equilibrium solution.

Case 1: u(θ) = 0.

Since u(θ) = 0, we have

Thus

So the eigenvalues of characteristic Equation (82) are

Here g(u) is a Heaviside function and g′(u) is Dirac delta function. That means g′(0) = +∞. So if b < 0 and c < 0 all the eigenvalues are negative, then equilibrium solution u(θ) = 0 is stable. Otherwise, if b > 0 or c > 0 there is at least one eigenvalue larger than 0, then equilibrium solution u(θ) = 0 is unstable.

In order to discuss the stability of other equilibrium solutions, we set $(V^T(\frac{1}{n}W)V)(V^{T}GV)={(g_{i,j})}_{n\times n,}$ thus the characteristic Equation (82) is given by

i.e.,

where

As n approaches infinity, Equation(87) become

So we can further determine the stability of other equilibrium solutions.

Case 2: $u(\theta )=\frac{b}{\pi }\text{cos }\theta ,b>0$.

When the equilibrium is $u(\theta )=\frac{b}{\pi }\text{cos }\theta (b>0)$, we have:

Use the same computing method, we obtain

Therefore characteristic (Equation 82) is

Thus the eigenvalues are

Because of b > 0 when parameters satisfy c < 2b all eigenvalues are negative, thus equilibrium solution $u(\theta )=\frac{b}{\pi }\text{cos }\theta (b>0)$ is stable.

Case 3: $u(\theta )=\frac{c}{\pi }\text{cos }2\theta ,c>0$.

When the equilibrium is $u(\theta )=\frac{c}{\pi }\text{cos }2\theta (c>0)$, by computation we get

Therefore the characteristic Equation (82) is

and the eigenvalues are:

All the eigenvalues are negative, so $u(\theta )=\frac{c}{\pi }\text{cos }2\theta (c>0)$ is stable.

Case 4: $u(\theta )=\frac{b}{\pi }\sqrt{\frac{c+b}{2c}}\text{cos }\theta \pm \frac{\sqrt{c^2-b^2}}{2\pi }\text{cos }2\theta ,0<b<c\le 2b$.

When the equilibrium is $u(\theta )=\frac{b}{\pi }\sqrt{\frac{c+b}{2c}}\text{cos }\theta +\frac{c}{\pi }\frac{\sqrt{c^2-b^2}}{2c}\text{cos }2\theta (0<b<c\le 2b)$, we obtain

Therefore characteristic equation is

For convenient we set $k=\frac{c}{b}$, the Equation (97) becomes

Thus the eigenvalues are

Obviously when parameters satisfy 0 < b < c ≤ 2b, i.e., 1 < k ≤ 2, all eigenvalues are negative. So the equilibrium solution $u(\theta )=\frac{b}{\pi }\sqrt{\frac{c+b}{2c}}\text{cos }\theta +\frac{c}{\pi }\frac{\sqrt{c^2-b^2}}{2c}\text{cos }2\theta$ is stable. By the same analysis method we also get that solution $u(\theta )=\frac{b}{\pi }\sqrt{\frac{c+b}{2c}}\text{cos }\theta -\frac{c}{\pi }\frac{\sqrt{c^2-b^2}}{2c}\text{cos }2\theta$ 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:

Since g₂₃ = g₃₂, we have

Since ${(b{g}_{22}-c{g}_{33})}^2+4bc{g}_{23}^2>0$, we know all eigenvalues of the four remaining equilibrium solutions are real numbers and the greatest eigenvalue is

where,

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(\theta )=\frac{b}{\pi }\sqrt{\frac{2c-b}{2c}}\text{cos}(\theta -{\theta }_0)+\frac{b}{2\pi }\text{sin }2(\theta -{\theta }_0)$ (63) we set the parameters as a = 0, b = 3, c = 2, and θ₀ = 0 to obtain $u(\theta )=\frac{3}{2\pi }(1+2\text{sin }\theta )\text{cos }\theta$. 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_{\theta }=\frac{3}{\pi }\text{cos}(\theta -{\theta }_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(\theta )=\frac{3}{\pi }\text{cos}(\theta -\frac{6.5\pi }{42})$. We can see that the system starts from an unstable equilibrium state but converges gradually to this stable equilibrium.

3.2.3. 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(\theta )=\frac{b}{\pi }\text{cos}(\theta -{\theta }_0)$;

(iii) When c > 0, system (32) has symmetric and two-peaked equilibrium solutions $u(\theta )=\frac{c}{\pi }\text{cos }2(\theta -{\theta }_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

where

The distribution of all possible equilibrium solutions in system (32) is shown in Figure 9.

3.2.4. 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 head-direction 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(\theta )=\frac{b}{\pi }\text{cos }\theta$ and the two-peaked attractor state $u(\theta )=\frac{b}{\pi }\text{cos }2\theta$, with parameters a = 0, b = 3, and c = 2. Figure 10C shows the shifting mechanism for the stable special two peaked attractor state $u(\theta )=\frac{b}{\pi }\sqrt{\frac{c+b}{2c}}\text{cos }\theta -\frac{\sqrt{c^2-b^2}}{2\pi }\text{cos }2\theta$, with parameters a = 0, b = 1 and c = 1.5.

4. 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)=\frac{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

i.e.,

the constant solution u(θ) = 0 is stable. Furthermore $b=\frac{8}{k}$ and $c=\frac{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₀ = cos θ, the two-peaked state u₀ = 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 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 head-direction 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.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References