Bump attractors and waves in networks of leaky integrate-and-fire neurons
BBUMP ATTRACTORS AND WAVES IN NETWORKS OF LEAKYINTEGRATE-AND-FIRE NEURONS
DANIELE AVITABILE ∗ , JOSHUA L. DAVIS † , AND
KYLE C. A. WEDGWOOD ‡ Abstract.
Bump attractors are wandering localised patterns observed in in vivo experimentsof spatially-extended neurobiological networks. They are important for the brain’s navigationalsystem and specific memory tasks. A bump attractor is characterised by a core in which neuronsfire frequently, while those away from the core do not fire. These structures have been found insimulations of spiking neural networks, but we do not yet have a mathematical understanding of theirexistence, because a rigorous analysis of the nonsmooth networks that support them is challenging.We uncover a relationship between bump attractors and travelling waves in a classical network ofexcitable, leaky integrate-and-fire neurons. This relationship bears strong similarities to the onebetween complex spatiotemporal patterns and waves at the onset of pipe turbulence. Waves in thespiking network are determined by a firing set, that is, the collection of times at which neuronsreach a threshold and fire as the wave propagates. We define and study analytical properties of the voltage mapping , an operator transforming a solution’s firing set into its spatiotemporal profile. Thisoperator allows us to construct localised travelling waves with an arbitrary number of spikes at thecore, and to study their linear stability. A homogeneous “laminar” state exists in the network, andit is linearly stable for all values of the principal control parameter. Sufficiently wide disturbancesto the homogeneous state elicit the bump attractor. We show that one can construct waves with aseemingly arbitrary number of spikes at the core; the higher the number of spikes, the slower thewave, and the more its profile resembles a stationary bump. As in the fluid-dynamical analogy, suchwaves coexist with the homogeneous state, are unstable, and the solution branches to which theybelong are disconnected from the laminar state; we provide evidence that the dynamics of the bumpattractor displays echoes of the unstable waves, which form its building blocks.
1. Introduction.
Understanding how networks of coupled, excitable units gen-erate collective patterns is a central question in the life sciences and, more generally,in applied mathematics. In particular, the study of network models is ingrained inneuroscience applications, as they provide a natural way to describe the interactionof neurons within a population, or of neural populations within the cortex. In thepast decades, a large body of work in mathematical neuroscience has addressed thedevelopment and analysis of neurobiological networks, with the view of studying theorigin of large-scale brain activity [35, 13, 21], and mapping single-cell and popula-tion parameters to experimental observations, including in vivo and in vitro corticalwaves [70, 47, 43], electroencephalogram recordings [77], and patterns in the visualcortex [17].This paper presents a novel mathematical characterisation of a prominent exampleof spatiotemporal pattern in neuroscience applications, and draws an analogy inspiredby recent progress in the fluid-dynamics literature on transition to turbulence in apipe [7]. We focus on the so-called bump attractor , a localised pattern of neuralactivity observable in experiments and numerical simulations of spatially-extended,neurobiological networks [68, 87]. Bump attractors have been associated to workingmemory , the temporary storage of information in the brain, and experimental evidencesupporting their existence has been found in the navigational systems of rats [52] and ∗ † Defence Science and Technology Laboratory, Cyber and Information Systems Division. ‡ Living Systems Institute, College of Engineering, Mathematics and Physical Sciences. In the neuroscience literature the term bump attractor refers sometimes to a network producinga localised pattern, as opposed to the pattern itself. Similarly, some authors use ring attractor for anetwork with ring topology, generating a localised activity bump. Here, we use these terms to referto patterns, following the standard convention in the dynamical system literature.1
D. AVITABILE, J. L. DAVIS, K. C. A. WEDGWOOD
Fig. 1 . (a) Bifurcation diagram of travelling waves in a continuous integrate-and-fire model.(b) Bump attractor in a discrete integrate-and-fire model with 5000 neurons (dots represent neu-ronal firing events). Model descriptions and parameters will be given later in section , Table ,and Figures and ). The bifurcation diagram in (a) shows selected branches of stable (blue) andunstable (grey) travelling waves, in the continuation parameter β , that is, the timescale at whichneuron process incoming currents. Waves are measured using their width ∆ and are indexed by thenumber of advected spikes. The profile of TW , a representative wave with spikes, is shown. Alarge number of waves (TW –TW in the picture, but many more unstable branches are omitted)coexist with the trivial homogeneous state, which is the only steady state in the model, and which isstable for all values of β . Thin waves are stable, with a small basin of attraction. Sufficiently large,localised disturbances of the homogenous state lead to the formation of a bump with a characteris-tic width: the bump in (b) is marked as (2) in (a). The region in parameter space where bumpsare observed is crowded with unstable travelling waves, with a large number of spikes, and a widthcomparable to the one of the bump. Branches of waves are detached from the homogeneous state;they originate at critical points called grazing points (blue dots in (a)); waves that are born stablebecome unstable at oscillatory bifurcations (grey dots in (a)). flies [51, 79], and in oculomotor responses of monkeys [85].In a bump attractor, the neural activity is localised around a particular positionin the network (see Figure 1(b)) which may encode, for instance, the animal’s headposition. Bumps are elicited by transient localised stimuli, such as visual cues at aspecific locations, but are sustained autonomously by the network once the stimulus isremoved (the network dynamics is attracted to the bump). These coherent structuresdisplay a characteristic wandering motion, and may exhibit discontinuous jumps ifthe impinging stimulus undergoes sudden spatial shifts [51]. Mathematical neuroscience has a long-standing fas-cination with localised bumps of activity. Neural field models, which represent thecortex as a continuum, were introduced in the 1970s, and spatially-localised solutionsto these models appeared already in seminal papers on the subject, by Wilson andCowan [84], and by Amari [1]. Since then, many authors have studied localised so-lutions in neural fields, addressed the derivation of neural field equations from firstprinciples, their relevance to a wide variety of neural phenomena, and their rigorousmathematical treatment. We refer the reader to [35, 13, 21] for exhaustive introduc-tions on this topic.
UMPS AND WAVES IN SPIKING NETWORKS
Spiking neural networks are intermediate, bottom-up models which couple neu-rons with idealised dynamics. The salient feature of spiking models is that the firingof a neuron is described as an event, and no attempt is made to model the temporalevolution of the membrane potential during and after the spike [49, 39, 13]. Spik-ing neural networks are specified by 3 main ingredients: (i) an ordinary differentialequation (ODE) for the membrane potential of each neuron; (ii) rules to define theoccurrence and effects of a spike; (iii) the network coupling.Since the introduction of the first single-cell spiking model by Lapicque [57], theso-called leaky integrate-and-fire model , more realistic variants have been proposed,and spiking neural networks have become a widely adopted tool in theoretical neu-roscience [78, 15, 39]. In specific spiking models, analytical progress has been madefor single neurons and spatially-independent networks using coordinate transforma-tions [34, 62], dimension reduction [58, 63], and probabilistic methods [26] (see also thereviews [72, 9]). Exact mean-field reductions, amenable to standard pattern-formationanalysis, have been derived in selected spatially-extended networks [53, 36, 16, 74], butgenerally the study of bumps in spiking models has been possible only with numericalsimulations [54, 18].The present paper investigates localised patterns supported in discrete and contin-uous networks of globally coupled leaky integrate-and-fire neurons. In direct numericalsimulations, we use a well-known discrete model, proposed by Laing and Chow [54],whose details will be given later. For now it will suffice to consider a cursory formu-lation of the model, simulated in Figure 1(b). The network describes the idealised,dimensionless voltage dynamics of n all-to-all coupled neurons, evenly-spaced in acortex with ring geometry,(1.1) ˙ v i = − v i + I i ( t ) + n (cid:88) j =1 S ij ( v j , β ) , i = 1 , . . . , n. The dynamics of the i th neuron’s membrane voltage is specified in terms of an Ohmicleakage current − v i , an external current I i ( t ), and voltage-dependent currents, re-ceived from other neurons via synaptic connections; the latter currents, indicated by S ij , have a characteristic time scale β , and are caused by v j crossing a fixed threshold D. AVITABILE, J. L. DAVIS, K. C. A. WEDGWOOD (when the j th neuron fires ). After a firing event, marked with a dot in Figure 1(b),the neuron’s voltage is instantly reset to a default value, from which it can evolveagain, following an ODE of type (1.1). Discrete and continuous networks of this typeare canonical models of neural activity, widely adopted in the mathematical neuro-science literature [82, 81, 30, 31, 14, 54, 32, 65, 19, 64, 62, 40]. It is now establishedthat such networks supports bump attractors and localised waves, but an explanationof the mathematical origins of the former is still lacking.This paper presents a new approach to the problem, and uncovers a novel bifur-cation structure for localised travelling waves of the network, shedding light onto thenature of the bump attractor. Our findings suggest an intriguing analogy betweenthe bump attractor in the integrate-and-fire network and the phenomenon of transi-tion to turbulence in a pipe. The analogy between the bifurcation scenarios of thesetwo problems is notable, and we use it here to summarise our results, highlightingsimilarities between the respective bifurcation structures and dynamical regimes. Stemming from the pioneering ex-periments of Reynolds [69], a large body of work in fluid dynamics has addressed howhigh-speed pipe flows transition from a laminar state, whose analytical expressionis known in closed form, to complex spatio-temporal patterns, characteristic of theturbulent regime (see [7] for a recent review). In this context, the Navier-Stokes equa-tions are studied as a deterministic dynamical system, subject to changes in Reynoldsnumber, the principal control parameter. Experiments and computer simulations in-dicate that the laminar state is stable to infinitesimal perturbations (linearly stable)up to large values of the control parameter (up to at least Reynolds number 10 innumerical computations) [73, 25, 61, 80, 60]. However, when a disturbance is ap-plied at sufficiently large Reynolds numbers, a transition to turbulence is observed,depending sensitively on the applied stimulus [25, 45]. Current opinions view the tran-sition as being determined by travelling wave solutions to the Navier-Stokes equations[75, 29, 37, 83, 66, 41]. These invariant states, whose spatial profiles display hallmarksof the turbulent transition, coexist with the laminar state at intermediate Reynoldsnumbers, are linearly unstable, and provide an intricate blueprint for the dynamics, inthat orbits may visit transiently these repelling solutions in phase space. Importantly,the waves lie on branches that are disconnected from the stable laminar state, andemerge at saddle-node bifurcations [37, 83]: this turbulence mechanism is thereforedifferent from other paradigmatic routes to chaos, involving the destabilisation of thelaminar state, and the progressive appearance of more complicated structures via acascade of instabilities [56, 46, 71]. In a series of recent papers addressing turbulencefrom a dynamical-system viewpoint, Barkley proposed an analogy between pipe flowsand excitable media, using the propagation of an electrical pulse along the axon ofa neuron as a metaphor for localised turbulence puffs [5, 6, 8, 7]. The present paperoffers a specular view, at a different scale: we are motivated by studying a canonical,complex neurobiological network of coupled excitable neurons, supporting localisedspatio-temporal chaos, and we find a compelling similarity between the bifurcationstructure of waves in this system, and the one of waves in the pipe turbulence.With reference to Figure 1, the principal control parameter of the problem is β , the timescale of synaptic currents: a low β gives small, persisting currents, while β → ∞ gives large instantaneous currents. A homogeneous steady state exists and islinearly stable for all values of β (∆ = 0 line in Figure 1(b)), but transient localisedstimuli trigger the bump attractor [54]. In the analogy, the homogeneous equilibrium UMPS AND WAVES IN SPIKING NETWORKS localised spikes , or having a non-localisedprofile [31, 12, 14, 32, 65]. The travelling waves of interest to us, however, have alocalised profile, and advect a large number of spikes, such as the one presented inFigure 1(a). These structures are not accessible with the current techniques, hence wedevelop here analytical and numerical tools to construct them. We define particulartype of solutions, which retain a fixed number of spikes in time; this class of solutionsis sufficiently general to incorporate travelling waves with an arbitrary, finite, numberof spikes, and small perturbations to them. We introduce the voltage mapping , a newoperator which formalises an idea previously used in the literature for spiking [31,12, 14, 32, 65, 3] and non-spiking networks [1, 35, 13, 21]. The voltage mappingis based on level sets describing firing events, and it allows efficient travelling waveconstructions and stability computations.Using the voltage mapping, we construct numerically waves with more than 200concurrent spikes. These waves are spatially localised, linearly unstable, and coexistwith the trivial (laminar) state (see Figure 1(a)). As in the turbulence analogy,the waves contain features of the bump attractor: they pack a seemingly arbitrarynumber of spikes within the width of a bump attractor, and they advect them at anarbitrarily slow speed, depending on β , and on the number of carried spikes. As in thefluid-dynamical analogy, waves are disconnected from the laminar state. Owing to theintrinsic non-smoothness of the network, the waves emerge primarily at grazing points(as opposed to saddle-node bifurcations seen in the fluid-dynamical analogy), whichare however seen in certain parameter regimes. In addition, we present numericalevidence that the dynamics of the bump attractor displays echoes of the unstablewaves which, as in the fluid-dynamics analogy, form building blocks for the localisedstructure. Also, the characteristic wandering of the bump attractor, whose excursionsbecome more prominent as β increases, is supported by this purely deterministic system.The paper is structured as follows: in section 2 we introduce the discrete model,characterise it as a non-smooth threshold network, and present numerical simulationsof bumps and waves; in section 3 we introduce the continuum model, the voltagemapping, and the construction of travelling waves; in section 4 we discuss travellingstability, we present numerical results in section 5, and we conclude in section 6.
2. Coherent structures in the discrete model.
We begin by introducingthe discrete model by Laing and Chow [54]. We characterise it as a piecewise-lineardynamical system, and we show numerical simulations of coherent structures. Animportant difference from the work by Laing and Chow is that we consider a de-terministic model, which we call the Discrete Integrate-and-Fire Model (DIFM). Weremark that the neurons considered here, taken in isolation, are in an excitable regime ,that is, they exhibit an all-or-none response, based on the input they receive. This isconsiderably different from the so-called oscillatory regime , in which neurons, whendecoupled from the network, display oscillations [82, 14, 62, 40].
The DIFM is a spatially-extended system of n identical integrate-and-fire neurons, posed on S = R / L Z , that is, a ring of period2 L . Neurons are indexed using the set N n = { , . . . , n } and occupy the discrete, D. AVITABILE, J. L. DAVIS, K. C. A. WEDGWOOD evenly spaced nodes x i = − L + 2 iL/n ∈ S , for i ∈ N n . Neurons are coupled viatheir synaptic connections, which are modelled by a continuous, bounded, even andexponentially decaying function w : S → R : the strength of the connections from the k th to the i th neuron depends solely on the distance | x i − x k | , measured around thering, hence we write it as W ik = w ( x i − x k ), for all i, k ∈ N n .To the i th neuron is associated a real-valued time-dependent voltage function v i ( t ), and the coherent structures of interest are generated when voltages { v i } attaina threshold value (when neurons fire ). The DIFM is formally written as follows:˙ v i ( t ) = I i ( t ) − v i ( t ) + 2 Ln (cid:88) k ∈ N n (cid:88) j ∈ N W ik α ( t − τ jk ) − (cid:88) j ∈ N δ ( t − τ ji ) , i ∈ N n , (2.1) v i (0) = v i , i ∈ N n . (2.2)At time τ ji , when the voltage v i reaches the value 1 from below for the j th time, afiring event occurs; a more precise definition of these spiking times will be given below.The formal evolution equation (2.1) expresses the modelling assumption that, whena neuron fires, its voltage is instantaneously reset to 0 (hence the Dirac delta), and aso-called post-synaptic current is received by all other neurons in the network, withintensity proportional to the strength of the synaptic connections. The time-evolutionof this current is modelled via the post-synaptic function α ( t ) = p ( t ) H ( t ), expressedas the product of a continuous potential function p and the Heaviside function H ,hence the post-synaptic current is zero before a spike.In this paper we present concrete calculations for(2.3) α ( t ) = β exp( − βt ) H ( t ) , w ( x ) = a exp( − b | x | ) − a exp( − b | x | ) , with β, a , a , b , b >
0, albeit the analytical and numerical framework presentedbelow is valid for more generic choices, subject to general assumptions which willbe made precise in subsection 3.2. The function α models exponentially-decayingcurrents with rate − β and initial value β , hence the limit β → ∞ approximatesinstantaneous currents. Currents with an exponential rise and decay are also used inliterature. The synaptic coupling function w is chosen so that connections are positive( excitatory ) on the lengthscale 1 /b , and negative ( inhibitory ) on the lengthscale 1 /b .In addition to the post-synaptic current, neurons are subject to an external stim-ulus I i ( t ). In certain time simulations, coherent structures will be elicited with theapplication of a transient, heterogeneous stimulus of the form(2.4) I i ( t ) = I + d H ( t − τ ext ) / cosh( d x i ) , i ∈ N n . Our investigation, however, concerns asymptotic states of the autonomous homoge-neous case I i ( t ) ≡ I , hence one should assume d = 0, unless stated otherwise. Adescription of model parameters and their nominal values can be found in Table 1. Laing and Chow studied and simulated a stochasticversion of the DIFM, using the Euler method and a first-order interpolation schemeto obtain the firing times [54]. We use here a different approach: in preparation forour analytical and numerical treatment of the problem, we write the formal model(2.1)–(2.2) as a system of 2 n piecewise-linear ODEs. To this end we introduce thesynaptic input variables(2.5) s i ( t ) = 2 Ln (cid:88) k ∈ N n (cid:88) j ∈ N W ik α ( t − τ jk ) , i ∈ N n UMPS AND WAVES IN SPIKING NETWORKS Table 1
Parameter descriptions and nominal values.
Parameter Symbol Value(s)Number of neurons n { } Domain half-width L { } Synaptic efficacy and time scale β [0,25]Synaptic excitation coefficient a a b b I τ ext d { } Time-dependent external input spatial scale d { } and combining (2.3) and (2.1) we obtain formally˙ v i ( t ) = I i ( t ) − v i ( t ) + s i ( t ) − (cid:88) j ∈ N δ ( t − τ ji )˙ s i ( t ) = − βs i ( t ) + 2 Lβn (cid:88) k ∈ N n (cid:88) j ∈ N W ik δ ( t − τ jk ) i ∈ N n . One way to define the associated non-smooth dynamical system is to express themodel as an impacting system, by partitioning the phase space R n via a switchingmanifold, on which a reset map is prescribed (see [27] and references therein for adiscussion on non-smooth and impacting systems). Here, we specify the dynamicsso as to expose the firing times { τ jk } , as opposed to the switching manifold: this isnatural in the mathematical neuroscience context, and it prepares our analysis of thecontinuum model. Since { τ jk } are the times at which orbits in R n reach the switchingmanifold, a translation between the two formalisms is possible.Following these considerations, we set τ i = 0 for all i ∈ N n , introduce the notation f ( · ± ) = lim µ → + f ( · ± µ ), and define firing times as follows (2.6) τ ji = inf (cid:8) t ∈ R : t > τ j − i , v i ( t − ) = 1 , ˙ v i ( t − ) > (cid:9) , i ∈ N n , j ∈ N . We arrange firing times in a monotonic increasing sequence { τ j k i k } qk =1 such that(2.7) (0 , T ] = (cid:91) k ∈ N q +1 (cid:0) τ j k − i k − , τ j k i k (cid:3) , τ j i < τ j i ≤ . . . ≤ τ j q i q < τ j q +1 i q +1 = T, for some time horizon T >
0, and obtain the desired set of 2 n piecewise-linear ODEs(2.8) ˙ v i = I i − v i + s i , ˙ s i = − βs i i ∈ N n , t ∈ (cid:91) k ∈ N q +1 (cid:0) τ j k − i k − , τ j k i k (cid:3) , Note that { τ i } i are not firing times, but auxiliary symbols for the definition of firing times (2.6).Indeed, since the sums in (2.1) run for j ∈ N , the { τ i } i are immaterial for the dynamics. D. AVITABILE, J. L. DAVIS, K. C. A. WEDGWOOD (a) v i ( t ) { ( x i , t ) : v i ( t ) = 1 } t x i s i ( t ) x i
UMPS AND WAVES IN SPIKING NETWORKS s i ( t ) v i ( t ) { ( x i , t ) : v i ( t ) = 1 } (b) s i v i . . x i t x i x i x i s i v i . . t x i x i
The DIFM supports standing andtravelling localised structures, as in the stochastic setting [54]. Bumps form robustlywhen we prescribe homogeneous initial conditions with a short transient stimulus(Equation (2.4) with τ ext = 2). Since I i ( t ) ≡ I for all t > τ ext , the structures observedover long-time intervals are solutions to a homogeneous, non-autonomous problem.As seen in Figure 2, the bump wanders when β is increased. In passing, we note Typically we set v i = u ∈ (0 , , s i = 0, for i ∈ N n , but the coherent structures discussed inthe paper can also be found with random, independent and identically distributed initial voltages,for instance v i ∼ U ([0 , U is the uniform distribution. D. AVITABILE, J. L. DAVIS, K. C. A. WEDGWOOD that this phenomenon is not due to stochastic effects, as studied in other contexts [50,48, 3], because the DIFM is deterministic. For sufficiently large β , the system exhibitsstable travelling structures: in Figure 3 we show two coexisting waves, found for β = 4 . d and intensity d of the transient stimulus.In each case we plot the voltage and synaptic profiles, and associated raster plots. Wenotice different firing patterns in the waves, involving 2 and 4 firings, respectively: thewave with 2 firings travels faster, and its voltage and synaptic profiles are narrower.We found coexisting waves with a greater number of firings and progressively lowerspeed, whose existence and bifurcation structure will be at the core of the followingsections. The patterns pre-sented so far are found in the DIFM with a finite number of neurons ( n = 80). At firstsight, the raster plots of the waves seem to indicate that neurons fire simultaneouslyin pairs (Figure 3(a)) or quartets (Figure 3(b)) as the structure travels across the net-work. A closer inspection of the instantaneous profiles v i ( t ) reveal that this is not thecase, as the threshold (red dashed line) is attained by a single neuron in Figure 3(a),and by two neurons in Figure 3(b): neurons in a raster pair fire alternatively at shorttimes from each other, whereas a quartet displays a more complex firing pattern.Hence, for finite n , the propagating structures displayed in Figure 3 are not strictlytravelling waves, in the sense that the profile is not stationary in the comoving frame;their dynamics is that of saltatory waves [22, 86, 3]. The saltatory nature of the waves,however, is an effect of the network size: as we increase n , the amplitude of temporaloscillations in the comoving frame scales as O ( n − ), and the spatio-temporal profileconverges to the one of a travelling wave as n → ∞ .In addition, the structure in Figure 2(a) is not a bump, in the sense that it is not aspatially heterogeneous steady state of the DIFM, because the pattern is sustained byfiring events. Indeed, the only equilibrium supported by the DIFM is the homogeneousstate v i ( t ) ≡ I , s i ( t ) ≡ i ∈ N n , which is linearly stable for all values of β , as can bededuced by inspecting system (2.8).It will be shown below that we can gain insight into the structure shown inFigure 2(a) (and its wandering) by constructing travelling waves and investigatingtheir stability in a continuum version of the DIFM.
3. Travelling waves in the continuum model.
As stated in section 2, theprofiles { v i ( t ) } i and { s i ( t ) } i in Figure 3 behave like travelling wave solutions as n → ∞ . Motivated by this observation, we study travelling waves in a continuum,translation-invariant version of the DIFM: we set d = 0 in the stimulus (2.4), considera continuum spatial domain, and pose the model on R as opposed to S , obtaining ∂ t v ( x, t ) = − v ( x, t ) + I + (cid:88) j ∈ N (cid:90) ∞−∞ w ( x − y ) α (cid:0) t − τ j ( y ) (cid:1) dy − (cid:88) j ∈ N δ (cid:0) t − τ j ( x ) (cid:1) , ( x, t ) ∈ R × R . (3.1)The formal evolution equation presented above, which we henceforth call the con-tinuous integrate-and-fire model (CIFM), has been proposed and studied by severalauthors in the mathematical neuroscience literature [31, 42, 10, 12, 65, 64]. In theCIFM, firing-time functions τ j ( x ) indicate that the neural patch at position x fires UMPS AND WAVES IN SPIKING NETWORKS j th time, and replace the discrete model’s firing times τ jk . A graph of thefiring functions replaces the raster plot in the discrete model, so that a travelling wavein the CIFM corresponding to the n → ∞ limit of the structure in Figure 3(a), forinstance, will involve 2 linear firing functions τ , τ , with τ ( x ) < τ ( x ) for all x ∈ R .The existence of travelling waves solutions in (3.1) with a single spike has beenstudied by Ermentrout [31] who presented various scalings of the wavespeed as afunction of control parameters. A general formalism for the construction and linearstability analysis of wavetrains (spatially-periodic travelling solutions) was introducedand analysed by Bressloff [12], who derived results in terms of Fourier series expan-sions. The construction of travelling waves with multiple spikes was later studiedby O¸san and coworkers [64], albeit stability for these states was not presented andcomputations were limited to a few spikes, for purely excitatory connectivity kernels.The common thread in the past literature on this topic is the idea that travelling waveconstruction and stability analysis rely entirely on knowledge of the firing function τ j (as in the DIFM, with firing times). A similar approach has been used effectivelyin Wilson-Cowan-Amari neural field equations, where it is often called interfacial dy-namics (see [ ? ] for the first study of this type, [20] for a recent review, and [23, 38],amongst others, for examples of spatio-temporal patterns analysis).Here we present a new treatment of travelling wave solutions that draws from thisidea; we introduce an operator, that we call the voltage mapping , with the followingaims: (i) Expressing a mapping between firing functions and solution profiles, withthe view of replacing the formal evolution equation (3.1) for travelling waves with m spikes (where m is arbitrary). (ii) Finding conditions for the linear stability of thesewaves. (iii) Using root-finding algorithms to compute travelling waves and study theirlinear stability. We will relate to existing literature in our discussion. Before analysing solutions to the CIFM, we discuss the notationused in this section. For fixed m ∈ N , we use | · | ∞ to denote, with a little abuse ofnotation, the ∞ -norm on C m . We will use | · | for the standard modulus in C . Wedenote by C ( X, Y ) the set of continuous functions from X to Y , and use C ( X ) when Y = R . We denote by B ( X ) the set of real-valued bounded functions defined on a set X , and by BC ( X ) the set of real-valued, bounded and continuous functions definedon X , respectively; both spaces are endowed with the supremum norm (cid:107) · (cid:107) ∞ . Wedenote by L ( R ) the pace of Lebesgue-integrable functions defined on R . Further, fora positive number η , we define the exponentially weighted space L η ( R ) = (cid:110) u : R → R : (cid:107) u (cid:107) L η = (cid:90) R e ηx | u ( x ) | dx < ∞ (cid:111) , which is a Banach space equipped with the norm (cid:107) · (cid:107) L η , and the following space ofexponentially growing functions C η ( R , C m ) = (cid:110) u ∈ C ( R , C m ) : (cid:107) u (cid:107) C m,η = sup x ∈ R e − η | x | | u ( x ) | ∞ < ∞ (cid:111) , which is a Banach space equipped with the norm (cid:107) · (cid:107) C m,η . We begin by discussing in what sense a voltage function v satisfies the CIFM The index j is used as a superscript in the firing times, but for notational convenience we useit as a subscript in the firing functions, so that τ j ( x k ) ≈ τ jk . D. AVITABILE, J. L. DAVIS, K. C. A. WEDGWOOD formal evolution equation (3.1). While we eschew the definition of the CIFM as a dy-namical system on a Banach space (a characterisation that is currently unavailable inthe literature), we note that progress can be made for voltage profiles with a constantand finite number of spikes for t ∈ R . This class of solutions is sufficiently large totreat travelling waves, and small perturbations to them.We make a few assumptions on the network coupling, and we restrict the type offiring functions and solutions of interest, as follows: Hypothesis w is an even func-tion in C ( R ) ∩ L η ( R ), for some η >
0. The post-synaptic function α : R → R ≥ canbe written as α ( t ) = p ( t ) H ( t ), where H is the Heaviside function, and p : R ≥ → R isa bounded and everywhere differentiable Lipschitz function, hence p, p (cid:48) ∈ B ( R ). Definition m -spike CIFM solution). Let m ∈ N and I ∈ R . A function v m : R → R is an m -spike CIFM solution if there exists τ = ( τ , . . . , τ m ) ∈ C ( R , R m ) such that τ < . . . < τ m on R and v m ( x, t ) = I + (cid:88) j ∈ N m (cid:90) t −∞ (cid:90) ∞−∞ exp( z − t ) w ( x − y ) α ( z − τ j ( y )) dy dz − (cid:88) j ∈ N m exp( τ j ( x ) − t ) H ( t − τ j ( x )) , ( x, t ) ∈ R (3.2) v m ( x, t ) = 1 , ( x, t ) ∈ F τ , (3.3) v m ( x, t ) < , ( x, t ) ∈ R \ F τ , (3.4) where F τ = (cid:91) j ∈ N m { ( x, t ) ∈ R : t = τ j ( x ) } . We call τ and F τ the firing functions and the firing set of v m , respectively. The definition above specifies how we interpret solutions to (3.1), and is composed ofthree ingredients: (i) Equation (3.2), which derives from integrating (3.1) on ( −∞ , t ),and expresses a mapping between the set of m firing functions τ and the voltageprofile; (ii) System (3.3), which couples the firing functions by imposing the thresholdcrossings; (iii) A further condition on v m , ensuring that the solution has exactly m spikes, attained at the firing set; this is necessary because, as we shall see below, it ispossible to find a set of m functions τ satisfying Equations (3.2)–(3.3), but exhibitinga number of threshold crossings greater than m .We now aim to characterise m -spike CIFM solutions by means of a voltage map-ping , which can be conveniently linearised around a firing set, and is a key tool toconstruct waves and analyse their stability. Inspecting (3.2) we note that the voltageprofile features two contributions, one from the (synaptic) coupling functions w and α , and one from reset conditions. This observation leads to the following definitions: Definition
Let u : R → R . We define thesynaptic operator, S , and the reset operator, R , by ( Su )( x, t ) = (cid:90) t −∞ (cid:90) ∞−∞ exp( z − t ) w ( x − y ) α ( z − u ( y )) dy dz, ( x, t ) ∈ R , (3.5) ( Ru )( x, t ) = − exp( u ( x ) − t ) H ( t − u ( x )) , ( x, t ) ∈ R . (3.6)These operators map univariate functions, such as a firing function, to bivariate func-tions, contributing to the spatio-temporal voltage profile. The following lemma shows UMPS AND WAVES IN SPIKING NETWORKS
Lemma
If Hypothesis holds, then for the operators S , R in Definition we have S : C ( R ) → BC ( R ) and R : C ( R ) → B ( R ) , respectively.Proof. The statement for R is immediate, whereas proving the continuity of Su on R requires some estimates for the improper integral. A proof is given in Appendix A.We can now define the voltage mapping as follows, combining S and R . Definition
Let m ∈ N , I ∈ R and τ ∈ C ( R , R m ) . The m -spike voltage mapping, V m , is the operator defined as (3.7) V m τ = I + (cid:88) j ∈ N m ( Sτ j + Rτ j ) , where S and R are given in Definition . By construction, the voltage operator characterises m -spike CIFM solutions, asthe following proposition shows. Proposition
Let m ∈ N , I ∈ R . An m -spike CIFM solution exists if, andonly if, there exists τ ∈ C ( R , R m ) such that V m τ = 1 , in F τ , (3.8) V m τ < , in R \ F τ (3.9) Proof.
The statement follows by setting v m ( x, t ) = ( V m τ )( x, t ) and applying thedefinition of the voltage mapping, Equation (3.7).In some cases it is useful to replace the threshold conditions (3.3) and (3.8) byequivalent conditions involving left limits of the voltage function and mapping, re-spectively, as specified by the following result. Corollary v m ). Under the hypotheses of Lemma , v m = 1 in F τ if, and only if, lim µ → + v m ( x, τ i ( x ) − µ ) = v m ( x, τ i ( x ) − ) = 1 for all ( i, x ) ∈ N m × R .Proof. Using Lemma 3.4 one can show that v m ( x, τ i ( x )) = v m ( x, τ i ( x ) − ), as wederive in Appendix B.Proposition 3.6 implies that the voltage of an m -spike solution can be computedfor any ( x, t ) ∈ R once the firing functions τ are known. The spatio-temporal profileof an m -spike solution is determined entirely by its firing functions. This aspect,which underlies the formal evolution equation (3.1) and the literature which analysesit, is a key part of what follows and, as we shall see below, it also suggests a naturalway to compute travelling waves, and determine their linear stability. A first step inthis direction is the definition of travelling waves via the voltage mapping. m ). Following Proposition 3.6, wecan capture travelling waves with m spikes (TW m ) using the voltage mapping, and aset of parallel firing functions. Henceforth, we will assume without loss of generalitythat the propagating speed of the wave is positive: for any wave with c >
0, thereexists a wave with speed − c , and the wave profiles related by the transformation x → − x .4 D. AVITABILE, J. L. DAVIS, K. C. A. WEDGWOOD c = 0 . . ⇠ . ⇠
00 11 cT cT cT cT cT . . . c = 0 . ⌫ ( ⇠ )
The statement follows from evaluating the operators S and R at the fir-ing functions { τ j ( x ) = x/c + T j } j ∈ N m , and operating a change of variables. SeeAppendix C.Proposition 3.9 shows that the travelling wave profile is completely determinedby the vector ( c, T ) ∈ R > × R m , that is, ( c, T ) is a vector of coarse variables for thetravelling wave. In the discrete model we introduced an auxiliary spatially-extendedvariable for the model, the synaptic input { s i ( t ) } i defined in (2.5). In the continuummodel, the corresponding variable is the function s m ( x, t ) = (cid:80) j ∈ N m ( Sτ j )( x, t ), whichin a TW m satisfies s m ( x, t ) = σ m ( ct − x ; c, T ), with(3.11) σ m ( ξ ; c, T ) = 1 c m (cid:88) j =1 (cid:90) ∞ w ( y − ξ + cT j ) p ( y/c ) dy. Proposition 3.9 suggests a simple way tocompute a TW m , by determining its m + 1 coarse variables ( c, T ), as a solution tothe following coarse problem : UMPS AND WAVES IN SPIKING NETWORKS Problem m ). Find ( c, T ) ∈ R > × R m such that T < · · · < T m and T = 0 , (3.12) ν m ( cT − i ; c, T ) = 1 , for i ∈ N m , (3.13) ν m ( ξ ; c, T ) < , on R \ ∪ j ∈ N m { cT − j } . (3.14)Equation (3.13) of the coarse problem imposes that the travelling wave profile crossesthe threshold 1 when ξ → cT − j . As expected, if ν m is a travelling wave profile, thenso is ν m ( ξ + ξ ) for any ξ ∈ R ; Equation (3.12) fixes the phase of the travelling wave,by imposing that the profile crosses threshold as ξ → − .If m = 1, Equations (3.12)–(3.13) of the coarse problem reduce to a compatibilitycondition for the speed c , c (cid:90) −∞ (cid:90) ∞ exp( s ) w (cid:0) c ( y − s ) (cid:1) p ( y ) dy ds = I − , which implicitly defines an existence curve for TW in the ( c , I )-plane. This result isin agreement with what was found in [64, 31]. Existence curves in other parametersare also possible, and are at the core of the numerical bifurcation analysis presentedin detail in the sections below.For m >
1, the coarse problem must be solved numerically. A simple solutionstrategy is to find a candidate solution using Newton’s method for the system of m +1transcendental equations (3.12)–(3.13), with ν m given by Proposition 3.9, and withinitial guesses estimated from direct simulation of the discrete model with large n ,or from a previously computed coarse vector. The candidate solution can then beevaluated at arbitrary ξ ∈ R , hence it is accepted if (3.13) holds on a spatial gridcovering [ − L, L ] ⊂ R , with L (cid:29)
1. In passing, we note that this procedure is consid-erably cheaper than a standard travelling wave computation for PDEs, which requiresthe solution of a boundary value problem, and hence a discretisation of differentialoperators on R . Depending on the particular choice of α and w , the profile ν m iseither written in closed form, as is the case for the choices (2.3), or approximatedusing standard quadrature rules.A concrete calculation is presented in Figure 4, where we show travelling waveprofiles and speeds of a TW and a TW . In passing, we note that the synapticprofile of a TW m at a given time is similar to a bump, but displays modulationsat the core (visible in Figure 4), as predicted by the Heaviside switches in (3.11).Travelling waves with a large number of spikes, such as these ones, have not beenaccessible to date. Remark ν m ( cT − j ) = 1 propagate with positive speed, and this does not contradict the numerical simulations in Figure 3,where solutions profiles with v m ( x, τ j ( x ) − ) = 1 propagate with negative speed. Thisis a consequence of choosing ξ = ct − x (as in [64]), hence initial conditions for the timesimulations are obtained by reflecting ν m about the y axis, since v m ( x,
0) = ν m ( − x ).
4. Wave Stability.
The time simulations in section 2 demonstrate that, forsufficiently large values of β , travelling waves with a variable number of spikes coexistand are stable. It is natural to ask whether these waves destabilise as β , or anyother control parameter of the model, is varied. An example of a prototypical waveinstability is presented in Figure 5 for TW : a travelling wave is computed solving6 D. AVITABILE, J. L. DAVIS, K. C. A. WEDGWOOD τ + ϕ
Assume Hypothesis , andlet ( c, T ) be the coarse variables of a TW m with firing functions τ . Further, let L bethe linear operator defined by Lϕ = (cid:0) ( Lϕ ) , . . . , ( Lϕ ) m (cid:1) , where ( Lϕ ) i = (cid:88) j ∈ N m ( ϕ i − ϕ j )1 j
To prove that L : C η ( R , C m ) → C η ( R , C m ) is well defined and bounded,one shows that for any ϕ ∈ C η ( R , C m ) there exists a positive constant κ m,η such that (cid:107) Lϕ (cid:107) C m,η ≤ κ m,η (cid:107) ϕ (cid:107) C m,η . A proof is obtained by obtaining preliminary estimates on ψ ij , and | ( Lϕ )( x ) | ∞ , and relies on the decaying properties of w , stated in Hypothe-sis 3.1. Part 2 is proved by setting u = τ + εϕ ∈ C η ( R , R m ), linearising the mapping u (cid:55)→ V m u − | F u around τ , and exploiting properties of m -spikes solutions outlined inthe subsection 3.2. A proof is given in Appendix D. Remark L depends on the coarse variables ( c, T ),albeit we omit this dependence for notational simplicity.We are now ready to define linear stability for a TW m , which we adapt from[12]. Intuitively, we compare the firing set F τ of a TW m with the firing set F τ + ϕ ofa perturbed m -spike solution with (cid:107) ϕ (cid:107) C m,η (cid:28)
1, for which ϕ satisfy (4.1) to leadingorder. If the sets F τ and F τ + ϕ are close around t = 0 and remain close for all positivetimes, we deem the wave linearly stable. With reference to Figure 5(c), we observethat, when TW m crosses the axis t = 0, each one of its firing functions τ i is perturbedby an amount ϕ i ( − cT i ). Roughly speaking, a TW m is linearly stable if ϕ i ( − cT i ) beingsmall implies that ϕ i ( x ) stays small for all x ∈ ( − cT i , ∞ ) and i ∈ N m . If a wave islinearly stable and all ϕ i decay to 0 as x → ∞ we say that the wave is asymptoticallylinearly stable. More precisely: Definition m ). A TW m with coarse varaibles ( c, T ) is linearly stable to perturbations ϕ if ϕ ∈ ker L , and for each ε > there exists δ = δ ( ε ) > , such that if | ϕ i ( − cT i ) | < δ , then | ϕ i ( x ) | < ε for all ( i, x ) ∈ N m × ( − cT i , ∞ ) .A TW m is asymptotically linearly stable to perturbations ϕ if is linearly stable toperturbations ϕ and | ϕ ( x ) | ∞ → as x → ∞ . We have seen that a TW m can be constructed by solving a nonlinear problem inthe unknowns ( c, T ). The following lemma, which is the central result of this section,establishes that linear stability of a TW m with respect to exponential perturbationsof the firing functions can also be determined by finding roots of a ( c, T )-dependent,complex-valued function.8 D. AVITABILE, J. L. DAVIS, K. C. A. WEDGWOOD
Lemma m stability). Assume Hypothesis , let ( c, T ) be coarse vari-ables of a TW m , and let D a,b = { z ∈ C : a ≤ Re z ≤ b } . Further, let E be thecomplex-valued function (4.2) E : D − η,η → C , z (cid:55)→ det[ D − M ( z )] , where M ∈ C m × m , D = diag( D , . . . , D m ) ∈ R m × m , are the matrices with elements M ij ( z ) = e T ji (cid:20) j
5. Bifurcation structure of travelling waves.
The pseudo-arclength con-tinuation routines developed in [67, 2] have been used to compute solutions to Prob-lem 3.10, continue waves in parameter space, and investigate their stability. A TW m isconstructed by solving Problem 3.10 in the corase variables ( c, T ) ∈ R > × R m , whichis sufficient to reconstruct the wave profile (3.10), and the corresponding synaptic pro-file (3.11); in addition, starting from a solution to Problem 3.10, the linear asymptoticstability of a TW m is determined by finding roots of the ( c, T )-dependent nonlinearfunction E defined in (4.2).Figure 6 shows the bifurcation structure of TW , which is common to most trav-elling waves found in the model. The simulations in Section 2 suggest to take thesynaptic timescale parameter β as the principal continuation parameter. We use thewavespeed c as solution measure. A branch of solutions originates from a grazing point(G, see below for a more detailed explanation) and it is initially stable, before desta-bilising at a sequence of oscillatory bifurcations (HB –HB ), as seen in Figure 6(a). Inpassing, we note that there exists a second, fully unstable, branch of TW solutionscharacterised by a slower speed and a smaller width. This branch, which we omitfrom the bifurcation diagrams for simplicity, also originates at a grazing point. UMPS AND WAVES IN SPIKING NETWORKS c
57, they are fully unstablefor this parameter set. The diagrams provide evidence that there exists unstablewaves with arbitrarily many spikes and vanishingly small speed. It seems thereforenatural to postulate a relationship between these waves and the structures found by
UMPS AND WAVES IN SPIKING NETWORKS m branch with m > ξ = 0by construction (see Problem 3.10), while its righmost spike is at ξ m = cT m , whichis therefore a proxy for the wave’s width . Figure 8(a) shows c and T m , computed atthe grazing points, as functions of m : we find c = O ( m ), T m = O ( m − ), thereforewe expect the sequence { ξ m } m ∈ N to converge to a finite value ξ ∗ as m → ∞ . Thesedata show that, as the wavespeed tends to zero, the growing number of spikes aredistributed in a fixed interval [0 , ξ ∗ ] resembling therefore a stationary bump of width ξ ∗ . This conclusion is also supported by an inspection of the distribution of spikes { ξ i,m } i as m → ∞ . Instead of looking directly at this distribution, we plot theinverse inter-spike time 1 / ( T i +1 − T i ), as a function of ξ i , for the grazing point ofTW . Laing and Chow [54] call this quantity the gain function , and propose itas an emerging firing rate function for the asynchronous bumps of the DIFM [54].Importantly, Laing and Chow show that the gain function for asynchronous bumpsresembles the one of a neural network rate equation, but it is an emergent propertyof the DIFM. Figure 8(b) shows that gain function by Laing and Chow is in excellentagreement with the one obtained for TW (the slowest wave we computed in theCIFM), confirming that, from a macroscopic viewpoint, asynchronous bumps couldbe understood in the limit m → ∞ of TW m .We have further investigated the bump attractor state, in relation to the TW m :the analysis of the continuum model, in the region of parameter space where thebump attractor is observed, predicts the coexistence of the trivial attracting solution v ( x, t ) ≡ I , with arbitrarily slow, unstable waves, whose spatial profile approximatesthe one of a bump. We simulated the DIFM with n = 5000, initialising the model froman unstable travelling wave of the CIFM, TW , and estimated the instantaneousspeed c ( t ) of the numerical DIFM solution at q time points { t k : k ∈ N q } , using a levelset of the synaptic profile and finite differences, as follows: z ( t ) = max { x ∈ S : s ( x, t ) = 0 . } , c k = ( z ( t k ) − z ( t k − )) / ( t k − t k − ) , k ∈ N q . A CIFM travelling wave solution corresponds to a constant c : when the DCLM solu-tion displays a wave for large n , the sequence { c k } k converges to a constant value, ifone disregards small oscillations due to the finite n , and which vanish as n → ∞ . Onthe other hand we expect that no differentiable function c ( t ) exists for a bump attrac-tor, albeit some information may be contained in the mean, c , standard deviation, σ c ,and extrema, c min , c max of the deterministic scalar c k (5.1) ¯ c = 1 q (cid:88) k ∈ N q c k , σ c = 1 q − (cid:88) k ∈ N q ( c k − ¯ c ) , c min = min k ∈ N q c k , c max = max k ∈ N q c k , which we computed for various values of β , and superimposed on the bifurcationdiagram of the CIFM model, in Figure 9(a): we plot ¯ c (purple dots) and two intervalestimators, [¯ c − σ c , ¯ c + σ c ] (dark purple shade) and [ c min , c max ] (light purple shade);we recall that the CFIM admits branches of waves with positive and negative speed,both plotted in the figure. The bump attractor is characterised by a zero-mean speed,and interval estimators which grow with β , indicating that the instantaneous speedundergoes larger fluctuations as the bump meanders. We observe that the interval Recall that c is also a function of m , but we omit this dependence for ease of notation D. AVITABILE, J. L. DAVIS, K. C. A. WEDGWOOD
Fig. 9 . (a): Mean instantaneous speed ( ¯ c in (5.1) , purple dots) and interval estimators ( [¯ c − σ c , ¯ c + σ c ] and [ c min , c max ] , dark and light purple shades, respectively) in direct simulations of theDIFM, superimposed on TW m branches of the CIFM (an inset of Figure (b), which has beenreflected by the c = 0 axis to signpost waves with negative speed). The bump attractor is characterisedby ¯ c ≈ , and speed fluctuations which grow with β , while echoing the oscillatory instabilities of thebranch (red segments). (b): exemplary solutions in (a) displaying an initial advection, followed bya bump attractor (1,2) or a stable wave (3). (c): Histograms of the instantaneous speed in selectedtime intervals, indicated by the lateral coloured bars in (b). We observe transitions through weaklyunstable waves (orange interval in (3), blue intervals in (1, 2)), corresponding to sharp peaks in thehistograms. UMPS AND WAVES IN SPIKING NETWORKS m branches, as β increases.This behaviour is robust, for low and medium values of β , albeit the fine detailsof the dynamics depend on initial conditions. Figure 9(b) shows 3 examples whoseestimated speeds appear also in Figure 9(a). The space-time plots display an initialadvection, followed by a bump attractor, or a stable travelling wave. To highlight thetransitions, we computed histograms of c k in selected time intervals, indicated by blue,orange, yellow, and purple bars in Figure 9(b). The histograms provide numericalevidence of transitions near weakly unstable waves (orange intervals/histograms inexample 3, blue intervals/histograms in examples 1 and 2), which manifest themselvesvia sharp peaks. In addition to the waves studied thus far, we foundby direct simulation waves whose firing functions are split into well-separated groups,that is, firing functions in the same group are closer to each another than they are tothose in other groups, see Figure 10. We call these structures composite waves , as theymay be formed via the interaction of travelling waves with various numbers of spikes.As in other non-smooth dynamical systems [44], we expect that these solutions havediscontinuities that are rearranged with respect to a TW m .For illustrative purposes, we denote a composite waves with k ∈ N groups byTW m + . . . + TW m k , where { m i } ki =1 is a sequence of positive integers specifyingthe number of spikes in each group. There are constraints for the groups, dictated bydynamical considerations: for instance a TW + TW cannot exist, because a TW ,taken in isolation, is faster than a TW . The construction of asymptotic profiles andcomputation of linear stability for composite waves follow in the same way as definedin Sect. 3 and Sect. 4.In Figure 10(a), we show a selection of of composite waves near the TW branch.Roughly speaking, the wave profile along each depicted branch comprises a TW asits leading group, followed by two additional spike group that collectively form acompound satisfying the travelling wave conditions (e.g. branch 1 combines a TW ,a TW and a TW ). The branches of composite waves are separate from each otherand from the previously computed TW m branches in Figure 7, however, all branchespossess a bifurcation structure similar to the one of the TW m discussed in the pastsection. Moreover, we see that the magnitude of the speed of the composite wave isbounded above by the magnitude of the speed of the group at the leading edge of thewave (the slowest wave, TW in this case).Direct numerical simulation highlights that composite waves can be formed fromthe interaction of multi-spike waves as shown in the left panel of Figure 10(b). Here wechoose an initial condition with well separated TW , TW and TW profiles. Initially,these separated structures travel with different speeds (TW being the fastest andTW the slowest, in line with what was found in Figure 7(a)). After a transient,the waves come closer and form a compound (the composite wave), with a commonintermediate speed. The dynamics of composite waves depend greatly on the initialconditions: in the right panel of Figure 10(b), we see that an initial condition inwhich a TW lies between another TW and a TW leads to the extinction of theintermediate wave resulting in a composite wave with a total of 4 spikes.Composite waves also result from the collision between waves and wanderingbumps (Figure 10(c), left panel). Here, we see a transition of two bump states into acomposite wave that is compounded with a pre-existing TW . The interaction withthe TW causes the left-most bump to visit the branches of travelling wave solutions6 D. AVITABILE, J. L. DAVIS, K. C. A. WEDGWOOD (a) TW TW TW t x −12 10 450 t x −12 10(b) TW -3 9 xt bumpbump bump xt TW bump bumpTW + TW TW + TW + TW (c) TW + TW + TW TW + TW TW TWTW
12 14.6 HB TW c β
11 22 33 4 45 56 6 Fig. 10 . (a) Bifurcation diagram of selected composite waves. The red curve is a TW branch,as computed in Figure . The blue curves are branches of composite waves, featuring an approximateTW at the front of the wave. The composite waves are slightly slower than TW . The diagramshows selected profiles at the first oscillatory bifurcation points. (b) Examples of composite wavesobtained via collisions of multi-spike waves. (c) Collisions between m -spike propagating structuresand wandering bumps generate composite waves (left) or bump repulsion (right), depending on initialconditions. Simulations in panels (b)–(c) have a lattice spacing of ∆ x = 2 L/n = 0 . . whereupon the combined state settles on a stable TW + TW branch. This process isrepeated for the right-most bump, giving rise to an overall TW + TW + TW wave.In the right panel of the Figure 10(c), we see that the same kind of collision can insteadresult in the wave packet transitioning to a wandering bump itself, highlighting thedependence of the formation of composite waves on initial conditions. In this scenario,the bump state does not visit a stable travelling wave branch and so only transientlyadopts a weakly unstable wave profile before returning to a bump attractor state.
6. Conclusions.
We have provided evidence that the relationship between bumpattractors and travelling waves in a classical network of excitable, leaky integrate-and-fire neurons bears strong similarities to the one between complex spatiotemporalpatterns and waves at the onset of pipe turbulence. We made analytical and numerical
UMPS AND WAVES IN SPIKING NETWORKS m -spike waves starting frommild solutions to the formal evolution equation (3.1), and derive a relatively simpleexpression for the wave profile (3.10). While this approach may be harder to carry outin more detailed spiking models, the general idea of a relationship between localisedwaves and bumps in spiking networks could be investigated, by direct simulations, inmore realistic networks (spiking or not).An important open question concerns the definition of (3.1) and, more gener-ally, of spatially-continuous spiking networks, as dynamical systems posed on func-tion spaces. This problem has been circumvented here by defining a suitable classof solutions, introducing the voltage mapping, and then providing proofs of its rele-vance to the construction and stability of multiple-spike waves. We believe that a fulldynamical-systems characterisation of similar models will be a key ingredient to un-cover further links between localised waves and bumps in complex, spatially-extendedthreshold networks. Acknowledgments.
We are grateful to Stephen Coombes, Predrag Cvitanovi´c,Gregory Faye, Joel Feinstein, John Gibson, Joost Hulshof, and Edgar Knobloch. forinsightful discussions.
Appendix A. Proof of Lemma 3.4.
Lemma
A.1 (Re-statement of Lemma 3.4).
If Hypothesis holds, then for theoperators S , R in Definition we have S : C ( R ) → BC ( R ) and R : C ( R ) → B ( R ) ,respectively.Proof. Fix u ∈ C ( R ). The real-valued function z (cid:55)→ exp( − z ) H ( z ) is bounded in R , hence Ru ∈ B ( R ). To prove the result on S , we define the functions ψ ( y, t ) = (cid:90) t ∞ e z − t α ( z − u ( y )) dz, s ( x, t ) = ( Su )( x, t ) = (cid:90) ∞−∞ w ( x − y ) ψ ( y, t ) dy, and set K α = (cid:107) α (cid:107) ∞ , K w = (cid:107) w (cid:107) L ( R ) , whose existence is guaranteed by Hypothe-sis 3.1. The function s is bounded on R because | ψ | ≤ K α on R , hence | s ( x, t ) | ≤ K α (cid:90) ∞−∞ | w ( x − y ) | dy ≤ K α K w . In order to prove the continuity of s , it is useful to first show that ψ is continuousin t on R , uniformly in y . This claim is proved by noting that, for any y, t, τ ∈ R with8 D. AVITABILE, J. L. DAVIS, K. C. A. WEDGWOOD t < τ we have | ψ ( y, τ ) − ψ ( y, t ) | = (cid:12)(cid:12)(cid:12)(cid:12) (cid:90) t −∞ [ e z − τ − e z − t ] α ( z − u ( y )) dz + (cid:90) τt e z − τ α ( z − u ( y )) dz (cid:12)(cid:12)(cid:12)(cid:12) ≤ K α (cid:90) t −∞ e z − t − e z − τ dz + K α ( τ − t )= K α (cid:0) − e − ( τ − t ) + τ − t (cid:1) , which, combined with a similar argument for τ < t , leads to(A.1) | ψ ( y, τ ) − ψ ( y, t ) | ≤ K α (cid:0) − e −| τ − t | + | τ − t | (cid:1) , for all y, t, τ ∈ R . We prove the continuity of s by showing that | s ( ξ, τ ) − s ( x, t ) | → ξ, τ ) → ( x, t ).We consider the following inequality(A.2) | s ( ξ, τ ) − s ( x, t ) | ≤ | s ( ξ, τ ) − s ( x, τ ) | + | s ( x, τ ) − s ( x, t ) | , and we note that the second term in the right-hand side of (A.2) can be made ar-bitrarily small as ( ξ, τ ) → ( x, t ), owing to (A.1). Therefore it suffices to show thatthe first term in the right-hand side of (A.2) can also be made arbitrarily small as( ξ, τ ) → ( x, t ), that is, we must show that for any ε > δ > | s ( ξ, τ ) − s ( x, τ ) | < ε . To prove this statement, we use the boundedness of ψ and achange of variables in the integral to obtain the estimate | s ( ξ, τ ) − s ( x, τ ) | ≤ (cid:90) ∞−∞ | w ( y − ξ ) − w ( y − x ) | | ψ ( y, τ ) | dy ≤ K α (cid:90) ∞−∞ | w ( y + x − ξ ) − w ( y ) | dy, therefore, for any X > | s ( ξ, τ ) − s ( x, τ ) | ≤ K α (cid:90) − X −∞ | w ( y + x − ξ ) − w ( y ) | dy + K α (cid:90) X − X | w ( y + x − ξ ) − w ( y ) | dy + K α (cid:90) ∞ X | w ( y + x − ξ ) − w ( y ) | dy := K α ( I + I + I ) . We bound I as follows I ≤ (cid:90) ∞ X | w ( y + x − ξ ) | dy + (cid:90) ∞ X | w ( y ) | dy = (cid:90) ∞ X + x − ξ | w ( y ) | dy + (cid:90) ∞ X | w ( y ) | dy ≤ (cid:90) ∞ X −| x − ξ | | w ( y ) | dy and a similar reasoning gives an identical bound for I , I ≤ (cid:90) ∞ X −| x − ξ | | w ( y ) | dy. UMPS AND WAVES IN SPIKING NETWORKS
X > | s ( ξ, τ ) − s ( x, τ ) | ≤ K α (cid:90) ∞ X −| x − ξ | | w ( y ) | dy + K α (cid:90) X − X | w ( y + x − ξ ) − w ( y ) | dy. We now fix ε, δ >
0. Since w ∈ L ( R ), we can pick X so that | x − ξ | < δ ⇒ K α (cid:90) ∞ X −| x − ξ | | w ( y ) | dy < ε . Furthermore, by continuity of w there exists δ > | x − ξ | < δ ⇒ K α (cid:90) X − X | w ( y + x − ξ ) − w ( y ) | dy ≤ ε , hence for any ε > δ = min( δ , δ ) > | s ( ξ, τ ) − s ( x, τ ) | < ε ,which implies the continuity of s . We conclude that S : C ( R ) → BC ( R ). Appendix B. Proof of Corollary 3.7.
Corollary
B.1 (Restatement of Corollary 3.7).
Under the hypotheses of Lemma , v m = 1 in F τ if, and only if, lim µ → + v m ( x, τ i ( x ) − µ ) = v m ( x, τ i ( x ) − ) = 1 for all ( i, x ) ∈ N m × R .Proof. The condition v m ( x, t ) = 1 for ( x, t ) ∈ F τ is equivalent to v m ( x, τ i ( x )) = 1 ( i, x ) ∈ N m × R . From the defintions of S and R , and the continuity Sτ i on R (see Lemma 3.4) wehave, for all ( i, x ) ∈ N m × R v m ( x, τ i ( x )) = I + (cid:88) j ∈ N m ( Sτ j )( x, τ i ( x )) − (cid:88) j
Proposition
C.1 (TW m profile). A TW m with speed c satisfies ( V m τ )( x, t ) = ν m ( ct − x ; c, T ) , and its ( c, T ) -dependent travelling wave profile ν m is given by (C.1) ν m ( ξ ; c, T ) = I − (cid:88) j ∈ N m exp (cid:18) − ξ − cT j c (cid:19) H (cid:18) ξ − cT j c (cid:19) + 1 c (cid:88) j ∈ N m (cid:90) ξ −∞ exp (cid:18) z − ξc (cid:19) (cid:90) ∞ w ( y − z + cT j ) p ( y/c ) dy dz. Proof.
We set τ j ( x ) = x/c + T j for j ∈ N n . The first sum in (C.1) is immediate,as ( Rτ j )( x, t ) = exp (cid:18) cT j − ct + xc (cid:19) H (cid:18) ct − x − cT j c (cid:19) D. AVITABILE, J. L. DAVIS, K. C. A. WEDGWOOD
The second sum is obtained as follows:( Sτ j )( x, t ) = (cid:90) t −∞ (cid:90) ∞−∞ exp( s − t ) w ( x − y ) α ( s − τ j ( y )) dy ds = (cid:90) t −∞ (cid:90) ∞ exp( s − t ) w ( x + y − cs + cT j ) p ( y/c ) dy ds = 1 c (cid:90) ct − x −∞ (cid:90) ∞ exp (cid:18) s + x − ctc (cid:19) w ( y − s + cT j ) p ( y/c ) dy ds. Appendix D. Proof of Lemma 4.1.
Lemma
D.1 (Restatement of Lemma 4.1).
Assume Hypothesis , and let ( c, T ) be the coarse variables of a TW m with firing functions τ . Further, let L be the linearoperator defined by Lϕ = (cid:0) ( Lϕ ) , . . . , ( Lϕ ) m (cid:1) , where ( Lϕ ) i = (cid:88) j ∈ N m ( ϕ i − ϕ j )1 j
1. By main hypothesis u and τ are firing functions of two distinct m -spike CIFM solutions. We claim that thisimplies Equation (4.1). Indeed, since u is a firing function, then V m u = 1 on F u , thatis(D.6) 1 = I + (cid:88) j ∈ N m (cid:16) ( Su i )( x, u j ( x )) + ( Ru j )( x, u i ( x )) (cid:17) , ( i, x ) ∈ N m × R . We obtain( Su j )( · , u i ) = (cid:90) ∞−∞ w ( · − y ) (cid:90) −∞ e z α ( z + u i − u j ( y )) dz dy = ( Sτ j )( · , τ i ) + ε (cid:90) ∞−∞ w ( · − y ) (cid:0) ϕ i − ϕ j ( y ) (cid:1) (cid:90) −∞ e z α (cid:48) ( z + τ i − τ j ( y )) dz dy + O ( ε ) , where we have denoted by α (cid:48) = p (cid:48) H + pδ the distributional derivative of α . We now2 D. AVITABILE, J. L. DAVIS, K. C. A. WEDGWOOD manipulate the integral in the previous equation as follows (cid:90) ∞−∞ w ( · − y ) (cid:0) ϕ i − ϕ j ( y ) (cid:1) (cid:90) −∞ e z α (cid:48) ( z + τ i − τ j ( y )) dz dy = (cid:90) ∞−∞ w ( · − y ) (cid:0) ϕ i − ϕ j ( y ) (cid:1) (cid:90) τ i − τ j ( y ) −∞ e z + τ i − τ j ( y ) α (cid:48) ( z ) dz dy = (cid:90) ∞−∞ w ( · − y ) (cid:0) ϕ i − ϕ j ( y ) (cid:1) (cid:90) τ i − τ j ( y ) −∞ e z + τ i − τ j ( y ) (cid:0) p (cid:48) ( z ) H ( z ) + p ( z ) δ ( z ) (cid:1) dz dy = (cid:90) τ − j ( τ i ) −∞ w ( · − y ) (cid:0) ϕ i − ϕ j ( y ) (cid:1) e τ j ( y ) − τ i (cid:16) p (0) + (cid:90) τ i − τ j ( y )0 e z p (cid:48) ( z ) dz (cid:17) dy , hence for all i, j ∈ N m we obtain( Su j )( · , u i ) = ( Sτ j )( · , τ i ) + εe T ji (cid:90) ∞−∞ e − y/c w ( y ) ψ ij ( y ) (cid:2) ϕ i − ϕ j ( · − y ) (cid:3) dy + O ( ε )For the reset operator, we obtain, for all i, j ∈ N m Ru j ( · , u − i ) = − lim κ → + exp( − u i + κ + u j ) H ( u i − κ − u j )= Rτ j ( · , τ − i ) + ε lim κ → + exp( − T ij + κ )( ϕ i − ϕ j ) H ( T ij ) + O ( ε )= Rτ j ( · , τ − i ) + ε exp( T ji )( ϕ i − ϕ j )1 j
Lemma
E.1 (Restatement of Lemma 4.4).
Assume Hypothesis , let ( c, T ) becoarse variables of a TW m , and let D a,b = { z ∈ C : a ≤ Re z ≤ b } . Further, let E bethe complex-valued function (E.1) E : D − η,η → C , z (cid:55)→ det[ D − M ( z )] , where M ∈ C m × m , D = diag( D , . . . , D m ) ∈ R m × m , are the matrices with elements M ij ( z ) = e T ji (cid:20) j
0. From part 1 we deduce thatthere exists Φ ∈ ker[ D − M ( λ )], such that ϕ ( x ) = Φ e λx + Φ ∗ e λ ∗ x ∈ ker L . We notethat | Φ | ∞ can be fixed to an arbitrary nonzero constant, and we bound ϕ i as follows(E.2) | ϕ i ( x ) | ≤ | Φ | ∞ e µx ≤ K | Φ | ∞ , K = 2 max j ∈ N m e − cT j µ , ( i, x ) ∈ N m × [ − cT i , ∞ ) . We now fix ε >
0. The bound (E.2) and the choices δ = ε and | Φ | ∞ < ε/K implythat TW m is linearly stable, according to Definition 4.3. Using again (E.2) we obtainlim x →∞ | ϕ ( x ) | ∞ ≤ | Φ | ∞ lim x →∞ e µx = 0therefore TW m is linearly asymptotically stable. Appendix F. Two-parameter continuation of bifurcations.
Grazingpoints are found generically as a secondary control parameter, say γ , is varied. It ispossible to perform a 2-parameter continuation of the grazing point in the ( β, γ )-planeby continuing in γ solutions the following problem: Problem
F.1 (Grazing point computation). Find ( c, T , . . . , T m , T G , β G ) ∈ R > × R m +2 such that T < · · · < T m < T G and T = 0 , (F.1) ν ([ cT i ] − ; c, T , . . . , T m , β G ) = 1 , i ∈ N m , (F.2) ν ( cT G ; c, T , . . . , T m , β G ) = 1 , (F.3) ν (cid:48) ( cT G ; c, T , . . . , T m , β G ) = 0 . (F.4)We note that we have exposed the dependence of ν and ν (cid:48) on β in the previousproblem. Unlike the numerical continuation presented in the main text, here β is afree parameter, which is determined also by Newton’s method.Similarly, we can trace loci of Hopf bifurcations in the ( γ, β )-plane, by solving Problem
F.2 (Computation of Hopf bifurcations). Find ( c, T , . . . , T m , β HB , ω HB ) ∈ R > × R m +2 such that T < · · · < T m and T = 0 , (F.5) ν ([ cT i ] − ; c, T , . . . , T m , β HB ) = 1 , i ∈ N m , (F.6) E ( iω HB ; c, T , . . . , T m , β HB ) = 0 . (F.7)4 D. AVITABILE, J. L. DAVIS, K. C. A. WEDGWOOD c β G (a) m = 1 m = 2 m = 4 m = 8 m = 166 0 .
78 0 .
82 0 .
86 0 . I m = 1 m = 16 (b) c m = 8 β = 2 β = 3 β = 4 β = 5 β = 5.4 − ξ (c) Fig. 11 . Bifurcation diagram showing wave speed branches c for TW m states with m =1 , , , , , with purely excitatory kernel, a = 2 , b = 5 , a = 0 . (a) External current is fixed at I = 0 . and solutions continued in β , here m = 1 , , , branches switch stability at a fold, but for m = 16 the branch terminates at a grazing point as seen in the lateral-inhibition case previously.(b) The synaptic processing rate is fixed at β = 4 and the waves continued in I . Note that in bothcases a greater m leads to larger speeds. (c) Synaptic profiles for TW solutions labelled (1–5) in(a). Smaller β values result in more localised profiles. REFERENCES[1]
S.-i. Amari , Dynamics of pattern formation in lateral-inhibition type neural fields , BiolologicalCybernetics, 27 (1977), pp. 77–87.[2]
D. Avitabile , Numerical computation of coherent structures in spatially-extended sys-tems , May 2020, https://doi.org/10.5281/zenodo.3821169, https://doi.org/10.5281/zenodo.3821169.[3]
D. Avitabile and K. C. A. Wedgwood , Macroscopic coherent structures in a stochasticneural network: from interface dynamics to coarse-grained bifurcation analysis , Journal ofMathematical Biology, (2017), pp. 1–44, https://doi.org/10.1007/s00285-016-1070-9.UMPS AND WAVES IN SPIKING NETWORKS [4] J. Baladron, D. Fasoli, O. Faugeras, and J. Touboul , Mean-field description and propa-gation of chaos in networks of hodgkin-huxley and fitzhugh-nagumo neurons , The Journalof Mathematical Neuroscience, 2 (2012), p. 10.[5]
D. Barkley , Simplifying the complexity of pipe flow , Physical Review E, 84 (2011), p. 016309.[6]
D. Barkley , Pipe flow as an excitable medium , Revista Cubana de F´ısica, 29 (2012), pp. 1–27.[7]
D. Barkley , Theoretical perspective on the route to turbulence in a pipe , Journal of FluidMechanics, 803 (2016).[8]
D. Barkley, B. Song, V. Mukund, G. Lemoult, M. Avila, and B. Hof , The rise of fullyturbulent flow , Nature, 526 (2015), pp. 550–553.[9]
C. Bick, M. Goodfellow, C. R. Laing, and E. A. Martens , Understanding the dynamics ofbiological and neural oscillator networks through mean-field reductions: a review , Journalof Mathematical Neuroscience, 10 (2020), pp. 1–43.[10]
P. C. Bressloff , Mean-field theory of globally coupled integrate-and-fire neural oscillatorswith dynamic synapses , Physical Review E, 60 (1999), p. 2160.[11]
P. C. Bressloff , Synaptically generated wave propagation in excitable neural media , PhysicalReview Letters, 82 (1999), pp. 2979–2982.[12]
P. C. Bressloff , Traveling waves and pulses in a one-dimensional network of excitableintegrate-and-fire neurons , Journal of Mathematical Biology, 40 (2000), pp. 169–198.[13]
P. C. Bressloff , Waves in neural media , Lecture Notes on Mathematical Modelling in theLife Sciences, Springer, New York, (2014).[14]
P. C. Bressloff and S. Coombes , Dynamics of strongly-coupled spiking neurons. , NeuralComputation, 12 (2000), pp. 91–129.[15]
A. N. Burkitt , A review of the integrate-and-fire neuron model: I. Homogeneoussynaptic input , Biological Cybernetics, 95 (2006), pp. 1–19, https://doi.org/10.1007/s00422-006-0068-6.[16] ´A. Byrne, D. Avitabile, and S. Coombes , Next-generation neural field model: the evolutionof synchrony within patterns and waves , Physical Review E, 99 (2019), p. 012313.[17]
M. Camperi and X.-J. Wang , A model of visuospatial working memory in prefrontal cortex:recurrent network and cellular bistability , Journal of Computational Neuroscience, 5 (1998),pp. 383–405.[18]
A. Compte, N. Brunel, P. S. Goldman-Rakic, and X.-J. Wang , Synaptic mechanismsand network dynamics underlying spatial working memory in a cortical network model ,Cerebral Cortex, 10 (2000), pp. 910–923.[19]
A. Compte, M. V. Sanchez-Vives, D. A. McCormick, and X.-J. Wang , Cellular and net-work mechanisms of slow oscillatory activity (¡ 1 Hz) and wave propagations in a corticalnetwork model , Journal of Neurophysiology, 89 (2003), pp. 2707–2725.[20]
S. Coombes, P. beim Graben, and R. Potthast , Tutorial on Neural Field Theory , in NeuralFields, Springer, Berlin, Heidelberg, Berlin, Heidelberg, 2014, pp. 1–43.[21]
S. Coombes, P. beim Graben, R. Potthast, and J. Wright , Neural fields: Theory andApplications , Springer, 2014.[22]
S. Coombes and P. C. Bressloff , Saltatory waves in the spike-diffuse-spike model of activedendritic spines , Physical Review Letters, 91 (2003), pp. 81–4.[23]
S. Coombes and C. Laing , Pulsating fronts in periodically modulated neural field models ,Physical Review E, 83 (2011), p. 011912.[24]
S. Coombes, H. Schmidt, and D. Avitabile , Spots: breathing, drifting and scattering in aneural field model , in Neural fields, Springer, Heidelberg, Berlin, Heidelberg, 2014, pp. 187–211.[25]
A. Darbyshire and T. Mullin , Transition to turbulence in constant-mass-flux pipe flow ,Journal of Fluid Mechanics, 289 (1995), pp. 83–114.[26]
F. Delarue, J. Inglis, S. Rubenthaler, E. Tanr´e, et al. , Global solvability of a networkedintegrate-and-fire model of mckean–vlasov type , The Annals of Applied Probability, 25(2015), pp. 2096–2133.[27]
M. di Bernardo, C. J. Budd, A. R. Champneys, P. Kowalczyk, A. B. Nordmark, G. O.Tost, and P. T. Piiroinen , Bifurcations in Nonsmooth Dynamical Systems , SIAM Re-view, 50 (2008), pp. 629–701.[28]
J. R. Dormand and P. J. Prince , A family of embedded Runge-Kutta formulae , Journal ofComputational and Applied Mathematics, 6 (1980), pp. 19–26, https://doi.org/10.1016/0771-050x(80)90013-3.[29]
B. Eckhardt, H. Faisst, S. A., and J. Schumacher , Turbulence transition in shear flows. ,in Advances in Turbulence IX: Proc. Ninth European Turbulence Conference, Barcelona,I. P. Castro, P. E. Hancock, and T. G. Thomas, eds., 2002, pp. 701–708.[30]
B. Ermentrout , Neural networks as spatio-temporal pattern-forming systems , Reports on D. AVITABILE, J. L. DAVIS, K. C. A. WEDGWOODProgress in Physics, 61 (1998), p. 353.[31]
B. Ermentrout , The analysis of synaptically generated traveling waves , Journal of Computa-tional Neuroscience, 5 (1998), pp. 191–208.[32]
B. Ermentrout, J. Rubin, and R. Os¸an , Regular traveling waves in a one-dimensional net-work of theta neurons , SIAM Journal on Applied Mathematics, 62 (2002), pp. 1197–1221.[33]
G. B. Ermentrout, S. E. Folias, and Z. P. Kilpatrick , Spatiotemporal pattern formation inneural fields with linear adaptation , in Neural Fields, Springer, Berlin, Heidelberg, Berlin,Heidelberg, 2014, pp. 119–151.[34]
G. B. Ermentrout and N. Kopell , Parabolic bursting in an excitable system coupled with aslow oscillation , SIAM Journal on Applied Mathematics, 46 (1986), pp. 233–253.[35]
G. B. Ermentrout and D. H. Terman , Mathematical Foundations of Neuroscience , vol. 35,Springer Science & Business Media, 2010.[36]
J. M. Esnaola-Acebes, A. Roxin, D. Avitabile, and E. Montbri´o , Synchrony-inducedmodes of oscillation of a neural field model , Physical Review E, 96 (2017), p. 052407.[37]
H. Faisst and B. Eckhardt , Traveling waves in pipe flow , Physical Review Letters, 91 (2003),p. 224502.[38]
S. E. Folias and P. C. Bressloff , Breathing pulses in an excitatory neural network , SIAMJournal on Applied Dynamical Systems, 3 (2004), pp. 378–407.[39]
W. Gerstner, W. M. Kistler, R. Naud, and L. Paninski , Neuronal dynamics: from singleneurons to networks and models of cognition , Cambridge University Press, 2014.[40]
W. Gerstner, J. L. van Hemmen, and J. D. Cowan , What matters in neuronal locking? ,dx.doi.org, 8 (2008), pp. 1653–1676.[41]
J. F. Gibson, J. Halcrow, and P. Cvitanovi´c , Equilibrium and travelling-wave solutions ofplane couette flow , Journal of Fluid Mechanics, 638 (2009), pp. 243–266.[42]
D. Golomb and G. B. Ermentrout , Continuous and lurching traveling pulses in neuronal net-works with delay and spatially decaying connectivity , Proceedings of the National Academyof Sciences, 96 (1999), pp. 13480–13485.[43]
L. R. Gonz´alez-Ram´ırez, O. J. Ahmed, S. S. Cash, C. E. Wayne, and M. A. Kramer , Abiologically constrained, mathematical model of cortical wave propagation preceding seizuretermination , PLoS Computational Biology, 11 (2015), pp. e1004065–34.[44]
A. Granados, L. Alseda, and M. Krupa , The period adding and incrementing bifurcations:from rotation theory to applications , SIAM Review, 59 (2017), pp. 225–292.[45]
B. Hof, A. Juel, and T. Mullin , Scaling of the turbulence transition threshold in a pipe ,Physical Review Letters, 91 (2003), p. 244502.[46]
E. Hopf , A mathematical example displaying features of turbulence , Communications on Pureand Applied Mathematics, 1 (1948), pp. 303–322.[47]
X. Huang, W. C. Troy, Q. Yang, H. Ma, C. R. Laing, S. J. Schiff, and J.-Y. Wu , Spiralwaves in disinhibited mammalian neocortex , J. Neurosci., 24 (2004), pp. 9897–9902.[48]
J. Inglis and J. MacLaurin , A general framework for stochastic traveling waves and patterns,with application to neural field equations , SIAM Journal on Applied Dynamical Systems,15 (2016), pp. 195–234.[49]
E. M. Izhikevich , Dynamical systems in neuroscience , MIT press, 2007.[50]
Z. P. Kilpatrick and B. Ermentrout , Wandering bumps in stochastic neural fields , SIAMJournal on Applied Dynamical Systems, 12 (2013), pp. 61–94.[51]
S. S. Kim, H. Rouault, S. Druckmann, and V. Jayaraman , Ring attractor dynamics in theDrosophila central brain , Science, 356 (2017), pp. 849–853.[52]
J. J. Knierim and K. Zhang , Attractor dynamics of spatially correlated neural activity in thelimbic system , Annual Review of Neuroscience, 35 (2012), pp. 267–285.[53]
C. R. Laing , Exact neural fields incorporating gap junctions , SIAM Journal on Applied Dy-namical Systems, 14 (2015), pp. 1899–1929.[54]
C. R. Laing and C. C. Chow , Stationary Bumps in Networks of Spiking Neurons , NeuralComputation, 13 (2001), pp. 1473–1494.[55]
C. R. Laing and O. Omel’chenko , Moving bumps in theta neuron networks , Chaos: AnInterdisciplinary Journal of Nonlinear Science, 30 (2020), p. 043117.[56]
L. D. Landau , On the problem of turbulence , in Doklady Akademii Nauk USSR, vol. 44, 1944,p. 311.[57]
L. Lapicque and M. Lapicque , Recherches quantitatives sur l’excitation ´electrique des nerfstrait´ee comme une polarisastion , Journal de Physiologie et de Pathologie G´en´eral, 9 (1907),pp. 620–635.[58]
T. B. Luke, E. Barreto, and P. So , Complete classification of the macroscopic behavior ofa heterogeneous network of theta neurons , Neural Computation, 25 (2013), pp. 3207–3234.[59]
J. MacLaurin and P. Robinson , Determination of effective brain connectivity from activity
UMPS AND WAVES IN SPIKING NETWORKS correlations , Physical Review E, 99 (2019), p. 042404.[60] P. Manneville , On the transition to turbulence of wall-bounded flows in general, and planecouette flow in particular , European Journal of Mechanics-B/Fluids, 49 (2015), pp. 345–362.[61]
A. Meseguer and L. N. Trefethen , Linearized pipe flow to Reynolds number , Journalof Computational Physics, 186 (2003), pp. 178–197.[62] R. E. Mirollo and S. H. Strogatz , Synchronization of pulse-coupled biological oscillators ,SIAM Journal on Applied Mathematics, 50 (2006), pp. 1645–1662.[63]
E. Montbri´o, D. Paz´o, and A. Roxin , Macroscopic description for networks of spiking neu-rons , Physical Review X, 5 (2015), p. 021028.[64]
R. Os¸an, R. Curtu, J. Rubin, and B. Ermentrout , Multiple-spike waves in a one-dimensional integrate-and-fire neural network , Journal of Mathematical Biology, 48 (2004),pp. 243–274.[65]
R. Os¸an and B. Ermentrout , The evolution of synaptically generated waves in one- andtwo-dimensional domains , Physica D: Nonlinear Phenomena, 163 (2002), pp. 217–235.[66]
C. C. Pringle, Y. Duguet, and R. R. Kerswell , Highly symmetric travelling waves inpipe flow , Philosophical Transactions of the Royal Society A: Mathematical, Physical andEngineering Sciences, 367 (2009), pp. 457–472.[67]
J. Rankin, D. Avitabile, J. Baladron, G. Faye, and D. J. Lloyd , Continuation of lo-calized coherent structures in nonlocal neural field equations , SIAM Journal on ScientificComputing, 36 (2014), pp. B70–B93.[68]
A. D. Redish, A. N. Elga, and D. S. Touretzky , A coupled attractor model of the rodenthead direction system , Network: Computation in Neural Systems, 7 (1996), pp. 671–685.[69]
O. Reynolds , An experimental investigation of the circumstances which determine whetherthe motion of water shall be direct or sinuous, and of the law of resistance in parallelchannels , Philosophical Transactions of the Royal Society of London, (1883), pp. 935–982.[70]
K. A. Richardson, S. J. Schiff, and B. J. Gluckman , Control of Traveling Waves in theMammalian Cortex , Physical Review Letters, 94 (2005), p. 028103.[71]
D. Ruelle and F. Takens , On the nature of turbulence , Les rencontres physiciens-math´ematiciens de Strasbourg-RCP25, 12 (1971), pp. 1–44.[72]
L. Sacerdote and M. T. Giraudo , Stochastic integrate and fire models: a review on mathe-matical methods and their applications , in Stochastic Biomathematical Models, Springer,2013, pp. 99–148.[73]
H. Salwen, F. W. Cotton, and C. E. Grosch , Linear stability of poiseuille flow in a circularpipe , Journal of Fluid Mechanics, 98 (1980), pp. 273–284.[74]
H. Schmidt and D. Avitabile , Bumps and oscillons in networks of spiking neurons , Chaos:An Interdisciplinary Journal of Nonlinear Science, 30 (2020), p. 033133.[75]
A. Schmiegel and B. Eckhardt , Fractal stability border in plane couette flow , Physical Re-view Letters, 79 (1997), p. 5250.[76]
L. F. Shampine and M. W. Reichelt , The MATLAB ODE suite , SIAM Journal on ScientificComputing, 18 (1997), pp. 1–22.[77]
M. L. Steyn-Ross, D. A. Steyn-Ross, J. W. Sleigh, and D. R. Whiting , Theoretical pre-dictions for spatial covariance of the electroencephalographic signal during the anesthetic-induced phase transition: increased correlation length and emergence of spatial self-organization , Physical Review E, 68 (2003), p. 021902.[78]
H. C. Tuckwell , Introduction to Theoretical Neurobiology. Volume 1: Linear Cable Theoryand Dendritic Structure , vol. 8, Cambridge University Press, 1988.[79]
D. Turner-Evans, S. Wegener, H. Rouault, R. Franconville, T. Wolff, J. D. Seelig,S. Druckmann, and V. Jayaraman , Angular velocity integration in a fly heading circuit ,eLife, 6 (2017), p. e23496.[80]
C. W. van Doorne and J. Westerweel , The flow structure of a puff , Philosophical Transac-tions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 367 (2009),pp. 489–507.[81]
C. van Vreeswijk , Partial synchronization in populations of pulse-coupled oscillators , PhysicalReview E, 54 (1996), pp. 5522–5537.[82]
C. van Vreeswijk, L. F. Abbott, and G. B. Ermentrout , When inhibition not excitationsynchronizes neural firing , Journal of Computational Neuroscience, 1 (1994), pp. 313–321.[83]
H. Wedin and R. R. Kerswell , Exact coherent structures in pipe flow: Travelling wavesolutions , Journal of Fluid Mechanics, 508 (2004), pp. 333–371.[84]
H. R. Wilson and J. D. Cowan , A mathematical theory of the functional dynamics of corticaland thalamic nervous tissue , Kybernetik, 13 (1973), pp. 55–80.[85]
K. Wimmer, D. Q. Nykamp, C. Constantinidis, and A. Compte , Bump attractor dynam- D. AVITABILE, J. L. DAVIS, K. C. A. WEDGWOOD ics in prefrontal cortex explains behavioral precision in spatial working memory , NatureNeuroscience, 17 (2014), pp. 431–439.[86]
M. Wolfrum, O. E. Omel’chenko, and J. Sieber , Regular and irregular patterns of self-localized excitation in arrays of coupled phase oscillators , Chaos: An InterdisciplinaryJournal of Nonlinear Science, 25 (2015), pp. 053113–8.[87]