Amélie Aussel

Laure Buhry

Radu Ranta

Amelie Aussel

Stability conditions of Hopfield ring networks with discontinuous piecewise-affine activation functions*

published or not. The documents may come

Stability conditions of Hopfield ring networks with discontinuous piecewise-affine activation functions

I. INTRODUCTION

In the hippocampus and the entorhinal cortex, some neurons, such as place cells, grid cells and head-direction cells, are capable of firing at different rates depending on an animal's position ( [START_REF] O'keefe | The Hippocampus as a Cognitive Map[END_REF]). When organized in a network, these cells function as a "cognitive map" of one's environment, providing the basis of spatial orientation skills, but are also involved in short-term memory, as they can encode spatial information and retain it for some time even when all external stimuli are removed. Different computational models have been developed to represent these structures, for examples see [START_REF] Mcnaughton | Path integration and the neural basis of the 'cognitive map[END_REF], mainly based on Hopfield-like neural networks ( [START_REF] Hopfield | Neural networks and physical systems with emergent collective computational abilities[END_REF]). In Hopfield neural networks, each neuron's activity s k is computed from the weighted sum of the activity of all the N neurons of the network, through an activation function φ. As the synaptic decay of biological neurons is known to have an exponential shape, we have the following differential equation:

ṡk = - 1 τ s k + φ( N -1 i=0 w ki s i + b k ) (1)
or, in a matrix format:

ṡ = - 1 τ s + φ(Ws + b) (2)
where W is called the (recurrent) synaptic connectivity matrix, τ is a synaptic time constant and b is the external *This work was not supported by any organization 1 Amelie Aussel is with the Centre de Recherche en Automatique de Nancy (CRAN-CNRS UMR 7039 University of Lorraine, FRANCE), and the Laboratoire Lorrain de Recherche en Informatique et ses applications (LORIA UMR 7503, University of Lorraine-INRIA-CNRS, F-54506).

amelie.aussel@loria.fr INRIA-CNRS, F-54506). laure.buhry@loria.fr 3 Radu Ranta is with the Centre de Recherche en Automatique de Nancy (CRAN-CNRS UMR 7039 University of Lorraine, FRANCE).

radu.ranta@univ-lorraine.fr input provided to the system (see [START_REF] Pineda | Dynamics and Architecture for Neural Computation[END_REF]). The term Ws + b can be interpreted as the total synaptic input received by the neurons, and the term φ(Ws + b) as their firing rates (the number of action potentials they emit in one second).

In the first studies of the latter model, the synaptic connectivity matrix W was symmetric, and the activation function φ was C 1 and strictly increasing (often a sigmoid function), which ensured the convergence of the system to a stable point for any stimulation b (as shown in [START_REF] Hopfield | Neurons with graded response have collective computational properties like those of two-states neurons[END_REF]). However, when some of these constraints were lifted in order to account for more realistic networks, the stability analysis of the model became a much more complicated task (for a review, see [START_REF] Zhang | A Comprehensive Review of Stability Analysis in Continuous-Time Recurrent Neural Networks[END_REF]). In particular, a few authors have studied Hopfield-like neural networks with non-C 1 activation functions φ, in order to establish the conditions for a unique stable point to exist (for example see [START_REF] Forti | Global Convergence of Neural Networks with Discontinuous Neuron Activations[END_REF]).

In [START_REF] Richard | Permitted and Forbidden Sets in Symmetric Threshold-Linear Networks[END_REF], a necessary and sufficient condition for the stability of a system with a non-C 1 activation functions φ was presented. More precisely, with φ(x) = max(x, 0) it was shown that the system converges to an equilibrium point for all constant b and all initial conditions, if and only if the matrix

( 1 τ I -W) is copositive; that is, x T ( 1 τ I -W)x > 0 for all nonnegative x, except x = 0.
However, it is not clear how to apply this reasoning for a more complex activation function φ. Also, the copositivity of a matrix is not something easy to ensure when tuning the connectivity matrix W. In addition, it provides no information about the form of the solution at the stable points.

Our point is that some other stability conditions can be obtained, which may be weaker than those developed in [START_REF] Richard | Permitted and Forbidden Sets in Symmetric Threshold-Linear Networks[END_REF], but are also easier to verify when choosing the synaptic connectivity matrix W. These conditions take advantage of the periodic structure that can be imposed on W, when considering place cells or head-direction cells.

Indeed, because networks made of these cells have to be able to encode a continuum of positions, the associated model should also present a continuum of stable points, and such a continuous attractor dynamics is often achieved by choosing a synaptic connectivity matrix with some kind of periodicity (see [START_REF] Mcnaughton | Path integration and the neural basis of the 'cognitive map[END_REF]). This structure is coherent with experimental results, as it has been observed that hippocampal neurons are highly interconnected and tend to form associative networks (see [START_REF] Le Duigou | Recurrent synapses and circuits in the CA3 region of the hippocampus: an associative network[END_REF]).

In Section II of the present article, we will detail the model we used, a ring-shaped, Hopfield-like neural network with a piecewise affine activation function φ, then we will present the results we obtained regarding its stability of the form of its stable states in Section III. Finally, in Section IV we will show how these results can apply on a simulation of a network with a gaussian synaptic connectivity matrix.

II. DEFINITION OF THE RING NETWORK MODEL

A. Network dynamics

We studied the stability of a network consisting of N identical neurons, which activity s followed the dynamics (2) with an activation function φ defined by:

φ(x) = αx + β if x ≥ 0 0 if x < 0 (3)
where α and β are positive real constants.

The stimulation b is chosen strictly positive and uniform on all neurons (∀k ∈ 0

; N -1 , b k = b).
This definition of φ ensures that neuronal firing rates are always positive, and it also represents the fact that neurons can only fire if the total input they receive is greater than a certain threshold. With α > 0 and β = 0, the function φ represents the behaviour of type-I neurons (capable of firing at frequencies close to zero when excited with stimulus just above their threshold), and with α > 0 and β > 0, it represents the behaviour of type-II neurons (that have a minimal firing rate).

It can be shown from equations (2) and (3) that, whatever the initial conditions, the network activity s will always become and remain non-negative after a certain amount of time, and therefore, for simplicity reason, we will only consider non-negative s.

With our definition of the activation function φ (equation (3)), each neuron can follow two distinct linear differential equations, depending on the input it receives:

• On the one hand, if a neuron k receives an input that is greater than or equal to zero, it will be called "excited" and follow the equation

: ṡk = -1 τ s k +α N -1 i=0 w ki s i + αb + β
• On the other hand, if it receives an input that is less than zero, it will be called "inhibited" and follow the equation: ṡk = -1 τ s k Therefore, at any instant t, our whole network can be described by one of 2 N different linear systems depending on which neurons are excited or inhibited. We can then consider our model as a switched system in the form:

ṡ(t) = - 1 τ I N + αP Σ(t) W s(t) + P Σ(t) (αb + β)1 (4)
where I N is the identity matrix of size N-by-N, 1 is the Nby-1 all-ones vector and P Σ(t) is a diagonal N-by-N matrix containing ones at indexes corresponding to the excited neurons at the instant t and zeros elsewhere.

B. The ring structure condition

In order to model structures such as head-direction cells, we considered ring-shaped, rotation invariant networks, which imposed some constraints on the synaptic connectivity matrix W.

Each neuron was assigned a position so that they were all uniformly distributed on a circle: for k in 0; N -1 the angular position of neuron k was θ k = -π + 2kπ N . Then, because all neurons were identical (same time constant τ , same activation function φ), the rotation invariance was obtained by connecting all neurons to their neighbours identically. In other words, the synaptic interaction between any two neurons does not depend on their position on the ring, but only on the angular distance between them (w ij = f (θ i -θ j )). This imposes that W is Toeplitz. So that there was no discontinuity for θ i -θ j = kπ, W was even chosen circulant (which means that any line of W can be obtained as a right-shift a the previous one).

We also chose to keep W symmetric, which means that the function f is even, and that the synaptic interaction associated with an angular distance is the same whether we go clockwise or anti-clockwise. Our synaptic connectivity matrix is therefore of the form:

W =          w 0 w 1 w 2 • • • w 2 w 1 w 1 w 0 w 1 • • • w 3 w 2 w 2 w 1 w 0 • • • w 4 w 3 . . . . . . w 2 w 3 w 4 • • • w 0 w 1 w 1 w 2 w 3 • • • w 1 w 0          (5)
The elements of W are all real values, but they can be either positive or negative. A positive coefficient represents an excitatory synapse and a negative one represents an inhibitory synapse.

C. Connectivity matrix and submatrices eigenvalues

In order to analyze the stability of the system in the next section, we had to determine the eigenvalues of the matrix A Σ = (-1 τ I N + αP Σ W) for all possible matrices P Σ , as defined in II.

Let us consider a matrix P Σ with k zeros along its diagonal (so that it corresponds to a system with k inhibited neurons), and let us consider

A Σ = (-1 τ I N + αP Σ W) and A 0 = (-1 τ I N + αW).
The matrix P Σ W is obtained by replacing by zeroes all the k lines of W that corresponds to inhibited neurons. Therefore, as it is also explained in [START_REF] Richard | Permitted and Forbidden Sets in Symmetric Threshold-Linear Networks[END_REF], by developing the determinant of A Σ along these lines, we find that A Σ has k eigenvalues equal to -1 τ , and the other are the eigenvalues of the matrix ÃΣ , obtained by removing the rows and columns of A 0 corresponding to inhibited neurons.

Because W is a real symmetric matrix, it is noteworthy that ÃΣ is also a real symmetric matrix and therefore all the eigenvalues of the system are real. Moreover, we can infer from [START_REF] Horn | Matrix analysis[END_REF] (Cauchy interlacing theorem) that the eigenvalues of ÃΣ all lie between the minimum and the maximum of the eigenvalues of A 0 .

As W is also circulant, its eigenvalues can be calculated as the discrete Fourier transform of its first row's coefficients: for all m in 0; N -1 , the corresponding eigenvalue λ m of W is equal to:

λ m = N -1 n=0 w n e 2mniπ N (6)
And the corresponding eigenvalue of A 0 is:

µ m = αλ m - 1 τ (7)
In the next section, we will now present our results regarding the stability of the system defined here, as well as the possible stable states it can reach.

III. STABILITY ANALYSIS A. Form of the stable states

All the stable states of the system presented in section II can be seen as the equilibrium point of one of the subsystems defined by a combination of excited and inhibited neurons Σ:

-

1 τ I N + αP Σ W s + P Σ (αb + β)1 = 0 (8)
The equilibrium points of the different subsystems can be easily obtained.

The equilibrium point s 0 of the subsystem with P Σ = 0 (all neurons inactive) is

s 0 = 0 (9)
The equilibrium point s cons of the subsystem with P Σ = I N (all neurons active) is a state on which all neurons share the same activity, due to the circulant property of W:

s cons = αb + β 1 τ -αλ 0 1 (10)
The equilibrium point s Σ of a system with P Σ = I N , P Σ = 0 (at least one active and one inactive neuron) verifies:

• For all inhibited neuron k, s Σ k = 0. • As for the excited neurons, sΣ = Ã-1 Σ (αb + β)1 where sΣ is the vector obtained by deleting all the elements corresponding to inhibited neurons from s Σ .

However, not all the equilibrium points of the different subsystems are also stable points for the whole system. For example, the state defined by s 0 = 0 is the equilibrium point of the subsystem when all neurons are inactive, but it is not a stable point of the whole system, as is this state ṡ0 = -1 τ s 0 + φ(Ws 0 + b) = (αb + β)1 = 0 (because b, α > 0, and β ≥ 0).

In fact, an equilibrium point s Σ of a subsystem defined by Σ is a stable point for the whole system if and only if it verifies:

∀k active in Σ, s Σ k > 0 and W k s Σ + b ≥ 0 (11) ∀k inactive in Σ, s Σ k = 0 and W k s Σ + b < 0 (12)
where W k denotes the k th line of W. In this section, we will focus on the two following forms of the stable points of the whole system: Definition 1: A stable state of the system is called a consensus if all the neurons share the same non-zero activity. A stable state which is not a consensus is called a bump.

An example of consensus and bump are shown on Fig. 1.

It can be easily checked that the only possible consensus of the system is the state s cons defined previously ((10)).

As ring network models are generally used to encode information in the form of localised activity around some angular position, we want to achieve an equilibrium point which is not a consensus.

B. Sufficient condition for stability

In [START_REF] Richard | Permitted and Forbidden Sets in Symmetric Threshold-Linear Networks[END_REF], a necessary and sufficient condition for the stability of the system was presented. More precisely, with our notations and with β = 0, it can be shown that:

Theorem 1: The following statements are equivalent: 1) All positive eigenvectors of all submatrices of -A 0 have positive eigenvalues.

2) The matrix -A 0 is copositive; that is, x T (-A 0 )x > 0 for all nonnegative x, except x = 0. 3) For all constant b and all initial conditions, the dynamics converge to an equilibrium point. However, we could not apply this reasoning with β > 0, and even with β = 0, the copositivity of the matrix -A 0 is not something easy to ensure when tuning the connectivity matrix W. What's more, it provides no information about the form of the solution at the stable points.

Then, let us first propose a sufficient condition for the stability of the system.

Theorem 2: If all the eigenvalues of W are strictly less than 1 ατ , then the system converges to an equilibrium point, for all constant b and all initial conditions. Under these conditions, the consensus s cons is an equilibrium point of the system if and only if λ 0 ≥ -b βτ . Proof: Let us suppose that the connectivity matrix W's eigenvalues verify: ∀m ∈ 0; N -1 , λ m < 1 ατ Then all the eigenvalues µ m of A 0 are strictly negative. This also implies that all the possible matrices A Σ = -1 τ I N + P Σ W defining our systems are negative definite. The function V (x) = x T x is then a common Lyapunov function for all our subsystems, as we have A T Σ + A Σ negative definite for all possible A Σ . This ensures that the system converges to an equilibrium point, for all constant b and all initial conditions.

Let us now prove that in that case, the consensus s cons is an equilibrium point of the system if and only if we also have λ 0 ≥ -b βτ . Let us consider a network that verifies ∀m ∈ 0; N -1 , λ m < 1 ατ , and the state s cons in which for all neuron, the neural activity is equal to:

s cons = αb + β 1 τ -αλ 0 > 0 (13)
We can compute for each neuron k (k ∈ 0; N -1 ):

W k s cons + b = N -1 i=0 w ik s i + b = s cons N -1 i=0 w ik + b
As W is circulant, we have ∀k ∈ 0; N -1 , N -1 i=0 w ik = λ 0 . Therefore:

W k s cons + b = s cons λ 0 + b = b + βτ λ 0 1 -ατ λ 0
As explained in III-A, this state is a stable state for the system if and only if W k s cons + b ≥ 0 for all k. Therefore the state s cons is a stable state of the system if and only if λ 0 ≥ -b βτ .

Summarizing this subsection, we have proposed a conservative sufficient for stability (λ i < 1/(ατ ), ∀ i) corresponding to region 1 in figure 2 (splitted in two subregions depending on the type of equilibrium point). Moreover, we have defined strict necessary and sufficient conditions for the existence of a stable consensus (region 1a).

The next subsections partially tighten the stability conditions for the other possible equilibirum (bump), both analytically (subsection C) and numerically (section IV).

C. Necessary condition for stability

In this section we will prove the following statement: Theorem 3: If λ 0 ≥ 1 ατ then at least one neuron's activity diverges to infinity, for all constant b and all initial conditions (we will say that the system diverges).

Such a network, with neural activity that diverges to infinity, is of course highly unrealistic, so the region λ 0 ≥ 1 ατ should be avoided when choosing the synaptic connectivity matrix.

Proof: From the definition of the activation φ (equation ( 3)), because α and β are positive, it is clear that: ∀x ∈ R, φ(x) ≥ αx Therefore, we can write the following inequality:

N -1 k=0 ṡk ≥ N -1 k=0 - 1 τ s k + α( N -1 i=0 w ki * s i + b k ) (14)
Also, because the synaptic connectivity matrix W in circulant, we can show that:

N -1 k=0 N -1 i=0 w ki * s i = N -1 i=0 w i N -1 k=0 s k = λ 0 N -1 k=0 s k
Then, we can deduce from (14) that:

N -1 k=0 ṡk ≥ - 1 τ + αλ 0 N -1 k=0 s k + N αb
As all the neuronal activities s k are positive, and the stimulation b is strictly positive, if the connectivity matrix W is such that λ 0 ≥ 1 ατ , at any given time we will have

N -1
k=0 ṡk > 0 which means that the system will diverge to infinity.

The eigenvalue λ 0 therefore plays a key role in the behaviour of the system, as it provides a necessary condition for stability of the system as well as necessary and sufficient condition for the existence of a consensus. As this eigenvalue is obtained from the sum of synaptic weights, it can be interpreted in terms of balance between excitation and inhibition in the network: stronger excitatory connections in the network (which increase λ 0 ) tend to make the network diverge, whereas stronger inhibitory connections (which decrease λ 0 ) tend to stabilize it and promote localized activity, as could be expected.

D. What about the remaining case ?

We will now focus on the case that was not addressed before, that is when W has at least one eigenvalue which is greater than 1 ατ and λ 0 ≤ 1 ατ . Under these conditions, we could not find any more specific results regarding the stability of the system than those stated in Theorem 1, but we addressed the issue of the form of the solution when the system converges to an equilibrium point.

If λ 0 ≥ -b βτ , as we also have λ 0 ≤ 1 ατ , then it can be easily shown that the consensus defined by s cons = αb+β 1 τ -αλ0 1 is an equilibrium point of the system (with a similar reasoning as in the proof for Theorem 2). However, as W has at least one eigenvalue that is strictly greater than 1 ατ , it implies that A 0 has at least one eigenvalue which is strictly greater than 0, and so this consensus is an unstable fixed point.

Therefore, when the system converges to an equilibrium point, if the initial state was not s = s cons with s cons as defined in [START_REF] O'keefe | The Hippocampus as a Cognitive Map[END_REF], or if λ 0 ≤ -b βτ then this equilibrium point must be a bump.

In order to verify our results, as well as tighten our stability conditions for the case λ 0 ≤ 1 ατ but ∃i = 0, λ i > 1 ατ , we focused on studying a network in which the function used to generate the synaptic connectivity matrix W was known, and performed numerical simulations. Note that plausible synaptic connectivity matrices will usually have one mode (connection weights decrease monotonically with the distance). In this case, λ 1 will be the highest eigenvalue (corresponding to the lowest non null frequency revealed by the Fourier transform), which makes it possible to analyze the stability of the system regarding this eigenvalue only.

More precisely, the synaptic connectivity matrix W was built from a gaussian function as follows:

∀i ∈ 0; N -1 , ∀j ∈ 0; N -1 , w ij = f (θ i -θ j ) + µ ( 15
)
with f : x -→ e -x 2 2σ 2 if x = 0 0 if x = 0
This connectivity represents a ring in which all neurons are excitatory, with an excitatory synaptic interaction f decreasing as a gaussian of width σ with distance. Because neurons usually don't send synapses to themselves, the excitatory synaptic interaction is equal to zero for i = j. The presence of inhibitory interneurons with a uniform action on the network is represented by introducing effective inhibition as a global shift of µ < 0 of the synaptic connectivity, as suggested in [START_REF] Gerstner | Neuronal Dynamics, From single neurons to networks and models of cognition[END_REF]. The connectivity matrix W can then be tuned by chosing only the two parameters σ and µ appropriately.

The parameters of the network were chosen as follows: N =1000, τ =10ms, α=2, β=10, b=1. The initial state of the network was chosen randomly, and the differential equations were solved using Euler method with a time step of 0.5ms.

All the simulations were performed using Matlab (MAT-LAB 8.1, The MathWorks Inc., Natick, MA, 2000).

With this synaptic connectivity matrix, for a large N , the eigenvalues of W can be approximated using the error function defined by erf(x) = 1 √ π

x -x e -t 2 dt. More precisely, we have:

λ 0 ≈ (N -1)µ + N σ √ 2π erf π σ √ 2 (16)
and, among all the other eigenvalues of W, the largest one is always:

λ 1 ≈ N σ 2 √ 2π e -σ 2 erf π σ √ 2 + i σ √ 2 +erf π σ √ 2 -i σ √ 2 + µ (17)
We could then define the behaviour of the system in the following regions, depending on σ and µ:

• Region 1a (λ 0 < 1 ατ , λ 1 < 1 ατ and λ 0 ≥ -b βτ ): the system converges to a consensus. • Region 1b (λ 0 < 1 ατ , λ 1 < 1 ατ and λ 0 < -b βτ ): the system converges to a bump.

• Region 2 (λ 0 ≥ 1 ατ ): the system diverges to infinity.

As for Region 3, characterized by λ 0 < 1 ατ and λ 1 ≥ 1 ατ , even though we could not obtain theoretical results, we used numerical methods to predict the stability of the network as well as the exact form of the solution.

To do this for a given matrix W, we considered all the different combinations of active and inactive neurons Σ on the ring (with regard to the rotation invariance of the network), such that ÃΣ was negative definite (with ÃΣ as defined in II-C). For each of these Σ, we computed the theoretical state s Σ in which the network would be if it were to stabilize with active and inactive neurons as defined by Σ:

• for each inactive neuron k, s Σ k = 0 • sΣ = Ã-1 Σ (αb + β)1
, with sΣ the vector obtained by selecting only the indices corresponding to active neurons in s Σ . We then checked the nullity of ṡΣ .

When no Σ could be found for which ṡΣ = 0, we concluded that the network defined by W was unstable (Region 3a), whereas on the other case, we listed the possible forms in which the network could stabilize (Region 3b).

In Region 3b, our observations were that only one form of solution could be stable for a given matrix W, and this form always corresponded to a Σ where all active neurons were grouped together on the ring, thus forming a single bump of activity on the ring.

With our simulations, we could then confirm our theoretical predictions and complement our numerical observations (as summed up in Fig. 2):

• In Region 1a, the system converges to a consensus.

• In Region 1b, the system converges to a bump.

• In Region 2, the system diverges to infinity.

• In Region 3a, the system diverges to infinity.

• In Region 3b, the system converges to a bump of the predicted form. In order to build a network that will converge to a bump, one should then choose the synaptic weight matrix parameters σ and µ corresponding to region 1b or 3b.

V. CONCLUSION

In this work, we studied the stability of ring neural networks with piecewise-affine, discontinuous activation functions. In comparison to existing literature, we obtained some simpler necessary conditions and sufficient conditions for convergence of the system depending on the eigenvalues of the synaptic connectivity matrix, and we were able to link them to the form of the solutions. When these theoretical conditions were not met, we were able to predict the behaviour of the system numerically on an example network. All of these results were then confirmed using simulations of a network with gaussian-like connectivity matrix.

Understanding the theoretical aspects of the behaviour of such systems could prove very useful in computational neurosciences, in order to build more complex models of place cells or head-direction cells, and easily tune them to resemble physiological networks.

As for future works, it would be interesting to see if the results of this study can be extended to activation functions with more than two affine regions (e.g., with some saturation), to piecewise non-linear activation functions, or to systems with more complex dynamics (e.g., with some noisy input, or with delays on the synaptic transmission). Another perspective would be to study this network in terms of memory, i.e., how long the network can maintain accurate localized activity when the stimulation is removed.

Fig. 1 .

 1 Fig. 1. Example of the two types of stable points for ring networks, achieved from the same initial conditions with two different synaptic connectivity matrices. Top plot: Initial conditions. Bottom plot: two stable states -consensus (red dashed line), and bump (blue solid line).

Fig. 2 .

 2 Fig. 2. Stability analysis of a neural network with a Gaussian connectivity matrix of parameters σ and µ. Red circled line: λ 0 = 1 ατ ; Blue dashed line: λ 1 = 1 ατ ; Magenta squared line: λ 0 = -b βτ ; Black line: Border of the region of convergence, as obtained in simulations; Green dash-dot line: Border of the region of convergence, as predicted numerically.

Laure Buhry is with the Laboratoire Lorrain de Recherche en Informatique et ses applications (LORIA UMR 7503, University of Lorraine-