Skip to main content

ORIGINAL RESEARCH article

Front. Comput. Neurosci., 23 January 2020
This article is part of the Research Topic From Neuronal Network to Artificial Neural Network: Structure, Function and Intelligence View all 7 articles

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

  • 1School of International Economics, China Foreign Affairs University, Beijing, China
  • 2Department of Biomedical Engineering, Johns Hopkins University School of Medicine, Baltimore, MD, United States

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.

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

FIGURE 1
www.frontiersin.org

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

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

τu(θ,t)t=u(θ,t)+w(θ,t)g(u(θ,t))+I(θ,t),   θ[0,2π),    (1)

where the convolution is defined by

w(θ,t)g(u(θ,t))=12π02πw(θφ,t)g(u(φ,t))dφ.    (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

τu(θ,t)t=u(θ,t)+12π02πw(θφ,t)g(u(φ,t))dφ                +I(θ,t),   θ[0,2π).    (3)

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 etτ of both sides of this equation. It is easy to get the solution of Equation (3) as follows:

u(θ,t)=u(θ,0)etτ+0tesτ(02πw(θφ,s)g(u(φ,s))dφ)ds2πτetτ              +0tesτI(θ,s)dsτetτ.    (4)

Obviously, as t → +∞ the first term u(θ,0)/etτ tends to zero. If w(θ, t) and g(u) are bounded function, then the second term is bounded. This is because:

| 02πw(θφ,s)g(u(φ,s))dφ |         02π|w(θφ,s)g(u(φ,s))|dφ2π M,    (5)

where M is the maximum value of the integrand. At the same time we have:

| 0tesτ(02πw(θφ,s)g(u(φ,s))dφ)ds2πτetτ |           M0tesτdsτetτM,(t+).    (6)

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:

w(θ,t)=n=+wn(t)einθ,      I(θ,t)=n=+In(t)einθ.    (7)

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:

u(θ,t)=n=+un(t)einθ,      g(u(θ,t))=n=+gn(t)einθ.    (8)

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

[w(θ,t)g(u(θ,t))]n(t)=wn(t)·gn(t),n=0,±1,±2,.    (9)

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

τdun(t)dt=un(t)+wn(t)·gn(t)+In(t),n=0,±1,±2,.    (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.,

w(θ,t)=n=m+mwn(t)einθ,    (11)

then when |n| > m the Fourier series coefficients of solution u(θ, t) is:

τdun(t)dt=un(t),    (12)

i.e.,

un(t)=un(0)etτ.    (13)

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

u(θ,t)t=0.    (14)

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:

u(θ)=12π02πw(θφ,t)g(u(φ))dφ+I(θ,t).    (15)

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

u(θ)=12π02πw(θφ)g(u(φ))dφ+I.    (16)

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:

u(θθ0)=12π02πw(θθ0φ)g(u(φ))dφ+I.    (17)

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

u(θθ0)=12πθ02π+θ0w(θθ0(ϕθ0))g(u(ϕθ0))d(ϕθ0)+I.    (18)

That is,

u(θθ0)=12π02πw(θϕ)g(u(ϕθ0))dϕ+I.    (19)

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.

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:

w(θ)=a0+a1 cos θ+b1 sin θ+a2 cos 2θ             +b2 sin 2θ++an cos nθ+bn sin nθ,    (20)

and then an equilibrium solution satisfies

u(θ)=12π02π(a0+a1 cos(θφ)+a2 cos 2(θφ)            ++an cos n(θφ))g(u(φ))dφ            +12π02π(b1 sin(θφ)+b2 sin 2(θφ)            ++bn sin n(θφ))g(u(φ))dφ.    (21)

Since

12π( 02πan cos n(θφ)g(u(φ))dφ           +02πbn sin n(θφ)g(u(φ))dφ )=An cos nθ         +Bn sin nθ,    (22)

where

An=12π(an02πcos nφg(u(φ))dφbn02πsin nφg(u(φ))dφ),Bn=12π(an02πsin nφg(u(φ))dφ+bn02πcos nφg(u(φ))dφ),    (23)

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

u(θ)=A+A1 cos θ+B1 sin θ+A2 cos 2θ+B2 sin 2θ          ++An cos θ+Bn sin θ.    (24)

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:

E=02π0g(u(θ,t))(g1(V)I)dVdθ      14π02π02πw(θφ)g(u(φ,t))g(u(θ,t))dφdθ.    (25)

Obviously, the energy function is bounded and the time derivative of Equation (25) is:

dE/dt=02π(g1(g(u(θ,t)))I)dg(u(θ,t))dtdθ                14π02π02πw(θφ,t)( g(u(φ,t))dg(u(θ,t))dt                   +g(u(θ,t))dg(u(φ,t))dt )dφdθ.    (26)

We find if w(θ) is even, i.e., w(θ − φ) = w(φ − θ), then

02π02πw(θφ,t)g(u(φ,t))dg(u(θ,t))dtdφdθ            =02π02πw(θφ,t)g(u(θ,t))dg(u(φ,t))dtdφdθ.    (27)

Thus its time derivative becomes:

dE/dt=02π(u(θ,t)I)dg(u(θ,t))dtdθ               12π02π02πw(θφ,t)g(u(φ,t))dg(u(θ,t)dtdφdθ               =02π[u(θ,t)I12π02πw(θ               φ,t)g(u(φ,t))dφ]dg(u(θ,t))dtdθ               =τ02πdg1(g(u(θ,t)))dg(u(θ,t))·(dg(u(θ,t))dt)d2θ    (28)

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., dg-1(V)dV>0. Since the time parameter τ > 0, we obtain dEdt0, where dEdt=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:

u(θ)=a0+b1 cos θ+b2 sin θ+c1 cos 2θ+c2 sin 2θ,    (29)

where,

a0=a02πg(u(φ))dφ,b1=b02πcos φ g(u(φ))dφ,      b2=b02πsin φ g(u(φ))dφ,c1=c02πcos 2φ g(u(φ))dφ,    c2=c02πsin 2φ g(u(φ))dφ.    (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:

g(u)=11+ek(uu0),    (31)

where k > 0 is the gain and u0 is the threshold. For convenience in this paper we set the threshold u0 = 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:

τu(θ,t)t=u(θ,t)+12π02π( a+b cos(θφ)                       +c cos 2(θφ) ) H(u(φ,t))dφ.    (32)

The equilibrium solution equation becomes:

u(θ)=12π02π(a+b cos(θφ)+c cos 2(θφ))H(u(φ))dφ.    (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.

FIGURE 2
www.frontiersin.org

Figure 2. The synaptic weights for parameters a = 0, b = 3 and c = 2. (A) Diagram of weight profile w(θ) = a + b cos θ + c cos 2θ. (B) The synaptic weight matrix of the ring network.

3.2.1. All Possible Equilibrium States

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

u(θ)=a0+b1 cos θ+b2 sin θ+c1 cos 2θ+c2 sin 2θ.    (34)

Using basic trigonometric formula, we rewrite this form as:

u(θ)=A+B cos(θθ1)+C cos 2(θθ2),    (35)

where A = a0, B=b12+b22, C=c12+c22, θ1=arctanb2b1 and 2θ2=arctanc2c1. 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 phase-shift θ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.

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

u(θ)=H(A)2π02π(a+b cos(θφ)+c cos 2(θφ))dφ,

i.e.,

u(θ)=aH(A).    (36)

Since a = 0, the constant A=a02π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:

u(θ)=B cos(θθ1)+C cos 2(θθ2),    (37)

where B and C are unknown constants and θ1 and θ2 are unknown angles.

Case 2: u(θ) = B cos(θ − θ0), 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:

u(θ)=12π02π(b cos(θφ)+c cos 2(θφ))H(B cos φ)dφ

i.e.,

u(θ)=12ππ2π2(b cos(θφ)+c cos 2(θφ))dφ=bπ cos θ.    (38)

So if b > 0, then u(θ)=bπcos(θ-θ0) is the equilibrium solution for any θ0; 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(θ)=bπcos(θ-π) in Figure 3A. Figure 3B shows numerical simulation with initial state u0 is 0.1bπ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.

FIGURE 3
www.frontiersin.org

Figure 3. Equilibrium state for parameters a = 0, b = 3, and c = 2. (A) The equilibrium solution u(θ)=bπcos(θ-π). (B) The time evolution when initial state is u0 is 0.1bπcos(θ-π) with additional small perturbations, where the blue curve is the equilibrium solution u(θ)=bπcos(θ-π). Parameters: number of head-direction cells N = 500, τ = 1, time step Δt = 0.1, and maximum steps 500.

Case 3: u(θ) = C cos 2(θ − θ0), 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(θ − θ0)(C > 0). Plugging u(θ) = C cos 2θ into (Equation 33) we find the following equality:

u(θ)=12π02π(b cos(θφ)+c cos 2(θφ))H(C cos 2φ)dφ                 =12ππ4π4(b cos(θφ)+c cos 2(θφ))dφ                 +12π3π45π4(b cos(θφ)+c cos 2(θφ))dφ                 =cπ cos 2θ.    (39)

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

FIGURE 4
www.frontiersin.org

Figure 4. Another equilibrium state for the same parameters as Figure 3. (A) The equilibrium solution u(θ)=cπcos(θ-π2). (B) The time evolution when the initial state u0 is 0.1cπcos 2(θπ2)) with additional small perturbations. The blue curve is the equilibrium solution u(θ)=cπcos(θ-π2).

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

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

B cos(θθ1)+C cos 2(θθ1(θ2θ1))=B cos ϕ+C cos 2(ϕϕ0),    (40)

where ϕ = θ − θ1 and ϕ0 = θ2 − θ1. Due to the property of equilibrium solutions we just need to find the basic equilibrium solutions

u(ϕ)=B cos ϕ+C cos 2(ϕϕ0).    (41)

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.

FIGURE 5
www.frontiersin.org

Figure 5. Equilibrium states for parameters a = 0, b = 1, and c = 1.5, with other parameters the same as Figure 3. (A) A sketch of u(θ) = B cos ϕ + C cos 2(ϕ − ϕ0), BC > 0. (B) The profile of synaptic weight function w(θ) = a + b cos θ + c cos 2θ. (C) The equilibrium solution u1(θ)=bπc+b2ccos θ+ c2-b22πcos 2θ. (D) The time evolution for initial state u0 = 0.1u1(θ) with additional small perturbations. The blue lines indicates the equilibrium solution. (E) The equilibrium solution u2(θ)=bπc+b2ccos θ-c2-b22πcos 2θ. (F) The time evolution for initial state u0 = 0.1u2(θ) with additional small perturbations. The blue lines indicates the equilibrium solution.

As a equilibrium solution, we plug u(ϕ) = B cos ϕ + C cos 2(ϕ − ϕ0) into Equation (33) to obtain

u(ϕ)=B cos ϕ+C cos  2(ϕϕ0)=12π02π(b cos(ϕφ)                 +c cos  2(ϕφ))H(B cos ϕ+C cos  2(ϕϕ0))dφ                 =12παα(b cos(ϕφ)+c cos  2(ϕφ))dφ,    (42)

i.e.,

u(ϕ)=bπsin α-α2cos (ϕ-α+α2)          +c2πsin (α-α) cos 2(ϕ-α+α2).    (43)

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:

B cosϕ±C cos 2ϕ=bπsin α cos ϕ+c2πsin 2α cos 2ϕ.    (44)

Comparing the corresponding coefficients we have B=bπsinα and ±C=c2πsin2α. Since α is a root of equation u(ϕ) = 0, next we should solve equation u(ϕ) = 0. We find sin α = 0, cos α = 0 and cos 2α=-bc. Because B > 0 and C > 0, we keep cos 2α=-bc and reject the others. If cos 2α=-bc is a valid root, we have |bc|<1, i.e., |b| < |c|. So we obtain sinα=c+b2c and cos α=±c-b2c. At the same time we have B=bπc+b2c and C=cπc2-b24c2. For the equations to hold there is only one positive domain in period 2π, we have BC > 0, i.e., bπc+b2ccπc2-b24c2>0, and the parameter sb and c must satisfy 2bc > 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:

u(θ)=bπc+b2c cos (θ-θ0)±c2-b22πcos 2(θ-θ0).    (45)

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 u0 = 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

u(ϕ)=B cos ϕ+C cos 2(ϕϕ0)=12π02π(b cos(ϕφ)              +c cos 2(ϕφ))H(B cos ϕ+C cos 2(ϕϕ0))dφ              =12παα(b cos(ϕφ)+c cos 2(ϕφ))dφ              +12πββ(b cos(ϕφ)+c cos 2(ϕφ))dφ              =b2π( cos α cos α+ cos β cos β) sin ϕ              b2π( sin α sin α+ sin β sin β) cos ϕ              +c4π( cos 2α cos 2α+ cos 2β cos 2β) sin 2ϕ              c4π( sin 2α sin 2α+ sin 2β sin 2β) cos 2ϕ.    (46)
FIGURE 6
www.frontiersin.org

Figure 6. (A) A sketch of u(θ) = B cos ϕ + C cos 2(ϕ − ϕ0), 0 < B < C. (B) Four solid lines represent the four roots of Equation (49), i.e., cos α, cos α′, cos β and cos β′, while two dash lines represent all possibilities of cos − cos α′ + cos β − cos β′. Parameter k = 0.8.

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.    (47)

For convenience let x = cos ϕ and k=BC we square both side of the Equation (47), then the Equation (47) becomes

(kx+(2x2-1) cos 2ϕ0)2=4(1-x2)x2(1- cos 2ϕ02),    (48)

i.e.,

4x4+4k cos 2ϕ0x3+(k2-4)x2-2k cos 2ϕ0x+cos 22ϕ0=0.    (49)

According to the four roots of Equation (49), we can find that if and only if 2ϕ0 is equals to 0,π2,π and 3π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ϕ0 equals 0, π2, π and 3π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ϕ0-axis are immovability. And the change of picture is successive. That means if and only if 2ϕ0 is equals to 0,π2,π and 3π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)0 = 0.

If 2ϕ0 = 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

u(ϕ)=B cos ϕ+C cos 2ϕ=bπ(sin α- sin β) cos ϕ              +c2π(sin 2α- sin 2β) cos 2ϕ.    (50)

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

{ bπ(sin α sin β) cos α+c2π(sin 2α sin 2β) cos 2α=0,bπ(sin α sin β) cos β+c2π(sin 2α sin 2β) cos 2β=0.    (51)

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

bπ(sin α- sin β)(cos α- cos β)            +c2π(sin 2α- sin 2β)(cos 2β- cos 2β)=0.    (52)

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

b(cos α- cos β)-c(2 sin α cos α-2 sin β cos β)(sin β+ sin α)=0.    (53)

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

2C cos 2ϕ+B cos ϕ-C=0.    (54)

So we have cos α·cos β=-12. Plug cos β=-12cos α into (Equation 53) and compute cos α. In order to keep meaningful during computational process we set 0 < b < 2c and bc12, then we have cos α=s-s2-12, where s=2-bc+2bc2. Therefore we have the following conclusion:

If 0 < b < 2c and bc12, then for any θ0 the system has equilibrium solutions which are

u(θ)=bπ(1m24m212m)cos(θθ0)+cπ(m1m2             +4m214m2)cos 2(θθ0),    (55)

where

m= cos α=2-bc+2bc-b2c2-2bc-2bc2bc+42bc2.    (56)

(2)0 = π.

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

If 0 < b < 2c and bc12, then for any θ0 the system has equilibrium solutions which are

u(θ)=bπ(4m2-12m-1-m2) cos (θ-θ0)    -cπ(4m2-14m2+m1-m2) cos 2(θ-θ0),    (57)

where

m= cos α=2-bc+2bc+b2c2-2bc-2bc2bc+42bc2.    (58)

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.

FIGURE 7
www.frontiersin.org

Figure 7. Equilibrium states for parameters a = 0, b = 3 and c = 2, which satisfy 0 < b < 2c and bc12. The network contains N = 1000 head-direction cells, τ = 1 and time step Δt = 0.1. (A) The equilibrium solution (55). (B) The time evolution of the equilibrium solution (55), where blue line represents the equilibrium solution (55). (C) The equilibrium solution (57). (D) The time evolution of equilibrium solution (57), where blue line represents the equilibrium solution (57).

(3) 2ϕ0=π2.

If 2ϕ0=π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 −α, π2, π + α and 3π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:

u(θ)=B cos θ+C sin 2θ=bπsin α cos θ+cπcos2α sin 2θ.    (59)

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

u(-α)=bπsin α cos α-cπcos2 α sin 2α=0    (60)

i.e.,

b sin α cos α-2c cos3α sin α=0.    (61)

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

u(θ)=bπ2c-b2c cos θ+b2πsin 2θ.    (62)

In a words, if 0 < b < 2c, then for any θ0 the system has equilibrium solutions which are:

u(θ)=bπ2c-b2ccos(θ-θ0)+b2πsin 2(θ-θ0).    (63)

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

FIGURE 8
www.frontiersin.org

Figure 8. Equilibrium states for parameters a = 0, b = 3, and c = 2, which satisfy 0 < b < 2c and bc12. The network contains N = 1000 head-direction cells, τ = 1 and time step Δt = 0.1. (A) The time evolution of equilibrium solution (63) conforms that the state is stationary. (B) The time evolution of equilibrium solution (63) with disturbances in the initial state indicates that the equilibrium is unstable. (C) Comparison chart, where the blue line represents the equilibrium solution (63), the red line represents the final state starting from (63) after 100, 000 time steps, which matches perfectly with the one-peaked stable equilibrium solution u(θ)=bπcos(θ-6.5π42) as indicated by the green stars.

(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

u(θ)=bπ2c-b2ccos(θ-θ0)-b2πsin 2(θ-θ0).    (64)

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.

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 1·2πn, 2·2πn, …, n·2πn. Therefore the head-direction neural network which contains n units is as follows

τt(u(θ1,t)u(θ2,t)u(θn,t))=-(u(θ1,t)u(θ2,t)u(θn,t))+1n(w(θ1-θ1)w(θ1-θ2)w(θ1-θn)w(θ2-θ1)w(θ2-θ2)w(θ2-θn)w(θn-θ1)w(θn-θ2)w(θn-θn))(g(u(θ1,t))g(u(θ2,t))g(u(θn,t))).    (65)

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:

g(u(θ,t))=g(u(θ)+ε(θ,t))=g(u(θ))+g(u(θ))ε(θ,t))+O(ε2).    (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

τEt=-(U+E)+1nW(G1+GE)+O(E2)    (67)

where

E=(ε(θ1,t)ε(θ2,t) ε(θn,t)),W=(w(θ1θ1)w(θ1θ2)w(θ1θn)w(θ2θ1)w(θ2θ2)w(θ2θn)w(θnθ1)w(θnθ2)w(θnθn)),U=(u(θ1)u(θ2)u(θn)),G1=(g(u(θ1))g(u(θ2))g(u(θn))),G=(g(u(θ1)) 000g(u(θ2))000g(u(θn))),    (68)

Since u(θ) is the equilibrium solution, we have U=1nWG1. Thus head-direction neural network (67) can be simplified as follows:

τEt=-E+1nWGE+O(E2)=(-I+1nWG)·E+O(E2),    (69)

and the linear part is

τEt=JE,    (70)

where J=-I+1nWG is the Jacobian matrix and I is n × n identity matrix. The general solution to this linear ordinary differential equations is:

E=c1eλ1tV1+c2eλ2tV2++cneλntVn,    (71)

where λ1, λ2, ⋯, λn are eigenvalues of the Jacobian matrix, V1, V2, ⋯, Vn are corresponding eigenvectors, and c1, c2, ⋯, cn 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 1nW. They are

Vj=1n(ρj0,ρj1,,ρjn-1)T,j=0,1,,n-1,    (72)

where ρj=exp(2πinj) (j = 0, 1, ⋯, n−1) are the n-th roots of unity and i is the imaginary unit. The corresponding eigenvalues6 are given by

λj=1n(w(02πn)ρj0+w(12πn)ρj1++w((n-1)2πn)ρjn-1)T,               j=0,1,,n-1.    (73)

Since

ρjn-k=exp(j(n-k)2πni)=exp(j-k2πni)=exp(jk2πni)¯              =ρjk¯,j=0,1,,n-1,    (74)

thus

ρjn-k+ρjk=ρjk¯+ρjk=2 cos (2kπnj),j=0,1,,n-1.    (75)

So we change the form of eigenvalues as follows:

λj=1n[w(02πn)cos(02πnj)+w(12πn)cos(12kπnj)            ++w((n1)2πn)cos((n1)2πnj)].    (76)

i.e.,

λj=1nk=0n-1w(k2πn) cos (k2πnj),j=0,1,,n-1.    (77)

As n tends to infinite the eigenvalues become:

λj=limn1nk=0n-1w(k2πn)cos(k2πnj)=12π02πw(φ)cos(jφ)dφ,             j=0,1,,n-1.    (78)

Because of the synaptic weight w(θ) = b cos θ + c cos 2θ, the eigenvalues of matrix 1nW are:

λ0=0,λ1=b2,λ2=c2,λ3=0,,λn-1=0.    (79)

Set V = (V0, V1, V2, ⋯, Vn−1), obviously V is an normal orthogonal matrix and the characteristic equation becomes

f(λ)=|VT(λ+1)IV-VT(1nWG)V|            =|(λ+1)I-(VT(1nW)V)(VTGV)|,    (80)

where

VT1nWV=(000000b200000c2000000000000)=K.    (81)

We have the characteristic equation as follows

f(λ)=|(λ+1)I-K(VTGV)|.    (82)

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

VTGV=VT(g(0)000g(0)000g(0))V=g(0)VTIV=g(0)I.    (83)

Thus

VT(1nW)V·VTGV=g(0)·I·(000000b200000c2000000000000)13g    =(000000g(0)b200000g(0)c2000000000000)    (84)

So the eigenvalues of characteristic Equation (82) are

λ1=-1,λ2=-1+g(0)b2,λ3=-1+g(0)c2,λ4==λn=-1.    (85)

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 (VT(1nW)V)(VTGV)=(gi,j)n×n, thus the characteristic Equation (82) is given by

f(λ)=(λ+1)n-2·|λ+1-b2g22λ+1-b2g23λ+1-c2g32λ+1-c2g33|

i.e.,

f(λ)=(λ+1)n2( (λ+1)2(b2g22+c2g33)(λ+1)               +bc4(g22g33g23g32) ),    (86)

where

g22=1n(g(u(θ1))(ρ10)2+g(u(θ2))(ρ11)2 ++g(u(θn))(ρ1n1)2),g33=1n(g(u(θ1))(ρ20)2+g(u(θ2))(ρ21)2 ++g(u(θn))(ρ2n1)2),g32=g32=1n(g(u(θ1))ρ10ρ20+g(u(θ2))ρ11ρ21 ++g(u(θn))ρ1n1ρ2n1).    (87)

As n approaches infinity, Equation(87) become

g22=12π02πg(u(θ))e2θidθ,g33=12π02πg(u(θ))e4θidθ,g23=g32=12π02πg(u(θ))e3θidθ.    (88)

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

Case 2: u(θ)=bπcos θ,b>0.

When the equilibrium is u(θ)=bπcos θ(b>0), we have:

g22=12π(0πg(u(θ))e2iθdθ+π2πg(u(θ))e2iθdθ)         =12π(bπbπg(u)e2iarccos(πbu)πb1(πbu)2du          + bπbπg(u)e2i(2πarccos(πbu))πb1(πbu)2du)=1b.    (89)

Use the same computing method, we obtain

g33=1b,     g32=g32=0.    (90)

Therefore characteristic (Equation 82) is

f(λ)=(λ+1)n2((λ+1)2+(12c2b)(λ+1)c4b).    (91)

Thus the eigenvalues are

λ1==λn2=1,λn1=32,λn=1+c2b.    (92)

Because of b > 0 when parameters satisfy c < 2b all eigenvalues are negative, thus equilibrium solution u(θ)=bπcos θ(b>0) is stable.

Case 3: u(θ)=cπcos 2θ,c>0.

When the equilibrium is u(θ)=cπcos 2θ(c>0), by computation we get

g22=0,    g33=1c,    g32=g32=0.    (93)

Therefore the characteristic Equation (82) is

f(λ)=(λ+1)n1((λ+1)+12),    (94)

and the eigenvalues are:

λ1==λn1=1,λn=32.    (95)

All the eigenvalues are negative, so u(θ)=cπcos 2θ(c>0) is stable.

Case 4: u(θ)=bπc+b2ccos θ±c2-b22πcos 2θ,0<b<c2b.

When the equilibrium is u(θ)=bπc+b2ccos θ+cπc2-b22ccos 2θ(0<b<c2b), we obtain

g22=2b2c2+bcb2 g33=4b22c2(2c2+bcb2)c,                            g32=g32=2(2b+c)cb2c2c2+bcb2.    (96)

Therefore characteristic equation is

f(λ)=(λ+1)n2( (λ+1)2+c2b22c2+bcb2(λ+1)               b2c2+bc32(2c2+bcb2)2 ).    (97)

For convenient we set k=cb, the Equation (97) becomes

f(λ)=(λ+1)n2( (λ+1)2+k212k2+k1(λ+1)               k2+k32(2k2+k1)2 ).    (98)

Thus the eigenvalues are

λ1==λn2=1,λn1,n=1+1+2k+k22k3±14k+4k2+2k37k4+4k5+4k62(13k+4k3).    (99)

Obviously when parameters satisfy 0 < b < c ≤ 2b, i.e., 1 < k ≤ 2, all eigenvalues are negative. So the equilibrium solution u(θ)=bπc+b2ccos θ+cπc2-b22ccos 2θ is stable. By the same analysis method we also get that solution u(θ)=bπc+b2ccos θ-cπc2-b22ccos 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

f1(λ)=(λ+1)2(b2g22+c2g33)(λ+1)+bc4(g22g33g23g32),

thus all the eigenvalues are:

λ1==λn2=1, λn1,n          =1+bg22+cg33±(bg22cg33)2+4bc(g22g33g23g32)4,

Since g23 = g32, we have

λn1,n=1+bg22+cg33±(bg22cg33)2+4bcg2324.

Since (bg22-cg33)2+4bcg232>0, we know all eigenvalues of the four remaining equilibrium solutions are real numbers and the greatest eigenvalue is

λmax=λn=1+bg22+cg33+(bg22cg33)2+4bcg2324,    (100)

where,

g22=12π02πg(u(θ))e2θidθ,     g33=12π02πg(u(θ))e4θidθ,g23=12π02πg(u(θ))e3θidθ.

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-b2ccos(θ-θ0)+b2πsin 2(θ-θ0) (63) we set the parameters as a = 0, b = 3, c = 2, and θ0 = 0 to obtain u(θ)=32π(1+2sin θ)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.

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(θ)=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

u(θ)=bπc+b2ccos(θθ0)±c2b22πcos 2(θθ0);    

(v) When 0 < b < 2c, the system (32) has four kinds of special equilibrium solutions:

u(θ)=bπ2c-b2c cos (θ-θ0)±b2πsin 2(θ-θ0)

and

u(θ)=±bπ(1m24m212m)cos(θθ0)±cπ(m1m2+4m214m2)cos 2(θθ0),    

where

m=2-bc+2bcb2c2-2bc-2bc2bc+42bc2.

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

FIGURE 9
www.frontiersin.org

Figure 9. Phase diagram for the distribution of all possible equilibrium solutions in head-direction ring network.

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:

w1(θ)=w(θ)+αdw(θ)dθ,    (101)

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(θ)=bπc+b2ccos θ-c2-b22πcos 2θ, with parameters a = 0, b = 1 and c = 1.5.

FIGURE 10
www.frontiersin.org

Figure 10. (A) Shifting of equilibrium solution u(θ)=bπcos(θ), with parameters a = 0, b = 3, c = 2; (B) Shifting of equilibrium solution u(θ)=cπcos(2θ), with parameters a = 0, b = 3, c = 2; (C) Shifting of equilibrium solution u(θ)=bπc+b2ccos θ-c2-b22πcos 2θ, with parameters a = 0, b = 1, 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)=11+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)b2<1 and g(0)c2<1,    (102)

i.e.,

b<8k and c<8k,    (103)

the constant solution u(θ) = 0 is stable. Furthermore b=8k and c=8k 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 u0 = cos θ, the two-peaked state u0 = 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.

FIGURE 11
www.frontiersin.org

Figure 11. Equilibrium states in the ring network with sigmoidal gain function g(u)=11+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
www.frontiersin.org

Figure 12. The phase diagram for the distribution of equilibrium states in the ring network with the sigmoidal gain function g(u)=11+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, u0 = cosθ and u0 = 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).

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

Amari, S. I. (1977). Dynamics of pattern formation in lateral-inhibition type neural fields. Biol. Cybernet. 27, 77–87.

PubMed Abstract | Google Scholar

Arleo, A., and Gerstner, W. (2001). Spatial orientation in navigating agents: modeling head-direction cells. Neurocomputing 38, 1059–1065. doi: 10.1016/S0925-2312(01)00572-0

CrossRef Full Text | Google Scholar

Ben-Yishai, R., Bar-Or, R. L., and Sompolinsky, H. (1995). Theory of orientation tuning in visual cortex. Proc. Natl. Acad. Sci. U.S.A. 92, 3844–3848.

PubMed Abstract | Google Scholar

Bressloff, P. C. (2012). Spatiotemporal dynamics of continuum neural fields. J. Phys. A Math. Theor. 45:033001. doi: 10.1088/1751-8113/45/3/033001

CrossRef Full Text | Google Scholar

Cohen, M. A., and Grossberg, S. (1983). Absolute stability of global pattern formation and parallel memory storage by competitive neural networks. IEEE Trans. Syst. Man Cybernet. 5, 815–826.

Google Scholar

Finkelstein, A., Derdikman, D., Rubin, A., Foerster, J. N., Las, L., and Ulanovsky, N. (2015). Three-dimensional head-direction coding in the bat brain. Nature 517, 159–164. doi: 10.1038/nature14031

PubMed Abstract | CrossRef Full Text | Google Scholar

Goodridge, J. P., and Touretzky, D. S. (2000). Modeling attractor deformation in the rodent head-direction system. J Neurophysiol. 83, 3402–3410. doi: 10.1152/jn.2000.83.6.3402

PubMed Abstract | CrossRef Full Text | Google Scholar

Gray, R. M. (2005). Toeplitz and circulant matrices: a review. Commun. Inform. Theory 2, 155–239. doi: 10.1561/0100000006

CrossRef Full Text | Google Scholar

Hahnloser, R. H., Sarpeshkar, R., Mahowald, M. A., Douglas, R. J., and Seung, H. S. (2000). Digital selection and analogue amplification coexist in a cortex-inspired silicon circuit. Nature 405:947. doi: 10.1038/35016072

PubMed Abstract | CrossRef Full Text | Google Scholar

Hopfield, J. J. (1984). Neurons with graded response have collective computational properties like those of two-state neurons. Proc. Natl. Acad. Sci. U.S.A 10, 3088–3092.

Google Scholar

Jacob, P.-Y., Casali, G., Spieser, L., Page, H., Overington, D., and Jeffery, K. (2017). An independent, landmark-dominated head-direction signal in dysgranular retrosplenial cortex. Nat. Neurosci. 20, 173–175. doi: 10.1038/nn.4465

PubMed Abstract | CrossRef Full Text | Google Scholar

Kakaria, K. S., and de Bivort, B. L. (2017). Ring attractor dynamics emerge from a spiking model of the entire protocerebral bridge. Front. Behav. Neurosci. 11:8. doi: 10.3389/fnbeh.2017.00008

PubMed Abstract | CrossRef Full Text | Google Scholar

Knierim, J. J., and Zhang, K. (2012). Attractor dynamics of spatially correlated neural activity in the limbic system. Ann. Rev. Neurosci. 35, 267–285. doi: 10.1146/annurev-neuro-062111-150351

PubMed Abstract | CrossRef Full Text | Google Scholar

Page, H. J., and Jeffery, K. J. (2018). Landmark-based updating of the head direction system by retrosplenial cortex: a computational model. Front. Cell. Neurosci. 12:191. doi: 10.3389/fncel.2018.00191

PubMed Abstract | CrossRef Full Text | Google Scholar

Pouget, A., Zhang, K., Deneve, S., and Latham, P. E. (1998). Statistically efficient estimation using population coding. Neural Comput. 10, 373–401.

PubMed Abstract | Google Scholar

Redish, A. D., Elga, A. N., and Touretzky, D. S. (1996). A coupled attractor model of the rodent head direction system. Network Comput. Neural Syst. 7, 671–685.

Google Scholar

Robertson, R. G., Rolls, E. T., Georges-François, P., and Panzeri, S. (1999). Head direction cells in the primate pre-subiculum. Hippocampus 9, 206–219.

PubMed Abstract | Google Scholar

Sharp, P. E., Blair, H. T., and Cho, J. (2001). The anatomical and computational basis of the rat head-direction cell signal. Trends Neurosci. 24, 289–294. doi: 10.1016/S0166-2236(00)01797-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Skaggs, W. E., Knierim, J. J., Kudrimoti, H. S., and McNaughton, B. L. (1995). A model of the neural basis of the rats sense of direction. Adv. Neural Inf. Process. Syst. 7, 173–180.

Google Scholar

Song, P., and Wang, X. J. (2005). Angular path integration by moving “hill of activity”: a spiking neuron model without recurrent excitation of the head-direction system. J. Neurosci. 25, 1002–1014. doi: 10.1523/JNEUROSCI.4172-04.2005

CrossRef Full Text | Google Scholar

Stringer, S. M., Trappenberg, T. P., Rolls, E. T., and Araujo, I. (2002). Self-organizing continuous attractor networks and path integration: one-dimensional models of head direction cells. Network Comput. Neural Syst. 13, 217–242. doi: 10.1080/net.13.2.217.242

PubMed Abstract | CrossRef Full Text | Google Scholar

Taube, J. S. (2007). The head direction signal: origins and sensory-motor integration. Annu. Rev. Neurosci. 30, 181–207. doi: 10.1146/annurev.neuro.29.051605.112854

PubMed Abstract | CrossRef Full Text | Google Scholar

Taube, J. S., Muller, R. U., and Ranck, J. B. (1990). Head-direction cells recorded from the postsubiculum in freely moving rats. I. Description and quantitative analysis. J. Neurosci. 10, 420–435.

PubMed Abstract | Google Scholar

Tee, G. J. (2007). Eigenvectors of block circulant and alternating circulant matrices. N. Zeal. J. Math. 36, 195–211.

Google Scholar

Wilson, H. R., and Cowan, J. D. (1972). Excitatory and inhibitory interactions in localized populations of model neurons. Biophys. J. 12, 1–24.

PubMed Abstract | Google Scholar

Xie, X., Hahnloser, R. H., and Seung, H. S. (2002). Double-ring network model of the head-direction system. Phys. Rev. E 66:041902. doi: 10.1103/PhysRevE.66.041902

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, K. (1996). Representation of spatial orientation by the intrinsic dynamics of the head-direction cell ensemble: a theory. J. Neurosci. 16, 2112–2126.

PubMed Abstract | Google Scholar

Zhang, W. H., Wang, H., Chen, A., Gu, Y., Lee, T. S., Wong, K. M., and Wu, S. (2019). Complementary congruent and opposite neurons achieve concurrent multisensory integration and segregation. eLife 8:e43753. doi: 10.7554/eLife.43753

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: continuous attractor network, ring network, neural field, dynamical system, Fourier analysis, stability

Citation: Wang C and Zhang K (2020) Equilibrium States and Their Stability in the Head-Direction Ring Network. Front. Comput. Neurosci. 13:96. doi: 10.3389/fncom.2019.00096

Received: 03 September 2019; Accepted: 23 December 2019;
Published: 23 January 2020.

Edited by:

Yilei Zhang, Nanyang Technological University, Singapore

Reviewed by:

Yuanyuan Mi, Chongqing University, China
Guangtao Zhai, Shanghai Jiao Tong University, China

Copyright © 2020 Wang and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Caixia Wang, wangcaixia@cfau.edu.cn

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.