Stability of the splay state in networks of pulse-coupled neurons.

We analytically investigate the stability of splay states in the networks of N globally pulse-coupled phase-like models of neurons. We develop a perturbative technique which allows determining the Floquet exponents for a generic velocity field and implement the method for a given pulse shape. We find that in the case of discontinuous velocity fields, the Floquet spectrum scales as 1/N2 and the stability is determined by the sign of the jump at the discontinuity. Altogether, the form of the spectrum depends on the pulse shape, but it is independent of the velocity field.PACS: 05.45.Xt, 84.35.+i, 87.19.lj.

1 Introduction

The first objective of (neural) network theory is the identification of asymptotic regimes. Previous research activity led to the discovery of fully- and partially-synchronised states, clusters and splay or asynchronous states in pulse-coupled networks [1-4]. It has also been made clear that ingredients such as disorder (the diversity of neurons and the structure of connections) are very important in determining the asymptotic behaviour, as well as the possible presence of delayed interactions and plasticity [5,6]. However, even if one restricts the analysis to identical, globally-coupled oscillators, there are very few theoretical results and they mostly concern fully-synchronised regime or specific types of neurons (e.g. the leaky integrate-and-fire model) [4,7,8].

In this paper, we develop a perturbative analysis for the stability of splay states (also known as antiphase states [9], ‘ponies on a merry-go-round’ [10], or rotating waves [11]) in ensembles of N globally pulse-coupled identical neurons. In a splay state, all the neurons follow the same periodic dynamics and their phases are evenly shifted. Accordingly, the phase, and potential, separation is of order . Splay states have been identified in experimental measurements performed on electronic circuits [11] and on multimode lasers [12]. Theoretical studies have been devoted to splay states in fully-coupled Ginzburg-Landau equations [13], Josephson arrays [14,15], laser models [16], traffic models [17], unidirectionally delay-coupled Stuart-Landau oscillators [18] and pulse-coupled neuronal networks [2]. In the latter context, splay states have been mainly investigated in leaky-integrate-and-fire (LIF) neurons [2,3,19-21], but some studies have been also devoted to the θ-neurons [22] and to more realistic neuronal models [23]. Finally, splay states are important in that they provide the simplest instance of asynchronous behaviour and can be thereby used as a testing ground for the stability of a more general class of dynamical regimes.

Our model neurons are characterised by a membrane potential u that is continuously driven by the velocity field from the resetting value toward the threshold (see the next section for a more precise definition). As threshold and resetting value can be identified with one another and thereby u interpreted as a phase, it will be customary to refer to the case as to that of a discontinuous velocity field. Additionally, we assume that the post-synaptic potential (PSP) has a stereotyped shape, the so-called α-pulse, that is characterised by identical rise and decay time [2]. The linear stability analysis reveals that the eigenvectors are characterised by different spatial frequencies (when moving from the neuron with the smallest to that one with the largest membrane potential). It is therefore convenient to use the frequency ϕ (scaled to the average phase separation ) to parametrise the Floquet spectrum. As already discussed in [21], there exist two components, namely short (SWs) and long (LWs) wavelengths. SWs vary on ‘microscopic’ scales, i.e. correspond to finite ϕ values: they are typically marginally stable in the thermodynamic limit (). LWs vary on scales of order , i.e. correspond to vanishing frequencies: they have been studied in the context of mean-field theory, i.e. by analysing a suitable functional equation for the probability distribution of the membrane potential u[2,24]. By developing an approach that is valid for arbitrary coupling strength and is perturbative in the inverse system-size , here we prove that the Floquet spectrum scales as and is proportional to . We are also able to determine the spectral shape and find it to be independent of the structure of the velocity field. The transition from SWs to LWs is signalled by the occurrence of a singularity in the spectrum for the frequency . In the crossover region (where k is large but small compared to N), we show that the exponents remain finite and coincide with those determined in the weak coupling limit by Abbott and Van Vreeswijk [2] with their mean field approach. This result is non-trivial, since it is not a priori obvious that the ‘macroscopic’ description discussed in [2] is fully contained in the ‘microscopic’ description derived in this paper, as they refer to two different levels of description.

More specifically, we first build a suitable event-driven map and expand it in powers of (a posteriori, we have verified that it is necessary to reach the fourth order). Afterwards, an expression of the splay state is determined: this task corresponds to finding a fixed point of the event-driven map in a suitably moving reference frame - analogously to what has previously been done in specific contexts [21,25,26]. In practise this task is carried out by first taking the continuum limit for various orders and then obtaining suitable differential equations. The solutions of such equations show that all finite-size corrections for both the period T and the membrane potential vanish up to the third order. Next, the stability analysis is carried out to determine the leading term of the Floquet spectrum. This task involves the introduction of a suitable Ansatz to decompose each eigenvector into the linear superposition of a slow and rapidly oscillating component. The following continuum limit shows that the two components satisfy an ordinary and a differential equation, respectively.

Altogether, the proof of our main result requires determining all terms up to the third order in the expansion of the splay state solution, while some third-order terms are not necessary for the tangent space analysis. Going beyond discontinuous fields would require extending our analysis to account for higher-order terms and this might not even be sufficient to characterise analytic velocity fields. In fact, previous numerical simulations [25] suggest that the Floquet exponents scale with higher powers of depends on which derivatives of are eventually discontinuous. Moreover, it is worth recalling that in the case of a strictly sinusoidal field, the theorem proved by Watanabe and Strogatz [27] implies that Floquet exponents ( for a splay solution) vanish exactly for any value of N.

From the analysis of the SW spectra, one can conclude that the splay state is stable in excitatory (inhibitory) networks whenever (). These conditions can be extended also to the crossover region, where our results coincide with those obtained in [2] (in the limit of a small coupling). Our analytical studies cannot, however, predict the behaviour of the LW component that may be responsible for the emergence of new collective solutions in excitatory networks [3,28]. The overall scenario is partially reminiscent of the stability of synchronous and clustered regimes that is determined by the sign of the first derivative of the velocity-field averaged on the interval (the latter problem has been investigated in excitatory pulse-coupled integrate-and-fire oscillators subject to δ-pulses [1,29]).

Section 2 is devoted to the introduction of the model and to a brief presentation of the main results, including an expression for the leading correction to the period for the LIF model, to provide evidence that it is typically of the fourth order. A general perturbative expression for the map is derived in Section 3, while Section 4 is devoted to deriving the splay state solution up to the third order in . The main result of the paper is discussed in Section 5, where the Floquet spectra are finally obtained. Section 6 contains some general remarks and a discussion of open problems. The technical details of some lengthy calculations have been confined in the appendices: Appendix A is devoted to the derivation of the splay state solution; Appendix B contains the derivation of the leading term (of order four) of the period T for the LIF model; Appendix C is concerned with the linear stability analysis.

2 Model and main results

We consider a network of N identical neurons (rotators) coupled via a mean-field term. The dynamics of the ith neuron writes as

where represents the membrane potential, is the forcing field and g is the coupling constant. When the membrane potential reaches the threshold value , a spike is sent to all neurons (see below for the relationship between the single spikes and the global forcing field ) and it is reset to . The resetting procedure is an approximate way to describe the discharge mechanism operating in real neurons. The function F represents a velocity field for the isolated neuron and it is assumed to be everywhere positive (thus ensuring that the neurons repetitively fire, since they are supra-threshold), while is the velocity field seen by the neuron i in the presence of a coupling with other neurons. While we consider both excitatory () and inhibitory networks (), it is easy to show that ℱ remains always positive to ensure the existence of splay states. For the simple choice

the model reduces to the well-known case of LIF neurons. The evolution of a membrane potential for a LIF suprathreshold neuron () is reported in Figure 1(a).

Fig. 1

(a) Temporal course of the membrane potential for a suprathreshold LIF neuron in the absence of synaptic stimuli. The pulses are exactly collocated in correspondence of the firing times. represents the reset value and the threshold value. (b) ‘Raster plot’ for a splay state in a network of neurons: neuron indices are reported as a function of firing times. The dots correspond to the times when the neurons emit an action potential. Indices are ordered according to the (first) time the neurons have reached the threshold.

The field E is the linear superposition of the pulses emitted in the past when the membrane potential of each single neuron had reached the threshold value. By following Ref. [2], we assume that the shape of a pulse emitted at time is given by , where is the pulse-width. This is equivalent to saying that the total field evolves according to the equation

where the sum in the rhs represents the source term due to the spikes emitted at times .

It is convenient to transform the continuous-time model into a discrete-time mapping. We do so by integrating the equations of motion from time to time (where is the time immediately after the nth pulse has been emitted). The resulting map for the field variables reads

where is the interspike time interval and, for the sake of simplicity, we have introduced the new variable .

In this paper we focus on a specific solution of the network dynamics, namely on splay states, which are asynchronous states, where all neurons fire periodically with the period T and two successive spike emissions occur at regular intervals . A typical ‘raster plot’ for this state is reported in Figure 1(b).

In the large-N limit, it is natural to consider as a smallness parameter and thereby to expand the evolution equations in powers of . In order to perform this expansion, the unique condition to require is the differentiability of the velocity field in the definition interval . The only exception is represented by the boundaries of the interval where discontinuity is allowed.

The first result of this paper is that under the assumption that the velocity field is differentiable at least four times, the dependence of the period T onto the size N is of order . In the specific case of LIF neurons, we show in Appendix B that the leading correction δT to the infinite size result is indeed of order and, more precisely,

where encodes the information on the pulse dynamics (see Eq. (71)). We did not dare to estimate the quartic contribution for generic velocity fields, not only because the algebra would be utterly complicated, but also since our main motivation is to determine the leading contributions in the stability analysis, and it turns out that it is sufficient to determine the splay state up to the third order.

The study of the stability requires determining the Floquet spectrum, i.e. the complex eigenvalues of a given periodic orbit of the period T. With reference to a system of size N and by following [21], the Floquet multipliers can be written as

where represents a zeroth-order phase, while and are the real and imaginary parts of the Floquet exponents, respectively.

Notice that since the total number of exponents is (the zero exponent has been removed by taking the Poincaré section), we are going to miss one of them. Furthermore, as shown in [21], the Nth and th exponents are associated to the field evolution and they will be not considered in our analysis.

In the following we prove that the leading term of the spectrum is

i.e. for discontinuous velocity fields, the real part of the spectrum scales as , while the imaginary part is of even higher order.

For continuous fields, it has been numerically observed that the scaling of the spectrum is at least [25]. In other words, the shape of the spectrum is universal, apart from a multiplicative factor that vanishes if and only if , i.e. for true phase rotators where coincides with .

In the limit , Eq. (7) is singular and seems to diverge. However, one should recall that , i.e. it is a discrete variable in a finite system. By expanding for , one finds that

This expression, which holds in the limit of (), can be compared with the results obtained in [2], in the small-coupling limit (), for sufficiently large k. The equivalence of the methods has two implications: (i) the crossover component is ‘universal’ as it is valid also for large coupling constants; (ii) the ‘macroscopic’ stability is fully contained in our ‘microscopic’ analysis. In fact, numerical studies reveal a perfect correspondence also for the LW component that is not amenable to an analytic treatment [25,30].

From Eqs. (7) and (8), it follows that the stability of the splay state can be inferred, for arbitrary coupling strength, from the sign of : in excitatory (inhibitory) networks, the state is stable whenever (). The same result was previously reported in the weak coupling limit in [2].a It is, however, necessary to point out that such condition(s) do not account for instabilities that can arise in the LW component. This is, e.g., the case of the onset of partial synchronisation via a supercritical Hopf bifurcation [3,28].

3 Event-driven map

By following Ref. [25,31], it is convenient to pass from a continuous to a discrete time evolution rule by introducing the event-driven map which connects the network configuration at subsequent spike emissions occurring at time and . The membrane-potential value just before the emission of the st spike can be obtained by formally integrating Eq. (1),

where the minus superscript means that the map construction has not yet been completed. This task is accomplished by ordering the membrane potentials from the largest () to the smallest () value and by passing to a comoving frame that advances with the firing neuron, i.e. by shifting the neuron index by one unit , where the first subscript indicates that the variable is determined at time . As a result, the evolution is described by an event-driven map that is composed of Eq. (4) together with the following recursive relation:

where . The time at which the st spike is emitted can be determined implicitly from Eq. (10) by setting since, by definition of the model, . In this reference frame the splay state corresponds to a fixed point of the map.

Now we perform a perturbative expansion of both terms and (see Eq. (9)) in powers of . This is justified by the smallness of δt which is . More precisely, the first integral appearing on the rhs of Eq. (9) is solved perturbatively by introducing a polynomial expansion of around ,

Explicit expressions for the time derivates of can be obtained from Eq. (1) and its time derivatives,

where one can further eliminate with the help of Eq. (3).

By inserting the expansion (11) into the expression of , expanding the function around and performing the trivial integrations, one obtains

where we have introduced the short-hand notation for (and analogously for ℱ).

In the case of , it is possible to obtain an exact expression for the integral, which can be then expanded in powers of as follows:

By finally assembling Eqs. (10), (12), (13), we obtain a perturbative expression for the evolution rule of the membrane potential,

4 Splay state solution

The splay state is a fixed point of the event-driven mapping with a constant interspike interval . Since the fixed-point solutions do not depend on the index n, they are denoted as

Substituting Eq. (15) into Eq. (4), one obtains

By introducing Eq. (16) in Eq. (10) and eliminating the n dependence, we obtain a recursive equation for the variable . The fixed-point solution corresponds to the ‘trajectory’ that, starting from , ends in and can be found by tuning the ‘parameter’ . The existence of one or more solutions is related to the dependence on T. Simple calculations reveal that , i.e. it is independent of T; moreover, since is the integral from time 0 to time of a positive defined function (), it is a monotonically increasing function of T which vanishes in zero. Accordingly, for any function F, the minimal value of is g (obtained for ), while the maximal value is unbounded from above. Therefore, there exists one and only one solution provided that .

The variables and depend explicitly on N, and it is natural to expect that the period T itself varies with N. By replacing the expansion of T in powers of the smallness parameter ,

in Eq. (16) we obtain perturbative expressions for , and , namely Eqs. (60), (61), (62) reported in Appendix A. As for the membrane potential, it is necessary to introduce the formal expansion

Substituting the expansions of T, , and into Eq. (14) (after having dropped the n dependence and expanded the ‘F’ functions), one obtains an equation for and , namely

where the variables are defined in Appendix A. Notice the dependence on variables is hidden in terms.

In the large-N limit, one can introduce the continuous spatial coordinate . In practise, this is tantamount to write

It is important to stress that the event-driven neuronal evolution in the comoving frame implies that , i.e. the first neuron will fire at the next step, and , i.e. the membrane potential of the last neuron has been just reset to zero. This implies that and , while for any .

Furthermore, by expanding around , one obtains

By inserting this expansion into Eq. (19), we obtain an equation that can be effectively split into terms of different order that will be analysed separately. Notice that by retaining terms of order h, it is possible to determine the original variables at order .

4.1 Zeroth-order approximation

By assembling first-order terms, we obtain the evolution equation for the zeroth-order membrane potential, namely

This equation is equal to the evolution equation of the membrane potential for a constant field E, with x playing the role of (inverse) time. Please notice that up to the first order, (see Eq. (61)). An implicit and formal solution of Eq. (22) is

where we have imposed the condition . However, there is the second condition to impose, namely . This second condition transforms itself in the equation defining the interspike time interval , when (i.e. in the thermodynamic limit)

This result is, so far, quite standard and could have been easily obtained by just assuming a constant field E in Eq. (1). If we introduce the formal relation in Eq. (22), we obtain

which can be easily integrated

giving the following relation (already derived in [28], by following a different approach)

where, for later convenience, we have introduced

and where, for the sake of simplicity, the prime denotes derivative with respect to the variable and the dependence of F and on has been dropped.

4.2 First-order approximation

By collecting the terms of order , one obtains

An explicit expression for the second derivative of appearing in Eq. (29) can be computed by deriving Eq. (22) with respect to x. This allows rewriting Eq. (29) in a simplified form, namely

By imposing , one obtains the general solution of Eq. (30),

where is defined by Eq. (28). The further condition to be satisfied, , implies and thereby we have , i.e. first-order corrections vanish both for the period and the membrane potential.

4.3 Second-order approximation

Second-order corrections can be estimated by assembling terms of order and by imposing the previously determined conditions and ,

Once evaluated from Eq. (22), the above ODE reduces to

which has the same structure as Eq. (30). Since one has also to impose the same boundary conditions as for the first order, namely , we can conclude that and, consequently, . Therefore, second-order corrections are absent too.

4.4 Third-order approximation

By assembling terms of order , once imposed that first- and second-order corrections vanish, one obtains

By replacing with its expression derived from Eq. (22), Eq. (33) takes the same form as in the two previous examined cases, namely

Therefore, we can safely conclude that third-order terms vanish too.

The LIF model can be solved exactly for any value of N, starting from the asymptotic value (). As shown in Appendix B, it turns out that the leading corrections are of the fourth order for both the period T and the membrane potential.

5 Linear stability analysis

The fixed-point analysis has revealed that the finite-size corrections to the stationary solutions are of order . Since such deviations do not affect the leading terms of the linear stability analysis (as it can be verified a posteriori) they will be simply neglected. Therefore, for the sake of simplicity, from now on and will be simply referred to as T and .

The evolution rule in the tangent space is obtained by differentiating Eq. (4) and then by expanding in powers of τ (this is equivalent to expanding in powers of , as the dependence of τ on N would only generate higher-order terms),

where the variables are polynomials of τ defined in Appendix C.1.

By further differentiating Eq. (14) around the fixed-point solution, one obtains

Finally, can be determined by differentiating Eq. (14) for

where the auxiliary variables are defined in Appendix C.1.

As usual, the eigenvalue problem can be solved by introducing the Ansatz,

where labels the eigenvalues, which must also be expanded as

where is, in principle, a complex number and, for the sake of simplicity, we have dropped its dependence on k. Finally, as already shown at the zeroth order, the eigenvalues correspond to a pure rotation (specified by ) with no expansion or contraction, i.e..

By inserting the above Ansätze in the map expression (35), (36), (37), (38), one obtains, after eliminating δP, δE and δτ, a closed equation for ,

that is the object of our investigation. The overline means that the function is evaluated in , corresponding to the infinite N limit.

5.1 Continuum limit

Similarly to the splay state estimation, it is convenient to take the continuum limit. However, at variance with the previous case, now one should take into account also the presence of fast scales associated to the ‘spatial’ dependence of .

Therefore, the correct Ansatz is slightly more complicated, and we have to separate slowly and rapidly oscillating terms,

where the complex exponential term accounts for the fast oscillations of the eigenvectors, while

are slowly varying variables. On the one hand, the existence of the slow component follows from the analogy with the real-space evolution. On the other hand, the presence of the rapidly oscillating terms , first noticed in Ref. [2] in the uncoupled limit, suggests the presence of the second slow field, namely . Anyway, the correctness and uniqueness of Ansatz (42) is ensured a posteriori by the consistency of the equations that are obtained for the various perturbation orders.

Next, we can finally introduce the continuous variable , as previously done in a real space (see Eq. (20)),

where . This allows expanding around , similarly to what has been done in Eq. (21). At variance with the computation of the fixed point, now there are also terms like and , whose computation requires a similar expansion, but around . By incorporating all the expansion terms within Eq. (41), we have finally an equation, where terms of different orders are naturally separated from one another. The calculations are summarised in Appendix C, and the final equation is (86). By separately treating the different orders, we obtain differential and ordinary equations for the Θ and Π variables. It turns out that it is necessary to consider in parallel different orders in the fast and slow terms to obtain Θ and Π to the same order. As a consequence, we will see that it is sufficient to expand up to order .

5.2 Zeroth-order approximation

By assembling terms of order in Eq. (86), multiplied by the fast oscillating factor , we obtain a first-order linear differential equation for , namely

where is the first-order correction to the Floquet exponent which should be determined. It is important to remind that the prime denotes derivative with respect to the variable , which has been simply redefined U, as previously mentioned. The solution is

where we made use of the definition (28) and is a suitable integration constant.

By assembling now the slow terms of the zeroth order and reminding the definition of , we find the following algebraic equation:

With the help of Eq. (46), we obtain

We can now impose the boundary condition . This implies that

From Eq. (27), we can conclude that , for (values outside this range give identical solutions). Since the vector was assumed to be characterised by the phase (see Eq. (42)), the phase factor would imply (through Eq. (46)) that has to be shifted by some amount, contrary to the initial assumption. Accordingly, the only solution consistent with the original Ansatz is , i.e., and from Eq. (46),

i.e. the eigenvectors are independent of the phase and equal to one another. In other words, the degeneracy has not been lifted.

5.3 First-order approximation

By assembling the fast terms of order and by setting , we find that satisfies the following first order differential equation:

whose solution is

where is an integration constant associated with the solution of the previous equation.

By collecting the slow terms of order in Eq. (86), one obtains the algebraic equation

whose solution is

By imposing the boundary condition , it is possible to evaluate ,

Again from Eq. (27) and using the same argument as in the previous section, we find that and thereby (from Eq. (51))

Altogether, we can conclude that the second-order correction to the Floquet exponent vanishes as well, and one cannot lift the degeneracy among the eigenvectors.

5.4 Second-order approximation

By assembling fast terms of order appearing in Eq. (86) and by setting , the following first-order differential equation for can be derived:

whose solution is

where is an integration constant associated with the solution of the previous differential equation.

Furthermore, by collecting the slow terms of order , we obtain the algebraic equation

By imposing that the above equation is satisfied for , it reduces to

Finally, by imposing the boundary condition , it is possible to determine ,

Accordingly, is real and depends on the difference between and , confirming the numerical findings in [25]. Therefore, the imaginary terms are smaller than .

In the specific example of a leaky integrate-and-fire neuron, the expression for reduces to

since, by using the equations that characterise LIF neurons, the following relation holds:

All in all, Eq. (57) generalises the expression found for the LIF model Eq. (58) [25].b

6 Conclusions

We have derived analytically the short-wavelength component of the Floquet spectrum of the splay solution in a fully-coupled network composed of generic suprathreshold pulse-coupled phase-like neurons in the large-N limit. Our analysis has revealed that, for discontinuous velocity fields, the spectrum scales as and the stability is controlled by the sign of the difference between the velocity at reset and at threshold. The shape of the spectrum is otherwise independent of the velocity field. It would be interesting to investigate the role of the pulse shape as well. As long it follows from a linear evolution equation, such as Eq. (3), it is indeed possible to replicate the analysis carried out in this paper. Numerical studies suggest that the scaling behaviour is truly universal, while the shape of the Floquet spectrum depends on the pulse shape [30]. It would be interesting to discover whether and which pulse shapes may give rise to SW instabilities.

Networks of LIF neurons coupled via δ-like pulses are characterised by a finite (in)stability of the whole spectrum [21]. The difference with the case of α-pulse is so strong that it cannot be reconciled even by taking the limit (zero pulse-width), indicating that the limits and zero pulse-width do not commute [21]. This reveals that even the development of general stability theory of the simple splay states requires some further progress.

Finally, notice that although our analytical approach is able to cover the entire SW component and the crossover region, it does not cover the truly long wavelengths which require going beyond a perturbative approach.

Appendix A: Fixed-point expansion (general case)

The expansion of the exact expressions (16) for and leads to

where we have reported also the expansion of that is necessary to pass from expression (14) to (19). Please notice that while the membrane potentials and the period are expanded up to , as in (18) and (17), here we limit the expansion to terms, since the field variables appearing in the event-driven map are integrated over an interspike-interval (see (9)).

To proceed further, we need also to introduce the expansions of the velocity field and of its derivatives,

where the overline means that the function is computed in , which corresponds to the infinite N limit.

By replacing the membrane potentials, the period, the self-consistent fields and the velocity field with their expansions, the event-driven map (14) can be formally rewritten for the splay state as (19) with the introduction of the following auxiliary variables:

Appendix B: Fixed-point expansion (LIF model)

In the case of the LIF neuron (see Eq. (2)), the fixed point of the event-driven map reads

where

Its solution is

By expanding Eq. (69) for and for a generic j, one can derive perturbative expressions for the period T and the membrane potential, respectively. Let us start by substituting the expressions (17), (60), (61), (62) in Eq. (68). This leads to the expansion

where

accounts for the dependence on the field dynamics. Now, with the help of Eqs. (17), (70) and expanding the exponential terms up to the fourth order, we obtain a closed equation for the interspike interval,

where

while identifies a term of order that is multiplied by . Since, while proceeding from lower- to higher-order terms, we find that (for ), it is not necessary to give the explicit expression of the functions as they do not contribute at all.

One can equivalently expand

where

while we do not provide explicit expressions for as they turn out to be irrelevant.

Now we are in the position to analyse the different orders.

B.1 Zeroth order

By assembling the terms of order one in Eq. (72), we obtain

This is an implicit definition of the asymptotic interspike time

Analogously, we can find an explicit equation for the membrane potential by assembling the terms of order one in Eq. (73)

In the thermodynamic limit, the solution for becomes

which coincides with Eq. (23) with .

B.2 From first to third order

By separately assembling the terms of order (for ) in Eq. (72), we obtain

which implies that since . Moreover, by assembling the terms of order in Eq. (73), we obtain

which thereby implies that first-, second- and third-order corrections vanish also for the membrane potential.

B.3 Fourth order

The order which reveals a different scenario is the fourth one. By assembling the terms of order in Eq. (72), we obtain

whose explicit expression is reported in Eq. (5). By analogously assembling the terms of order in Eq. (73), we obtain

which becomes, in the thermodynamic limit,

Appendix C: Expansion in tangent space around the fixed point

C.1 Introduction

The auxiliary variables required to complete the definition of the tangent-space evolution rule (35), (36), (37), (38) are as follows:

where the dependence of τ on n has been dropped, since we are considering a linearisation around the splay state.

In order to find the Floquet eigenvalue , one should substitute the Ansätze (39) into Eqs. (35), (36). This allows to find explicit expressions for δP and δE as a function of , τ and δτ, namely

where we have introduced the shorthand notation

By substituting δP and δE, as given by (82) and (83), into Eq. (38), we can express δτ directly in terms of

where we exploited the equality , which follows from the fact that in the thermodynamic limit .

By inserting the expressions in Eqs. (82), (83), (85) into Eq. (37), we find a single equation for the eigenvalues and eigenvectors,

By now substituting the expansion (40) and retaining the leading terms, we obtain Eq. (41).

C.2 limit

Once the continuous variables (44) have been introduced, it is necessary to estimate and by expanding such variables around zero. By inserting the resulting expansion for into the expressions for and , we obtain, respectively,

An analogous procedure for leads to

where , , are defined according to the following equations:

We now expand up to the order , thus neglecting higher orders because they contribute to the definition of Π variable and we need terms at least of order (one order lower than needed to define Θ). By inserting the Ansatz (42) and the previous expansions in Eq. (41), we finally obtain a closed equation for the eigenvalues and eigenvectors,

where we have introduced the shorthand notation ℬ in order to characterise a term of order , whose explicit expression is not necessary, since it turns out to contribute to the definition of the Π variable, and it is, therefore, one order beyond what we need. Moreover, notice that the terms appearing within round brackets in the rhs of the above equation can be shown to be zero due to exact algebraic cancellations that emerge from the solution of the equation order by order. Finally,

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

SO carried out the analytical calculations and drafted the manuscript. AT participated in the design of the study and in the drawing up of the article. AP conceived the study and participated in its design and coordination. All authors read and approved the final manuscript.