B Tounsi

Y Hayel

D Quadri

L Brotcorne

T Jimenez

Sensitivity analysis of stochastic user equilibrium and its application to delivery services pricing

Keywords: Stochastic user equilibrium, sensitivity analysis, bi-level programming, service network design

Introduction

In e-commerce the delivery of products is a crucial part for the success of a shop. An efficient delivery system should offer various services and predict customers' behavior. The latter are influenced by the price of a delivery service, but also by its quality (congestion effect induced by customers' choices). In this study, we introduce a bilevel model to optimize a delivery system. At the upper level, the provider controls services' tariffs. At the lower level, users react by choosing their delivery service according to an utility function which incorporates the provider's tariffs and the congestion effect. We model the customers' reaction using Stochastic User Equilibrium (SUE). We also present a sensitivity analysis for the SUE that gives the explicit expression of the derivatives of customers distribution with respect to services' tariffs. Based on a local search that exploit the derivatives information, a new heuristic algorithm for the bilevel delivery services pricing problem is developed and compared to others existing methods.

The delivery system

For the delivery of ordered products, customers choose among a set J of services. We consider two types of services: delivery at home (DH) and delivery at relay station (RS). The disutility of service j perceived by a customer depends on tariff A j set by the provider, and on the quality of service or the congestion effect induced by customers' choices. Denoting by λ the customers arrival rate, and by p j the fraction of customers choosing service j ∈ I, the general form of the disutility

c i of service i is c i (p i ) = A i + α i f i (p i )
, where f is the congestion function of service i and α i a monetary conversion coefficient .

-DH service is modeled using an M/D/1 queue with a First-In-First-Out (FIFO) service discipline. The congestion function is given by the service delay

f (p i , λ, d) = d + λp i d 2 2(1 -λp i d)
.

-RS service is modeled using a M/M/K/K queue where K is the capacity of the relay. The congestion function is given by the Erlang's formula

f (p i , λ, K) = (λpi) K K! K j=0 (λpi) j j!
A first way to model customers' behavior is the Deterministic User Equilibrium (DUE) also known as the Wardrop equilibrium. Here, customers have full rationality and select the service with the smallest cost. In the real world, the full rationality assumption appears to be very strong (lack of information, individual preference, ...). It is more accurate to consider that users can make errors in their decision. Thus we use Stochastic User Equilibrium (SUE) to model user's decision process. The Logit based model (L-SUE) is often used for this issue. In L-SUE, the error term follows a Gumble distribution. In this work we compare LSUE to a Nested-Logit based SUE (N-SUE), more appropriate to the case study. In the latter, correlated services are grouped in nests. For our application, DH services form a nest, and RS ones form another nest. We consider that, at first, a customer selects a nest, and then he selects the delivery service. In the two cases, The SUE is computed using the Method of Successive Averages (Sheffi 1995, [START_REF] Bekhor | Application of cross-nested logit route choice model in stochastic user equilibrium traffic assignment[END_REF]).

The control problem

Knowing that customers follow a stochastic user equilibrium, the provider can control services parameters in a way that optimize some criteria. This yields to a Mathematical Program with an Equilibrium Constraint (MPEC). The provider (leader) controls service tariffs A i , delivery capacities D i and relay capacities K i . Denoting by u the leader's variables, u = (A i , D i , K i ), and by p the user equilibrium, p(u) = (p i (u)) i∈I . (MPEC) can be formulated as follows:

(M P EC)    max F (u, p) (1) s.t L i ≤ u i ≤ U i ∀i ∈ I (2) p = p ⋆ (u) (3)
F (u, p) = i∈I λA i p i -a i D i -b i K i . F (u, p) = i∈I p i (u)c i (u, p).
Where constraint (3) describes that p is a solution for equilibrium problem parameterized by u. The leader's objectif function can be maximizing profit F (u, p) or minimizing the total cost F (u, p).

The resulting program is in general non convex, and hence it is difficult to find a global optimum [START_REF] Yang | An algorithm for the inflow control problem on urban freeway networks with user-optimal flows[END_REF]. Sensitivity Analysis based heuristic (SABH), has been used to find efficient solution for such MPEC in the field of traffic control [START_REF] Yang | An algorithm for the inflow control problem on urban freeway networks with user-optimal flows[END_REF][START_REF] Friesz | Sensitivity analysis based heuristic algorithms for mathematical programs with variational inequality constraints[END_REF]). The key point of SABH lies in how to evaluate the user equilibrium (follower variables) at a given point of services' tariffs (leader variables). Sensitivity analysis gives explicit expression of the derivatives of the formers with respect to the latters. Having these informations, a descent direction (the gradient of F (u, p)) is build and used to update leader variables [START_REF] Friesz | Sensitivity analysis based heuristic algorithms for mathematical programs with variational inequality constraints[END_REF]). This process is iterated until no improvement of the leader's objective value is provided.

In this work we apply this method for the resolution of a bilevel program with a stochastic user equilibrium (MPSEC), in the case where the provider controls only services' tariffs. We also, introduce a new way to find the descent direction based on a local search that exploits sensitivity informations. Numerical experimentations show better results for the new descent direction regarding the quality of the solution and the execution time.

Further works will treat the case where the provider controls tariffs and capacities of services. This leads to a challenging MPSEC where a part of the upper level variables are integer.

Conclusion and perspectives

This work combines Bilevel Programming and Stochastic User Equilibrium (SUE) for the optimization of a delivery services system. A new heuristic algorithm based on sensitivity analysis is presented. Closely to the delivery system, futur works could consider heterogenous customers in a multi-class model.


A Le Yaouanc

D Melikhov

V Morénas

L Oliver

O Pène

J.-C Rayna

P V Landshoff

Duality in semileptonic inclusive B-decays in potential models: regular versus singular potentials

Making use of the nonrelativistic potential model for the description of mesons, and working in the Shifman-Voloshin limit, we compare the integrated rate Γ(B → X c lν) calculated as a sum of the individual decay rates to the quantum-mechanical analog of the OPE. In the case of a potential regular at the origin, we find a well-defined duality violation, which is, however, exponentially small. It corresponds to the charm resonances kinematically forbidden in the decay process, but apparently picked up by the OPE. For singular potentials, we do not obtain a full OPE series, but only a limited Taylor expansion, since the coefficients become infinite beyond some order. In this case, we do not find an indication of duality violation: the difference is smaller than the last term of the limited expansion. This emphasizes that the case of singular potentials, which may be relevant for QCD, deserves further study.

The theoretical framework based on the Operator Product Expansion (OPE) determines in QCD the heavy meson inclusive decay rate as series in inverse powers of the heavy quark mass, with the coefficients proportional to the meson matrix elements of the local operators of increasing dimensions [1,2]. The calculation is based on representing the decay rate as the contour integral in the complex q 0 -plane. The OPE makes the contour integrals easily calculable term by term and provides the decay rate as a 1/m Q series.

There are, however, potentially dangerous points in this calculation: E-mail address: melikhov@thphys.uni-heidelberg.de (D. Melikhov).

(i) the OPE series is at best asymptotically convergent even for large absolute values of the complex q 0 ;

(ii) the integration contour for the decay rate contains a segment near the physical region, where the OPE cannot be justified [1].

This might lead to the violation of duality for the decay rate, i.e., to the difference between the OPEcalculated decay rate and the result of summing the individual decay rates of the opened channels. This issue was also discussed by N. Isgur [3].

In this Letter we discuss the semileptonic decay rate in the small velocity limit and use the nonrelativistic potential model for the description of mesons. We perform a short-time expansion in operators of increasing dimensions which we call OPE and which has indeed some common features (but also important dif-ferences) with the OPE expansion in the field theory. We consider the two cases: regular confining potentials 1 and singular potentials. Note that, in the past, the nonrelativistic potential model has been widely used in the discussion of QCD sum rules [4,5]. The question of the distinction between regular (harmonic) and singular (Coulomb) potentials has been emphasised by Novikov et al. [4] (in the context of heavy quark QCD sum rules).

We work in the Shifman-Voloshin (SV) limit Λ ≪ δm = m bm c ≪ m c ,m b . In this limit both the amplitude and the decay rate can be formally obtained as a double expansion in 1/m c and 1/δm. We consider lowest orders in 1/m c ,u pt o1 /m 2 c ,a n dall orders in 1/δm. Note that this involves terms of much higher order than usually done when one expands in 1/m Q with m b /m c ,o ra sw e l lδm/m Q , fixed. Our double expansion allows on the contrary to go much further in 1/δm, and this might allow to display subtle duality violations.

For the regular potential we obtain the full 1/δm expansion, which is only asymptotic to the physical width expanded to the same order in 1/m c .T h e difference 2 is of order δm/m 2 c exp(-δm/Λ),w h i c h means exponentially small duality violation.

For the singular potential we do not obtain the full 1/δm expansion: following the same procedure as for the regular potential leads to infinite coefficients beyond some order in 1/δm. In this case, we find that the truncated expansion satisfies duality up to this order.

We consider the inclusive semileptonic decay B → X c lν in the SV limit and treat mesons as nonrelativistic bound states of spinless quarks in a confining potential (a detailed calculation is given in [6]). This model maximally simplifies both constructing the OPE series and calculating the sum of the exclusive channels. For the sake of argument we consider the case of leptons coupled to hadrons through the scalar current. In this case the leptonic tensor is reduced to a scalar function L(q 2 ). The amplitude T depends on the two variables, and we choose them as q 0 and q 2 in 1 A regular potential is a potential which is an analytic function of r at r = 0. For example, the potential V(r)≃| r| falls out of this class. 2 As we shall see this expansion contains only a finite number of nonzero terms.

Fig. 1. Singularities of the amplitude T(q 0 , q 2 ) in the complex q 0 -plane. Circles are poles, corresponding to low-lying charm states, and the cross marks the location of the pole in the free b → c quark transition.

the B-rest frame:

T q 0 , q 2 = 1 i dx e -iqx B|T J(x),J + (0) |B (1) = X | B|J |X(-q) | 2 M B -E X (-q) -q 0 .
The sum in (1) runs over all hadron states with the appropriate quantum numbers. The states are normalized as follows p| p ′ =(2π)3 δ( pp ′ ),a n d E X (-q) is the energy of the state X with the total 3-momentumq. At fixed q 2 , T(q 0 , q 2 ) has a cut in the complex q 0 -plane along the real axis for q 0 <M B -M Dq 2 /2(m c + m d ), see Fig. 1.

A part of this cut for | q| <q 0 <M B -M Dq 2 /2(m c + m d ) corresponds to the decay process. The decay rate can be represented as the contour integral in the complex q 0 -plane over the contour C( q 2 ) (Fig. 1)

(2) Γ(B→ X c lν) = d q 2 | q| C( q 2 ) dq 0 2πi L q 2 T q 0 , q 2 .
The contour C( q 2 ) selects at any given q 2 only states kinematically allowed in the decay B → X c lν.I ti s tightly attached to the points P ± with the coordinates (| q|, ±i0), otherwise it can be freely deformed in the region where the function T 0 (q 0 , q 2 ) is analytic. The amplitude can be expanded in a series

(3) T q 0 , q 2 = i c i q 0 , q 2 B| O i |B ,
where O i are operators of increasing dimensions and c i (q 0 , q 2 ) are the c-number coefficients. Introducing the expansion (3) into (2) gives the integrated rate as an OPE series. 3

The model

Let us proceed along the lines of Ref. [6]. We treat the leptonic part relativistically, but for the description of mesons as bound states of spinless quarks use the nonrelativistic potential model with a confining potential. We consider the decay in the B-rest frame. The Hamiltonian of the b q system at rest has the form

H bd = m b + m d + ĥbd , ĥbd = k 2 /2m b + k 2 /2m d + V bd (r), such that ( ĥbd -ǫ B )|B =0, ( H bd -M B )|B =0, M B = m b + m d + ǫ B .
The Hamiltonian of the c q system produced in the semileptonic b → clν decay reads

H cd ( q) = m c + m d + ( k + q) 2 /2m c + k 2 /2m d + V cd (r).
The eigenstates of this Hamiltonian are |D n ( q) such that

H cd ( q) -E D n ( q) D n ( q) = 0,
where

E D n ( q) = M D n + q 2 /2(m c + m d ), M D n = m c + m d + ǫ D n .
The Q q potential can be expanded as follows:

V Qq = V 0 + V 1 /2m Q + V 2 /2m 2 Q +•••.
contributions of distances below the scale 1/µ are referred to the Wilson coefficients, while contributions of distances above this scale are referred to the matrix elements. As a result, both the Wilson coefficients and the matrix elements acquire the µ-dependence.

In the potential model, the T -product of the two currents is also expanded in a series of operators of increasing dimensions, but the resulting c-number coefficients and the matrix elements in Eq. ( 3) are scale independent.

Sum rules

The relationship between the sum over the individual channels and the meson matrix elements of the operators is established by the sum rules. Let us introduce δ n ( q) through the relation (δ

n ( q) = ǫ D n -ǫ B - q 2 m d /2m c (m c + m d )) (4) M B -q 0 -E D n ( q) = δm -q 0 -q 2 /2m c -δ n ( q).
The δ n ( q) is the eigenvalue of the operator δH( q)

(5)

M B -H cd ( q) = δm -q 2 /2m c -δH( q)
with |D n ( q) the corresponding eigenstates. The sum rules are obtained by inserting the full system of the eigenstates |D n ( q) into B| (δH( q)) i |B :

(6) B| δH( q) i |B = ∞ n=0 F n ( q) 2 δ n ( q) i ,
where F n ( q) = B|D n ( q) is the B → D n transition form factor. This relation represents the sum over all c d resonances in terms of the B-meson matrix element of the operators (δH ( q)) i . For the potential regular at the origin r = 0 the sum over n is convergent for any i, whereas for the singular potential both sides of Eq. ( 6) are convergent for small i and diverge for large i.

At the moment we proceed formally and discuss this problem in more detail in Section 5.

Duality relation for the amplitude

Making use of the sum rules (6), we represent the amplitude as a sum of the operators:

T(q 0 , q) (7) = ∞ n=0 |F n ( q)| 2 M B -q 0 -E n ( q) (8) = 1 δm -q 2 2m c -q 0 n=0 i=0 |F n ( q)| 2 (δ n ( q)) i δm -q 2 2m c -q 0 i (9) = 1 δm -q 2 2m c -q 0 i=0 B|(δH ( q)) i | B δm -q 2 2m c -q 0 i .
This expression is the duality relation for the amplitude: the sum (7) runs over the infinite number of the charm resonances, and the sum (8) runs over the infinite number of the operators of the increasing dimensions (the OPE series). In fact, the location of singularities in the complex q 0 -plane in the series ( 7) and ( 9) is quite different: in [START_REF] Shifman | Talk at the Workshop Continuous Advances in QCD[END_REF] it is an infinite set of single poles at the different locations corresponding to different charm resonances, and in (8) it is an infinite set of poles of the increasing order at the same point. However, this set of equations is only a formal one; in fact, ( 7) is a summable series leading to a finite result in all cases; on the other hand, the situation of Eq. ( 9) is more subtle. In the singular case, the coefficients are infinite beyond some order, and one must accordingly truncate the series. In the regular case, the Eq. ( 9) is only an asymptotic series: notice that the geometric sum over i in Eq. ( 8) has a domain of convergence which is repelled to infinity with n.

Let us illustrate it with a simple example: assume that F 2 n ≃ e -n and E n ≃ n. Then the analog of the above equations takes the form ( 10)

∞ n=0 e -n z -n = n=0 e -n z i=0 n z i ≃ 1 z i=0 i! z i .
The last step is obtained by changing the order of summation and using the relation ∞ n=0 e -n n i ≃ i! The series (10) in i is only asymptotic and not even Borel summable.

Such a factorial divergence appears in the example of M.A. Shifman [START_REF] Shifman | Talk at the Workshop Continuous Advances in QCD[END_REF]. In this example, the residues are constant with the excitation number. This would imply that the direct method of our Letter will give formally infinite coefficients for the expansion, which means that it fails completely. In our example, on the contrary, the residues are rapidly decreasing, and all the coefficients are finite. The interest of our example is to show that nevertheless a factorial divergence is also appearing.

From the amplitude T under the form Eq. ( 7) or Eq. ( 9), respectively, by integration over the same contour C, we can obtain either the width as a sum over the exclusive final states, or as the OPE series. The expression ( 9) is an accurate approximation to [START_REF] Shifman | Talk at the Workshop Continuous Advances in QCD[END_REF] only when q 0 is far from the singularities of T(q 0 , q). The contour C can be deformed away from the singularities except near its fixed end points. When integrating over q 0 this is a possible source of discrepancy, i.e., of duality violation. Consequently, we are now going to estimate the integral of expression [START_REF] Shifman | Talk at the Workshop Continuous Advances in QCD[END_REF], i.e. the sum over the exclusive channels, and the integral of expression (9), i.e., the OPE prediction, and compare both results.

The OPE calculation of the decay rate

Let us first proceed with the amplitude in the form ( 9) and obtain the OPE expression for the decay rate. We consider the leptonic tensor of the general form L(q 2 ) = (q 2 ) N . For technical reasons, it is convenient to isolate h bd in the expression for δH( q) as follows ( 11)

δH( q) = h bd -ǫ B + k q m c + 1 m c - 1 m b k 2 + V 1 2 .
Substituting ( 11) in ( 9) and performing the necessary integrations gives a series in 1/m c [6]

Γ OPE (B → X c lν) Γ(b→ clν) = 1 + B| k 2 |B 2m 2 c -(2N + 3) B|V 1 |B 2m 2 c ( 12
)
+ 2N+3 i=1 (-1) i C i+2 2N+5 2N + 5 B| O i |B m 2 c δm i + O Λ 2 δm m 3 c , with C i n = n! i!(n-i)! and O i = k(h bd -ǫ B ) i k.
A n important feature of the OPE series (12) is that the leading-order term reproduces the free-quark decay rate, and the first correction emerges only in the 1/m 2 c order (cf. [1,2]).

Summation of the exclusive channels

Now let us sum the rates of the exclusive channels. The B → D n transition form factors have the form [6]

F 2 0 ( q) = 1 -ρ 2 0 q 2 /m 2 c + O q 4 /m 4 c + O δm 2 β 2 /m 4 c , F 2 n ( q) = ρ 2 n q 2 /m 2 c + O q 4 /m 4 c + O δm 2 β 2 /m 4 c .
Since | q| δm in the decay region, these expressions allow calculating the decay rate to the accuracy δm 2 /m 2 c . Explicitly, we obtain [6]:

Γ(B→ D 0 lν) Γ(b→ clν) = 1 - 3ρ 2 0 2N + 5 δm 2 m 2 c + 3 2 m d 1 + m d /m c δm m 2 c -(2N + 3) B| k 2 + V 1 |B 2m 2 c , Γ(B→ D n lν) Γ(b→ clν) = 3ρ 2 n 2N + 5 δm 2 m 2 c -3 ρ 2 n ∆ n δm m 2 c (13) + 1 m 2 c 2N+5 i=2 (-1) i C i 2N+5 2N + 5 (3ρ 2 n ∆ i n ) δm i-2 , where ∆ n = ǫ D n -ǫ D 0 .
The main contribution is given by the B → D 0 transition. Excited states contribute only starting from the (δm) 2 /m 2 c order in the SV limit. Notice that each of the exclusive rates contains terms of the order δm 2 /m 2 c and Λδm/m 2 c which are absent in the OPE series.

Summing over all opened exclusive channels gives

Γ(B→ X c lν) Γ(b→ clν) = 1 - δm 2 m 2 c 3 2N + 5 ρ 2 0 - n max n=1 ρ 2 n + 3 δm m 2 c m d /2 1 + m d /m c - n max n=1 ρ 2 n ∆ n -(2N + 3) × B| k 2 + V 1 |B 2m 2 c (14) + 2N+5 i=2 (-1) i C i 2N+5 2N + 5 3 n max n=1 ρ 2 n ∆ i n m 2 c δm i-2 .
The sum over the charm resonances is truncated at n max , which is the total number of the resonance levels opened at q 2 = 0. For the confining potential and in the SV limit n max is found from the relation ∆ n max ≃ δm.

Check of duality for regular potentials

The transition radii in the expression [START_REF] Schwinger | Particles, Sources and Fields[END_REF] are not independent and related to each other through the sum rules. These sum rules can be obtained from (6). Expanding both sides of (6) in powers of 1/m Q and taking the linear q 2 term gives the set of the sum rules [6]: for i = 0 one finds the Bjorken sum rule [START_REF] Bjorken | Invited talk at Les Rencontres de Physique de la Vallée d'Aoste[END_REF], for i = 1 -the Voloshin sum rule [9], for i 2 -higher moment sum rules:

i = 0: ∞ n=1 ρ 2 n = ρ 2 0 , i = 1: ∞ n=1 ρ 2 n ∆ n = m d /2 1 + m d /m c , (15) i 2: ∞ n=1 ρ 2 n ∆ i n = 1 3 B| k(h bd -ǫ B ) i-2 k|B .
Using these relations to rewrite the OPE result (12) as the sum over hadronic resonances, the difference between the OPE and the exclusive sum (the dualityviolating contribution) explicitly reads

δ Γ ≡ Γ OPE (B → X c lν) -Γ(B→ X c lν) Γ(b→ clν) = 3 δm 2 m 2 c 2N+5 i=0 (-1) i C i 2N+5 (2N + 5)δm i n>n max ρ 2 n (∆ n ) i + O Λ 2 δm m 3 c = δm 2 m 2 c 3 2N + 5 n>n max ρ 2 n 1 - ∆ n δm 2N+5 (16) + O Λ 2 δm m 3 c .
Quite remarkably, δ Γ happens to be equal to the sum of the extrapolated widths for charm states beyond the kinematical limit. A similar expression is found in QCD 2 [12]. Clearly, the duality-violating effect is connected with the charm states forbidden kinematically in the decay process. Notice that Γ OPE (B → X c lν) -Γ(B → X c lν) < 0, because ∆ n >δm for n>n max ,and2N + 5 is odd.

To estimate the size of the duality-violation effects, the behavior of the transition radii and the relation between ∆ n and n max , which will be given by the behaviour of the excitation energies also at large n, are needed. For quite a general form of the confining potential we can write the following relations for ∆ n at large n∆ n ΛCn a for n>n max and ∆ n max = ΛC(n max ) a ≃ δm, with C and a some positive numbers. In particular, this estimate is valid for the confining potentials with a power behavior at large r.T h i s estimate for ∆ n is only depending on the behaviour of the potential at large distances.

The behavior of the radii ρ 2 n at large n are then connected with the finiteness of the r.h.s. of the sum rules (15): for a potential regular at r = 0, the matrix elements in the r.h.s. of the sum rules are finite for any i, which means that the radii ρ 2 n are decreasing with n faster than any power. Essentially this means that ρ 2 n ≃ exp(-n), and, therefore, the duality-violating effect in the decay rate in (16) is of order δ Γ ≃ δm 2 /m 2 c exp(-δm/Λ). One of such examples, the harmonic oscillator potential, is discussed in [10].

Singular potentials

However, if the potential is singular at r = 0, the situation changes dramatically. First, only a few first number of the matrix elements B| O i |B are finite. 4 We can try to proceed along the same lines but then have to truncate the series in 1/δm at the last finite term. We want to estimate the difference between this truncated series and the exclusive sum.

Let us illustrate this considering a potential with a Coulomb behavior at small r, V ≃α/r,a n d confining at large r.T h e n B| k(h bdǫ B ) i k|B are finite for i 1, but diverge starting from i = 2. We then find that 4 The appearance of infinite coefficients in the OPE series is probably due to a breakdown of the power series expansion, for instance by fractional powers or logarithms of m Q as seems to be the case in the pure Coulomb case [11]. For similar phenomena in a perturbation expansion, see [13].

In particular, it has been noted by Novikov et al. [4] that, in the case of the Coulomb potential, half integral powers of the expansion parameter occur in the calculation of the moments of the ratio R. However, in our case, preliminary work shows that there are not only half integer powers of Λ/δm, but also logarithms of this quantity. In fact, we note in this respect that the full Schwinger formula for the ratio R contains also logarithms in addition to half integer powers [START_REF] Schwinger | Particles, Sources and Fields[END_REF]. .

More generally, if the above matrix element begins to diverge for some value i = K + 1, the formulas are to be replaced by:

(19) ρ 2 n 1 n 1+ε 1 n a K+2 , (20) δ Γ ≃ Λ 2 m 2 c Λ δm (K+ε/a)
.

Notice that this δ Γ is smaller than the last retained term in the OPE series which is of order

Λ 2 m 2 c Λ δm K .
Therefore, the 'duality violation' is just smaller than the last retained term as for the asymptotic series. This means in fact that there is no indication of duality violation at this computable order. This is independent of a, therefore of the large distance behavior of the potential.

Conclusion

Summarizing our results, the amplitude T(q 0 , q 2 ) (the T -product, Eq. ( 1)) can be expanded in inverse powers of δmq 2 /2m cq 0 , the so-called OPE expansion. Exact duality would mean that the OPE series was convergent and equal to T . Actually, this is not exactly the case. Even in the favourable case of the regular potentials (at r = 0), the OPE series is not convergent, it is only asymptotic to the actual T .F o r singular potentials, the coefficients are simply infinite beyond a certain order.

Besides these problems concerning the amplitude, additional problems appear for the expansion of the width, which is given by a contour integral of T in the q 0 complex plane: the OPE expansion is accurate far from the singularities in q 0 , while the contour has fixed end points in the complex plane close to the singularities (Fig. 1). In view of this situation, we have computed explicitly the difference between the OPE and the actual width. For singular potentials, the series must be truncated, and the difference is found smaller than the last retained term.

As to the perspectives opened by this work, we must first emphasize that singular potentials seem more interesting than regular ones. Indeed, in QCD the effective quark potential is singular, a smoothed Coulomb singularity. Moreover, in QCD 2 , one can suspect some similarity with a linear potential | r|, which is also singular at the origin in the sense of this Letter. For a singular potential, we have seen that the entire series must be truncated at some order, because the coefficients become eventually infinite. We think that such infinite coefficients in an entire series expansion correspond to the fact that the correct expansion is not entire but must include fractional powers and/or logarithms in the expansion parameter, i.e., δm.I nQ C D , one can argue that the operator matrix elements are finite due to renormalisation, but nevertheless the coefficients still contain logarithms of heavy masses. In the nonrelativistic case, the object of the present Letter, the method which has been does not lead to definite conclusions as regards duality for singular potentials: namely, to the order we are able to calculate in this Letter, we find that there is no duality violation, but this leaves open the question of duality violation at some higher order. 5 To proceed further, one would to the order 1/m 4 Q and it may be believed that duality has been fully demonstrated in higher orders [12]. However, a comment is in order here. In [12], it was shown that the matrix element of the leading operator B| QQ|B is dual to the sum of the widths of the full tower of resonances. Therefore, one can suspect that there is a difference between the actual width and the OPE, that is of higher order 1/m 5 Q , corresponding to the extrapolated width of the kinematically forbidden states. This difference, however, has the same order 1/m 5 Q as the matrix elements of the higher dimension operators [12]. It was then assumed that both quantities are dual to each other, but the corresponding OPE coefficients were not calculated and we have not found where this assumption was demonstrated.

have to devise new methods to obtain the above conjectured generalized expansions.

The OPE series in the potential model has an important distinctions from the Wilsonian scheme in the field theory where

In the context of QCD 2 , one has demonstrated duality up

Acknowledgements

We are grateful to G. Korchemsky, B. Stech for discussions, and specially to N. Uraltsev for detailed correspondence. D.M. was supported by BMBF under project 05 HT 9 HVA3. Laboratoire de Physique Théorique is Unité Mixte de Recherche CNRS-UMR8627.


Pascal Chossat

Maciej Krupa

Heteroclinic cycles in Hopfield networks

Keywords: heteroclinic cycles, Hopfield networks, learning rule, network architecture AMS classification: 34C37, 37N25, 68T05, 92B20

Learning or memory formation are associated with the strengthening of the synaptic connections between neurons according to a pattern reflected by the input. According to this theory a retained memory sequence is associated to a dynamic pattern of the associated neural circuit. In this work we consider a class of network neuron models, known as Hopfield networks, with a learning rule which consists of transforming an information string to a coupling pattern. Within this class of models we study dynamic patterns, known as robust heteroclinic cycles, and establish a tight connection between their existence and the structure of the coupling.

Introduction

A simplest example of a heteroclinic cycle is a sequence of saddle type equilibrium points joined in a circle by connecting orbits. Generically heteroclinic cycles are not robust under perturbations of the system (changes of parameters), but for special classes of systems they may occur robustly, typically due to the presence of invariant hyperplanes. Examples of special structures leading to robust heteroclinic cycles are symmetry, the existence of invariant planes corresponding to extinction of some species in Lotka-Volterra systems or existence of synchrony subspaces in coupled cell systems. Heteroclinic networks are a generalization of heteroclinic cycles to sets of equilibria with more complicated connection structure. More generally, heteroclinic sets may consist of invariant sets connecting periodic or more complicated saddle type dynamics. The study of heteroclinic cycles was motivated by examples in fluid mechanics (systems with symmetry) [START_REF] Chossat | Methods in Equivariant Bifurcation and Dynamical Systems[END_REF] and in population biology [START_REF] Hofbauer | Evolutionary Games and Population Dynamics[END_REF]. See [START_REF] Krupa | Robust heteroclinic cycles[END_REF] for an introduction to the subject.

More recently, Rabinovich and co-workers have proposed applications of robust heteroclinic cycles in neuroscience, see [START_REF] Rabinovich | Dynamical Principles in Neuroscience[END_REF] for an early review. Among the contexts proposed in [START_REF] Rabinovich | Dynamical Principles in Neuroscience[END_REF] where heteroclinic dynamics could be relevant were central pattern generators (CPGs) and memory formation. These two applications were validated by some more detailed biological studies [START_REF] Szucs | Robust Microcircuit Synchronization by Inhibitory Connections[END_REF] [START_REF] Bick | Dynamical origin of the effective storage capacity in the brain's working memory[END_REF]. The CPGs are circuits controlling the motoric function, and are known to support a variety of complex oscillations corresponding to different movements of the body. As shown in this paper, by slightly changing the coupling structure in the model one can obtain a variety of heteroclinic cycles and thereby complex periodic solutions. The idea of the memory application is similar -modifications of the coupling, arising from the action of the input, lead to the occurrence of periodic orbits, existing near heteroclinic cycles whose properties reflect the structure of the input.

The focus of this work is to study Hopfield networks, which are the simplest models of memory circuits, with the goal of investigating the presence of heteroclinic cycles. Hopfield introduced thirty years ago [START_REF] Hopfield | Neural networks and physical systems with emergent collective computational abilities[END_REF] this model for learning sequences and associative memory in neural networks, in their simplest possible form. In the continuous time version, for each neuron i in a network of N neurons, the activity is modeled by the following equation (activity model):

ui = -u i - N i=1 J ij g(u j ) + I i , i = 1, . . . N. (1)
where u j is the activity variable (membrane potential) of neuron j and I i is a constant external input on neuron i. The function g : R → (0, 1) is strictly increasing and invertible, for example a sigmoid. The quantity v j = g(u j ) is the firing rate of neuron j, that is the time rate of spikes which are emitted by the neuron. Classically the function g(u) = (1 + e -u ) -1 is used. The coupling coefficients J ij define a N × N matrix J called the connectivity matrix. A positive (resp. negative) coefficient corresponds to an inhibitory (resp. excitatory) input from j to i. When J = 0 all neurons have the same state of rest u i = I i . When coupling is switched on, other equilibria may exist depending on the coefficients J ij . It is often assumed that J is a symmetric matrix, implying that the dynamics of the network always converges to an equilibrium. Each equilibrium is defined by a sequence of values (u 1 , . . . , u N ) called a pattern. Depending upon the inputs I i , one or another equilibrium will be reached, a process that is interpreted as retrieving a pattern which has been earlier memorized through a tuning of the coupling coefficients in the network (Hebbian rule).

The assumption that J is symmetric is unnecessary to storing information, and besides experiments have shown that dynamical patterns are often present in neural circuits and seem to play an important role in various aspects, in particular generating periodic oscillations in CPG. In a recent work Chuan et al. [START_REF] Zhang | Storing cycles in Hopfield-type networks with pseudo inverse learning rule: Admissibility and network topology[END_REF] developed a method of converting strings of information, in the form of sequences of vectors with entries ±1, into coupling matrices so the Hopfield network with the resulting coupling architecture would have storing cycles, i.e. periodic orbits carrying the information of the underlying sequence of vectors. The defining feature of a storing cycle associated to such sequence is that it visits the vicinity of each of its vector, preserving their order in the sequence. In another study Chuan et al. [START_REF] Zhang | Storing cycles in Hopfield-type networks with pseudoinverse learning rule: retrievability and bifurcation analysis[END_REF] used Hopf bifurcation analysis to find storing cycles.

A natural observation is that heteroclinic cycles between equilibria given by the elements of the sequence provide a natural approximation to storing cycles. However the simulations of [START_REF] Zhang | Storing cycles in Hopfield-type networks with pseudo inverse learning rule: Admissibility and network topology[END_REF] and [START_REF] Zhang | Storing cycles in Hopfield-type networks with pseudoinverse learning rule: retrievability and bifurcation analysis[END_REF] gave no evidence of the existence of heteroclinic cycles.

In this work we show that after a small modification the systems studied in [START_REF] Zhang | Storing cycles in Hopfield-type networks with pseudo inverse learning rule: Admissibility and network topology[END_REF] support heteroclinic cycles. Our approach draws on the work of Fukai and Tanaka [START_REF] Fukai | A Simple Neural Network Exhibiting Selective Activation of Neuronal Ensembles: From Winner-Take-All to Winners-Share-All[END_REF], who observed that by of replacing a non-differentiable term in the firing rate equations by a constant one obtains a Lotka-Volterra system, which supports robust heteroclinic cycles [START_REF] Hofbauer | Evolutionary Games and Population Dynamics[END_REF]. This approach was subsequently used by [START_REF] Rabinovich | Dynamical Principles in Neuroscience[END_REF] and [START_REF] Ashwin | Criteria for robustness of heteroclinic cycles in neural microcircuits[END_REF] in their study robust heteroclinic cycles in firing rate models. In this work we continue the approach of [START_REF] Fukai | A Simple Neural Network Exhibiting Selective Activation of Neuronal Ensembles: From Winner-Take-All to Winners-Share-All[END_REF], introducing some refinements to their approximation of the firing rate equations. We point out that the original firing rate equations cannot support hateroclinic cycles due to the presence of non-smooth terms and introduce two methods of regularizing the equations. When the systems studied in [START_REF] Zhang | Storing cycles in Hopfield-type networks with pseudo inverse learning rule: Admissibility and network topology[END_REF] are modified using either of our approaches heteroclinic cycles do exist and there is a direct correspondence between the input string/vector sequence/coupling structure and the resulting heteroclinic cycle. In this work we carry out a detailed study of this correspondence.

Hopfield networks

Storing cycles and network architecture

System (1) is often transformed to the firing rate formulation, by letting the firing rates x i = g(u i ) be the dependent variables. In this section we make the same choice of g as the authors of [START_REF] Zhang | Storing cycles in Hopfield-type networks with pseudo inverse learning rule: Admissibility and network topology[END_REF], namely g(u) = tanh(λu), λ is a parameter controlling the steepness of g.

(2)

In Section 2.2, where we review some of the work of [START_REF] Ashwin | Criteria for robustness of heteroclinic cycles in neural microcircuits[END_REF], [START_REF] Bick | Dynamical origin of the effective storage capacity in the brain's working memory[END_REF] and [START_REF] Fukai | A Simple Neural Network Exhibiting Selective Activation of Neuronal Ensembles: From Winner-Take-All to Winners-Share-All[END_REF], we make a brief switch to a different but equivalent choice of g used by these authors. System (1) transformed to the firing rate variables with g given by (2) has the form

ẋj = (1 -x 2 j ) (λ(Jx) j -f (x j )) , x = (x 1 , . . . , x n ) ∈ [-1, 1] n (3)
where

f (x) = g -1 (x) = arctanh(x) = 1 2 ln 1 + x 1 -x ( 4
)
and J is the coupling matrix. We further decompose J as follows:

J = c 0 I + c 1 J 1 ( 5
)
where I is the identity matrix, c 0 and c 1 are non negative coefficients and c 1 = 1 -c 0 . Provided that λc 0 > 1 the equation λc 0 β = arctanh(β) has a couple of non zero solutions ±β λ with 0 < β λ < 1. Therefore when c 1 = 0, any vector of the form β λ (ξ 1 , . . . , ξ n ) with ξ j = ±1 is a stable steady-state of (3). If we think of vectors of the form (ξ 1 , . . . , ξ n ) as information strings in a neural network, then the above steady-states represent stored memory states. However it is well-known that memory states need not be steady (see [START_REF] Gencic | Parallel processing in neural systems and computers[END_REF] and references therein). If c 1 > 0 the steady-states may become unstable or even disappear, but nevertheless information may still be dynamically stored. We now explain the idea of information storage by means of limit cycles of (3) (storing cycles), as explored in [START_REF] Zhang | Storing cycles in Hopfield-type networks with pseudo inverse learning rule: Admissibility and network topology[END_REF] and then we introduce our idea to use robust heteroclinic cycles instead. The basic question adressed in [START_REF] Zhang | Storing cycles in Hopfield-type networks with pseudo inverse learning rule: Admissibility and network topology[END_REF] is the following: given an information string, can it be stored by a Hopfield network in the form of dynamic information, more specifically a limit cycle? Concretely, the information is given in the form of a string of binary n-vectors (with components equal to ±1). The learning rule, consistent with Hebbian learning, is an algorithm specifying how the information string structures the coupling matrix c 0 I + c 1 J (we forget from now on the subscript 1 in J 1 ). This learning rule will be described in detail in Section 3.1. The main research question of [START_REF] Zhang | Storing cycles in Hopfield-type networks with pseudo inverse learning rule: Admissibility and network topology[END_REF] is whether the system with the coupling structure resulting from applying the learning rule supports stable limit cycles that code the original information string in the sense that the periodic orbit passes through the quadrants of R n corresponding to the elements of the information string, following its order.

In this article we focus on a different version of such encoding by the dynamics, choosing a robust heteroclinic cycle as the invariant object encoding the information string. The condition we impose is that the cycle should connect equilibrium points located at vertices of the cube [-1, 1] n corresponding to the elements of the information string, following its order. This is a rather natural condition, yet the first obstacle we must overcome is that with f as given by ( 4) the RHS of (3) is not C 1 on the cube [-1, 1] n , so that heteroclinic cycles cannot exist. We discuss this problem in more detail and propose a solution in the next section, which also relates to the work of [START_REF] Ashwin | Criteria for robustness of heteroclinic cycles in neural microcircuits[END_REF], [START_REF] Bick | Dynamical origin of the effective storage capacity in the brain's working memory[END_REF] and [START_REF] Fukai | A Simple Neural Network Exhibiting Selective Activation of Neuronal Ensembles: From Winner-Take-All to Winners-Share-All[END_REF].

The Lokta-Volterra approximation to Hopfield equations

The articles [START_REF] Ashwin | Criteria for robustness of heteroclinic cycles in neural microcircuits[END_REF], [START_REF] Bick | Dynamical origin of the effective storage capacity in the brain's working memory[END_REF] and [START_REF] Fukai | A Simple Neural Network Exhibiting Selective Activation of Neuronal Ensembles: From Winner-Take-All to Winners-Share-All[END_REF] consider the question of the existence of robust heteroclinic cycles in the firing rate version of (1) and show that such cycles exist for a Lotka-Volterra approximation of the system. In this section we use a different combination of g and f = g -1 consistent with choice made in these articles. Specifically we will use the functions:

g(u) = 1 1 + e -u and f (x) = ln x 1 -x . (6)
The coefficients J ij are assumed to be all positive so that the synaptic couplings are all of inhibitory type. We define the firing rate by x j = g(λ -1 u j ) and transform (1) to the firing rate formulation. After applying a time rescaling we obtain the following system.

ẋi = x i (1 -x i ) -λf (x i ) - n i=1 J ij x j + I , i = 1, . . . n. (7)
Note that system [START_REF] Ermentrout | Mathematical Foundations of Neuroscience[END_REF] is well defined and continuous on the cube [0, 1] n , but it is not smooth on the faces, with the term

x i (1 -x i )f (x i ) (8)
being the source of non-smootheness. As we are interested in heteroclinic cyles that lie on the edges, with connections in the faces, this becomes a problem for the existence and stability of the cycle. Since our purpose is merely to illustrate the problem of the lack of smoothness we restrict our attention to the simplest case n = 3. Then [START_REF] Ermentrout | Mathematical Foundations of Neuroscience[END_REF] has the form

ẋ1 = x 1 (1 -x 1 ) (-λf (x 1 ) -J 11 x 1 -J 12 x 2 -J 13 x 3 + I) ẋ2 = x 2 (1 -x 2 ) (-λf (x 2 ) -J 21 x 1 -J 22 x 2 -J 23 x 3 + I) ẋ3 = x 3 (1 -x 3 ) (-λf (x 3 ) -J 31 x 1 -J 32 x 2 -J 33 x 3 + I) . (9)
For simplicity we assume J jj = 1. The goal is to construct heteroclinic cycles connecting equilibria of the form:

(ρ, 0, 0), (0, ρ, 0) and (0, 0, ρ) where -λf (ρ) -ρ + I = 0.

The Jacobian matrix at such equilibria is given as follows:

  -ρ(1 -ρ) -λ -J 12 (1 -ρ)ρ -J 13 (1 -ρ)ρ 0 -λf (0) -J 21 ρ + I -λ 0 0 0 -λf (0) -J 31 ρ + I -λ   .
Since f blows up at 0 the Jacobian is undefined. If the term ( 8) is neglected in the RHS of each equation of ( 9) then a heteroclinic cycle can be easily found, with

J =   1 1.25 0 0.875 1 1.25 3 0.625 1   (10)
giving an example [START_REF] Ashwin | Criteria for robustness of heteroclinic cycles in neural microcircuits[END_REF].

Regularization

We propose two approaches to regularize the function f . First approach, which we use in this paper, is to replace f defined in (6) by its Taylor polynomial at x = 0.5. We denote such Taylor polynomial of degree q as f q . Note that the sequence {f q } diverges at 0 and 1 when q → ∞, but converges uniformly to f on any compact subinterval in (0, 1).

Another approach is to replace f by

f ε (x) = ln x + ε 1 + ε -x , ( 11
)
where ε is a small parameter. The function f ε is well defined on the interval [0, 1], yet its properties are similar to f , in particular its derivative at 0 and 1 equals 1/(ε(1 + ε), thus is very large. If we replace f by f q or f ε in (9) then, depending on the relative size of λ and q or ε, any of the three possibilities can arise:

1. the cycle does not exist, 2. the cycle exists and is unstable, 3. the cycle exists and is stable.

Fig. 1 shows simulations of ( 9) with f replaced by f ε . The matrix J is as given in [START_REF] Hofbauer | Evolutionary Games and Population Dynamics[END_REF], I = 0.8 and λ = 0.01. The value of ε is varied showing an example of each of the possible cases.

3 Hopfield networks with coupling given by the learning rule of [START_REF] Personnaz | Collective computational properties of neural networks: new learning mechanisms[END_REF] 3

.1 The equations and network architecture

We now return to the formulation (3) with f replaced by its q-th order polynomial expansion f q at x = 0 as described in Section 2.3. The equation now reads

ẋ = (I -diag x • x) (λ(c 0 x + c 1 Jx) -f q (x)) , x ∈ [-1, 1] n , (12)
with f q (x) = (f q (x 1 ), . . . , f q (x n )) and

f q (x) = x + x 3 3 + • • • + x q q .
The power series f ∞ has a radius of convergence equal to 1. It follows that given any interval (-1 + ǫ, 1 -ǫ), the approximation of f by f q can be as good as we wish provided that q is large enough. We now give a formal description of the information string and introduce the learning rule. A binary pattern (or simply a pattern) is a vector ξ of binary states of n neurons:

ξ = (x 1 , . . . , x n ) t with x j = ±1. Let Σ = ξ 1 , ξ 2 , . . . , ξ p (13)
be a sequence of p patterns. This matrix is called a cycle if there exists a connectivity matrix J such that the corresponding network of n neurons visits sequentially and cyclically the patterns defined by Σ. In other words each column ξ j can be associated with a state of the system such that the signs of the cell variables are equal to the signs of the corresponding components of ξ j . We shall always assume p ≥ n.

Let P be the matrix of the cyclic permutation

(x 1 , x 2 , . . . , x p-1 , x p ) → (x 2 , x 3 , . . . , x p , x 1 ). The cycle Σ is called admissible if there exists J such that JΣ = ΣP, (14)
has a solution [START_REF] Personnaz | Collective computational properties of neural networks: new learning mechanisms[END_REF]. This relation expresses a necessary condition for the network (3) to possess a solution that periodically takes the signs defined by the patterns ξ 1 , . . . , ξ p . Note that, if Σ is admissible then a solution exists in the form where Σ + is the Moore-Penrose pseudo-inverse of Σ, and if Σ has full rank it is unique. A cycle Σ is called simple [START_REF] Zhang | Storing cycles in Hopfield-type networks with pseudo inverse learning rule: Admissibility and network topology[END_REF] if there exists a vector η ∈ R p such that each row of Σ equals ηP l j , for some p > l j ≥ 0. We define

J = ΣP Σ + (15)
W η = span {ηP j : j = 0, . . . , p -1}. (16)
J =        0 1 0 . . . 0 0 0 0 1 . . . 0 0 . . . . . . . . . . . . . . . 0 0 0 . . . 0 a 0 a 1 a 1 . . . a n-2 a n-1        . ( 18
)
where a 0 , . . . , a n-1 are rational coefficients. If Σ has full rank then J is uniquely defined and a 0 = 0 (in this case the cycle is minimal in the sense of [START_REF] Zhang | Storing cycles in Hopfield-type networks with pseudo inverse learning rule: Admissibility and network topology[END_REF]).

Proof. By construction we can write

Σ =      η ηP . . . ηP n-1      and ΣP =      ηP ηP 2 . . . ηP n     
.

By admissibility ΣP = JΣ and moreover ηP n must be a linear combination of the ηP j 's. Hence (18) follows. The a j 's are rational because the vectors ηP j have integer coordinates. If Rank(Σ) = n then J is non singular, hence a 0 = 0.

Example 1. We consider Σ as follows, with p = 6:

Σ =   1 1 1 -1 -1 -1 1 1 -1 -1 -1 1 1 -1 -1 -1 1 1   , (19)
Let η = (1, 1, 1, -1, -1, -1). Note that the rows of Σ are η, ηP and ηP 2 . Note also that ηP 3 = -η. It follows that the rows of ΣP are ηP , ηP 2 and -η, i.e. the second, the third and the negative of the first row of Σ. Hence  

0 1 0 0 0 1 -1 0 0   Σ = ΣP. ( 20
)
Since the rows of Σ are independent the matrix ΣΣ T is invertible. Hence [START_REF] Rabinovich | Dynamical Principles in Neuroscience[END_REF] has a unique solution which, by (20), must be given by:

J =   0 1 0 0 0 1 -1 0 0   . ( 21
)
Note that Σ + = Σ T (ΣΣ T ) -1 and J satisfies [START_REF] Rabinovich | Information flow dynamics in the brain[END_REF]. This matrix provides a simple example of heteroclinic cycle, which we illustrate in Fig. 2: the reader can check on this numerical simulation that indeed trajectories follow the pattern defined by Σ. Observe that the trajectory closely follows the edges of the cube connecting the equilibria in the pattern. The existence of the cycle follows from Theorem 2 in Section 4.1. A direct proof is also easy to obtain.

Classification of simple consecutive cycles

Suppose that p is fixed and note that every consecutive cycle is uniquely determined by the choice of η and n. If n ≥ p then such a cycle is always admissible. If n < p then there are only very special choices of η and n such that the cycle is admissible. In this section we will address the question of finding the conditions on η so that there exists an n < p such that the cycle determined by η and n is admissible. In order to avoid confusion with prime numbers we will, throughout this section, use the letter m instead of p to denote the dimension of η. We will return to the original notation of [START_REF] Zhang | Storing cycles in Hopfield-type networks with pseudo inverse learning rule: Admissibility and network topology[END_REF] in the subsequent sections. Consider a simple cycle as defined in Section 3.1, with η corresponding to the first row of Σ. If Σ is admissible the there exists J in the form given by (18) such that ( 14) is satisfied. Let (a 0 , a 1 , . . . , a n-1 ) be the last row of J (see ( 18)) and let

ψ(x) = x n -a n-1 x n-1 -. . . -a 1 x -a 0 . ( 22
)
It follows from ( 14) that η ∈ ker(ψ(P )). In this section we use the following result:

Theorem 1. Let V be a non-trivial invariant subspace of the action of P on R m . Then there exists a polynomial φ(x), which is a divisor of x m -1, such that V = ker(φ(P )). Moreover, for any φ the inclusion V ⊂ ker( φ(P )) holds if and only if φ is a divisor of φ.

Theorem 1 follows from some classical results of algebra, which we will review in the appendix, thereby providing the proof. We now state two corollaries of Theorem 1 which we will use to characterize the possible choices of η for which dim(W η ) < m.

Corollary 1. If n = dim (W η ) < m then W η = ker(ψ(P )) and ψ is a divisor of x m -1.

Proof. It is easy to see that W η ⊂ ker(ψ(P )). We will prove that the opposite inclusion holds and that ψ is a divisor of x m -1. By Theorem 1 there exists φ a divisor of x m -1 such that ker(φ(P )) = W η and φ divides ψ. Suppose that φ is a proper divisor of ψ. Then

n = dim (W η ) ≤ deg(φ) < n,
which is a contradiction. It follows that ψ = φ. Hence the corollary holds.

For φ a minimal divisor of x m -1 (over Q) let ψ = (x m -1)/φ and let W φ = ker(ψ). The following result leads to a characterization of ηs such that n = dim(W η ) < m.

Corollary 2. If n = dim(W η ) < m then η ∈ W φ , for some φ a minimal divisor of x m -1 (over Q).
Proof By Corollary 1 there exists ψ, a divisor of x m -1, such that W η = ker( ψ). By unique decomposition into prime factors over Q there exists a minimal divisor φ of x m -1 which divides ψ. Let ψ = (x m -1)/φ. Clearly ψ divides ψ. It follows from Theorem 1 that η ∈ ker(ψ) = W φ .

In the remainder of this section we will derive the conditions on η needed for η ∈ W φ for some (the simplest) choices of φ, where φ is a minimal divisor of x n -1. We begin by recalling the decomposition of x n -1 into irreducible polynomials over Q. For a positive integer k let

φ k (x) = k-1 i=0 x i .
The polynomials φ p , where p is a prime number, are irreducible over Q. For a prime number p and a non-negative integer j we define

φ p,j (x) = p-1 i=0 x ip j .
Note that φ p,1 = φ p . The polynomials φ p,j (x) are the irreducible factors over Q of the polynomial φ p k .

Suppose m = m 1 m 2 . . . m l , with m j = p k j j , j = 1, . . . , l, and let

Φ m (x) = φ m l j=1 φ m j (x).
The polynomial Φ m (x) is called the cyclotomic polynomial of degree m and is irreducible. It now follows that the decomposition of x m -1 into irreducible factors over Q given by

x m -1 = (x -1)Φ m (x) l j=1   k j -1 i=0 φ p j ,i (x)   . ( 23
)
All the possible factors of x m -1 over Q are products of the irreducible factors appearing in (23), hence all the possible choices of W φ are obtained that way. As announced above we now describe some of the spaces W φ by simple conditions on the components of η.

Proposition 2. If φ = φ p j , for some prime number p j then W φ consists of vectors η = (b 0 , . . . , b m-1 ) satisfying

m/p l -1 i=0 b ip j = m/p l -1 i=0 b ip j +1 = . . . m/p l -1 i=0 b ip j +p j -1 . ( 24
)
Proof we use the following identity:

ψ(x) = (x -1)φ m/p j (x p j ).
Hence

Ψ(P )η = (P -I)   m/p l -1 i=0 b ip l , m/p l -1 i=0 b ip l +1 , . . . , m/p l -1 i=0 b ip l +m-1   . (25)
(The indexing of the components of b in (25) must be understood modulo m.) It follows that the RHS of ( 25) is equal to the 0 vector if (24) holds.

We now state the condition on η for φ = Φ m (x). We begin with the following elementary lemma (the proof is left to the reader).

Lemma 1. If k divides m then ker(P k -I) = {η : there exists v ∈ {-1, 1} k such that η = (v, v, . . . , v)}. ( 26
) Proposition 3. Suppose φ = Φ m (x). Then η ∈ ker(P m 1 -I) + . . . + ker(P m l -I). ( 27
)
Proof Note that

ψ(x) = (x m -1) Φ m (x) = (x -1) l j=1 φ m j (x).
Further note that, for each j, (x -1)φ m j (x) = x m j -1. Hence, for each j, ker(P m j -I) ⊂ ker(ψ(P )).

Moreover, for j = j ′ ker(P m j -I) ∩ ker(P m j ′ -I) = span{1 m }.

The result follows.

Remark 1. Since the coordinates of η are ±1, it follows that η must be contained in one of the spaces ker(P m l -I).

Remark 2. The conditions for the other minimal factors of x m -1 are slightly more complicated and we will not state them here. They are, however, not hard to derive.

Example 2. We consider η = (b 0 , b 1 , . . . , b 14 ) = (1, 1, 1, 1, 1, 1, 1 -1, -1, -1, 1, 1, 1, 1, 1) with m = 15 = 3 • 5.
Note that η satisfies (24) with p = 3. Hence, by Proposition 2, η ∈ W φ = ker(ψ) with φ = 1+x+x 2 ; ψ = (x-1)(1+x 3 +x 6 +x 9 +x 12 ) = -1+x-x 3 +x 4 -x 6 +x 7 -x 9 +x 10 -x 12 +x 13 Hence ηP 13 = η -ηP + ηP 3 -ηP 4 + ηP 6 -ηP 7 + ηP 9 -ηP 10 + ηP 12 .

(28)

Note that the last row of ΣP is equal to the LHS of (28). If we define J as in (18) with (a 0 , . . . , a 12 ) = (1, -1, 0, 1, -1, 0, 1, -1, 0, 1, -1, 0, 1) then the RHS of (28) equals JΣ. Hence, by (28), identity [START_REF] Rabinovich | Dynamical Principles in Neuroscience[END_REF] holds with J as specified. Arguing as in Example 3 we conclude that Σ generated by η and n = 5 is admissible with J whose last row equals (1, 0, 0, 0, 0).

Example 4. An interesting class of admissible cycles exists for m even with n = m/2. Note that in this case x m -1 = (x n + 1)(x n -1), i.e. x n + 1 divides x m -1. Further note that

ker (P n + I) = {η ∈ R m = (ν, -ν), ν ∈ R n }. ( 30
)
Let Σ be a cycle constructed with some η ∈ ker (P n + I), n = m/2. Then Σ is admissible. Moreover, by a similar argument as in Example 2 we conclude that the last row of J equals (-1, 0, . . . , 0). In particular Example 1 of Section 3.1 is a special case of this construction. This type of admissible cycle is called antisymmetric in [START_REF] Zhang | Storing cycles in Hopfield-type networks with pseudo inverse learning rule: Admissibility and network topology[END_REF].

Example 5. Since (1 + x + . . . + x m-1 )(x -1) = x m -1 the space ker(I + P . . . + P m-1 ) corresponds to vectors η for which Σ with n = m -1 is admissible. Moreover ker(I + P . . . + P m-1 ) = {η :

η i = 0}.
Note that for η's whose entries are ±1 this means that the number of coordinates equal to 1 is the same as the number of coordinates equal to -1. Hence m must be even. In this case the last row of J is (-1, -1, . . . , -1).

Consecutive Hopfield cycles and their heteroclinic cycles

We now come to the study of heteroclinic cycles for admissible consecutive simple cycles governed by equation ( 3), hence with J as in (18). Then the equation reads as a system

ẋ1 = (1 -x 2 1 ) (λc 0 x 1 + λc 1 x 2 -f q (x 1 )) ẋ2 = (1 -x 2 2 ) (λc 0 x 2 + λc 1 x 3 -f q (x 2 )) . . . ẋn = (1 -x 2 n ) (λc 0 x n + λc 1 (a 1 x 1 + • • • + a n x n ) -f q (x n )) (31)
Following [START_REF] Zhang | Storing cycles in Hopfield-type networks with pseudo inverse learning rule: Admissibility and network topology[END_REF], we also assume that the two coefficients which control the relative contributions of J 0 and J to each neuron satisfy

(H) 0 ≤ c 0 < 1 and c 0 + c 1 = 1.
We aim at studying the existence and stability of heteroclinic cycles connecting vertex equilibria, i.e. equilibria with entries ±1, for this system. By construction, the edges, faces and simplices of the hypercube {(±1, . . . , ±1)} are invariant under the dynamics of 31. Let ξ = (x 1 , . . . , x n ) be a vertex equilibrium: x k = ±1 for all k. Linearizing (31) at ξ leads to a system of equations uk = σ k u k where: we can express the eigenvalues as follows:

σ k = 2(f q (1) -λ) if x k x k+1 = 1, k < n σ k = 2(f q (1) -λ(c 0 -c 1 )) if x k x k+1 = -1, k < n σ n = 2 f q (1) -λ(c 0 + c 1 x n n j=1 a j x j ) , (32)
Note that under the above conditions on c 0 and c 1 , which we assume from now on, a necessary and sufficient condition for the existence of negative and positive eigenvalues with k < n is that

λ(c 0 -c 1 ) < f q (1) < λ (33)
This is always possible to realize since |c 0 -c 1 | < 1. Then σ k < 0 if x k x k+1 = 1 and > 0 otherwise.

Remark 3. The equation (12) (hence (31)) is invariant by the symmetry S : x → -x. Therefore any time a cycle Σ admits a heteroclinic cycle, the cycle -Σ admits the opposite heteroclinic cycle obtained by applying S.

In the following we shall always consider Σ's up to this symmetry.

Heteroclinic edge cycles Definition 1. A heteroclinic cycle is called an "edge cycle" if it connects a cyclic sequence of vertex

equilibria through heteroclinic orbits lying on the edges of the hypercube {(±1, . . . , ±1)}. We also request that the unstable manifold at each equilibrium in the cycle has dimension 1 (therefore is contained in an edge).

The condition about the unstable manifolds is necessary for asymptotic stability of the edge cycles. If σ - k and σ + k denote respectively the contracting and expanding eigenvalues along the heteroclinic trajectories, the edge cycle is asymptotically stable if (see [START_REF] Krupa | Robust heteroclinic cycles[END_REF])

|Πσ - k | > Πσ + k (34)
The example 1 provides a simple case of an asymptotically stable edge cycle, see Fig. 2. We show below that all asymptotically stable edge cycles have the same simple structure.

Theorem 2. Let hypothesis (H) hold. An edge cycle exists for (31) if and only if condition (33) holds as well as the following:

λ(c 0 + c 1 (a 1 + • • • + a n ) < f q (1) < λ(c 0 + c 1 (-a 1 + • • • -a n-1 + a n ).
This edge cycle connects the sequence of 2n equilibria

(1, 1, . . . , 1) → (1, 1, . . . , -1) → . . . (-1, -1, . . . , -1) → • • • → (-1, 1, . . . , 1). ( 35
)
Proof. Let ξ = (x 1 , . . . , x n ) be an equilibrium in the cycle. Note first that according to (32), in order to have one unique positive eigenvalue σ k with k < n, the following must be true: (33) holds and (i) all x j with j ≤ k have the same sign, (ii) x k x k+1 = -1 and (iii) x j has the sign of x k+1 for k + 1 < j. Let ξ ′ be the next equilibrium in the cycle, then we must have x ′ j = x j for all j = k and x ′ k = -x k . Observe that we then have σ ′ k < 0. It is straightforward to check that under (33), there is no equilibrium point lying on the edge joining ξ to ξ ′ and therefore that a heteroclinic connection exists on this edge. Now let's assume that the positive eigenvalue is σ n . Then all x j 's, j = 1, . . . n, must be equal and the condition σ n > 0 can be written

f q (1) -λ(c 0 + c 1 (a 1 + • • • + a n ) > 0. Also we request x ′ n = -x n and σ ′ n < 0, which can be written f q (1) -λ(c 0 + c 1 (-a 1 + • • • -a n-1 + a n ) < 0.
As for the case k < n we can check that if these inequalities are satisfied a heteroclinic orbit joins ξ to ξ ′ .

From the above construction we deduce that the edge cycle must connect the equilibria in the sequence (35).

Corollary 3. Edge cycles are in one-to-one correspondance with connectivity matrices (18) with a 1 = -1 and a j = 0 for j > 1. Moreover, under hypothesis (H), they are asymptotically stable iff λ < c 0 /f q (1).

Proof. It follows from the above theorem that the matrix Σ for an edge cycle is defined by

Σ = (η, ηP, . . . , ηP n-1 ) T
where η is the vector (v, -v) with v = (1, . . . , 1) (n times). Note that ηP n = -η. It follows that the rows of ΣP are ηP, . . . , ηP n-1 and -η. Hence J with a 1 = -1 and a j = 0 is solution of

JΣ = ΣP.
Since the rows of Σ are independant the matrix ΣΣ T is invertible, hence the solution is unique.

It is straightforward to check that (34) is true in this case iff λ < c 0 /f q (1).

Note that Example 1 provides the simplest case of an edge cycle.

Heteroclinic non-edge cycles

We have seen in the previous section that in order for a vertex equilibrium ξ = (x 1 , . . . , x n ) of (31) to have a unique positive eigenvalue, a necessary condition was that a change of sign in the sequence of coordinates x j occurs at most once. The sign of σ n is a special case, it depends on the coefficients a 1 , . . . , a n . Suppose now that ξ has two positive eigenvalues, along directions x k and x l . The corresponding two dimensional unstable manifold lies in the face defined by the fixed coordinates x j when j = k, l. Assuming k < l < n, the four vertices on this face are ξ, ξ ′ = (x 1 , . . . , -x k , . . . , x n ), ξ ′′ = (x 1 , . . . , -x l , . . . , x n ) and ξ = (x 1 , . . . , -x k , . . . , -x l , . . . , x n ). The question which we address now is whether there can exist stable heteroclinic cycles which involve saddle-sink connecting trajectories from ξ to ξ. This situation can of course be generalized to more than two unstable eigenvalues, if there are more than two switches of sign in the x j 's. Let us first look at an example in low dimension. In all the following we assume hypothesis (H) holds.

Example 6. Consider 3 neurons (n = 3) and 4 equilibria such that η = (1, 1, -1, -1). Defining P as before and η = (1, 1, -1, -1), we build Σ to form a consecutive cycle with n = 3 and p = 4: {η, ηP, ηP 2 }. Hence

Σ =   1 1 -1 -1 1 -1 -1 1 -1 -1 1 1   . ( 36
)
Clearly the third row is the opposite of the first one, hence this matrix has rank 2. Nevertheless the cycle is admissible because rank(Σ) = rank(η) where rank(η) is the rank of the matrix [η T , (ηP ) T , (ηP 2 ) T , (ηP 3 ) T ] T (Theorem 2 of [START_REF] Zhang | Storing cycles in Hopfield-type networks with pseudo inverse learning rule: Admissibility and network topology[END_REF]). Since ηP 3 = -η -ηP -ηP 2 , it follows that

J =   0 1 0 0 0 1 -1 -1 -1   . ( 37
)
Note that this example illustrates the criterion derived in Section 3.2, Example 5. The equations read

ẋ1 = (1 -x 2 1 ) (λc 0 x 1 + λc 1 x 2 -f q (x 1 )) ẋ2 = (1 -x 2 2 ) (λc 0 x 2 + λc 1 x 3 -f q (x 2 )) ẋ3 = (1 -x 2 3 ) (λc 0 x 3 -λc 1 (x 1 + x 2 + x 3 ) -f q (x 1 )) (38)
Numerical simulations exhibit a heteroclinic cycle for (3) with J given above, see Fig. 3. Observe that after short transient time x 1 and x 3 are opposite and move (in opposite directions) while x 2 is fixed at ±1. This indicates that a heteroclinic orbit (if it exists) connects opposite vertices in the faces x 2 = ±1. Now if we set x 3 = -x 1 in (39) with x 2 = ±1, we see that the first and third equations are identical. Therefore the diagonal axis joining the vertices (1, ±1, -1) to (-1, ±1, 1) is flow-invariant.

Moreover the eigenvalues at opposite vertices along each of these axes have opposite signs as in the previous sections, showing that a saddle-sink connection exists on these diagonal axes. A more detailed calculation shows that on each of these faces, the dynamics looks like in Fig. 4. Let us explain why the dynamics aligns itself on the diagonal connection (in black in Fig. 4). The unstable eigenvalue at the point (-1, -1, -1) is given by σ n (n = 3) in (32): σ 3 = 2λ(-c 0 + 3c 1 ) + 2f q (1).

We now come back to the general problem.

Lemma 2. Let ξ = (x 1 , . . . , x k , . . . , x l , . . . , x n ) be a vertex equilibrium for (31). If there are m switches of sign in the sequence x 1 , . . . , x n , then ξ has either m or m + 1 unstable directions. The latter case occurs if σ n > 0 in (32). The unstable manifold of ξ is contained in the hyperface generated by its unstable eigenvectors. Proof. One already knows that if j < n, then σ j > 0 iff x j x j+1 < 0. However when j = n, the sign of the eigenvalue depends upon the coefficients a 1 , . . . , a n . The last claim is straightforward from (31).

Lemma 3. An equilibrium ξ possesses a 2-dimensional unstable manifold if and only if the two conditions are satisfied: (i) either the sequence of coordinates in ξ undergoes two switches of signs and σ n < 0, or one switch of signs and σ n > 0; (ii) if x k and x l are the unstable directions, then x k x l < 0.

Proof. Condition (i) is clear from Lemma 2. If condition (ii) is not satisfied, then an additional change of sign must occur somewhere between x k and x l and therefore an additional positive eigenvalue must exist.

The next lemma characterizes when σ n > 0 when ξ is a column vector of a consecutive cycle Σ.

Lemma 4. Let ξ = (x 1 , . . . , x n ) be an equilibrium in a consecutive cycle. We write ξ ′ = Jξ = (x ′ 1 , . . . , x ′ n ) T where J is the connectivity matrix (18). Then under the condition (33), σ n > 0 if and only if x n x ′ n < 0.

Proof. Recall that σ n = -2 λ(c 0 + c 1 n j=1 a j x j ) -f q (1) . By construction n j=1 a j x j = x ′ n , therefore by (32),

σ n = -2 (λ(c 0 + c 1 x n x ′ n ) -f q (1)
). The claim follows by (33).

Lemma 5. Let ξ = (x 1 , . . . , x n ) be a vertex equilibrium in a simple consecutive cycle with connectivity matrix J and ξ ′ = Jξ = (x ′ 1 , . . . , x ′ n ). Suppose that ξ possesses a 2-dimensional unstable manifold along directions x k and x l , which by construction implies x ′ k = -x k and x ′ l = -x l . Then a heteroclinic

u = (1 -u 2 ) λc 0 u + λc 1 x ′ k -f q (u) v = (1 -v 2 ) λc 0 v + λc 1 x ′ l -f q (v) (39)
Set v = -u. Then the two above equations are identical because by Lemma 3 we have x k x l = -1, which also implies x ′ k x ′ l = -1. The saddle-sink connection along this segment follows from the same analysis as in the "edge" case. Suppose now that l = k + 1. Then in ξ ′ the coordinates of indices k, k + 1 have opposite signs, which implies that the eigenvalue σ ′ k of ξ ′ is positive. Therefore ξ ′ is a saddle or a source in the face joining ξ to ξ ′ , which proves that no saddle-sink connection ξ → ξ ′ can exist. This lemma can be generalized in a straightforward way to more that two unstable eigenvalues. Definition 2. Let Σ = (ξ 1 , . . . , ξ p ) be a simple, admissible consecutive cycle. Σ has adjacent switches if in one column (at least), the sign of the entries change two or more times consecutively.

The following theorem summarizes the previous results. Theorem 3. Let hypothesis (H) hold. For the admissible simple consecutive cycle Σ = (ξ 1 , . . . , ξ p ) with Hopfield equations (31), the equilibria ξ 1 , . . . , ξ p are connected by a robust heteroclinic cycle if and only if: (i) condition (33) is satisfied; (ii) Σ has no adjacent switches. The heteroclinic connections ξ i → ξ i+1 lie either along the corresponding edge of the unit cube in phase space, or inside the corresponding face, the dimension of which is equal to the number q of switches of sign of coordinates from ξ i to ξ i+1 . In the latter case these connections form a q-dimensional manifold, and in this manifold one of them is the diagonal segment joining ξ i to ξ i+1 .

The example 1 above illustrates this theorem for a network of three neurons. The first column in (36) contains one switch of sign but the second column contains 2 non adjacent switches. The resulting dynamics close to the heteroclinic cycle is shown in Fig. 3. The next example also concerns a network with three neurons, however it is a counter-example to existence of a heteroclinic cycle.

Example 7. Let's take Σ in the following form:

Σ =   -1 1 1 1 1 -1 1 -1 1   . (40)
This matrix defines a simple minimal consecutive cycle: it is 3 × 3 and invertible. Since P 3 = Id, it is easy to find that

J =   0 1 0 0 0 1 1 0 0   . ( 41
)
Observe that in Σ the third row has two adjacent switches of sign. The numerical simulation shows a dynamics converging to the equilibrium (1, 1, 1) (Fig. 5).

We can't rule out the possibility that condition (ii) in Theorem 3 is not satisfied, but the network still possesses a heteroclinic cycle. However in this case the cycle will be different from the one defined by Σ. Two different Σ's can give the same connectivity matrix. Next is an example of this kind. By Example 5 of Section 3.2 the connectivity matrix is

Example 8. This example shows that the assumption of Theorem 3 concerning the absence of adjacent switches is essential. Lets consider the following minimal consecutive cycle, which was introduced in [4] (Example 6).

Σ =       1 1 -1 1 -1 -1 1 -1 1 -1 -1 1 -1 1 -1 -1 1 1 1 -1 -1 1 1 -1 -1 -1 1 1 -1 1       . (42)
J =       0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 -1 -1 -1 -1 -1       . ( 43
)
Observe that Σ possesses two adjacent sign switches in the first column. The simulation of the dynamics in this case shows a heteroclinic cycle, however not the one which would correspond to the cycle formed by the columns of Σ (Fig. [START_REF] Dummit | Abstract Algebra 3rd Edition[END_REF]), which is consistent with Theorem 3. In the figure we see that x 3 moves first from +1 to -1 while x 1 = x 2 = +1 and x 4 = x 5 = -1, then x 1 moves from +1 to -1 and simultaneously x 5 moves from -1 to +1, then x 2 and x 4 do the same, and the process repeats itself. The corresponding cycle is given by the following matrix (see also Section 4.3.2):

Σ =       1 1 1 -1 -1 -1 1 1 -1 -1 -1 1 1 -1 -1 -1 1 1 -1 -1 -1 1 1 1 -1 -1 1 1 1 -1       . ( 44
)
In the following example the rank of Σ is not maximal. Example 9. Let's take

Σ =     1 1 -1 -1 -1 1 1 -1 -1 -1 1 1 -1 -1 -1 1 1 1 -1 -1 1 1 1 -1     . (45)
Observe that the last row is opposite to the first one, which we call η, and rank(Σ) = 3. The cycle is admissible and since ηP 4 = -ηP , we have that

J =     0 1 0 0 0 0 1 0 0 0 0 1 0 -1 0 0     . (46)
A numerical simulation is shown in Fig. 7.

Some classes of simple consecutive cycles with non-edge heteroclinic cycles 4.3.1 Simple consecutive cycles with p = n

Therefore Σ is a square matrix. By construction, if η denotes the first row of Σ, then Σ = η T , (ηP ) T , . . . , (ηP n-1 ) T . Then it follows from Theorem 2 in [START_REF] Zhang | Storing cycles in Hopfield-type networks with pseudo inverse learning rule: Admissibility and network topology[END_REF] that the cycle is admissible. Moreover ηP n = η, which implies

Observe that for n = 3 adjacent switches always exist in this case (see Example 2). However for all square simple consecutive cycles of a given dimension n > 3 and such that no adjacent switches occurs, we can conclude that several non-edge heteroclinic cycles can coexist and their number increases with n.

Example 10. n = 4. The only square consecutive cycle with non adjacent switches is generated by η = (1, 1, -1, -1). Time series of the heteroclinic cycle shown in Fig. 8.

Example 11. n = 6. Then the following square consecutive cycles have no adjacent switches:

η 1 = (1, 1, 1, -1, -1, -1) and η 2 = (1, 1, 1, 1, -1, - 1
). The first heteroclinic cycle has connections on three different faces while the second cycle has connections on six different faces. Example 2 shows the construction of J for η satisfying (24) with p = 15 and n = 13. We leave the obvious generalization of this construction to the reader. Figure 9 shows a heteroclinic cycle obtained for η = (1, 1, 1, 1, 1, 1, 1, -1, -1, -1, 1, 1, 1, 1, 1), p = 15 and n = 13. For η ∈ ker(P n -I), as stipulated by Proposition 3 and Remark 1 the heteroclinic cycles which arise are the same as for η truncated to a single block. Multiple blocks simply correspond to repeated passages through the same heteroclinic cycle. For example, if we use, as in Example 2 of Section 3.2, η = (1, 1, 1, -1, -1, 1, 1, 1, -1, -1, 1, 1, 1, -1, -1), with n = 5, we obtain a heteroclinic cycle in R 5 corresponding to η = (1, 1, 1, -1, -1). A triple passage through this cycle gives η as stated above.

Conclusion

In this work we have studied robust heteroclinic cycles in Hopfield networks with coupling given by by the learning rule of [START_REF] Personnaz | Collective computational properties of neural networks: new learning mechanisms[END_REF]. We gave an extensive classification of heteroclinic cycles couplings of a simple consecutive type. In particular we established a tight relation between the structure of the coupling and the heteroclinic cycles supported by the resulting network. This work is a part of the general program of establishing connections between heteroclinic/homoclinic dynamics and neural processing (see [START_REF] Rabinovich | Dynamical Principles in Neuroscience[END_REF] for an outline of this program).

An interesting direction for continuing this work is to determine if a correspondence between cycles in the coupling and robust heteroclinic cycles carries over to more realistic settings. As a first attempt of such a generalization we intend to introduce delays, as delays arise naturally in neural coupling and can play a functional role [START_REF] Ermentrout | Mathematical Foundations of Neuroscience[END_REF]. Other generalizations include considering systems with noise, more realistic models of neurons, or generalizations of the learning rule. Theorem 5. If T is a cyclic linear operator on R n and v ∈ R n is a cyclic generator, then the invariant subspaces of R n for T are in one-to-one correspondence with the pairs (f, g) such that f g = f T . The space corresponding to such a pair (f, g) is [g(T )(v)] = Im(g(T )) = ker(f (T )).

In particular, the minimal invariant subspaces correspond to the pairs (f, g) where f is an irreducible factor of f T in R[x] and the maximal invariant subspaces correspond to the pairs (f, g) where g is an irreducible factor of f T in R[x] .

Note that Theorem 5 implies Theorem 4.

Figure 1 :

 1 Figure1: Simulations with J is as given in[START_REF] Hofbauer | Evolutionary Games and Population Dynamics[END_REF], I = 0.8 and λ = 0.01 illustrate the three possible cases for the regularized system corresponding to three different values of ε. Panel A shows a heteroclinic cycle, with ε = 0.12, panel B shows a periodic orbit close to an unstable cycle, with ε = 0.07, panel C shows the dynamics attracted to an equilibrium for the case when no cycle exists, with ε = 0.02.

Figure 2 :

 2 Figure 2: A trajectory of (3) with J given by (21) and with c 0 = 0.6, λ = 8. Initial conditions close to (1, 1, 1).

Example 3 .

 3 We consider η = (b 0 , b 1 , . . . , b 14 ) = (1, 1, 1, -1, -1, 1, 1, 1, -1, -1, 1, 1, 1,-1, -1) with m = 15 = 3 • 5. Note that η ∈ ker(P 5 -I), or, in other words, ηP 5 = η.(29)

Figure 3 :

 3 Figure 3: Simulations of (3) with c 0 = 0.6, λ = 2.5, random initial conditions close to origin. Colour code: blue= x 1 , red= x 2 , green= x3.

Figure 4 :

 4 Figure 4: Sketch of phase portrait on the face x 2 = -1 for (39). The cross indicates an equilibrium lying on the edge off vertices.

Figure 5 :

 5 Figure 5: Time series for the network with connectivity matrix (41), c 0 = 0.6, λ = 2, initial condition close to the first equilibrium (-1, 1, 1). Colour code: blue= x 1 , red= x 2 , green= x 3 .

Figure 6 :

 6 Figure 6: Time series (after transients) for the network with connectivity matrix (43), c 0 = 0.6, λ = 2.5, random initial conditions. Colour code: blue= x 1 , red= x 2 , green= x 3 , black= x 4 , yellow= x 5 .

Figure 7 :

 7 Figure 7: Time series (after transients) for the network with connectivity matrix (46), c 0 = 0.6, λ = 2.5, random initial conditions. Colour code: blue= x 1 , red= x 2 , green= x 3 , black= x 4 .

Figure 8 :

 8 Figure 8: Time series (after transients) for the consecutive cycle with p = n = 4, c 0 = 0.6, λ = 3.4, random initial conditions. Colour code: blue= x 1 , red= x 2 , green= x 3 , black= x 4 .

4. 3 . 3

 33 Simple consecutive cycles with n < p odd given by Propositions 2 and 3

Figure 9 :

 9 Figure 9: Time series (after transients) for the consecutive cycle with p = 15, n = 13, c 0 = 0.6, λ = 3.4, random initial conditions. Top picture = time series of x 1 , . . . , x 6 , bottom picture = time series of x 7 , . . . , x 13 . Color code: blue= x 1 , x 7 , red= x 2 , x 8 , green= x 3 , x 9 , black= x 4 , x 10 , yellow= x 5 , x 11 , cyan= x 6 , x 12 , magenta= x 13 .

  By Theorem 2 in[START_REF] Zhang | Storing cycles in Hopfield-type networks with pseudo inverse learning rule: Admissibility and network topology[END_REF] a simple cycle is admissible if and only if dim W η = Rank(Σ). If in addition we can write Σ = η t , (ηP ) t , . . . , (ηP n-1 ) t t (17) then the simple cycle is called consecutive. The following proposition is essentially contained in Sec.

5.1.1 of [4]. Proposition 1. If a simple consecutive cycle is admissible, then there exists a J satisfying (14) of the form

Acknowledgement

This work was partially supported by the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement no. 269921 (BrainScaleS), no. 318723 (Mathemacs), and by the ERC advanced grant NerVi no. 227747.

A Appendix on invariant spaces for linear actions

The material presented here is standard in algebra and is related to the rational or Frobenius normal form for matrices, see for example [START_REF] Dummit | Abstract Algebra 3rd Edition[END_REF]. We have not found a reference that presents the necessary results in a concise fashion, thus we feel that there is a need for this appendix.

Fix a linear transformation T : R n → R n . We are interested in understanding the (T -)invariant subspaces of R n , that is, those subspaces W with T (W ) ⊆ W , and in particular in identifying the maximal proper invariant subspaces. Here we will argue that invariant subspaces are very closely linked to the action of the polynomial ring R[x] on the linear transformations on R n given by (f, T ) → f (T ). Rather than studying specifically the permutation P we consider a more general context of cyclic transformations. A transformation T is cyclic if there exists a vector v ∈ R n such that R n = span{T j v, j = 0, 1, . . .}.

We will prove the following result.

Theorem 4. Let T R n → R n be a cyclic transformation. There exists an polynomial f T of degree n (the minimal polynomial) with the following property. The invariant spaces for the action of T on R n are in one to one correspondence with non-trivial factors of the polynomial f T . If f is such a factor then ker(f (T )) is the corresponding invariant space.

It will be clear from the arguments below that f P = x n -1. Hence Theorem 1 is a direct consequence of Theorem 4.

Generalities on the correspondence between invariant subspaces and polynomials

For each invariant subspace W ⊆ R n , we get a map

We denote the monic generator of this ideal by f W . Note that W = ker(f W (T )). Note also that for W = R n , the polynomial f R n is the minimal polynomial of T , and we will denote it by f T . A few observations in this setting will be useful.

Proof. This is a consequence of the fact that, for any h ∈ R[x], the linear transformations h(T ) and f (T ) commute so that:

Proof. To see this, first we can write

For the statement about images we have

and thus f (T ) and f ′ (T ) have equal images. The proof for kernels is similar.

Proof. Suppose v ∈ ker(f (T )) ∩ ker(g(T ). We write 1 = hf + kg and thus

We refer to [v] T as the cyclic subspace generated by v. Clearly, every minimal subspace must be of this form. We will prove in Lemma 10 that every invariant subspace has this form. Now let v ∈ R n and, again from the action of the polynomial ring, we obtain a map

Again, this map is linear, and its kernel is a principal ideal of R[x], the monic generator of which we will denote by f v . The image of this map is

Proof. The first statement follows as if f (T ) annihilates v then it also annihilates g(T )(v) for any

The second statement follows since, by the first isomorphism theorem, we have

In the proof of the next lemma, we need the fact that R[x] is a Unique Factorization Domain (UFD), which means that each non-zero polynomial f ∈ R[x] may be written as f = af n 1 1 . . . f n k k where a is a real number and each f i is an irreducible divisor of f which is also prime (that is, for all

Lemma 10. Let W be an invariant subspace. Then there is w ∈ W with f w = f W .

Proof. To see this, first note that f W = lcm(f w | w ∈ W ). (lcm denotes the least common multiple). Thus, for any irreducible divisor f of f W and for k the largest power of f that divides f W , there must be a w ∈ W so that f k divides f w . Now taking g = f w /f k , we see by Lemma 9 that f w ′ = f k where w ′ = g(T )(w). Doing this for each irreducible divisor of f W , the sum of the resulting w ′ s is the required element by Corollary 4, since f w ′ = f k implies w ′ ∈ ker(f k (T )).

Cyclic transformations

Recall that a linear transformation T

Proof. This follows from Lemmas 8 and 9: Take v ∈ R n with [v] T = R n . Then f v = f T is the minimal polynomial and thus f = f w for w = g(T )(v) and dim([w] T ) = deg(f ). Finally note that [w] T = span(T i (g(T )(v)) = g(T )(T i (v)) | i ∈ N) = g(T )(R n ) = ℑ(g(T )).

Lemma 12. If T is cyclic and f g = f T , then ℑ(g(T ))) = ker(f (T )).

Proof. For any w with w = g(T )(w ′ ) we have f (T )(w) = f (T )(g(T )(w ′ )) = 0, so ℑ(g(T )) ⊆ ker(f (T )) and we have n = deg(f ) + deg(g) = dim(ℑ(g(T ))) + dim(ℑ(f (T )))

≤ dim(ker(f (T ))) + dim(ℑ(f (T ))) = n.

Thus ℑ(g(T ))) = ker(f (T )).

The following result now follows.