Biological Cybernetics (2023) 117:259–274 https://doi.org/10.1007/s00422-023-00969-6 ORIG INAL ART ICLE Periodic solutions in next generation neural field models Carlo R. Laing1 ·Oleh E. Omel’chenko2 Received: 24 January 2023 / Accepted: 12 July 2023 / Published online: 3 August 2023 © The Author(s) 2023 Abstract 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. Keywords Neural field model · Riccati equation · Theta neuron · Ott/Antonsen · Self-consistency 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 fir- ing of many neurons. Such oscillations have relevance to rhythmic movement (Lindén et al. 2022), epilepsy (Jiruska et al. 2013; Netoff and Schiff 2002), schizophrenia (Uhlhaas and Singer 2010) neural communication (Reyner-Parra and Huguet 2022) and EEG/MEGmodelling (Byrne et al. 2022), among others. In different networks, oscillations may arise from mechanisms such as synaptic delays (Devalle et al. 2018), the interaction of excitatory and inhibitory popu- lations (Börgers and Kopell 2005; Schmidt and Avitabile 2020), having sufficient connectivity in an inhibitory network (di Volo and Torcini 2018), or the finite width of synaptic pulses emitted by neurons (Ratas and Pyragas 2016). The Communicated by Benjamin Lindner. B Oleh E. Omel’chenko oleh.omelchenko@uni-potsdam.de Carlo R. Laing c.r.laing@massey.ac.nz 1 School of Mathematical and Computational Sciences, Massey University, Private Bag 102-904 NSMC, Auckland, New Zealand 2 Institute of Physics and Astronomy, University of Potsdam, Karl-Liebknecht-Str. 24/25, 14476 Potsdam, Germany 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 (Ermentrout and Kopell 1986), Winfree oscillators (Ariaratnam and Stro- gatz 2001) and quadratic integrate-and-fire (QIF) neurons (Latham et al. 2000) are some of the simplest. These three types of model neurons have the advantage that their mathe- matical form often allows infinite networks of such neurons to be analysed exactly using the Ott-Antonsen method (Ott and Antonsen 2008, 2009). 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 con- nection with modelling head direction (Zhang 1996) and working memory (Funahashi et al. 1989; Wimmer et al. 2014), for example, and been studied by others (Esnaola- Acebes et al. 2017; Laing and Chow 2001; Kilpatrick and Ermentrout 2013). We consider a network of theta neu- rons. The theta neuron is a minimal model for a neuron which undergoes a saddle-node-on-invariant-circle (SNIC) bifurcation as a parameter is varied (Ermentrout and Kopell 1986). The theta neuron is exactly equivalent to a quadratic integrate-and-fire (QIF) neuron, under the assumption of infi- nite firing threshold and reset values (Montbrió et al. 2015; Avitabile et al. 2022; Devalle et al. 2017). The coupling in the network is nonlocal synaptic coupling, implemented using a spatial convolution with a translationally-invariant coupling kernel. 123 http://crossmark.crossref.org/dialog/?doi=10.1007/s00422-023-00969-6&domain=pdf http://orcid.org/0000-0002-6086-2978 http://orcid.org/0000-0003-0526-1878 260 Biological Cybernetics (2023) 117:259–274 We studied such a model in the past (Omel’chenko and Laing 2022), concentrating on describing spatially-uniform states and also stationary “bump” states in which there is a spatially-localised group of active neurons while the remain- der of the network is quiescent. We determined the existence and stability of such states and also found regions of parame- ter 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 (Omel’chenko 2023, 2022) which is significantly faster than the standard approach. This technique was successfully applied to describe trav- eling and breathing chimera states in nonlocally coupled phase oscillators (Omel’chenko 2023, 2022). 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 for- mulas, 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 thatWinfree oscillators are also treatable by the proposed technique. The structure of the paper is as follows. In Sect. 1 we present the discrete network model and its continuum-level description. The bulk of the paper is in Sect. 2,wherewe show how to describe periodic states in a self-consistent way, and show numerical results from implementing our algorithms. Other models are considered in Sect. 3 and we end in Sect. 4. The Appendix contains useful results regarding the complex Riccati equation. 1 Theta neuron networkmodel Themodelwe consider first is that inOmel’chenko andLaing (2022), Laing (2014a) and Laing and Omel’chenko (2020), which we briefly present here. The discrete network consists of N synaptically coupled theta neurons described by dθ j dt = 1 − cos θ j + (1 + cos θ j )(η j + κ I j ), j = 1, . . . , N , (1) 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 j th neuron is κ I j where I j (t) = 2π N N∑ k=1 K jk Pn(θk(t)) (2) where Pn(θ) = an(1 − cos θ)n is a pulsatile function with a maximum at θ = π (when the neuron fires) and an is chosen according to the normalization condition ∫ 2π 0 Pn(θ)dθ = 2π. Increasing n makes Pn(θ) “sharper” and more pulse-like. The excitability parameters η j are chosen from a Lorentzian distribution with mean η0 and width γ > 0 g(η) = γ π 1 (η − η0) 2 + γ 2 . 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 K (x) = 1 2π (1 + A cos x) (3) 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 level — periodic in time; see Fig. 3 in Omel’chenko and Laing (2022). 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 betweenneurons onopposite 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 (Pinto and Ermentrout 2001; Esnaola-Acebes et al. 2017). Using the Ott/Antonsen ansatz (Ott and Antonsen 2008, 2009), one can show that in the limit N → ∞, the long-term dynamics of the network (1) can be described by ∂z ∂t = (iη0 − γ )(1 + z)2 − i(1 − z)2 2 +κ i(1 + z)2 2 KHn(z), (4) 123 Biological Cybernetics (2023) 117:259–274 261 where (Kϕ)(x) = ∫ 2π 0 K (x − y)ϕ(y)dy (5) is the convolution of K and ϕ and Hn(z) = an ⎡ ⎣C0 + n∑ q=1 Cq ( zq + zq ) ⎤ ⎦ , where Cq = n∑ k=0 k∑ m=0 δk−2m,q(−1)kn! 2k(n − k)!m!(k − m)! . For our computations we set n = 2, so that H2(z) = (2/3)[3/2 − (z + z̄) + (z2 + z̄2)/4] 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 (Pietras 2023). Equation (4) is an integro-differential equation for a complex-valued function z(x, t), where x ∈ [0, 2π ] is posi- tion on the ring. z(x, t) is a local order parameter and can be thought of as the average of eiθ for neurons in a small neigh- bourhood 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 θ (Laing 2014a). 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 (Laing 2015; Montbrió et al. 2015) that π−1Re W is the flux through θ = π or the instantaneous firing rate of neurons at position x and time t . Similarly, if Vj = tan(θ j/2) is the voltage of the j th QIF neuron, the mean voltage at position x and time t is given by ImW . 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 D = {z ∈ C : |z| < 1} is the unit disc in the complex plane. Equations of the form (4)–(5) are sometimes referred to as next generation neural field models (Byrne et al. 2019; Coombes and Byrne 2019) as they have the form of a neural field model (an integro-differential equation for a macro- scopic quantity such as “activity” (Laing et al. 2002; Amari 1977)) but are derived exactly from a network like (1), rather than being of a phenomenological nature. Fig. 1 A typical periodic solution of Eq. (4). a arg (z(x, t)). b |z(x, t)|. c A realization of this solution in a network of N = 4096 theta neurons described by Eq. (1). θ j is shown in color. Parameters: A = −5, η0 = −0.7, κ = 1, γ = 0.01 2 Periodic states In this paper we focus on states with periodically oscillat- ing macroscopic dynamics. For the mean field Eq. (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). Figure 1c 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 (Laing and Troy 2003). This is another example of next generation neural field mod- els showing more complex time-dependent behaviour than classical ones (Laing and Omel’chenko 2020). 123 262 Biological Cybernetics (2023) 117:259–274 A straight-forward way to study a periodic orbit like that in Fig. 1 would be discretise Eq. (4) on a spatially-uniform grid and approximate the convolution using matrix/vector multiplication or otherwise, resulting in a large set of cou- pled ordinary differential equations. The periodic solution of these could then be studied using standard techniques (Laing 2014a), 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 Omel’chenko (2023, 2022), which allows us to perform the same calcula- tions with only∼ N operations. The main ingredients of this method are explained in Sects. 2.1 and 2.2. They include the description of the properties of complex Riccati equation and the derivation of the self-consistency equation for periodic solutions of Eq. (4). In Sect. 2.3 we explain how the self- consistency equation can be solved in the case of coupling function (3). Then in Sect. 2.4 we report some numerical results obtained with our method. In addition, in Sect. 2.5 we perform a rigorous linear stability analysis of periodic solutions of Eq. (4) by considering the spectrum of the corre- spondingmonodromy operator. Finally, in Sect. 2.6, we show how the mean field Eq. (4) can be used to predict the average firing rate distribution in the neural network (1). 2.1 Periodic complex Riccati equation andMöbius transformation By the time rescaling u(x, t) = z(x, t/ω) we can rewrite Eq. (4) in the form ω ∂u ∂t = (iη0 − γ )(1 + u)2 − i(1 − u)2 2 +κ i(1 + u)2 2 KHn(u), (6) such that the above T -periodic solution of Eq. (4) corre- sponds to a 2π -periodic solution ofEq. (6). Then, dividing (6) byω and reordering the termswe can see that Eq. (6) is equiv- alent to a complex Riccati equation ∂u ∂t = i ( W (x, t) + ζ − 1 2ω ) + 2i ( W (x, t) + ζ + 1 2ω ) u +i ( W (x, t) + ζ − 1 2ω ) u2, (7) with the 2π -periodic in t coefficient W (x, t) = κ 2ω KHn(u) (8) and ζ = η0 + iγ 2ω . (9) In Omel’chenko (2023) 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 func- tion W (x, t), parameters ζ ∈ Cup = {z ∈ C : Im z > 0} and ω > 0, for every fixed x ∈ [0, 2π ] Eq. (7) has a unique sta- ble 2π -periodic solutionU (x, t) that lies entirely in the open unit discD. Denoting the corresponding solution operator by U : Cper([0, 2π ];R) × Cup × (0,∞) → Cper([0, 2π ];D), we can write the 2π -periodic solution of interest as U (x, t) = U ( W (x, t), η0 + iγ 2ω ,ω ) . (10) Note thatCper([0, 2π ];R) here denotes the space of all real- valued continuous 2π -periodic functions, while the notation Cper([0, 2π ];D) stands for the space of all complex con- tinuous 2π -periodic functions with values in the open unit discD. Importantly, the variable x appears in formula (10) as a parameter so that the function W (x, ·) ∈ Cper([0, 2π ];R) with a fixed x is considered as the first argument of the oper- ator U . As for the operator U , although it is not explicitly given, its value can be calculated without resource-demanding iter- ative methods by solving exactly four initial value problems for Eq. (7). The rationale for this approach can be found in (Omel’chenko 2023, Section 4) and is repeated for com- pleteness 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) andW (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 functionsw j (t), we calculate the corresponding functions u j (t) by performing the following four steps. (i) We (somewhat arbitrarily) choose three initial condi- tions u1j (0) = −0.95, u2j (0) = 0, u3j (0) = 0.95 and solve Eq. (7) with these initial conditions to obtain solutions ukj (t), j = 1, . . . N ; k = 1, 2, 3. In Omel’chenko (2023) (see also Proposition 2 in the Appendix) it is shown that for γ > 0 these solutions lie in the open unit disc D. (ii) Since at each point in space the Poincaré map of Eq. (7)with 2π -periodic coefficients coincideswith aMöbius mapM(u) (seeOmel’chenko 2023 for detail), we can use the relations ukj (2π) = M j (ukj (0)), k = 1, 2, 3, to reconstruct these maps M j , j = 1, . . . N . The corresponding formulas are given in Omel’chenko (2023, Section 4). (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 equa- tion, therefore in general it has two solutions in the complex 123 Biological Cybernetics (2023) 117:259–274 263 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 self- consistency equation, the solution of which allows us to determine a 2π -periodic solution of Eq. (6). 2.2 Self-consistency equation Supposing that Eq. (6) has a 2π -periodic solution, then using formula (8) we can calculate the corresponding func- tion 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 W (x, t) = κ 2ω KHn ( U ( W (x, t), η0 + iγ 2ω ,ω )) , (11) obtained by inserting (10) into (8). In the following we con- sider 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 (1984), Strogatz (2000) and Shima and Kuramoto (2004), but such approaches have always focused on steady states, whereas we consider peri- odic solutions here. Note that from a computational point of view, the self- consistency Eq. (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 har- monics (see Sect. 2.3). 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 Sect. 2.5. Proposition 1 Let the pair (W (x, t), ω) be a solution of the self-consistency Eq. (11) and let U (x, t) be defined by (10). Then ∣∣∣∣exp (∫ 2π 0 M(x, t)dt )∣∣∣∣ < 1 where M(x, t) = i ω [ (κKHn(U ) + η0 + iγ )(1 +U (x, t)) + 1 −U (x, t)] . (12) Proof For every fixed x ∈ [0, 2π ], the function U (x, t) yields a stable 2π -periodic solution of the complex Riccati Eq. (7). The linearization of Eq. (7) around this solution reads dv dt = M(x, t)v, where M(x, t) = 2i ( W (x, t) + ζ + 1 2ω ) +2i ( W (x, t) + ζ − 1 2ω ) U (x, t). 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 sta- ble 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. 2.3 Numerical implementation Equation (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 ∫ 2π 0 dx ∫ 2π 0 W (x, t) sin(2t)dt = 0. (13) In the following we focus on the case of the cosine cou- pling (3). It is straight-forward to show that (KHn)(x) can be written as a linear combination of 1, cos x and sin x . How- ever, 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, tak- ing into account that the functionW (x, t) is real, we seek an approximate solution of the system (11), (13) using a Fourier- Galerkin ansatz W (x, t) = 2F∑ m=0 (vm + wm sin x)ψm(t) (14) 123 264 Biological Cybernetics (2023) 117:259–274 where vm andwm are real coefficients andψm(t) are trigono- metric basis functions ψ0(t) = 1, ψm(t) = √ 2 cos(nt), if m = 2n with n ∈ N, ψm(t) = √ 2 sin(nt), if m = 2n − 1 with n ∈ N. 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. Using the scalar product 〈u, v〉 = 1 (2π)2 ∫ 2π 0 dx ∫ 2π 0 u(x, t)v(x, t)dt we project Eq. (11) on different spatio-temporal Fourier modes to obtain the system vm = κ 2ω 〈 Hn ( U ( W (x, t), η0 + iγ 2ω ,ω )) , ψm(t) 〉 , (15) wm = κA 2ω 〈 Hn ( U ( W (x, t), η0 + iγ 2ω ,ω )) , ψm(t) sin x〉 , (16) for m = 0, 1, . . . , 2F . Equations (15) and (16), together with (13), are a set of 2(2F +1)+1 = 4F +3 equations for the 4F + 3 unknowns v0, v1, . . . , v2F , w0, w1, . . . , w2F , ω, which must be solved simultaneously. We solve them using Newton’s method and find convergence within 3 or 4 itera- tions. In simple terms, suppose we have somewhat accurate esti- mates of v0, v1, . . . , v2F , w0, w1, . . . , w2F , ω. These can be inserted into (14) to calculate the functionW (x, t). Then one can calculate a periodic solution of Eq. (7) with the speci- fied W (x, t) by formula (10) and finally calculate Hn(U) and insert this into (15) and (16) and perform the projec- tions. We want the difference between the new values of v0, v1, . . . , v2F , w0, w1, . . . , w2F 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 stan- dard way (Laing 2014b). Note that these calculations involve discretising the spatial domain with N points. However, the number of unknowns (4F + 3) is significantly less than N . Fig. 2 Period of the type of solution shown in Fig. 1 as a function of η0 (solid curve). The circles show valuesmeasured from direct simulations of Eq. (4). Other parameters: A = −5, κ = 1, γ = 0.01 2.4 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, andbecomeunstable to a wandering pattern at η0 is increased through approxi- mately −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 diagramswhere parameter ranges of two patterns with differ- ent symmetries are separated by heteroclinic or homoclinic bifurcations were found for non-locally coupled Kuramoto- 123 Biological Cybernetics (2023) 117:259–274 265 Fig. 3 Another periodic solution of Eq. (4). a arg (z(x, t)). b |z(x, t)|. Parameters: A = −5, η0 = −2.5, κ = 1, γ = 0.01 Fig. 4 Period of the type of solution shown in Fig. 3 as a function of η0 (solid curve). The circles show valuesmeasured from direct simulations of Eq. (4). Other parameters: A = −5, κ = 1, γ = 0.01 type phase oscillators (Omel’chenko 2020) and seem to be a general mechanism which, however, needs additional inves- tigation. 2.5 Stability of breathing bumps Given a T -periodic solution a(x, t) of Eq. (4), we can per- form its linear stability analysis, using the approach proposed in Omel’chenko (2022). Before doing this, we write Hn(z) = anC0 + 2Re[Dn(z)] where Dn(z) = an n∑ q=1 Cqz q , to emphasise that Hn(z) is always real. Now, we insert the ansatz z(x, t) = a(x, t)+v(x, t) intoEq. (4) and linearize the resulting equationwith respect to small perturbations v(x, t). Thus, we obtain a linear integro-differential equation ∂v ∂t = μ(x, t)v + κ i(1 + a(x, t))2 2 K ( D′ n(a)v + D′ n(a)v ) , (17) where μ(x, t) = [i(η0 + κKHn(a)) − γ ](1 + a(x, t)) +i(1 − a(x, t)) (18) and D′ n(z) = d dz Dn(z) = an n∑ q=1 qCqz q−1. Note thatEq. (17) coincideswithEq. (5.1) fromOmel’chenko and Laing (2022), except that the coefficients a(x, t) and μ(x, t) are now time-dependent. Since Eq. (17) contains the complex-conjugated term v, it is convenient to consider this equation along with its complex-conjugate ∂v ∂t = μ(x, t)v − κ i(1 + a(x, t))2 2 K ( D′ n(a)v + D′ n(a)v ) . This pair of equations can be written in the operator form dV dt = A(t)V + B(t)V , (19) where V (t) = (v1(t), v2(t))T is a function with values in Cper([0, 2π ];C2), and A(t)V = ( μ(x, t) 0 0 μ(x, t) ) ( v1 v2 ) , 123 266 Biological Cybernetics (2023) 117:259–274 and B(t)V = iκ 2 ( (1 + a(x, t))2 0 0 −(1 + a(x, t))2 ) × ( K ( D′ n(a)v1 ) K ( D′ n(a)v2 ) ) . (20) For every fixed t the operatorsA(t) and B(t) are linear oper- ators from Cper([0, 2π ];C2) 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 Daleckiı̆ and Kreı̆n (2002, 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 E(t) = exp (∫ t 0 (A(t ′) + B(t ′))dt ′ ) . The analysis of Eq. (19) in the case when A(t) is a matrix multiplication operator and B(t) is an integral operator sim- ilar to (20) has been performed in (Omel’chenko 2022, Section 4). Repeating the same arguments we can demon- strate 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 σess = { exp (∫ T 0 μ(x, t)dt ) : x ∈ [0, 2π ] } ∪ {c.c.} (21) (ii) the discrete spectrum σdisc that consists of finitely many isolated eigenvalues λ, which can be found using a characteristic integral equation, as explained in (Omel’chenko 2022, Section 4). Note that if a(x, t) is obtained by solving the self- consistency Eq. (11) and hence it satisfies a(x, t/ω) = U (x, t) = U ( W (x, t), η0 + iγ 2ω ,ω ) , where (W (x, t), ω) is a solution of Eq. (11), then we can use Proposition 1 and formula (18) to show ∣∣∣∣exp (∫ T 0 μ(x, t)dt )∣∣∣∣ < 1 for all x ∈ [0, 2π ]. 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 Sect. 2. For both calculations the spatial domain has been discretized using 512 evenly spaced points In this case, the essential spectrum σess lies in the open unit disc D and therefore it cannot contribute to any linear insta- bility of the zero solution of Eq. (19). To illustrate the usefulness of formula (21), in Fig. 5a we plot the essential spectrum for the periodic solution shown in Fig. 1. In Fig. 5b 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 ordi- nary 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 calcu- lating the discrete spectrum by the method of Omel’chenko (2022, Section 4) is numerically difficult, so we do not do that here. 2.6 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 fk = 1 2π 〈 dθk dt 〉 T , 123 Biological Cybernetics (2023) 117:259–274 267 where the angled brackets 〈·〉T indicate a long-time average. In the case of large N , we can also consider the average firing rate f (x) = 1 #{k : |xk − x | < π/ √ N } ∑ |xk−x |<π/ √ N fk, (22) where xk = 2πk/N is the spatial positions of the kth neu- ron 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 fk 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 predic- tion 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 dθk dt = Re { 1 − eiθk + (ηk + κ Ik)(1 + eiθk ) } . (23) We recall that in deriving (4) from (1) we introduce a probability distribution ρ(θ, x, η, t) which satisfies a conti- nuity equation (Laing 2015; Omel’chenko et al. 2014; Laing 2014a). At a given time t , ρ(θ, x, η, t)dθdηdx is the proba- bility that a neuronwith 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)eiθdθ = z(x, t), (24) ∫ ∞ −∞ dη ∫ 2π 0 ηρ(θ, x, η, t)eiθdθ = (η0 + iγ )z(x, t), (25) which are obtained by a standard contour integration in the complex plane (Ott and Antonsen 2008). The relation (24) has been already used to calculate the continuum limit analog of (2) I (x, t) = ∫ 2π 0 K (x − y)Hn(z(y, t))dy = KHn(z). Inserting this instead of Ik(t) into Eq. (23) and replacing the time and index averaging in (22) with the corresponding averaging over the probability density, we obtain f (x) = lim τ→∞ 1 2πτ ∫ τ 0 ∫ ∞ −∞ ∫ 2π 0 Re { 1 − eiθ + (η + κ I (x, t))(1 + eiθ ) } ρ(θ, x, η, t)dθ dη dt . 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 showvaluesmeasured from direct simulations of Eq. (1). For the discrete simulation, N = 214 neu- rons were used and the average frequency profile, { f j }, 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 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 f (x) = lim τ→∞ 1 2πτ ∫ τ 0 Re {1 − z(x, t) +(η0 + iγ + κ I (x, t))(1 + z(x, t))} dt . 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 f (x) = 1 2πT ∫ T 0 Re {1 − a(x, t) +(η0 + iγ + κ I (x, t))(1 + a(x, t))} dt . (26) (Note that with simple time rescaling formula (26) can be rewritten in terms of a 2π -periodic solution of the complex Riccati Eq. (7), u(x, t) = a(x, t/ω).) The expression (26) is different from the firing rate expression given in Sect. 1, 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. 3 Other models We now demonstrate how the approach presented above can be applied to various other neural models. 123 268 Biological Cybernetics (2023) 117:259–274 3.1 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 (Roxin et al. 2005; Coombes and Laing 2009; Atay and Hutt 2006; Devalle et al. 2018). 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 I j (t) = 2π N N∑ k=1 K jk Pn(θk(t − τ)). (27) The mean field equation is now ∂z ∂t = (iη0 − γ )(1 + z)2 − i(1 − z)2 2 +κ i(1 + z)2 2 KHn(z(x, t − τ)). (28) We can write this equation in the same form as Eq. (7) but now W (x, t) = κ 2ω KHn(u(x, t − ωτ)), (29) and the corresponding self-consistency equation is also time- delayed W (x, t) = κ 2ω KHn ( U ( W (x, t − ωτ), η0 + iγ 2ω ,ω )) . We expandW (x, t) as in (14), and givenW (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 Hn(U (x, t)) in the scalar products, we need to use Hn(U (x, t − ωτ)). SinceU (x, t) is 2π -periodic in time, we can evaluate it at any time using just its values for t ∈ [0, 2π ]. Specifically, U (x, t − ωτ) = { U (x, 2π + t − ωτ), 0 ≤ t ≤ ωτ, U (x, t − ωτ), ωτ < t ≤ 2π. (30) Note that this approach would also be applicable if one had a distribution of delays (Lee et al. 2009; Laing and Longtin 2003) or even state-dependent delays (Keane et al. 2019). 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. 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 simula- tions of Eq. (28). Other parameters: A = −5, η0 = −2, κ = 1, γ = 0.1 3.2 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 approx- imation to a two-population model (Esnaola-Acebes et al. 2017). Here we consider a two population model which sup- ports a travelling wave. The mean field equations are ∂u ∂t = (iηu − γ )(1 + u)2 − i(1 − u)2 2 + i(1 + u)2 2 [weeKHn(u) − weiKHn(v)] , (31) ∂v ∂t = (iηv − γ )(1 + v)2 − i(1 − v)2 2 + i(1 + v)2 2 [wieKHn(u) − wiiKHn(v)] (32) 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: K (x) = 1 2π (1 + cos x) and there are four connection strengths within and between populations:wee,wei,wie andwii. Similar models have been 123 Biological Cybernetics (2023) 117:259–274 269 Fig. 8 A modulated travelling wave solution of Eqs. (31)–(32). a |u|; b |v|. Other parameters: η0 = 0.1, ηv = 0.1, ε = 0.01, γ = 0.03, wee = 1, wei = 0.7, wie = 0.3, wii = 0.1 studied in Blomquist et al. (2005), Pinto and Ermentrout (2001). For some parameter values, such a system supports a trav- elling 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 Omel’chenko (2023). However, here we consider a slightly different case: that where the mean drive to the excitatory population, ηu , is spatially modulated. We thus write ηu = η0 + ε sin x . For small |ε| the travelling wave persists, but not with a con- stant 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 ∂ ũ ∂t = i ( weeWu − weiWv + ζu − 1 2ω ) +2i ( weeWu − weiWv + ζu + 1 2ω ) ũ +i ( weeWu − weiWv + ζu − 1 2ω ) ũ2, (33) ∂ṽ ∂t = i ( wieWu − wiiWv + ζv − 1 2ω ) +2i ( wieWu − wiiWv + ζv + 1 2ω ) ṽ +i ( wieWu − wiiWv + ζv − 1 2ω ) ṽ2, (34) where ũ(x, t) ≡ u(x, t/ω), ṽ(x, t) ≡ v(x, t/ω), Wu(x, t) = 1 2ω KHn(u), (35) Wv(x, t) = 1 2ω KHn(v) (36) and ζu = ηu + iγ 2ω = η0 + ε sin x + iγ 2ω , ζv = ηv + iγ 2ω . In the same way as above, we can derive self-consistency equations of the form (11) for Wu(x, t) and Wv(x, t). Doing so, we obtain a system of two coupled equations Wu(x, t) = κ 2ω KHn ( U ( weeWu(x, t) − weiWv(x, t), η0 + ε sin x + iγ 2ω ,ω )) , Wu(x, t) = κ 2ω KHn ( U ( wieWu(x, t) − wiiWv(x, t), ηv + iγ 2ω ,ω )) . 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 Wu(x, t) = 2F∑ m=0 (vum + wu m sin x + zum cos x)ψm(t) (37) Wv(x, t) = 2F∑ m=0 (vv m + wv m sin x + zvm cos x)ψm(t) (38) These equations contain 6(2F + 1) unknowns and we find them (and ω) in the same way as above, by projecting the self-consistency equations forWu(x, t) andWv(x, t)onto the different spatio-temporal Fourier modes to obtain equations similar to (15)–(16). 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 param- eter values the model (31)–(32) can show periodic solutions which do not travel, like those shown in Sect. 2. Likewise, the model in Sect. 1 can support travelling waves for κ = 2. 3.3 Winfree oscillators One of the first models of interacting oscillators studied is the Winfree model (Ariaratnam and Strogatz 2001; Laing et al. 123 270 Biological Cybernetics (2023) 117:259–274 Fig. 9 Period, T , of a modulated travelling wave solution of Eqs. (31)– (32) as a function of heterogeneity strength ε. Circles are from direct simulation of Eqs. (31)–(32). Other parameters are as in Fig. 8 2021; Pazó and Montbrió 2014; Gallego et al. 2017). We consider a spatially-extended network of Winfree oscillators whose dynamics are given by dθ j dt = ω j + ε 2πQ(θ j ) N N∑ k=1 K jk P(θk) 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 functionwith 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 (Laing 2017) ∂z ∂t = ε 2 √ π KĤ(z) + (iω0 − �)z − ε 2 √ π z2KĤ(z), (39) where the integral operator K is again defined by (5) and Ĥ(z) = (2/3)[3/2 + z + z̄ + (z2 + z̄2)/4]. A typical periodic solution of Eq. (39) for the choice K (x) = 0.1 + 0.3 cos x, 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 W (x, t) = ε 2 √ πω KĤ(u), Fig. 10 A periodic solution of Eq. (39). a arg (z(x, t)). b |z(x, t)|. Parameters: ω0 = 1, � = 0.1, ε = 2.1 then Eq. (39) can be recast as a complex Riccati equation ∂u ∂t = W (x, t) + i ω0 + i� ω u − W (x, t)u2. (40) From the Proposition 2 in Omel’chenko (2023), it follows that for every ω,� > 0, ω0 ∈ R and for every real-valued, 2π -periodic in t functionW (x, t), Eq. (40) has a 2π -periodic solution lying in the open unit disc D. Denoting the corre- sponding solution operator Û : Cper([0, 2π ];R) × Cup → Cper([0, 2π ];D), weeasily obtain a self-consistency equation for periodic solu- tions of Eq. (40) W (x, t) = ε 2 √ πω KĤ ( Û ( W (x, t), ω0 + i� ω )) . (41) Since the Poincarémap of Eq. (40) coincideswith theMöbius transformation, we can again use the calculation scheme of Sect. 2.1 to find the value of operator Û . 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 variedwe obtain Fig. 11. As shown by the circles (from direct simulations of Eq. (39)) this solution is not always stable. 123 Biological Cybernetics (2023) 117:259–274 271 Fig. 11 Period, T , of periodic solutions of Eq. (39) of the form shown in Fig. 10. Circles are from direct simulations of Eq. (39). Parameters: ω0 = 1, � = 0.1 As ε is decreased the solution loses stability to a uniformly travelling wave, and the period of this wave is not plotted. 4 Discussion We considered time-periodic solutions of the Eqs. (4)–(5), which exactly describe the asymptotic dynamics of the net- work (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 numeri- cally solve the self-consistency equation in a computationally efficient manner. We showed the results of numerically con- tinuing several types of periodic solutions as a parameter was varied. We derived equations governing the stability of such peri- odic 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 quanti- ties already calculated in the self-consistency equation. We finished in Sect. 3 by demonstrating the application of our approach to several other models involving delays, two pop- ulations 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 dif- ferences. Other systems for which our approach should be applicable include two-dimensional networks which support moving or “breathing” solutions (Bataille-Gonzalez et al. 2021); however the coupling functionwould 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 (Segneri et al. 2020; Reyner-Parra and Huguet 2022; Schmidt et al. 2018). Acknowledgements TheworkofO.E.O.was supportedby theDeutsche Forschungsgemeinschaft under Grant No. OM 99/2-2. Author Contributions Conceptualization: CRL, OEO. Formal analy- sis: CRL OEO. Funding acquisition: OEO. Investigation: CRL, OEO. Methodology: CRL, OEO. Software: CRL. Validation: CRL, OEO. Visualization: CRL. Writing – original draft: CRL, OEO. Writing – review & editing: CRL, OEO. Funding Open Access funding enabled and organized by Projekt DEAL. Code availability Code for calculating any of the results presented here is available on reasonable request from the authors. Declarations Conflict of interest The authors have no competing interests to declare that are relevant to the content of this article. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adap- tation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indi- cate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, youwill need to obtain permission directly from the copy- right holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Appendix Let us consider a complex Riccati equation dz dt = c0(t) + c1(t)z + c2(t)z 2 (42) with 2π -periodic complex-valued coefficients c0(t), c1(t) and c2(t). In this section we prove a statement which is a modified version of Proposition 3 fromOmel’chenko (2023). 123 http://creativecommons.org/licenses/by/4.0/ http://creativecommons.org/licenses/by/4.0/ 272 Biological Cybernetics (2023) 117:259–274 Proposition 2 Suppose that there is c∗ > 0 such that Re (c0(t)z + c1(t) + c2(t)z) ≤ −c∗Re (z + 1) (43) 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 transfor- mation. 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 (Campos 1997;Wilczyński 2008). So we only need to show that such a Möbius transformationM(z) maps all points of the unit cir- cle |z| = 1 into the open unit discD (but not on its boundary). Then we can repeat the arguments of the proof of Proposi- tion 2 from Omel’chenko (2023). Suppose the opposite. Then Eq. (42) has a solution z∗(t) such that |z∗(0)| = |z∗(2π)| = 1. Due to the inequality (43) and the Proposition 1 fromOmel’chenko (2023) this solution satisfies |z∗(t)| ≤ 1 for all t ∈ [0, 2π ]. Moreover, since z = −1 is not a fixed point of Eq. (42), we can always choose t∗ ∈ (0, 2π) such that z∗(t) �= −1 for t ∈ [t∗, 2π). Then the Mean Value Theorem implies 0 ≤ |z∗(2π)|2 − |z∗(t∗)|2 = 2π d|z∗|2 dt (t∗∗) (44) for some t∗∗ ∈ (t∗, 2π). On the other hand, due to our assumptions, we have d|z∗|2 dt = 2Re ( c0(t)z∗(t) + c1(t) + c2(t)z∗(t) ) ≤ −2c∗Re (z∗(t) + 1) and therefore d|z∗|2 dt (t∗∗) ≤ −2c∗Re (z∗(t∗∗) + 1) < 0. This is a contradiction to (44), which completes the proof. 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 Eq. (7) from the main text is equivalent to Eq. (42) with c0(t) = c2(t) = i ( W (x, t) + ζ − 1 2ω ) and c1(t) = 2i ( W (x, t) + ζ + 1 2ω ) , where W (x, t) is a real-valued function, ω is a positive con- stant, and ζ = (η0 + iγ )/(2ω) with η0 ∈ R and γ > 0. In this case, we have Re (c0(t)z + c1(t) + c2(t)z) = −γ ω Re (z + 1) and therefore inequality (43) is satisfied. On the other hand, we have c0(t) + c1(t)(−1) + c2(t)(−1)2 = −2i/ω �= 0 and therefore z = −1 is not a fixed point of the corresponding Riccati equation. Thus, all of the conclusions of Proposition 2 hold true for Eq. (7). Remark 3 If the conditions of Proposition 2 are fulfilled, then the stable solution of Eq. (42) can be computed in the fol- lowing way. (i) One solves Eq. (42) on the interval t ∈ (0, 2π ] with three different initial conditions z(0) = zk ∈ D, k = 1, 2, 3, and obtains three solutions Zk(t). Since each zk lies in the open unit disc D this automatically implies |Zk(t)| < 1 for all t ∈ (0, 2π ]. (ii) One denotes wk = Zk(2π). Then, due to the prop- erties of Poincaré map one has wk = M(zk), k = 1, 2, 3, where M(z) is a Möbius transformation representing this map. The above three relations can be used to reconstruct the map M(z), namely M(z) = az + b cz + d where a = det ⎛ ⎝ z1w1 w1 1 z2w2 w2 1 z3w3 w3 1 ⎞ ⎠ , b = det ⎛ ⎝ z1w1 z1 w1 z2w2 z2 w2 z3w3 z3 w3 ⎞ ⎠ , c = det ⎛ ⎝ z1 w1 1 z2 w2 1 z3 w3 1 ⎞ ⎠ , d = det ⎛ ⎝ z1w1 z1 1 z2w2 z2 1 z3w3 z3 1 ⎞ ⎠ . (iii) Once the map M(z) is known, one can find its fixed points by solving the quadratic equation cz2 + dz − az − b = 0. This yields two roots z± = a − d ± √ (a − d)2 + 4bc 2c . 123 Biological Cybernetics (2023) 117:259–274 273 (iv) Choosing from the roots z+ and z− the one that lies in the unit discD, one obtains the initial condition that deter- mines 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é mapM(z) is strongly contracting so that |w1 − w2| + |w3 − w2| < 10−8, 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 (w1 + w2 + w3)/3. The above steps (i)–(v) can be understood as a constructive definition of the solution operator of Eq. (42),which for every admissible choice of 2π -periodic coefficients c1(t), c2(t) and c3(t) yields the corresponding stable 2π -periodic solution of Eq. (42). More detailed justification of this definition can be found in Omel’chenko (2023, Section 4). The algorithm in Sect. 2.1 consists of applying the procedure above at every point on the spatial grid (in parallel). References Amari SI (1977) Dynamics of pattern formation in lateral-inhibition type neural fields. Biol Cybern 27(2):77–87 Ariaratnam JT, Strogatz SH (2001) Phase diagram for the Winfree model of coupled nonlinear oscillators. Phys Rev Lett 86(19):4278 Atay FM, Hutt A (2006) Neural fields with distributed transmission speeds and long-range feedback delays. SIAM J Appl Dyn Syst 5(4):670–698 AvitabileD,DesrochesM,ErmentroutGB (2022)Cross-scale excitabil- ity in networks of quadratic integrate-and-fire neurons. PLoS Comput Biol 18(10):e1010569 Bataille-Gonzalez M, Clerc MG, Omel’chenko OE (2021) Moving spi- ral wave chimeras. Phys Rev E 104:L022203 Blomquist P,Wyller J, Einevoll GT (2005) Localized activity patterns in two-population neuronal networks. Physica D Nonlinear Phenom 206(3):180–212 Börgers C, Kopell N (2005) Effects of noisy drive on rhythms in networks of excitatory and inhibitory neurons. Neural Comput 17(3):557–608 Byrne A, Avitabile D, Coombes S (2019) Next-generation neural field model: the evolution of synchronywithin patterns andwaves. Phys Rev E 99:012313 Byrne Á, Ross J, Nicks R, Coombes S (2022) Mean-field models for EEG/MEG: from oscillations to waves. Brain Topogr 35(1):36–53 Campos J (1997) Möbius transformations and periodic solutions of complex Riccati equations. Bull Lond Math Soc 29:205–215 Coombes S, Laing C (2009) Delays in activity-based neural networks. Philos Trans R Soc A Math Phys Eng Sci 367(1891):1117–1129 Coombes S, Byrne Á (2019) Next generation neural mass models. In: Nonlinear dynamics in computational neuroscience, pp 1–16. Springer Daleckiı̆ JL,Kreı̆nMG (2002) Stability of solutions of differential equa- tions in Banach space, vol 43. American Mathematical Society, New York Devalle F, Roxin A, Montbrió E (2017) Firing rate equations require a spike synchrony mechanism to correctly describe fast oscillations in inhibitory networks. PLoS Comput Biol 13(12):e1005881 Devalle F, Montbrió E, Pazó D (2018) Dynamics of a large system of spiking neurons with synaptic delay. Phys Rev E 98(4):042214 di Volo M, Torcini A (2018) Transition from asynchronous to oscil- latory dynamics in balanced spiking networks with instantaneous synapses. Phys Rev Lett 121:128301 Ermentrout GB, Kopell N (1986) Parabolic bursting in an excitable system coupled with a slow oscillation. SIAM J Appl Math 46(2):233–253 Esnaola-Acebes JM, Roxin A, Avitabile D, Montbrió E (2017) Synchrony-induced modes of oscillation of a neural field model. Phys Rev E 96:052407 Funahashi S, Bruce CJ, Goldman-Rakic PS (1989) Mnemonic coding of visual space in the monkey’s dorsolateral prefrontal cortex. J Neurophysiol 61(2):331–349 Gallego R,Montbrió E, Pazó D (2017) Synchronization scenarios in the Winfree model of coupled oscillators. Phys Rev E 96(4):042208 Jiruska P, de Curtis M, Jefferys JGR, Schevon CA, Schiff SJ, Schindler K (2013) Synchronization and desynchronization in epilepsy: con- troversies and hypotheses. J Physiol 591(4):787–797 Keane A, Krauskopf B, Dijkstra HA (2019) The effect of state depen- dence in a delay differential equationmodel for the el niño southern oscillation. Philos Trans R Soc A 377(2153):20180121 Kilpatrick ZP, Ermentrout B (2013) Wandering bumps in stochastic neural fields. SIAM J Appl Dyn Syst 12(1):61–94 Kuramoto Y (1984) Chemical oscillations, waves, and turbulence. Springer, Berlin Laing CR (2017) Phase oscillator network models of brain dynamics. In:Moustafa A (ed) ComputationalModels of Brain and Behavior, chap. 37, pp 505–517. Wiley-Blackwell, Hoboken, NJ Laing CR (2014) Derivation of a neural field model from a network of theta neurons. Phys Rev E 90(1):010901 Laing CR (2014) Numerical bifurcation theory for high-dimensional neural models. J Math Neurosci 4(1):13 Laing CR (2015) Exact neural fields incorporating gap junctions. SIAM J Appl Dyn Syst 14(4):1899–1929 Laing C, Chow C (2001) Stationary bumps in networks of spiking neu- rons. Neural Comput 13(7):1473–1494 Laing CR, Longtin A (2003) Dynamics of deterministic and stochas- tic paired excitatory-inhibitory delayed feedback. Neural Comput 15(12):2779–2822 Laing CR, Omel’chenko O (2020) Moving bumps in theta neuron net- works. Chaos Interdiscip J Nonlinear Sci 30(4):043117 Laing CR, Troy W (2003) PDE methods for nonlocal models. SIAM J Appl Dyn Syst 2(3):487–516 Laing CR, TroyWC, Gutkin B, Ermentrout GB (2002) Multiple bumps in a neuronal model of working memory. SIAM J Appl Math 63(1):62–97 LaingCR, Bläsche C,Means S (2021) Dynamics of structured networks of Winfree oscillators. Front Syst Neurosci 15:631377 Latham P, Richmond B, Nelson P, Nirenberg S (2000) Intrinsic dynam- ics in neuronal networks. I. Theory. J Neurophysiol 83(2):808–827 Lee WS, Ott E, Antonsen TM (2009) Large coupled oscillator systems with heterogeneous interaction delays. Phys Rev Lett 103:044101 Lindén H, Petersen PC, Vestergaard M, Berg RW (2022) Movement is governed by rotational neural dynamics in spinal motor networks. Nature 610:526–531 Montbrió E, Pazó D, Roxin A (2015) Macroscopic description for net- works of spiking neurons. Phys Rev X 5:021028 Netoff TI, Schiff SJ (2002) Decreased neuronal synchronization during experimental seizures. J Neurosci 22(16):7297–7307 Omel’chenko OE (2020) Nonstationary coherence–incoherence pat- terns in nonlocally coupled heterogeneous phase oscillators. Chaos Interdiscip J Nonlinear Sci 30(4):043103 123 274 Biological Cybernetics (2023) 117:259–274 Omel’chenko O (2022)Mathematical framework for breathing chimera states. J Nonlinear Sci 32(2):1–34 Omel’chenko OE (2023) Periodic orbits in the Ott–Antonsen manifold. Nonlinearity 36:845–861 Omel’chenko O, Laing CR (2022) Collective states in a ring network of theta neurons. Proc R Soc A 478(2259):20210817 Omel’chenko O, Wolfrum M, Laing CR (2014) Partially coherent twisted states in arrays of coupled phase oscillators. Chaos 24:023102 Ott E, Antonsen TM (2008) Low dimensional behavior of large systems of globally coupled oscillators. Chaos 18(3):037113 Ott E, Antonsen TM (2009) Long time evolution of phase oscillator systems. Chaos 19(2):023117 Pazó D, Montbrió E (2014) Low-dimensional dynamics of populations of pulse-coupled oscillators. Phys Rev X 4:011009 Pietras B (2023) Pulse shape and voltage-dependent synchronization in spiking neuron networks. arXiv:2304.09813 Pinto DJ, Ermentrout GB (2001) Spatially structured activity in synap- tically coupled neuronal networks: II. Lateral inhibition and standing pulses. SIAM J Appl Math 62(1):226–243 Ratas I, Pyragas K (2016) Macroscopic self-oscillations and aging transition in a network of synaptically coupled quadratic integrate- and-fire neurons. Phys Rev E 94:032215 Reyner-Parra D, Huguet G (2022) Phase-locking patterns underlying effective communication in exact firing rate models of neural net- works. PLoS Comput Biol 18(5):e1009342 Roxin A, Brunel N, Hansel D (2005) Role of delays in shaping spa- tiotemporal dynamics of neuronal activity in large networks. Phys Rev Lett 94(23):238103 Schmidt H, Avitabile D (2020) Bumps and oscillons in networks of spiking neurons. Chaos 30:033133 Schmidt H, Avitabile D, Montbrió E, Roxin A (2018) Network mecha- nisms underlying the role of oscillations in cognitive tasks. PLoS Comput Biol 14(9):e1006430 Segneri M, Bi H, Olmi S, Torcini A (2020) Theta-nested gamma oscil- lations in next generation neural mass models. Front Comput Neurosci 14:47 Shima S, Kuramoto Y (2004) Rotating spiral waves with phase- randomized core in nonlocally coupled oscillators. Phys Rev E 69(3):036213 Strogatz S (2000) From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators. Physica D 143(1–4):1–20 Uhlhaas PJ, Singer W (2010) Abnormal neural oscillations and syn- chrony in schizophrenia. Nature Rev Neurosci 11(2):100–113 Wilczyński P (2008) Planar nonautonomous polynomial equations: the Riccati equation. J Differ Equ 244:1304–1328 Wimmer K, Nykamp DQ, Constantinidis C, Compte A (2014) Bump attractor dynamics in prefrontal cortex explains behavioral preci- sion in spatial working memory. Nature Neurosci 17(3):431–439 Zhang K (1996) Representation of spatial orientation by the intrinsic dynamics of the head-direction cell ensemble: a theory. J Neurosci 16(6):2112–2126 Publisher’s Note Springer Nature remains neutral with regard to juris- dictional claims in published maps and institutional affiliations. 123 http://arxiv.org/abs/2304.09813 Periodic solutions in next generation neural field models Abstract Introduction 1 Theta neuron network model 2 Periodic states 2.1 Periodic complex Riccati equation and Möbius transformation 2.2 Self-consistency equation 2.3 Numerical implementation 2.4 Results 2.5 Stability of breathing bumps 2.6 Formula for firing rates 3 Other models 3.1 Delays 3.2 Two populations 3.3 Winfree oscillators 4 Discussion Acknowledgements Appendix References