Periodic solutions in next generation neural field models

We consider a next generation neural field model which describes the dynamics of a network of theta neurons on a ring. For some parameters the network supports stable time-periodic solutions. Using the fact that the dynamics at each spatial location are described by a complex-valued Riccati equation we derive a self-consistency equation that such periodic solutions must satisfy. We determine the stability of these solutions, and present numerical results to illustrate the usefulness of this technique. The generality of this approach is demonstrated through its application to several other systems involving delays, two-population architecture and networks of Winfree oscillators.


INTRODUCTION
The collective behaviour of large networks of neurons is a topic of ongoing interest.One of the simplest forms of behaviour is a periodic oscillation, which manifests itself as a macroscopic rhythm created by the synchronous firing of many neurons.Such oscillations have relevance to rhythmic movement [36], epilepsy [20,38], schizophrenia [57] neural communication [50] and EEG/MEG modelling [9], among others.In different networks, oscillations may arise from mechanisms such as synaptic delays [14], the interaction of excitatory and inhibitory populations [7,52], having sufficient connectivity in an inhibitory network [58], or the finite width of synaptic pulses emitted by neurons [49].The modelling and simulation of such networks is essential in order to investigate their dynamics.
Among the many types of model neurons used when studying networks of neurons, theta neurons [16], Winfree oscillators [2] and quadratic integrate-and-fire (QIF) neurons [34] are some of the simplest.These three types of model neurons have the advantage that their mathematical form often allows infinite networks of such neurons to be analysed exactly using the Ott-Antonsen method [44,45].We continue along those lines in this paper.
Here we largely consider a spatially-extended network of neurons, in which the neurons can be thought of as lying on a ring.Such ring networks have been studied in connection with modelling head direction [61] and working memory [18,60], for example, and been studied by others [17,22,24].We consider a network of theta neurons.The theta neuron is a minimal model for a neuron which undergoes a saddle-node-on-invariantcircle (SNIC) bifurcation as a parameter is varied [16].The theta neuron is exactly equivalent to a quadratic integrate-and-fire (QIF) neuron, under the assumption of infinite firing threshold and reset values [4,15,37].The coupling in the network is nonlocal synaptic coupling, implemented using a spatial convolution with a translationally-invariant coupling kernel.
We studied such a model in the past [40], concentrating on describing spatially-uniform states and also stationary "bump" states in which there is a spatially-localised group of active neurons while the remainder of the network is quiescent.We determined the existence and stability of such states and also found regions of parameter space in which neither of these types of states were stable.Instead, we sometimes found solutions which were periodic in time.In this paper we study such periodic solutions using a recently-developed technique [39,43] which is significantly faster than the standard approach.
This technique was successfully applied to describe traveling and breathing chimera states in nonlocally coupled phase oscillators [39,43].In this paper, we generalize this technique to neural models and illustrate its possibilities with several examples.We start with a ring network of theta neurons and consider periodic bump states there.We describe a continuation algorithm, perform linear stability analysis of such states and derive some useful formulas, for example, for average firing rates.We show that the same approach works in the presence of delays and for two-population models.Finally, we show that Winfree oscillators are also treatable by the proposed technique.
The structure of the paper is as follows.In Sec.I we present the discrete network model and its continuumlevel description.The bulk of the paper is in Sec II, where we show how to describe periodic states in a selfconsistent way, and show numerical results from implementing our algorithms.Other models are considered in Sec.III and we end in Sec.IV.The Appendix contains useful results regarding the complex Riccati equation.

I. THETA NEURON NETWORK MODEL
The model we consider first is that in [25,31,40], which we briefly present here.The discrete network consists of N synaptically coupled theta neurons described by where each θ j ∈ [0, 2π] is an angular variable.The constant κ is the overall strength of coupling within the network, and the current entering the jth neuron is κI j where where is a pulsatile function with a maximum at θ = π (when the neuron fires) and a n is chosen according to the normalization condition Increasing n makes P n (θ) "sharper" and more pulselike.The excitability parameters η j are chosen from a Lorentzian distribution with mean η 0 and width γ > 0 The connectivity within the network is given by the weights K jk which are defined by K jk = K(2π(j − k)/N ) where the coupling kernel is for some constant A. Note that the form of coupling implies that the neurons are equally-spaced around a ring, with periodic boundary conditions.Such a network can support solutions which are -at a macroscopic levelperiodic in time; see Fig. 3 in [40].Such solutions are unlikely to be true periodic solutions of (1), since for a typical realisation of the η j , one or more neurons will have extreme values of this parameter, resulting in them not frequency-locking to other neurons.
Note that the model we study here has only one neuron at each spatial position, and for A > 0 the connections between nearby neurons are more positive than those between distant neurons.For A > 1 neurons on opposite sides of the domain inhibit one another, as the connections between them are negative, giving a "Mexican-hat" connectivity.For A < 0 connections between neurons on opposite sides of the domain are more positive than those between nearby neurons.Such a model with one population of neurons and connections of mixed sign can be thought of as an approximation of a network with populations of both excitatory and inhibitory neurons [17,48].
Using the Ott/Antonsen ansatz [44,45], one can show that in the limit N → ∞, the long-term dynamics of the network (1) can be described by where is the convolution of K and ϕ and where For our computations we set n = 2, so that where overline indicates the complex conjugate.Periodic solutions like those studied below were also found with n = 5, for example, so the choice of n is not critical.The wider question of the effects of pulse shape and duration is an interesting one [47].Eq. ( 4) is an integro-differential equation for a complex-valued function z(x, t), where x ∈ [0, 2π] is position on the ring.z(x, t) is a local order parameter and can be thought of as the average of e iθ for neurons in a small neighbourhood of position x.The magnitude of z is a measure of how synchronised the neurons are, whereas its argument gives the most likely value of θ [25].Using the equivalence of theta and QIF neurons, one can also provide a relevant biological interpretation of z.Namely, defining W ≡ (1 − z)/(1 + z), one can show [27,37] that π −1 Re W is the flux through θ = π or the instantaneous firing rate of neurons at position x and time t.Similarly, if V j = tan(θ j /2) is the voltage of the jth QIF neuron, the mean voltage at position x and time t is given by Im W .Note that physically relevant solutions of Eq. ( 4) must assume values |z| ≤ 1.In other words, we are interested only in solutions z ∈ D, where is the unit disc in the complex plane.
Equations of the form (4)-( 5) are sometimes referred to as next generation neural field models [8,11] as they have the form of a neural field model (an integro-differential equation for a macroscopic quantity such as "activity" [1,33]) but are derived exactly from a network like (1), rather than being of a phenomenological nature.

II. PERIODIC STATES
In this paper we focus on states with periodically oscillating macroscopic dynamics.For the mean field equation (4) these correspond to periodic solutions z = a(x, t) where a(x, t + T ) = a(x, t) for some T > 0. The frequency corresponding to T is denoted by ω = 2π/T .An example of a periodic solution is shown in Fig. 1, panels (a) and (b).Fig. 1 (c) shows a realization of the same solution in the discrete network (1).We note that this solution has a spatio-temporal symmetry: it is invariant under a time shift of half a period followed by a reflection about x = π.However, we do not make use of this symmetry in the following calculations.
The appearance of periodic solutions in this model is in contrast with classical one-population neural field models for which they do not seem to occur [32].This is another example of next generation neural field models showing more complex time-dependent behaviour than classical ones [31].
A straight-forward way to study a periodic orbit like that in Fig. 1 would be discretise Eq. ( 4) on a spatiallyuniform grid and approximate the convolution using matrix/vector multiplication or otherwise, resulting in a large set of coupled ordinary differential equations.The periodic solution of these could then be studied using standard techniques [25], but note that the computational complexity of this would typically scale as ∼ N 2 , where N is the number of spatial points used in the grid.Instead we propose here an alternative method based on the ideas from [39,43], which allows us to perform the same calculations with only ∼ N operations.The main ingredients of this method are explained in Sections II A and II B. They include the description of the properties of complex Riccati equation and the derivation of the selfconsistency equation for periodic solutions of Eq. ( 4).In Section II C we explain how the self-consistency equation can be solved in the case of coupling function (3).Then in Section II D we report some numerical results obtained with our method.In addition, in Section II E we perform a rigorous linear stability analysis of periodic solutions of Eq. ( 4) by considering the spectrum of the corresponding monodromy operator.Finally, in Section II F, we show how the mean field equation ( 4) can be used to predict the average firing rate distribution in the neural network (1).

A. Periodic complex Riccati equation and Möbius transformation
By the time rescaling u(x, t) = z(x, t/ω) we can rewrite Eq. (4) in the form such that the above T -periodic solution of Eq. ( 4) corresponds to a 2π-periodic solution of Eq. ( 6).Then, dividing (6) by ω and reordering the terms we can see that Eq. ( 6) is equivalent to a complex Riccati equation with the 2π-periodic in t coefficient and In [43] it was shown (see also Proposition 2 and Remark 2 in the Appendix for additional details) that independent of the choice of the real-valued periodic function W (x, t), parameters ζ ∈ C up = {z ∈ C : Im z > 0} and ω > 0, for every fixed x ∈ [0, 2π] Eq. ( 7) has a unique stable 2πperiodic solution U (x, t) that lies entirely in the open unit disc D. Denoting the corresponding solution operator by we can write the 2π-periodic solution of interest as Note that C per ([0, 2π]; R) here denotes the space of all real-valued continuous 2π-periodic functions, while the notation C per ([0, 2π]; D) stands for the space of all complex continuous 2π-periodic functions with values in the open unit disc D. Importantly, the variable x appears in formula (10) as a parameter so that the function W (x, •) ∈ C per ([0, 2π]; R) with a fixed x is considered as the first argument of the operator U.
As for the operator U, although it is not explicitly given, its value can be calculated without resourcedemanding iterative methods by solving exactly four initial value problems for Eq. ( 7).The rationale for this approach can be found in [43,Section 4] and is repeated for completeness in Remark 3 in Appendix.Below we describe its concrete implementation in the case of formula (10).
We assume that the spatial domain [0, 2π] is discretised with N points, x j , j = 1, 2, . . .N and that the functions U (x, t) and W (x, t) are replaced with their grid counterparts u j (t) = U (x j , t) and w j (t) = W (x j , t), respectively.Then, given a set of functions w j (t), we calculate the corresponding functions u j (t) by performing the following four steps.
(iii) Now that the Möbius maps M j (u) are known, their fixed points can be found by solving u * j = M j (u * j ) for each j.This equation is equivalent to a complex quadratic equation, therefore in general it has two solutions in the complex plane.For γ > 0, only one of these solutions lies in the unit disc D.
(iv) Using the fixed points u * j ∈ D of the Möbius maps M j (u) as an initial condition in Eq. ( 7) (i.e.setting u(x j , 0) = u * j ) and integrating (7) for a fourth time to t = 2π we obtain the grid counterpart of a 2π-periodic solution, U (x, t), that lies entirely in the unit disc D.
We now use this result to show how to derive a selfconsistency equation, the solution of which allows us to determine a 2π-periodic solution of Eq. ( 6).

B. Self-consistency equation
Supposing that Eq. ( 6) has a 2π-periodic solution, then using formula (8) we can calculate the corresponding function W (x, t).On the other hand, using formula (10) we can recover u(x, t).Then the new and the old expressions of u(x, t) will agree with each other if and only if the function W (x, t) satisfies a self-consistency equation obtained by inserting ( 10) into (8).In the following we consider Eq. ( 11) as a separate equation which must be solved with respect to W (x, t) and ω.Note that the unknown field W (x, t) has a problem-specific meaning: It is proportional to the current entering a neuron at position x at time t due to the activity of all other neurons in the network.The use of self-consistency arguments to study infinite networks of oscillators goes back to Kuramoto [23,55,56], but such approaches have always focused on steady states, whereas we consider periodic solutions here.Note that from a computational point of view, the selfconsistency equation ( 11) has many advantages.It allows us to reduce the dimensionality of the problem at least in the case of special coupling kernels with finite number of Fourier harmonics (see Section II C).Moreover, its structure is convenient for parallelization, since the computations of operator U at different points x are performed independently.Finally, the main efficiency is due to the fact that the computation of U is performed in non-iterative way.
In the next proposition, we will prove some properties of the solutions of Eq. ( 11), which will be used later in Section II E.
Proposition 1 Let the pair (W (x, t), ω) be a solution of the self-consistency equation (11) and let U (x, t) be defined by (10).Then Proof: For every fixed x ∈ [0, 2π], the function U (x, t) yields a stable 2π-periodic solution of the complex Riccati equation (7).The linearization of Eq. ( 7) around this solution reads where Moreover, using ( 8) and ( 9), we can show that the above expression determines a function identical to the function M (x, t) in (12).Recalling that U (x, t) is not only a stable but also an asymptotically stable solution of Eq. ( 7), see Remark 1 in Appendix, we conclude that the corresponding Floquet multiplier lies in the open unit disc D. This ends the proof.

C. Numerical implementation
Eq. ( 11) describes a periodic orbit, and since Eq. ( 7) is autonomous we need to append a pinning condition in order to select a specific solution of Eq. (11).For a solution of the type shown in Fig. 1 we choose In the following we focus on the case of the cosine coupling (3).It is straight-forward to show that (KH n )(x) can be written as a linear combination of 1, cos x and sin x.However, the system is translationally invariant in x, and we can eliminate this invariance from Eq. ( 11) by restricting this equation to its invariant subspace Span{1, sin x}.Then, taking into account that the function W (x, t) is real, we seek an approximate solution of the system (11), ( 13) using a Fourier-Galerkin ansatz where v m and w m are real coefficients and ψ m (t) are trigonometric basis functions Our typical choice of the number of harmonics in ( 14) is F = 10.To exactly represent W (x, t) in ( 14) would require an infinite number of terms in the series, so using a finite value of F introduces an approximation in our calculations.However, the excellent agreement between our calculations with F = 10 and those from full simulations of (4) (shown below) indicate that such an approximation is justified.
In simple terms, suppose we have somewhat accurate estimates of v 0 , v 1 , . . ., v 2F , w 0 , w 1 , . . ., w 2F , ω.These can be inserted into (14) to calculate the function W (x, t).Then one can calculate a periodic solution of Eq. ( 7) with the specified W (x, t) by formula (10) and finally calculate H n (U) and insert this into ( 15) and ( 16) and perform the projections.We want the difference between the new values of v 0 , v 1 , . . ., v 2F , w 0 , w 1 , . . ., w 2F and the initial values of them to be zero, and for (13) to hold.This determines the 4F + 3 equations we need to solve.The solutions of these equations can be followed as a parameter is varied in the standard way [26].Note that these calculations involve discretising the spatial domain with N points.However, the number of unknowns (4F + 3) is significantly less than N .

D. Results
The results of following the solution shown in Fig. 1 as η 0 is varied are shown in Fig. 2, along with values measured from direct simulations of Eq. ( 4).The period seems to become arbitrarily large as η 0 approaches −2.32, and the solution approaches a heteroclinic connection, spending more and more time near two symmetric states which are mapped to one another under the transformation x → −x.The solution becomes unstable at η 0 ≈ −0.5 through a subcritical torus (or secondary Hopf) bifurcation.(This was determined by finding the Floquet multipliers of the periodic solution in the conventional way; results not shown.)For these calculations we used N = 256 spatial points.
For more negative values of η 0 than those shown in Fig. 2, another type of periodic solution is stable: see Fig. 3.Such a solution does not have the spatio-temporal symmetry of the solution shown in Fig. 1.However, we can follow it in just the same way as the parameter η 0 is varied, and we obtain the results shown in Fig. 4.This periodic orbit appears to be destroyed in a supercritical Hopf bifurcation as η 0 is decreased through approximately −3.2, and become unstable to a wandering pattern at η 0 is increased through approximately −2.34.
Note that the left asymptote in Fig. 2 coincides with the right asymptote in Fig. 4. On the other hand, we note that two patterns shown in Figs. 1 and 3 have different spatiotemporal symmetries, therefore due to topological reasons they cannot continuously transform into each other.Similar bifurcation diagrams where parameter ranges of two patterns with different symmetries are separated by heteroclinic or homoclinic bifurcations were found for non-locally coupled Kuramoto-type phase oscillators [42] and seem to be a general mechanism which, however, needs additional investigation.

E. Stability of breathing bumps
Given a T -periodic solution a(x, t) of Eq. ( 4), we can perform its linear stability analysis, using the approach proposed in [39].Before doing this, we write where to emphasise that H n (z) is always real.Now, we insert the ansatz z(x, t) = a(x, t) + v(x, t) into Eq.( 4)  and linearize the resulting equation with respect to small perturbations v(x, t).Thus, we obtain a linear integro-differential equation and qC q z q−1 .
Note that Eq. ( 17) coincides with the Eq.(5.1) from [40], except that the coefficients a(x, t) and µ(x, t) are now time-dependent.Since Eq. ( 17) contains the complexconjugated term v, it is convenient to consider this equation along with its complex-conjugate This pair of equations can be written in the operator form where and For every fixed t the operators A(t) and B(t) are linear operators from C per ([0, 2π]; C 2 ) into itself.Moreover, they both depend continuously on t and thus their norms are uniformly bounded for all t ∈ [0, T ].
Recall that the question of linear stability of a(x, t) in Eq. ( 4) is equivalent to the question of linear stability of the zero solution in Eq. ( 17), and hence to the question of linear stability of the zero solution in Eq. (19).Moreover, using the general theory of periodic differential equations in Banach spaces, see [13, Chapter V], the last question can be reduced to the analysis of the spectrum of the monodromy operator E(T ) defined by the operator exponent The analysis of Eq. ( 19) in the case when A(t) is a matrix multiplication operator and B(t) is an integral operator similar to (20) has been performed in [39,Section 4].Repeating the same arguments we can demonstrate that the spectrum of the monodromy operator E(T ) is bounded and symmetric with respect to the real axis of the complex plane.Moreover, it consists of two qualitatively different parts: (i) the essential spectrum, which is given by the formula (ii) the discrete spectrum σ disc that consists of finitely many isolated eigenvalues λ, which can be found using a characteristic integral equation, as explained in [39, Section 4].
Note that if a(x, t) is obtained by solving the selfconsistency equation (11) and hence it satisfies where (W (x, t), ω) is a solution of Eq. ( 11), then we can use Proposition 1 and formula (18) to show exp In this case, the essential spectrum σ ess lies in the open unit disc D and therefore it cannot contribute to any linear instability of the zero solution of Eq. (19).
To illustrate the usefulness of formula (21), in Fig. 5 (a) we plot the essential spectrum for the periodic solution shown in Fig. 1.In Fig. 5 (b) we show the Floquet multipliers of the same periodic solution, where we have found the solution and its stability in the conventional way, of discretizing the domain and finding a periodic solution of a large set of coupled ordinary differential equations.In panel (b) we see several real Floquet multipliers that do not appear in panel (a); these are presumably part of the discrete spectrum.Note that calculating the discrete spectrum by the method of [39, Section 4] is numerically difficult, so we do not do that here.
F. Formula for firing rates One quantity of interest in a network of model neurons such as (1) is their firing rate.The firing rate of the kth neuron is defined by where the angled brackets • T indicate a long-time average.In the case of large N , we can also consider the average firing rate where x k = 2πk/N is the spatial positions of the kth neuron and the averaging takes place over all neurons in the (π/ √ N )-vicinity of the point x ∈ [0, 2π].Note that while the individual firing rates f k are usually randomly distributed due to the randomness of the excitability parameters η k , the average firing rate f (x) converges to a continuous (and even smooth) function for N → ∞.Moreover, the exact prediction of the limit function f (x) can be given, using only the corresponding solution z(x, t) of Eq. (4).To show this, we write Eq. ( 1) as We recall that in deriving ( 4) from ( 1) we introduce a probability distribution ρ(θ, x, η, t) which satisfies a continuity equation [25,27,41].At a given time t, ρ(θ, x, η, t)dθdηdx is the probability that a neuron with a position in [x, x + dx] and intrinsic drive in [η, η + dη] has its phase in [θ, θ + dθ].Moreover, in the case of the Lorentzian distribution of parameters η k , the probability distribution ρ(θ, x, η, t) satisfies the relations ∞ −∞ dη 2π 0 ηρ(θ, x, η, t)e iθ dθ = (η 0 + iγ)z(x, t), (25) which are obtained by a standard contour integration in the complex plane [44].The relation (24) has been already used to calculate the continuum limit analog of ( 2) Inserting this instead of I k (t) into Eq.( 23) and replacing the time and index averaging in (22) with the corresponding averaging over the probability denisty, we obtain The two inner integrals in the above formula can be simplified using the relations ( 24), ( 25) and the standard normalization condition for ρ(θ, x, η, t).Thus we obtain Moreover, if z = a(x, t) is a T -periodic solution of Eq. ( 4), then the long-time average is the same as an average over one period.So, in the periodic case, the continuum limit average firing rate equals (Note that with simple time rescaling formula ( 26) can be rewritten in terms of a 2π-periodic solution of the complex Riccati equation ( 7), u(x, t) = a(x, t/ω).)The expression ( 26) is different from the firing rate expression given in Sec.I, but both are equally valid.
Results for a pattern like that shown in Fig. 1 are given in Fig. 6, where we show both f (x) (from ( 26)) and values extracted from a long simulation of Eq. ( 1).The agreement is very good.

III. OTHER MODELS
We now demonstrate how the approach presented above can be applied to various other neural models.
FIG. 6.Average firing rate for a pattern like that shown in Fig. 1.The curve shows f (x) as given by ( 26).The dots show values measured from direct simulations of Eq. ( 1).For the discrete simulation, N = 2 14 neurons were used and the average frequency profile, {fj}, j = 1, 2 . . .N , was convolved with a spatial Gaussian filter of standard deviation 0.01 before plotting.For clarity, not all points are shown.Other parameters: A = −5, η0 = −0.9,κ = 1, γ = 0.01.

A. Delays
Delays in neural systems are ubiquitous due to the finite velocity at which action potentials propagate as well as to both dendritic and synaptic processing [3,12,14,51].Here we assume that all I j (t) are delayed by a fixed amount τ , i.e. we have Eq.( 1) but we replace (2) by The mean field equation is now We can write this equation in the same form as Eq. ( 7) but now and the corresponding self-consistency equation is also time-delayed We expand W (x, t) as in (14), and given W (x, t), we find the relevant 2π-periodic solution of Eq. ( 7) as above.The only difference comes in evaluating the projections (15) and ( 16).Instead of using H n (U (x, t)) in the scalar products, we need to use H n (U (x, t − ωτ )).For a periodic solution, the maximum and minimum over one period of this quantity is plotted.The black horizontal line corresponds to the steady state which is stable at τ = 2.5.Circles: measured from direct simulations of Eq. ( 28).Other parameters: Since U (x, t) is 2π-periodic in time, we can evaluate it at any time using just its values for t ∈ Note that this approach would also be applicable if one had a distribution of delays [30,35] or even statedependent delays [21].
As an example, we show in Fig. 7 the results of varying the delay τ on a solution of the form shown in Fig. 3. Increasing τ leads to the destruction of the periodic solution in an apparent supercritical Hopf bifurcation.

B. Two populations
Neurons fall into two major categories: excitatory and inhibitory.A model consisting of a single population with a coupling function of the form (3) is often used as an approximation to a two-population model [17].Here we consider a two population model which supports a travelling wave.The mean field equations are where u(x, t) is the complex-valued order parameter for the excitatory population and v(x, t) is that for the inhibitory population.The non-negative connectivity kernel between and within populations is the same: and there are four connection strengths within and between populations: w ee , w ei , w ie and w ii .Similar models have been studied in [6,48].
For some parameter values, such a system supports a travelling wave with a constant profile.Such a wave can be found very efficiently using the techniques discussed here, and that was done for a travelling chimera in [43].However, here we consider a slightly different case: that where the mean drive to the excitatory population, η u , is spatially modulated.We thus write For small | | the travelling wave persists, but not with a constant profile.An example is shown in Fig. 8.Note that such a solution is periodic in time.
By rescaling time we can write Eqs. ( 31)-(32) as where ũ(x, t) ≡ u(x, t/ω), ṽ(x, t) ≡ v(x, t/ω), and In the same way as above, we can derive self-consistency equations of the form (11) for W u (x, t) and W v (x, t).Doing so, we obtain a system of two coupled equations One difference between this model and the ones studied above is that the solution cannot be shifted by a constant amount in x to ensure that it is always even (or odd) about a particular point in the domain.Thus we need to write These equations contain 6(2F + 1) unknowns and we find them (and ω) in the same way as above, by projecting the self-consistency equations for W u (x, t) and W v (x, t) onto the different spatio-temporal Fourier modes to obtain equations similar to ( 15)-( 16).31)- (32).Other parameters are as in Fig. 8.
The results of varying the heterogeneity strength are shown in Fig. 9. Increasing heterogeneity decreases the period of oscillation, and eventually the travelling wave appears to be destroyed in a saddle-node bifurcation.
We conclude this section by noting that for some parameter values the model ( 31)- (32) can show periodic solutions which do not travel, like those shown in Sec.II.Likewise, the model in Sec.I can support travelling waves for κ = 2.

C. Winfree oscillators
One of the first models of interacting oscillators studied is the Winfree model [2,19,29,46].We consider a spatially-extended network of Winfree oscillators whose dynamics are given by where K jk = K(2π|j − k|/N ) for some 2π-periodic coupling function K, Q(θ) = − sin θ/ √ π and P (θ) = (2/3)(1 + cos θ) 2 is a pulsatile function with its peak at θ = 0.The ω j are randomly chosen from a Lorentzian with centre ω 0 and width ∆ and is the overall coupling strength.
In the limit N → ∞, using the Ott/Antonsen ansatz, one finds that the network is described by the equation [28] where the integral operator K is again defined by (5) and A typical periodic solution of Eq. ( 39) for the choice is shown in Fig. 10.If we rescale time by the frequency of periodic solution ω > 0, defining u(x, t) = z(x, t/ω), and denote then Eq. ( 39) can be recast as a complex Riccati equation From the Proposition 2 in [43], it follows that for every ω, ∆ > 0, ω 0 ∈ R and for every real-valued, 2π-periodic in t function W (x, t), Eq. ( 40) has a 2π-periodic solution lying in the open unit disc D. Denoting the corresponding solution operator we easily obtain a self-consistency equation for periodic solutions of Eq. ( 40) Since the Poincaré map of Eq. ( 40) coincides with the Möbius transformation, we can again use the calculation scheme of Sec.II A to find the value of operator U. Thus we can solve Eq. ( 41) numerically for the real-valued field W (x, t) and frequency ω.
Following the periodic solution shown in Fig. 10 as is varied we obtain Fig. 11.As shown by the circles (from direct simulations of Eq. ( 39)) this solution is not always stable.As is decreased the solution loses stability to a uniformly travelling wave, and the period of this wave is not plotted.

IV. DISCUSSION
We considered time-periodic solutions of the Eqs.( 4)-( 5), which exactly describe the asymptotic dynamics of the network (1) in the limit of N → ∞.At every point in space, Eq. ( 4) is a Riccati equation and we used this to derive a self-consistency equation that every periodic solution of Eq. (4) must satisfy.The Poincaré map of the Riccati equation is a Möbius map and we can determine this map at every point in space using just three numerical solutions of the Riccati equation.Knowing the Möbius map enables us to numerically solve the selfconsistency equation in a computationally efficient manner.We showed the results of numerically continuing several types of periodic solutions as a parameter was varied.We derived equations governing the stability of such periodic solutions, but solving these equations is numerically challenging.We also derived the expression for the mean firing rate of neurons in a network in terms of the quantities already calculated in the self-consistency equation.We finished in Sec.III by demonstrating the application of our approach to several other models involving delays, two populations of neurons, and a network of Winfree oscillators.
Our approach relies critically on the mathematical form of the continuum-level equations (they can be written as a Riccati equation) which are derived using the Ott/Antonsen ansatz, valid only for phase oscillators whose dynamics and coupling involve sinusoidal functions of phases or phase differences.Other systems for which our approach should be applicable include twodimensional networks which support moving or "breathing" solutions [5]; however the coupling function would have to be of the form such that the integral equivalent to (5) could be written exactly using a small set of spatial basis functions.Another application would be to any system which is periodically forced in time and responds in a periodic way [50,53,54].

Let us consider a complex Riccati equation
with 2π-periodic complex-valued coefficients c 0 (t), c 1 (t) and c 2 (t).In this section we prove a statement which is a modified version of Proposition 3 from [43].
Proposition 2 Suppose that there is c * > 0 such that for all z ∈ D and 0 ≤ t ≤ 2π, and that z = −1 is not a fixed point of Eq. ( 42).Then, the Poincaré map of Eq. ( 42) is described by a hyperbolic or loxodromic Möbius transformation.Moreover, the stable fixed point of this map lies in the open unit disc D, while the unstable fixed point lies in the complementary domain Ĉ\D, where Ĉ = C ∪ {∞} is the extended complex plane.For Eq. ( 42) this means that it has exactly one stable 2πperiodic solution and this solution satisfies |z(t)| < 1 for all 0 ≤ t ≤ 2π.
Proof: The fact that the Poincaré map of Eq. ( 42) is described by a Möbius transformation was shown elsewhere [10,59].So we only need to show that such a Möbius transformation M(z) maps all points of the unit circle |z| = 1 into the open unit disc D (but not on its boundary).Then we can repeat the arguments of the proof of Proposition 2 from [43].
Remark 1 If the conditions of Proposition 2 are fulfilled, then the stable solution of Eq. ( 42) is also asymptotically stable.This result follows from the properties of stable fixed points of hyperbolic and loxodromic Möbius transformations.
Remark 2 For every fixed x the equation ( 7) from the main text is equivalent to Eq. ( 42) with (iv) Choosing from the roots z + and z − the one that lies in the unit disc D, one obtains the initial condition that determines the periodic solution of interest.The latter can be computed by solving Eq. (42) with this initial condition.
(v) Sometime it may happen that the Poincaré map M(z) is strongly contracting so that where the value 10 −8 is chosen through experience.In this case, the calculations in steps (ii) and (iii) become inaccurate.Then the initial condition of the periodic solution of interest is approximately given by the average (w 1 + w 2 + w 3 )/3.
The above steps (i)-(v) can be understood as a constructive definition of the solution operator of Eq. ( 42), which for every admissible choise of 2π-periodic coefficients c 1 (t), c 2 (t) and c 3 (t) yields the corresponding stable 2π-periodic solution of Eq. (42).More detailed justification of this definition can be found in [43,Section 4].The algorithm in Sec.II A consists of applying the procedure above at every point on the spatial grid (in parallel).

FIG. 5 .
FIG. 5. (a)The essential spectrum given by(21) for the periodic solution shown in Fig.1.(b) Floquet multipliers of the same periodic solution found using the technique explained at the start of Sec.II.For both calculations the spatial domain has been discretized using 512 evenly spaced points.

FIG. 7 .
FIG.7.The vertical axis relates to averaging W (x, t) over x.For a periodic solution, the maximum and minimum over one period of this quantity is plotted.The black horizontal line corresponds to the steady state which is stable at τ = 2.5.Circles: measured from direct simulations of Eq. (28).Other parameters: A = −5, η0 = −2, κ = 1, γ = 0.1.