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

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 1 / N . 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 F ( u ) from the resetting value u = 0 toward the threshold u = 1 (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 F ( 1 ) ≠ F ( 0 ) 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 1 / α [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 1 / N ) 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 ( N → ∞ ). LWs vary on scales of order O ( 1 ) , 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 1 / N , here we prove that the Floquet spectrum scales as 1 / N 2 and is proportional to F ( 1 ) − F ( 0 ) . 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 ϕ → 0 . In the crossover region ϕ = k / N (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 1 / N (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 1 / N 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 1 / N depends on which derivatives of F ( u ) 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 N − 3 Floquet exponents ( N − 2 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 F ( 0 ) > F ( 1 ) ( F ( 0 ) < F ( 1 ) ). 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 d F / d u of the velocity-field averaged on the interval [ 0 , 1 ] (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 1 / N . 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.

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

where u i ( t ) represents the membrane potential, E ( t ) is the forcing field and g is the coupling constant. When the membrane potential reaches the threshold value u i ( t ) = u th = 1 , a spike is sent to all neurons (see below for the relationship between the single spikes and the global forcing field E ( t ) ) and it is reset to u i ( t ) = u R = 0 . 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 F i is the velocity field seen by the neuron i in the presence of a coupling with other neurons. While we consider both excitatory ( g > 0 ) and inhibitory networks ( g < 0 ), 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 ( a > 1 ) is reported in Figure 1(a).

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 t = 0 is given by E s ( t ) = α 2 t N e − α t , where 1 / α 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 t n < t .

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 t n to time t n + 1 (where t n is the time immediately after the nth pulse has been emitted). The resulting map for the field variables reads

where τ n = t n + 1 − t n is the interspike time interval and, for the sake of simplicity, we have introduced the new variable P : = α E + E ˙ .

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 τ n ≡ T / N . A typical ‘raster plot’ for this state is reported in Figure 1(b).

In the large-N limit, it is natural to consider 1 / N as a smallness parameter and thereby to expand the evolution equations in powers of 1 / N . In order to perform this expansion, the unique condition to require is the differentiability of the velocity field F ( u ) in the definition interval ] 0 , 1 [ . 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 F ( u ) is differentiable at least four times, the dependence of the period T onto the size N is of order o ( 1 / N 3 ) . 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 O ( 1 / N 4 ) and, more precisely,

where K ( α ) 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 ϕ k represents a zeroth-order phase, while λ k and ω k are the real and imaginary parts of the Floquet exponents, respectively.

Notice that since the total number of exponents is N + 1 (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 ( N + 1 ) 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 1 / N 2 , 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 O ( 1 / N 4 ) [25]. In other words, the shape of the spectrum is universal, apart from a multiplicative factor that vanishes if and only if F ( 1 ) = F ( 0 ) , i.e. for true phase rotators where u = 0 coincides with u = 1 .

In the limit ϕ k → 0 , Eq. (7) is singular and λ k seems to diverge. However, one should recall that ϕ k = 2 π k / N , i.e. it is a discrete variable in a finite system. By expanding for k ≪ N , one finds that

This expression, which holds in the limit of ( 1 ≪ k ≪ N ), can be compared with the results obtained in [2], in the small-coupling limit ( | g | ≪ 1 ), 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 F ( 1 ) − F ( 0 ) : in excitatory (inhibitory) networks, the state is stable whenever F ( 0 ) > F ( 1 ) ( F ( 0 ) < F ( 1 ) ). 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].

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 t n and t n + 1 . The membrane-potential value u i − ( t n + 1 ) just before the emission of the ( n + 1 ) 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 ( j = 1 ) to the smallest ( j = N ) value and by passing to a comoving frame that advances with the firing neuron, i.e. by shifting the neuron index by one unit u n + 1 , j − 1 = u j − ( t n + 1 ) , where the first subscript indicates that the variable is determined at time t n + 1 . 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 τ n = t n + 1 − t n . The time at which the ( n + 1 ) st spike is emitted can be determined implicitly from Eq. (10) by setting j = 1 since, by definition of the model, u n , 0 ≡ 1 . In this reference frame the splay state corresponds to a fixed point of the map.

Now we perform a perturbative expansion of both terms A 1 and A 2 (see Eq. (9)) in powers of δ t = t − t n . This is justified by the smallness of δt which is O ( 1 / N ) . More precisely, the first integral appearing on the rhs of Eq. (9) is solved perturbatively by introducing a polynomial expansion of u i ( t ) around t = t n ,

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

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

By inserting the expansion (11) into the expression of A 1 , expanding the function F ( u ) around u n , j and performing the trivial integrations, one obtains

where we have introduced the short-hand notation F n , j for F ( u n , j ) (and analogously for ℱ).

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

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

The splay state is a fixed point of the event-driven mapping with a constant interspike interval τ = T / N . 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 u ˜ j . The fixed-point solution corresponds to the ‘trajectory’ that, starting from u ˜ N = 0 , ends in u ˜ 0 = 1 and can be found by tuning the ‘parameter’ T / N . The existence of one or more solutions is related to the dependence on T. Simple calculations reveal that A 2 = g / N , i.e. it is independent of T; moreover, since A 1 is the integral from time 0 to time T / N of a positive defined function ( F > 0 ), it is a monotonically increasing function of T which vanishes in zero. Accordingly, for any function F, the minimal value of u 0 is g (obtained for T = 0 ), while the maximal value is unbounded from above. Therefore, there exists one and only one solution provided that g < 1 .

The variables E ˜ and P ˜ 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 1 / N ,

in Eq. (16) we obtain perturbative expressions for P ˜ , E ˜ and E ˜ ˙ , 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, u ˜ j , E ˜ and P ˜ into Eq. (14) (after having dropped the n dependence and expanded the ‘F’ functions), one obtains an equation for u ˜ j ( h ) and T ( h ) , namely

where the Q variables are defined in Appendix A. Notice the dependence on T ( h ) variables is hidden in Q terms.

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

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

Furthermore, by expanding U ( h ) ( x ) around x = j / N , 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 h − 1 .

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, E ˜ = 1 / T ( 0 ) (see Eq. (61)). An implicit and formal solution of Eq. (22) is

where we have imposed the condition U ( 0 ) ( 1 ) = 0 . However, there is the second condition to impose, namely U ( 0 ) ( 0 ) = 1 . This second condition transforms itself in the equation defining the interspike time interval T ( 0 ) , when N → ∞ (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 F ′ [ U ( 0 ) ( x ) ] = d F ( U ( 0 ) ) d U ( 0 ) 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 U ( 0 ) and the dependence of F and F ′ on U ( 0 ) has been dropped.

4.2 First-order approximation

By collecting the terms of order 1 / N 2 , one obtains

An explicit expression for the second derivative of U ( 0 ) ( x ) 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 U ( 1 ) ( 1 ) = 0 , one obtains the general solution of Eq. (30),

where H ( x ) is defined by Eq. (28). The further condition to be satisfied, U ( 1 ) ( 0 ) = 0 , implies T ( 1 ) = 0 and thereby we have U ( 1 ) ( x ) ≡ 0 , 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 1 / N 3 and by imposing the previously determined conditions T ( 1 ) = 0 and U ( 1 ) ( x ) = 0 ,

Once evaluated d 3 U ( 0 ) / d x 3 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 U ( 2 ) ( 0 ) = U ( 2 ) ( 1 ) = 0 , we can conclude that T ( 2 ) = 0 and, consequently, U ( 2 ) ( x ) ≡ 0 . Therefore, second-order corrections are absent too.

4.4 Third-order approximation

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

By replacing d 4 U ( 0 ) / d x 4 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 ( N → ∞ ). 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.

The fixed-point analysis has revealed that the finite-size corrections to the stationary solutions are of order o ( 1 / N 3 ) . 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 T ( 0 ) and u ˜ j ( 0 ) will be simply referred to as T and u ˜ j .

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 1 / N , as the dependence of τ on N would only generate higher-order terms),

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

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

Finally, δ τ n can be determined by differentiating Eq. (14) for j = 1

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

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

where μ k labels the eigenvalues, which must also be expanded as

where Γ ( h ) 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 ϕ k ) with no expansion or contraction, i.e. Γ ( 0 ) = 0 .

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 δ u j ,

that is the object of our investigation. The overline means that the function is evaluated in u ˜ j ( 0 ) , 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 ϕ k .

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 π j follows from the analogy with the real-space evolution. On the other hand, the presence of the rapidly oscillating terms e i ϕ k j , first noticed in Ref. [2] in the uncoupled limit, suggests the presence of the second slow field, namely ϑ j . 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 x = j / N , as previously done in a real space (see Eq. (20)),

where h = 0 , … , 3 . This allows expanding δ u j − 1 = π j − 1 + ϑ j − 1 e i ϕ k ( j − 1 ) around x = j / N , similarly to what has been done in Eq. (21). At variance with the computation of the fixed point, now there are also terms like U ( 1 / N ) and δ U ( 1 / N ) , whose computation requires a similar expansion, but around x = 0 . 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 δ U ( 1 / N ) up to order O ( 1 / N 3 ) .

5.2 Zeroth-order approximation

By assembling terms of order O ( 1 / N ) in Eq. (86), multiplied by the fast oscillating factor e i ϕ k j , we obtain a first-order linear differential equation for Θ ( 0 ) , namely

where Γ ( 1 ) 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 U ( 0 ) , which has been simply redefined U, as previously mentioned. The solution is

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

By assembling now the slow terms of the zeroth order and reminding the definition of F ( U ( x ) ) , we find the following algebraic equation:

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

We can now impose the boundary condition δ U ( 0 ) ( x = 1 ) = Θ ( 0 ) ( 1 ) + Π ( 0 ) ( 1 ) = 0 . This implies that

From Eq. (27), we can conclude that Γ ( 1 ) = 2 π i m , for 0 ≤ m < N (values outside this range give identical solutions). Since the vector δ u j was assumed to be characterised by the phase j ϕ j (see Eq. (42)), the phase factor Γ ( 1 ) = 2 π i m would imply (through Eq. (46)) that ϕ k has to be shifted by some amount, contrary to the initial assumption. Accordingly, the only solution consistent with the original Ansatz is m = 0 , i.e. Γ ( 1 ) = 0 , and from Eq. (46),

i.e. the eigenvectors are independent of the phase ϕ k 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 1 / N 2 and by setting Γ ( 1 ) = 0 , we find that Θ ( 1 ) satisfies the following first order differential equation:

whose solution is

where K ( 1 ) is an integration constant associated with the solution of the previous equation.

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

whose solution is

By imposing the boundary condition δ U ( 1 ) ( x = 1 ) = Θ ( 1 ) ( 1 ) + Π ( 1 ) ( 1 ) = 0 , it is possible to evaluate Γ ( 2 ) ,

Again from Eq. (27) and using the same argument as in the previous section, we find that Γ ( 2 ) = 0 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 1 / N 3 appearing in Eq. (86) and by setting Γ ( 1 ) = Γ ( 2 ) = 0 , the following first-order differential equation for Θ ( 2 ) can be derived:

whose solution is

where K ( 2 ) is an integration constant associated with the solution of the previous differential equation.

Furthermore, by collecting the slow terms of order 1 / N 2 , we obtain the algebraic equation

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

Finally, by imposing the boundary condition δ U ( 2 ) ( x = 1 ) = Θ ( 2 ) ( 1 ) + Π ( 2 ) ( 1 ) = 0 , it is possible to determine Γ ( 3 ) ,

Accordingly, Γ ( 3 ) is real and depends on the difference between F ( U ( x = 1 ) ) ≡ F ( 0 ) and F ( U ( x = 0 ) ) ≡ F ( 1 ) , confirming the numerical findings in [25]. Therefore, the imaginary terms ω k are smaller than 1 / N 2 .

In the specific example of a leaky integrate-and-fire neuron, the expression for Γ ( 3 ) 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