Macroscopic coherent structures in a stochastic neural network: from interface dynamics to coarse-grained bifurcation analysis

We study coarse pattern formation in a cellular automaton modelling a spatially-extended stochastic neural network. The model, originally proposed by Gong and Robinson [36], is known to support stationary and travelling bumps of localised activity. We pose the model on a ring and study the existence and stability of these patterns in various limits using a combination of analytical and numerical techniques. In a purely deterministic version of the model, posed on a continuum, we construct bumps and travelling waves analytically using standard interface methods from neural fields theory. In a stochastic version with Heaviside firing rate, we construct approximate analytical probability mass functions associated with bumps and travelling waves. In the full stochastic model posed on a discrete lattice, where a coarse analytic description is unavailable, we compute patterns and their linear stability using equation-free methods. The lifting procedure used in the coarse time-stepper is informed by the analysis in the deterministic and stochastic limits. In all settings, we identify the synaptic profile as a mesoscopic variable, and the width of the corresponding activity set as a macroscopic variable. Stationary and travelling bumps have similar meso- and macroscopic profiles, but different microscopic structure, hence we propose lifting operators which use microscopic motifs to disambiguate between them. We provide numerical evidence that waves are supported by a combination of high synaptic gain and long refractory times, while meandering bumps are elicited by short refractory times.


Introduction
In the past decades, single-neuron recordings have been complemented by multineuronal experimental techniques, which have provided quantitative evidence that the cells forming the nervous systems are coupled both structurally (Braitenberg and Schüz 1998) and functionally (for a recent review, see Yuste (2015) and references therein). An important question in neuroscience concerns the relationship between electrical activity at the level of individual neurons and the emerging spatio-temporal coherent structures observed experimentally using local field potential recordings (Einevoll et al. 2013), functional magnetic resonance imaging (Heuvel and Hulshoff Pol 2010) and electroencephalography (Nunez and Srinivasan 2006).
There exist a wide variety of models describing activity at the level of an individual neuron (Izhikevich 2007;Ermentrout and Terman 2010), and major research efforts in theoretical and computational neuroscience are directed towards coupling neurons in large-dimensional neural networks, whose behaviour is studied mainly via direct numerical simulations (Izhikevich and Edelman 2008;Fairhall and Sompolinsky 2014).
Recent studies have analysed neural fields with additive noise (Hutt et al. 2008;Faugeras et al. 2009;Kuehn and Riedler 2014), multiplicative noise (Bressloff and Webber 2012), or noisy firing thresholds (Brackley and Turner 2007), albeit these models are still mostly phenomenological. Even though several papers derive continuum neural fields from microscopic models of coupled neurons (Jirsa and Haken 1997;Bressloff 2009Bressloff , 2010Baladron et al. 2012), the development of a rigorous theory of multi-scale brain models is an active area of research.
Numerical studies of networks based on realistic neural biophysical models rely almost entirely on brute-force Monte Carlo simulations (for a very recent, remarkable example, we refer the reader to (Markram et al. 2015)). With this direct numerical simulation approach, the stochastic evolution of each neuron in the network is monitored, resulting in huge computational costs, both in terms of computing time and memory. From this point of view, multi-scale numerical techniques for neural networks present interesting open problems.
When few clusters of neurons with similar properties form in the network, a significant reduction in computational costs can be obtained by population density methods (Omurtag et al. 2000;Haskell et al. 2001), which evolve probability density functions of neural subpopulations, as opposed to single neuron trajectories. This coarse-graining technique is particularly effective when the underlying microscopic neuronal model has a low-dimensional state space (such as the leaky integrate-and-fire model) but its performance degrades for more realistic biophysical models. Developments of the population density method involve analytically derived moment closure approximations (Cai et al. 2004;Ly and Tranchina 2007). Both Monte Carlo simulations and population density methods give access only to stable asymptotic states, which may form only after long-transient simulations.
An alternative approach is offered by equation-free (Kevrekidis et al. 2003;Kevrekidis and Samaey 2009) and heterogeneous multiscale methods (Weinan and Engquist 2003;Weinan et al. 2007), which implement multiple-scale simulations using an on-the-fly numerical closure approximations. Equation-free methods, in particular, are of interest in computational neuroscience as they accelerate macroscopic simulations and allow the computation of unstable macroscopic states. In addition, with equation-free methods, it is possible to perform coarse-grained bifurcation analysis using standard numerical bifurcation techniques for time-steppers (Tuckerman and Barkley 2000).
The equation-free framework (Kevrekidis et al. 2003;Kevrekidis and Samaey 2009) assumes the existence of a closed coarse model in terms of a few macroscopic state variables. The model closure is enforced numerically, rather than analytically, using a coarse time-stepper: a computational procedure which takes advantage of knowledge of the microscopic dynamics to time-step an approximated macroscopic evolution equation. A single coarse time step from time t 0 to time t 1 is composed of three stages: (i) lifting, that is, the creation of microscopic initial conditions that are compatible with the macroscopic states at time t 0 ; (ii) evolution, the use of independent realisations of the microscopic model over a time interval [t 0 , t 1 ]; (iii) restriction, that is, the estimation of the macroscopic state at time t 1 using the realisations of the microscopic model.
While equation-free methods have been employed in various contexts (see Kevrekidis and Samaey (2009) and references therein) and in particular in neuroscience applications (Laing 2006;Laing et al. 2007Laing et al. , 2010; Siettos 2011, 2012;Laing and Kevrekidis 2015), there are still open questions, mainly related to how noise propagates through the coarse time stepper. A key aspect of every equationfree implementation is the lifting step. The underlying lifting operator, which maps a macroscopic state to a set of microscopic states, is generally non-unique, and lifting choices have a considerable impact on the convergence properties of the resulting numerical scheme . Even though the choice of coarse variables can be automatised using data-mining techniques, as shown in several papers by Laing, Kevrekidis and co-workers (Laing 2006;Laing et al. 2007Laing et al. , 2010, the lifting step is inherently problem dependent. The present paper explores the possibility of using techniques from neural field theory to inform the coarse-grained bifurcation analysis of discrete neural networks. A successful strategy in analysing neural fields is to replace the models' sigmoidal firing rate functions with Heaviside distributions (Bressloff 2012(Bressloff , 2014. Using this strategy, it is possible to relate macroscopic observables, such as bump widths or wave speeds, to biophysical parameters, such as firing rate thresholds. Under this hypothesis, a macroscopic variable suggests itself, as the state of the system can be constructed entirely via the loci of points in physical space where the neural activity attains the firing-rate threshold value. In addition, there exists a closed (albeit implicit) evolution equation for such interfaces (Coombes et al. 2012).
In this study, we show how the insight gained in the Heaviside limit may be used to perform coarse-grained bifurcation analysis of neural networks, even in cases where the network does not evolve according to an integro-differential equation. As an illustrative example, we consider a spatially-extended neural network in the form of a discrete time Markov chain with discrete ternary state space, posed on a lattice. The model is an existing cellular automaton proposed by Gong and Robinson (2012), and it has been related to neuroscience in the context of relevant spatio-temporal activity patterns that are observed in cortical tissue. In spite of its simplicity, the model possesses sufficient complexity to support rich dynamical behaviour akin to that produced by neural fields. In particular, it explicitly includes refractoriness and is one of the simplest models capable of generating propagating activity in the form of travelling waves. An important feature of this model is that the microscopic transition probabilities depend on the local properties of the tissue, as well as on the global synaptic profile across the network. The latter has a convolution structure typical of neural field models, which we exploit to use interface dynamics and define a suitable lifting strategy.
We initially study the model in simplifying limits in which an analytical (or semianalytical) treatment is possible. In these cases, we construct bump and wave solutions and compute their stability. This analysis follows the standard Amari framework, but is here applied directly to the cellular automaton. We then derive the corresponding lifting operators, which highlight a critical importance of the microscopic structure of solutions: one of the main results of our analysis is that, since macroscopic stationary and travelling bumps coexist and have nearly identical macroscopic profiles, a standard lifting is unable to distinguish between them, thereby preventing coarse numerical continuation. These solutions, however, possess different microstructures, which are captured by our analysis and subsequently by our lifting operators. This allows us to compute separate solution branches, in which we vary several model parameters, including those associated with the noise processes.
The manuscript is arranged as follows: In Sect. 2 we outline the model. In Sect. 3, we simulate the model and identify the macroscopic profiles in which we are interested, together with the coarse variables that describe them. In Sect. 4, we define a deterministic version of the full model and lay down the framework for analysing it. In Sects. 5 and 6, we respectively construct bump and wave solutions under the deterministic approximation and compute the stability of these solutions. In Sect. 7, we define and construct travelling waves relaxing the deterministic limit. In Sects. 8.1 and 8.2, we provide the lifting steps for use in the equation-free algorithm for the bump and wave respectively. In Sect. 9, we briefly outline the continuation algorithm and in Sect. 10, we show the results of applying this continuation to our system. Finally, in Sect. 11, we make some concluding remarks.

State variables for continuum and discrete tissues
In this section, we present a modification of a model originally proposed by Gong and Robinson (2012). We consider a one-dimensional neural tissue X ⊂ R. At each discrete time step t ∈ Z, a neuron at position x ∈ X may be in one of three states: a refractory state (henceforth denoted as −1), a quiescent state (0) or a spiking state (1). Our state variable is thus a function u : X × Z → U, where U = {−1, 0, 1}. We pose the model on a continuum tissue S = R/2LZ or on a discrete tissue featuring N + 1 evenly spaced neurons, We will often alternate between the discrete and the continuum setting, hence we will use a unified notation for these cases. We use the symbol X to refer to either S or S N , depending on the context. Also, we use u( · , t) to indicate the state variable in both the discrete and the continuum case: u( · , t) will denote a step function defined on S in the continuum case and a vector in U N with components u(x i , t) in the discrete case. Similarly, we write X u(x) dx to indicate S u(x) dx or 2L/N N j=0 u(x j ).

Model definition
We use the term stochastic model when the Markov chain model described below is posed on S N . An example of a state supported by the stochastic model is given in Fig. 1a.
In the model, neurons are coupled via a translation-invariant synaptic kernel, that is, we assume the connection strength between two neurons to be dependent on their mutual Euclidean distance. In particular, we prescribe that short range connections are excitatory, whilst long-range connections are inhibitory. To model this coupling, we use a standard Mexican hat function, and denote by W its periodic extension. . 1 a Example of a state u(x) ∈ U N and corresponding synaptic profile J (u)(x) ∈ R N in a stochastic network of 1024 neurons. b Schematic of the transition kernel for the network (see also Eqs. (5)- (7)). The conditional probability of the local variable u(x i , t + 1) depends on the global state of the network at time t, via the function q = f • J , as seen in (7) In order to describe the dynamics of the model, it is useful to partition the tissue X into the 3 pullback sets so that we can write, for instance, X u 1 (t) to denote the set of neurons that are firing at time t (and similarly for X u −1 and X u 0 ). Where it is unambiguous, we shall simply write X k or X k (t) in place of X u k (t). The synaptic input to a cell at position x i is given by a weighted sum of inputs from all firing cells. Using the synaptic kernel (1) and the partition (2), the synaptic input is then modelled as where κ ∈ R + is the synaptic gain, which is common for all neurons and 1 X is the indicator function of a set X .
Remark 1 (Dependence of J on u) Since X 1 depends on the state variable u, so does the synaptic input (3). Where necessary, we will write J (u)(x, t) to highlight the dependence on u. We refer the reader to Fig. 1 for a concrete example of synaptic profile.
The firing probability associated to a quiescent neuron is linked to the synaptic input via the firing rate function whose steepness and threshold are denoted by the positive real numbers β and h, respectively. We are now ready to describe the evolution of the stochastic model, which is a discrete-time Markov process with finite state space U N and transition probabilities specified as follows: for each x i ∈ S N and t ∈ Z where p ∈ (0, 1]. We give a schematic representation of the transitions of each neuron in the network in Fig. 1b. We remark that conditional probability of the local variable u(x i , t + 1) depends on the global state of the network at time t, via the function f • J . The model described by (1)- (7), complemented by initial conditions, defines a stochastic evolution map that we will formally denote as is a vector of control parameters.

Remark 2 (Microscopic, mesoscopic and macroscopic descriptions)
We will henceforth use the terms "microscopic", "mesoscopic" and "macroscopic" to refer to different state variables or model descriptions. Examples of these three state variables appear together in Figs. 2, 3 and 4 in Sect. 3, and we introduce them briefly here: Microscopic level Model (8) will be referred to as microscopic model and its solutions at a fixed time t as microscopic states. We will use these terms also when p = 1 and β → ∞, that is, when the evolution Eq. (8) is deterministic. Mesoscopic level In Remark 1, we associated to each microscopic state u a corresponding synaptic profile J , which is smooth, even when the tissue is discrete. We will not seek an evolution equation for the variable J , as the corresponding dynamical system would not reprent a reduction of the microscopic one. However, we will use J to bridge between the microscopic and macroscopic model descriptions; we therefore refer to J as a mesoscopic variable (or mesoscopic state). Macroscopic level Much of the present paper aims to show that, for the model under consideration, there exists a high-level model description, in the spirit of interfacial dynamics for neural fields (Bressloff 2012;Coombes et al. 2012;Bressloff 2014). The state variables for this level are points on the tissue where J (u)(x, t) attains the firing rate threshold h. We will denote these threshold crossings as {ξ i (t)} and we will discuss (reduced) evolution equations in terms of ξ i (t). The variables {ξ(t)} are therefore referred to as macroscopic variables and the corresponding evolution equations as the macroscopic model.  Table 1 3

Microscopic states observed via direct simulation
In this section, we introduce a few coherent states supported by the stochastic model. The main aim of the section is to show examples of bumps, multiple bumps and travelling waves, whose existence and stability properties will be studied in the following sections. In addition, we give a first characterisation of the macroscopic variables of interest and link them to the microscopic structure observed numerically.

Bumps
In a suitable region of parameter space, the microscopic model supports bump solutions (Qi and Gong 2015) in which the microscopic variable u(x, t) is nonzero only in a localised region of the tissue. In this active region, neurons attain all values in U. In Fig. 2, we show a time simulation of the microscopic model with N = 1024 neurons. At each time t, neurons are in the refractory (blue), quiescent (green) or spiking (yellow) state. We prescribe the initial condition by setting u(  Table 1 of a localised region, in which u(x i , 0) are sampled randomly from U. After a short transient, a stochastic microscopic bump is formed. As expected due to the stochastic nature of the system (Kilpatrick and Ermentrout 2013), the active region wanders while remaining localised. A space-time section of the active region reveals a characteristic random microstructure (see Fig. 2a). By plotting J (x, t), we see that the active region is well approximated by the portion of the tissue X ≥ = [ξ 1 , ξ 2 ] where J lies above the threshold h. A quantitative comparison between J (x, 50) and u(x, 50) is made in Fig. 2a. We interpret J as a mesoscopic variable associated with the bump, and ξ 1 and ξ 2 as corresponding macroscopic variables (see also Remark 2).

Multiple-bumps solutions
Solutions with multiple bumps are also observed by direct simulation, as shown in Fig. 3. The microstructure of these patterns resembles the one found in the single  Table 1 bump case (see Fig. 3a). At the mesoscopic level, the set for which J lies above the threshold h is now a union of disjoint intervals [ξ 1 , ξ 2 ], . . . , [ξ 7 , ξ 8 ]. The number of bumps of the pattern depends on the width of the tissue; the experiment of Fig. 3 is carried out on a domain twice as large as that of Fig. 2. The examples of bump and multiple-bump solutions reported in these figures are obtained for different values of the main control parameter κ (see Table 1), however, these states coexist in a suitable region of parameter space, as will be shown below.

Travelling waves
Further simulation shows that the model also supports coherent states in the form of stochastic travelling waves. In two spatial dimensions, the system is known to support travelling spots (Gong and Robinson 2012;Qi and Gong 2015). In Fig. 4, we show a time simulation of the stochastic model with initial condition Table 1 Parameter values for which the stochastic model supports a bump (Fig. 2), a multiple-bump solution (Fig. 3) and a travelling wave (Fig. 4). The value ∞ for the parameter β indicates that a Heaviside firing rate has been used in place of the sigmoidal function (4) u( In passing, we note that the state of the network at each discrete time t is defined entirely by the partition {X k } of the tissue; we shall often use this characterisation in the reminder of the paper. In the direct simulation of Fig. 4, the active region moves to the right and, after just 4 iterations, a travelling wave emerges. The microscopic variable, u(x, t), displays stochastic fluctuations which disappear at the level of the mesoscopic variable, J (x, t), giving rise to a seemingly deterministic travelling wave. A closer inspection ( Fig. 4a) reveals that the state can still be described in terms of the active region [ξ 1 , ξ 2 ] where J is above h. However, the travelling wave has a different microstructure with respect to the bump. Proceeding from right to left, we observe: 1. A region of the tissue ahead of the wave, x ∈ (ξ 2 , π), where the neurons are in the quiescent state 0 with high probability. 2. An active region x ∈ [ξ 1 , ξ 2 ], split in three subintervals, each of approximate width (ξ 2 − ξ 1 )/3, where u attains with high probability the values 0, 1 and −1 respectively. 3. A region at the back of the wave, x ∈ [−π, ξ 1 ), where neurons are either quiescent or refractory. We note that u = 0 with high probability as x → −π whereas, as x → ξ 1 , neurons are increasingly likely to be refractory, with u = −1.
A further observation of the space-time plot of u in Fig. 4b reveals a remarkably simple advection mechanism of the travelling wave, which can be fully understood in terms of the transition kernel of Fig. 1b upon noticing that, for sufficiently large β, q i = f (J (u))(x i ) ≈ 0 everywhere except in the active region, where q i ≈ 1. In Fig. 5, we show how the transition kernel simplifies inside and outside the active region and provide a schematic of the advection mechanism. For an idealised travelling wave profile at time t, we depict 3 subintervals partitioning the active region (shaded), together with 2 adjacent intervals outside the active region. Each interval is then mapped to another interval, following the simplified transition rules sketched in Fig. 5a: Schematic of the advection mechanism for the travelling wave state. Shaded areas pertain to the active region [ξ 1 (t), ξ 2 (t)], non-shaded areas to the inactive region X\[ξ 1 (t), ξ 2 (t)]. a In the active (inactive) region, q i = f (J (u))(x i ) ≈ 1 (q i ≈ 0), hence the transition kernel (5)-(7) can be simplified as shown. b At time t the travelling wave has a profile similar to the one in Fig. 4, which we represent in the proximity of the active region. We depict 5 intervals of equal width, 3 of which form a partition of [ξ 1 (t), ξ 2 (t)]. Each interval is mapped to another interval at time t + 1, following the transition rules sketched in (a). In one discrete step, the wave progresses with positive speed: so that J (x, t + 1) is a translation of J (x, t) 1. At the front of the wave, to the right of ξ 2 (t), neurons in the quiescent state 0 remain at 0 (rules for x / ∈ [ξ 1 , ξ 2 ]). 2. Inside the active region, to the left of ξ 2 (t), we follow the rules for in a clockwise manner: neurons in the quiescent state 0 spike, hence their state variable becomes 1; similarly, spiking neurons become refractory. Of the neurons in the refractory state, those being the ones nearest ξ 1 (t), a proportion p become quiescent, while the remaining ones remain refractory. The former remain at 0 whilst, of the latter, a proportion p transition into state 0, with the rest remaining at −1 (rules for x / ∈ [ξ 1 , ξ 2 ]). From this argument, we see that the proportion of refractory neurons in the back of the wave must decrease as ξ → −π .
The resulting mesoscopic variable J (x, t + 1) is a spatial translation by (ξ 2 (t) − ξ 1 (t))/3 of J (x, t). We remark that the approximate transition rules of Fig. 5a are valid also in the case of a bump, albeit the corresponding microstructure does not allow the advection mechanism described above.

Macroscopic variables
The computations of the previous sections suggest that, beyond the mesoscopic variable, J (x), coarser macroscopic variables are available to describe the observed patterns. In analogy with what is typically found in neural fields with Heaviside firing rate (Amari 1977;Bressloff 2014;Coombes and Owen 2004), the scalars {ξ i } defining the active region X ≥ = ∪ i [ξ 2i−1 , ξ 2i ], where J is above h, seem plausible macroscopic variables. This is evidenced not only by Figs. 2, 3 and 4, but also from the schematic in Fig. 5b, where the interval [ξ 1 (t), ξ 2 (t)] is mapped to a new interval [ξ 1 (t + 1), ξ 2 (t + 1)] of the same width. To explore this further, we extract the widths from the data in Figs. 2, 3 and 4, and plot the widths as a function of t. In all cases, we observe a brief transient, after which Δ i (t) relaxes towards a coarse equilibrium, though fluctuations seem larger for the bump and multiple bump when compared with those for the wave. In the multiple bump case, we also notice that all intervals have approximately the same asymptotic width (see Fig. 6b).

Deterministic model
We now introduce a deterministic version of the stochastic model considered in Sect. 2.2, which is suitable for carrying out analytical calculations. We make the following assumptions: 1. Continuum neural tissue. We consider the limit of infinitely many neurons and pose the model on S. 2. Deterministic transitions. We assume p = 1, which implies a deterministic transition from refractory states to quiescent ones (see Eq. (5)), and β → ∞, which induces a Heaviside firing rate f (I ) = (I − h) and hence a deterministic transition from quiescent states to spiking ones given sufficiently high input (see Eqs. 4,6).
In addition to the pullback sets X −1 , X 0 , and X 1 defined in (2), we will partition the tissue into active and inactive regions In the deterministic model, the transitions (5)-(7) are then replaced by the following rule We stress that the right-hand side of the equation above depends on u(x, t), since the partitions {X −1 , X 0 , X 1 } and {X < , X ≥ } do so (see Remark 1). As we shall see, it is sometimes useful to refer to the induced mapping of the pullback sets Henceforth, we will use the term deterministic model and formally write for (10), where the partition {X k } k∈U is defined by (2) and the active and inactive sets X ≥ , X < by (9).
Schematic of the analytical construction of a bump. A microscopic state whose partition comprises 3m + 2 strips is considered. The microscopic state, which is not an equilibrium of the deterministic system, has a characteristic width η 2 − η 1 , which differs from the width ξ m 2 − ξ m 1 of the mesoscopic bump J m . If we let m → ∞ while keeping η 2 − η 1 constant, then J m tends towards a mesoscopic bump J b and ξ m i → η i (see Proposition 1)

Macroscopic bump solution of the deterministic model
We now proceed to construct a bump solution of the deterministic model presented in Sect. 4. In order to do so, we consider a microscopic state with a regular structure, resulting in a partition, {X m k } k , with 3m + 2 strips (see Fig. 7) and then study the limit m → ∞.

Bump construction
Starting from two points η 1 , η 2 ∈ S, with η 1 < η 2 , we construct 3m intervals as follows We then consider states u m (x) = k∈U k1 X m k (x), with partitions given by and activity set X ≥ = [ξ m 1 , ξ m 2 ]. We note that, in addition to the 3m strips that form the active region of the bump, we also need two additional strips in the inactive region to form a partition of S. In general, {ξ m i } i = {η i } i , as illustrated in Fig. 7. Applying (10) or (11), we find Φ d (u m ) = u m , hence u m are not equilibria of the deterministic model. However, these states help us defining a macroscopic bump as a fixed point of a suitably defined map using the associated mesoscopic synaptic profile where we have highlighted the dependence of X m 1 on η 1 , η 2 . The proposition below shows that there is a well defined limit, J b , of the mesoscopic profile as m → ∞. We also have that ξ m i → η i as m → ∞ and that the threshold crossings of the activity set are roots of a simple nonlinear function.
Proposition 1 (Bump construction) Let W be the periodic extension of the synaptic kernel (1) and let h, κ ∈ R + . Further, let {A m i } 3m i=1 , X m 1 and J m be defined as in (13), (14) and (15), respectively, and let J b : S 3 → R be defined as The following results hold Proof We fix η 1 < η 2 and consider the 2L-periodic continuous mapping x → J b (x, η 1 , η 2 ), defined on S. We aim to prove that J m → J b uniformly in S. We pose Since W is continuous on the compact set S, it is also uniformly continuous on S.
Hence, there exists a modulus of continuity ω of W : We use ω to estimate |I m 1 (x) − I m 0 (x)| as follows: We have then |I m uniformly for all x ∈ S and η 1 , η 2 ∈ S with η 1 < η 2 , that is, result 1 holds true.
By hypothesis J b (0, 0, Δ) = h and, using a change of variables under the integral and the fact that W is even, it can be shown that J b (Δ, 0, Δ) = h, which proves result 2.

Corollary 1 (Bump symmetries) Let Δ be defined as in Proposition
Proof The assertion is obtained using a change of variables in the integral defining J b and noting that W is even.
The results above show that, ξ m i → η i as m → ∞, hence we lose the distinction between width of the microscopic pattern, η 2 − η 1 , and width of the mesoscopic Fig. 7, the factor 1/3 appearing in the expression for J b confirms that, in the limit of infinitely many strips, only a third of the intervals {A m j } j contribute to the integral. In addition, the formula for J b is useful for practical computations as it allows us to determine the width, Δ, of the bump. Fig. 7. In particular, Proposition 1 can be extended to a more general case of permuted intervals. More precisely, if we consider permutations, σ j , of the index sets {3 j −2, 3 j −1, 3 j} for j = 1, . . . , m and construct partitions

Remark 3 (Permuting intervals A m i ) A bump can also be found if the partition {X m k } of u m is less regular than the one depicted in
then the resulting J m converges uniformly to J b as m → ∞. The proof of this result follows closely the one of Proposition 1 and is omitted here for simplicity.

Bump stability
Once a bump has been constructed, its stability can be studied by employing standard techniques used to analyse neural field models (Bressloff 2012). We consider the map and study the implicit evolution The motivation for studying this evolution comes from Proposition 1, according to which the macroscopic bump ξ * = (0, Δ) is an equilibrium of (17), that is, b (ξ * , ξ * ) = 0. To determine coarse linear stability, we study how small perturbations of ξ * evolve according to the implicit rule (17). We set ξ(t) = ξ * + ε ξ(t), for 0 < 1 with ξ i = O(1) and expand (17) around (ξ * , ξ * ), retaining terms up to order ε, Using the classical ansatzξ (t) = λ t v, with λ ∈ C and v ∈ S 2 , we obtain the eigenvalue problem with eigenvalues and eigenvectors given by As expected, we find an eigenvalue with absolute value equal to 1, corresponding to a pure translational eigenvector. The remaining eigenvalue, corresponding to a compression/extension eigenvector, determines the stability of the macroscopic bump. The parameters A i , B i in Eq. (1) are such that W has a global maximum at x = 0, with W (0) > 0. Hence, the eigenvalues are finite real numbers and the pattern is stable if W (Δ) < 0. We will present concrete bump computations in Sect. 10.

Multi-bump solutions
The discussion in the previous section can be extended to the case of solutions featuring multiple bumps. For simplicity, we will discuss here solutions with 2 bumps, but the case of k bumps follows straightforwardly. The starting point is a microscopic structure similar to (14), with two disjoint intervals [η 1 , η 2 ), [η 3 , η 4 ) ⊂ S each subdivided into 3m subintervals. We form the vector η as m → ∞ uniformly in the variable x for all η 1 , . . . , η 4 ∈ S with η 1 < . . . < η 4 . In the expression above, J m and J b are the same functions used in Sect. 5.1 for the single bump.
In analogy with what was done for the single bump, we consider the mapping defined by Multi-bump solutions can then be studied as in Sect. 5. We present here the results for a multi-bump for L = π with threshold crossings given by where Δ satisfies A quick calculation leads to the eigenvalue problem x J Fig. 8 Stable mesoscopic multi-bump obtained for the deterministic model. We also plot the corresponding macroscopic bump ξ * (Eqs. 19-20) and coarse eigenvectors. Parameters are κ = 30, h = 0.9, p = 1, β → ∞, with other parameters as in Table 1 where . The real symmetric matrix in Eq. (21) has eigenvalues and eigenvectors given by As expected, we have one neutral translational mode. If the remaining 3 eigenvalues lie in the unit circle, the multi-bump solution is stable. A depiction of this multi-bump, with corresponding eigenmodes can be found in Fig. 8. We remark that the multi-bump presented here was constructed imposing particular symmetries (the pattern is even; bumps all have the same widths). The system may in principle support more generic bumps, but their construction and stability analysis can be carried out in a similar fashion.

Travelling waves in the deterministic model
Travelling waves in the deterministic model can also be studied via threshold crossings, and we perform this study in the present section. We seek a measurable function u tw : S → U and a constant c ∈ R such that almost everywhere in S and for all t ∈ Z. We recall that, in general, a state u(x, t) is completely defined by its partition, {X tw k (t)}. Consequently, Eq. (22) expresses that a travelling wave has a fixed profile u tw , whose partition, {X tw k }, does not depend on time. A travelling wave (u tw , c) satisfies almost everywhere the condition where Φ d is the deterministic evolution operator (12) and the shift operator is defined by σ x : u( · ) → u( · − x). The existence of a travelling wave is now an immediate consequence of the symmetries of W, as shown in the following proposition. An important difference with respect to the bump is that analytical expressions can be found for both microscopic and mesoscopic profiles, as opposed to Proposition 1, which concerns only the mesoscopic profile.
with partition Proof The assertion can be verified directly. We have hence the activity set for u tw is X tw Numerical simulations of the deterministic model confirm the existence of the mesoscopic travelling wave u tw in a suitable region of parameter space, as will be shown in Sect. 10. The main difference between u tw and the stochastic waves observed  9 Numerical investigation of the linear stability of the travelling wave of the deterministic system, subject to perturbations in the wake of the wave. We iterate the map Φ d starting from a perturbed state u tw + εũ, where u tw is the mesoscopic wave profile of Proposition 2, travelling with speed Δ, and εũ is non-zero only in two intervals of width 0.01 in the wake of the wave. We plot σ −tΔ Φ d (u tw + εũ) and the corresponding macroscopic profile as a function of t and we annotate the width of one of the perturbations. a For κ = 38, the wave is stable. b for sufficiently small κ, the wave becomes unstable in Fig. 4 is in the wake of the wave, where the former features quiescent neurons and the latter a mixture of quiescent and refractory neurons.

Travelling wave stability
As we will show in Sect. 10, waves can be found for sufficiently large values of the gain parameter κ. However, when this parameter is below a critical value, we observe that waves destabilise at their tail. This type of instability is presented in the numerical experiment of Fig. 9. Here, we iterate the dynamical system where u tw is the profile of Proposition 2, travelling with speed Δ, and the perturbation ε u tw is non-zero only in two intervals of width 0.01. We deem the travelling wave stable if u(z, t) → u tw (z) as t → ∞. For κ sufficiently large, the perturbations decay, as witnessed by their decreasing width in Fig. 9a. For κ = 33, the perturbations grow and the wave destabilises.
To analyse the behaviour of Fig. 9, we shall derive the evolution equation for a relevant class of perturbations to u tw . This class may be regarded as a generalisation of the perturbation applied in this figure and is sufficient to capture the instabilities observed in numerical simulations. We seek solutions to (23) with initial condition u(z, t) = k k1 X k (t) (z) with time-dependent partitions and activity set X ≥ (t) = [ξ 1 (t), ξ 2 (t)]. In passing, we note that for δ i = 0, the partition above coincides with {X tw k } in Proposition 2, hence this partition can be used as perturbation of u tw . Inserting the ansatz for u(ξ, t) into (23), we obtain a nonlinear implicit evolution equation, δ(t + 1), δ(t) = 0, for the vector δ(t) as follows (see Fig. 10) We note that the map above is valid under the assumption δ 3 (t) < δ 4 (t), which preserve the number of intervals of the original partition. As in Kilpatrick and Bressloff (2010), we note that this prevents us from looking at oscillatory evolution of δ(t). We set δ i (t) = ελ t v i , retain terms up to first order and obtain an eigenvalue problem for the matrix δ 5 (t + 1) = δ 6 (t) δ 6 (t + 1) = δ 7 (t)

Fig. 10
Visualisation of one iteration of the system (23): a perturbed travelling wave (top) is first transformed by Φ d using the rules (11) (centre) and then shifted back by an amount Δ (bottom). This gives rise to an implicit evolution equation δ(t + 1), δ(t) = 0 for the threshold crossing points of the perturbed wave, as detailed in the text . Once again, we have an eigenvalue on the unit circle, corresponding to a neutrally stable translation mode. If all other eigenvalues are within the unit circle, then the wave is linearly stable. Concrete calculations will be presented in Sect. 10.

Approximate probability mass functions for the Markov chain model
We have thus far analysed coherent states of a deterministic limit of the Markov chain model, and we now move to the more challenging stochastic setting. More precisely, we return to the original model (8) and find approximate mass functions for the coherent structures presented in Sect. 3 (see Figs. 2, 3, 4). These approximations will be used in the lifting procedure of the equation-free framework.
The stochastic model is a Markov chain whose 3 N -by-3 N transition kernel has entries specified by (1). It is useful to examine the evolution of the probability mass function for the state of a neuron at position x i in the network, μ k (x i , t) = Pr u(x i , t) = k , k ∈ U, which evolves according to We recall that f is the sigmoidal firing rate and that J is a deterministic function of the random vector, u(x, t) ∈ U N , via the pullback set X u 1 (t): As a consequence, the evolution equation for μ(x i , t) is non-local, in that J (x i , t) depends on the microscopic state of the whole network. We now introduce an approximate evolution equation, obtained by posing the problem on a continuum tissue S and by substituting J (x, t) by its expected value where μ : and In passing, we note that the evolution Eq. (25) is deterministic. We are interested in two types of solutions to (25): for all x ∈ S and t ∈ Z. 2. A travelling wave solution, that is, a mapping μ tw and a real number c such that μ(x, t) = μ tw (x − ct) for all x ∈ S and t ∈ Z.

Approximate probability mass function for bumps
We observe that, posing μ(y, t) = μ b (y) in (25), we have Motivated by the simulations in Sect. 3 and by Proposition 1, we seek a solution to (25) in the limit β → ∞, with We conclude that, for each x ∈ [0, Δ] (respectively x ∈ S\[0, Δ]), μ b (x) is the right · 1 -unit eigenvector corresponding to the eigenvalue 1 of the stochastic matrix Q ≥ (respectively Q < ). We find . 11 Comparison between the probability mass function μ b , as computed by (28)-(29), and the observed distribution μ of the stochastic model. a We compute the vector (μ b ) k , k ∈ U in each strip using (30) and visualise the distribution using vertically juxtaposed color bars, with height proportional to the values (μ b ) k , as shown in the legend. b A long simulation of the stochastic model supporting a stochastic bump u(x, t) for t ∈ [0, T ], where T = 10 5 . At each time t > 10 (allowing for initial transients to decay), we compute ξ 1 (t), ξ 2 (t), Δ(t) and then produce histograms for the random profile u(x − ξ 1 (t) − Δ(t)/2, t). c In the deterministic limit, the value of Δ is determined by (29), hence we have a Dirac distribution. d the distribution of Δ obtained in the Markov chain model. Parameters are as in Table 1 and, by imposing the threshold condition E[J ](Δ) = h, we obtain a compatibility condition for Δ, We note that if p = 1 we have where J b is the profile for the mesoscopic bump found in Proposition 1, as expected. In Fig. 11a, we plot μ b (x) as predicted by (28)-(29), for p = 0.7, κ = 30, h = 0.9. At each x, we visualise (μ b ) k for each k ∈ U using vertically juxtaposed color bars, with height proportional to the values (μ b ) k , as shown in the legend. For a qualitative comparison with direct simulations, we refer the reader to the microscopic profile u(x, 50) shown in the right panel of Fig. 2a: the comparison suggests that each u(x i , 50) is distributed according to μ b (x i ).
We also compared quantitatively the approximate distribution μ b with the distribution, μ(x, t), obtained via Monte Carlo samples of the full system (24). The distributions are obtained from a long-time simulation of the stochastic model supporting a microscopic bump u(x, t) for t ∈ [0, T ], with T = 10 5 . At each discrete time t, we compute the mesoscopic profile, J (u) (x, t), the corresponding threshold crossings and width: ξ 1 (t), ξ 2 (t), Δ(t) and then produce histograms for the random profile u(x − ξ 1 (t) − Δ(t) /2, t). The instantaneous shift applied to the profile is necessary to pin the wandering bump.
We note a discrepancy between the analytically computed histograms, in which we observe a sharp transition between the region x ∈ [0, Δ] and x ∈ S\[0, Δ], and the numerically computed ones, in which this transition is smoother. This discrepancy arises because Δ(t) oscillates around an average value Δ predicted by (29); the approximate evolution Eq. (25) does not account for these oscillations. This is visible in the histograms of Fig. 11c, d, as well as in the direct numerical simulation Fig. 6a.

Approximate probability mass function for travelling waves
We now follow a similar strategy to approximate the probability mass function for travelling waves. We pose μ( Proposition 2 provides us with a deterministic travelling wave with speed c = Δ. The parameter Δ is also connected to the mesoscopic wave profile, which has threshold crossings ξ 1 = −2Δ and ξ 2 = Δ. Hence, we seek for a solution to (25) in the limit where c is unknown. For simplicity, we pose the problem on a large domain whose size is commensurate with c, that is S = cT /R, where T is an even integer much greater than 1.
We obtain To make further analytical progress, it is useful to partition the domain S = cT /R in strips of width c, and impose that the wave returns back to its original position after T iterations, σ cT μ tw (z) = μ tw (z), while satisfying the compatibility condition h = E[J ](c). This leads to the system (30) With reference to system (30) we note that: 1. R(z, c) is constant within each strip I j , hence the probability mass function, μ tw (z), is also constant in each strip, that is, μ tw (z) = i ρ i 1 I i (c) (z) for some unknown vector (ρ −T /2 , . . . , ρ T /2 ) ∈ S 3T . 2. Each R j is a product of T 3-by-3 stochastic matrices, each equal to Q < or Q ≥ .
Furthermore, the matrices {R j } are computable. For instance, for the strip I −1 we have Consequently, μ tw (z) can be determined by solving the following problem in the unknown (ρ −T /2 , . . . , ρ T /2 , c) ∈ S 3T × R: Before presenting a quantitative comparison between the numerically determined distribution, μ tw (z), and that obtained via direct time simulations, we make a few efficiency considerations. In the following sections, it will become apparent that sampling the distribution μ tw (z) for various values of control parameters, such as h or κ, is a recurrent task, at the core of the coarse bifurcation analysis: each linear and nonlinear function evaluation within the continuation algorithm requires sampling μ tw (z), and hence solving the large nonlinear problem (31).
With little effort, however, we can obtain an accurate approximation to μ tw , with considerable computational savings. The inspiration comes once again from the analytical wave of Proposition 2. We notice that only the last equation of system (31) is nonlinear; the last equation is also the only one which couples {ρ j } with c. When p = 1 the wave speed is known as β → ∞, N → ∞ and p = 1 corresponds to the deterministic limit, hence E[J ](z) = J tw (z), which implies c = Δ and (ρ −1 ) 1 = 1. The stochastic waves observed in direct simulations for p = 1, however, display c ≈ (ξ 2 − ξ 1 )/3 = Δ and μ ≈ 1 in the strip where J achieves a local maximum (see, for instance Fig. 4, for which p = 0.4).
The considerations above lead us to the following scheme to approximate μ tw : (i) set c = Δ and remove the last equation in (31); (ii) solve T decoupled 3-by-3 eigenvalue Fig. 12 Similarly to Fig. 11, we compare the approximated probability mass function μ tw , and the observed distribution μ of the stochastic model. a The probability mass function is approximated using the numerical scheme outlined in the main text for the solution of (31); the strip I −1 is indicated for reference. b A set of 9 × 10 5 realisations of the stochastic model for a travelling wave are run for t ∈ [0, T ], where T = 1000. For each realisation s, we calculate the final threshold crossings ξ s 1 (T ), ξ s 2 (T ), and then compute histograms of u s (x − ξ s 2 (T ), T ). We stress that the strips in (a) are induced by our numerical procedure, while the ones in (b) emerge from the data. The agreement is excellent and is preserved across a vast region of parameter space (not shown). Parameters are as in Table 1 problems to find ρ i . Furthermore, if p remains fixed in the coarse bifurcation analysis, ρ i can be pre-computed and step (ii) can be skipped.
In Fig. 12a, we report the approximate μ tw found with the numerical procedure described above. An inspection of the microscopic profile u(x, 45) in the right panel of Fig. 4a shows that this profile is compatible with μ tw . We also compared quantitatively the approximate distribution with the distribution, μ(x, t), obtained with Monte Carlo samples of the full system (24). The distributions are obtained from M samples {u s (x, t)} M s=1 of the stochastic model for a travelling wave for t ∈ [0, T ]. For each sample s, we compute the thresholds, ξ s 1 (T ), ξ s 2 (T ), of the corresponding J (u s )(x, T ) and then produce histograms for u s (x − ξ s 2 (T ), T ). This shifting, whose results are reported in Fig. 12b, does not enforce any constant value for the velocity, hence it allows us to test the numerical approximation μ tw . The agreement between the two distributions is excellent: we stress that, while the strips in Fig. 12a are enforced by our approximation, the ones in Fig. 12b emerge from the data. We note a slight discrepancy, in that μ tw (−3Δ) ≈ 0, while the other distribution shows a small nonzero probability attributed to the firing state at ξ = −3Δ. Despite this minor disagreement, the differences between the approximated and observed distributions remain small across all parameter regimes of note and the approximations even retain their accuracy as β is decreased (not shown).

Coarse time-stepper
As mentioned in the introduction, equation-free methods allow us to compute macroscopic states in cases in which a macroscopic evolution equation is not available in closed form (Kevrekidis et al. 2003;Kevrekidis and Samaey 2009). To understand the general idea behind the equation-free framework, we initially discuss an example taken from one of the previous sections, where an evolution equation does exist in closed form.
In Sect. 5, we described bumps in a deterministic limit of the Markov chain model. In this description, we singled out a microscopic state (the function u m (x) with partition (14)) and a corresponding mesoscopic state (the function J m (x)), both sketched in Fig. 7. Proposition 1 shows that there exists a well defined mesoscopic limit profile, J b , which is determined (up to translations in x) by its threshold crossings ξ 1 = 0, ξ 2 = Δ. This suggests a characterisation of the bump in terms of the macroscopic vector (ξ 1 , ξ 2 ) or, once translation invariance is factored out, in terms of the macroscopic bump width, Δ. Even though the microscopic state u m is not an equilibrium of the deterministic system, the macroscopic state (0, Δ) is a fixed point of the evolution Eq. (17), whose evolution operator is known in closed form, owing to Proposition 1. It is then possible to compute Δ as a root of an explicitly available nonlinear equation.
We now aim to use equation-free methods to compute macroscopic equilibria in cases where we do not have an explicit evolution equation, but only a numerical procedure to approximate . As mentioned in the introduction, the evolution equation is approximated using a coarse time-stepper, which maps the macroscopic state at time t 0 to the macroscopic state at time t 1 using three stages: lifting, evolution, restriction. The specification of these stages (the lifting in particular) typically requires some closure assumptions, which are enforced numerically. In our case, we use the analysis of the previous sections for this purpose. In the following section, we discuss the coarse time-stepper for bumps and travelling waves. The multi-bump case is a straightforward extension of the single bump case.

Coarse time-stepper for bumps
The macroscopic variables for the bump are the threshold crossings {ξ i } of the mesoscopic profile J . The lifting operator for the bump takes as arguments {ξ i } and returns a set of microscopic profiles compatible with these threshold crossings: If β → ∞, u s (x) are samples of the analytical probability mass function μ b (x + Δ/2), where μ b is given by (28) with Δ = ξ 2 − ξ 1 . In this limit, a solution branch may also be traced by plotting (29).
If β is finite, we either extract samples from the approximate probability mass function μ b used above, or we extract samples u s (x) satisfying the following properties (see Proposition 1 and Remark 3): A more precise description of the latter sampling is given in Algorithm 1. As mentioned in the introduction, lifting operators are not unique and we have given above two possible examples of lifting. In our computations, we favour the second sampling method. The mesoscopic profiles, J , generated using this approach are wellmatched to E[J ] produced by the analytically derived probability mass functions (28). Numerical experiments demonstrate that this method is better than the first possible lifting choice at continuing unstable branches. This is most likely due to the fact that the latter method slightly overestimates the probability of neurons within the bump to be in the spiking state, and underestimates that of them being in the refractory state and this helps mitigate the problems encountered when finding unstable states caused by the combination of the finite size of the network and non-smooth characteristics of the model (when β is high).
The evolution operator is given by where ϕ T denotes T compositions of the microscopic evolution operator (8) and the dependence on the control parameter, γ , is omitted for simplicity. For the restriction operator, we compute the average activity set of the profiles. More specifically, we set : Profiles u 1 (x), . . . , u M (x) Comments : The profiles u s (x) are assumed to be symmetric around x = (ξ 1 + ξ 2 )/2. The operator round rounds a real number to a computational grid with stepsize δx = 2L/N . Pseudocode: We also point out that the computation stops if the two sets above are empty, whereupon, we set ξ s 1 = ξ s 2 = 0. The coarse time-stepper for bumps is then given by where the dependence on parameter γ has been omitted.

Coarse time-stepper for travelling waves
In Sect. 7.1, we showed that the probability mass function, μ tw (z), of a coarse travelling wave can be approximated numerically using the travelling wave of the deterministic model, by solving a simple set of eigenvalue problems. It is therefore natural to use μ tw in the lifting procedure for the travelling wave. In analogy with what was done for the bump, our coarse variables (ξ 1 , ξ 2 ) are the boundaries of the activity set associated with the coarse wave, X ≥ = [ξ 1 , ξ 2 ]. We then set where {u s (x i )} s are M independent samples of the probability mass functions μ tw (x i ), with c = (ξ 2 − ξ 1 )/3. The restriction operator for travelling waves is the same as used for the bump. The coarse time-stepper for travelling waves, Φ tw , is then obtained as in (32), with L b replaced by L tw .

Root finding and pseudo-arclength continuation
Once the coarse time-steppers, Φ b and Φ tw , have been defined, it is possible to use Newton's method and pseudo-arclength continuation to compute coarse states, continue them in one of the control parameters and assess their coarse linear stability. In this section, we will indicate dependence upon a single parameter γ ∈ R, implying that this can be any of the control parameters in (8).
For bumps, we continue in γ the nonlinear problem F b (ξ ; γ ) = 0, where A vector ξ such that F b (ξ ; γ ) = 0 corresponds to a coarse bump with activity set X ≥ = [0, ξ 2 ] and width ξ 2 , occurring for the parameter value γ , that is, we eliminated the translation invariance associated with the problem by imposing ξ 1 = 0. In passing, we note that it is possible to hardwire the condition ξ 1 = 0 directly in F b and proceed to solve an equivalent 1-dimensional system. Here, we retain the 2-dimensional formulation with the explicit condition ξ 1 = 0, as this makes the exposition simpler.
During continuation, the explicitly unavailable Jacobians are approximated using the first-order forward finite-difference formulas The finite difference formula for D ξ Φ b also defines the Jacobian operator used to compute stability: for a given solution ξ * of (34), we study the associated eigenvalue For coarse travelling waves, we define A solution (ξ, c) to the problem F tw (ξ, c; γ ) = 0 corresponds to a coarse travelling wave with activity set X ≥ = [0, ξ 2 ] and speed ξ 2 /3, that is, we eliminated the translation invariance and imposed a speed c in accordance with the lifting procedure L tw . As for the bump we can, in principle, solve an equivalent 1-dimensional coarse problem.

Numerical results
We begin by testing the numerical properties of the coarse time-stepper, the Jacobianvector products and the Newton solver used for our computations. In Fig. 14a, we evaluate the Jacobian-vector product of the coarse time stepper with p = 1, β → ∞ for bumps (waves) evaluated at a coarse bump (wave), in the direction ε ξ , where 0 < ε 1 and ξ is a random vector with norm 1. Since this coarse time stepper corresponds to the deterministic case, we expect the norm of the Jacobian-vector product to be an O(ε), as confirmed by the numerical experiment. In Fig. 14b, we repeat the experiment in the stochastic setting ( p = 0.4), for the travelling wave case with various number of realisations. As expected, the norm of the Jacobianvector action follows the O(ε) curve for sufficiently large ε: the more realisations are employed, the more accurately the O(ε) curve is followed.
We then proceed to verify directly the convergence history of the damped Newton solver. In Fig. 15a, we use a damping factor 0.5 and show the residual of the problem as a function of the number of iterations, showing that the method converges quickly to a solution. At first sight, it is surprising that the achievable tolerance of the problem does not change when the number of realisations increases. A second experiment, however, reported in Fig. 15b, shows that this behaviour is caused by the low system size: when we increase N from 2 7 to 2 9 , the achievable tolerance decreases by one order of magnitude. Gong and Robinson (2012), and Qi and Gong (2015) found wandering spots and propagating ensembles using direct numerical simulations on the plane. Here, we perform a numerical bifurcation analysis with various control parameters for the structures found in Sect. 3 on a one-dimensional domain.

Numerical Bifurcation analysis
In Fig. 16a, we vary the primary control parameter κ, the gain of the convolution term, therefore, we study existence and stability of the bumps and the travelling pulse For sufficiently high κ, these states coexist and are stable in a large region of parameter space. We stress that spatially homogeneous mesoscopic states J (x) ≡ J * , with 0 = J * or J * > h are also supported by the model, but are not treated here. Interestingly, the three solution branches are disconnected, hence the bump analysed in this study does not derive from an instability of the trivial state. A narrow unstable bump Δ Bifurcation diagrams for bumps (B), multibumps (MB) and travelling waves (TW) using κ as bifurcation parameter parameter. a Using the analytical results, we see that bump, multi-bump and travelling wave solutions coexist and are stable for sufficiently high κ (see main text for details). b The solution branches found using the equation-free methods agree with the analytical results. Parameters as in Table 1 except h = p = 1.0, β → ∞ the branch displays an asymptote (not shown) as the bump widens indefinitely. On a finite domain, like the one reported in the figure, there is a maximum achievable width of the bump, due to boundary effects. The travelling wave is also initially unstable, but does not stabilise at the saddle node bifurcation. Instead, the wave becomes stable at κ ≈ 33, confirming the numerical simulations reported in Fig. 9.
In Fig. 16b, we repeat the continuation for the same parameter values, but on a finite network, using the coarse time-steppers outlined in Sects. 8.1 and 8.2. The numerical procedure returns results in line with the continuum case, even at the presence of the noise induced by the finite size. The branches terminate for large κ and low Δ: this can be explained by noting that, if J (x) ≡ 0, then the system attains the trivial resting state u(x) ≡ 0 immediately, as no neuron can fire; on a continuum network, Δ can be arbitrarily small, hence the branch can be followed for arbitrarily large κ; on a discrete network, there is a minimal value of Δ that can be represented with a finite grid.
We now consider continuation of solutions in the stochastic model. In Fig. 17, we vary the transition probability, p, from the refractory to quiescent state. In panel (a), we show analytical results, given by solving (28)-(29), whilst panel (b) shows results found using the equation-free method. We find qualitatively similar diagrams in both cases, though we note some quantitative differences, owing to the finite size of the network and the finiteness of β: at the presence of noise, the stationary solutions exist for a wider region of parameter space (compare the folds in Fig. 17a, b); a similar situation arises, is also valid for the travelling wave branches.
The analytical curves of Fig. 17a do not contain any stability information, which are instead available in the equation-free calculations of Fig. 17b, confirming that bump and multi-bump destabilise at a saddle-node bifurcation, whereas the travelling wave becomes unstable to perturbations in the wake, if p is too large. The lower branch of the travelling wave is present in the analytical results, but not in the numerical ones, The solution branches found using the equation-free method agree qualitatively with the analytical results, and we can use the method to infer stability. For full details, please refer to the text. Parameters as in Table 1 except κ = 20.0, h = 0.9, with β → ∞ for (a) and β = 20.0 for (b) Fig. 18 Bifurcation in the control parameter β. For a large range of values, we observe very little change in Δ as β is varied. Parameters as in Table 1  as this branch is not captured by our lifting strategy: when we lift a travelling wave for very low values of Δ, we have that J < h for all x ∈ S N and the network attains the trivial state u(x) ≡ 0 in 1 or 2 time steps, thereby the coarse time stepper becomes ineffective, as the integration time T can not be reduced to 0. Gong and Robinson (2012); Qi and Gong (2015) found that refractoriness is a key component to generating propagating activity in the network. The bifurcation diagram presented here confirm this, as we recognise 3 regimes: for high p (low refractory time) the system supports stationary bumps, as the wave is unstable; for intermediate p, travelling and stationary bumps coexist and are stable, while for low p (high refractory time) the system selects the travelling wave.
In Fig. 18, we perform the same computation now varying β, which governs the sensitivity of the transition from quiescence to spiking. Here, we see that the wave and both bump solutions are stable for a wide range of β values and furthermore, that Bifurcation diagram for bumps in a heterogeneous network. To generate this figure, we replaced the coupling function with W (x, y) = W (x − y)(1 + W 0 cos(y/s)), with W 0 = 0.01, s = 0.5. We observe the snaking phenomenon in the approximate interval κ ∈ [38, 52]. The branches moving upwards and to the right are stable, whereas those moving to the left are unstable. The images on the right, obtained via direct simulation, depict the solution profiles on the labelled part of the branches. We note the similarity of the mesoscopic profiles within the middle of the bump. The continuation was performed for the continumm, deterministic model with parameters are κ = 30, h = 0.9 these states are largely insensitive to variations in this parameter, implying that the Heaviside limit is a good approximation for the network in this region of parameter space.
Finally, we apply the framework presented in the previous sections to study heterogeneous networks. We modulate the synaptic kernel using a harmonic function, as studied in Avitabile and Schmidt (2015) for a neural field. As in Avitabile and Schmidt (2015), the heterogeneity promotes the formation of a hierarchy of stable coexisting localised bumps, with varying width, arranged in a classical snaking bifurcation diagram as presented in Fig. 19. A detailed study of this bifurcation structure, while possible, is outside the scope of the present paper.

Discussion
In this article, we have used a combination of analytical and numerical techniques to study pattern formation in a Markov chain neural network model. Whilst simple in nature, the model exhibits rich dynamical behaviour, which is often observed in more realistic neural networks. In particular, spatio-temporal patterns in the form of bumps have been linked to working memory (Funahashi et al. 1989;Colby et al. 1995;Goldman-Rakic 1995), whilst travelling waves are thought to be important for plasticity (Bennett 2015) and memory consolidation (Massimini et al. 2004;Rasch and Born 2013). Overall, our results reinforce the findings of Gong and Robinson 2012, namely that refractoriness is key to generating propagating activity: we have shown analytically and numerically that waves are supported by a combination of high gains in the synaptic input and moderate to long refractory times. For high gains and short refractory times, the network supports localised, meandering bumps of activity.
The analysis presented in this manuscript highlights the multiscale nature of the model by showing how evolution on a microscopic level gives rise to emergent behaviour at meso-and macroscopic levels. In particular, we established a link between descriptions of the model at multiple spatial scales: the identified coarse spatiotemporal patterns have typified and recognisable motifs at the microscopic level, which we exploit to compute macroscopic patterns and their stability.
To connect our micro-and macroscopic variables, we take advantage of interface approaches, which are typically applied to continuum networks. A notable exception is offered by Chow and Coombes (2006), who consider a network based upon the lighthouse model (Haken 2000a, b). In a similar vein to our approach, they show how analysis of the discrete network can be facilitated by considering a continuum approximation and derive threshold equations to define bump solutions. This analysis also highlights that perturbations to the microscopic state, specifically the phase arrangement within the bump, can alter the dynamics of the bump edges.
Chow and Coombes found that wandering bump solutions in the lighthouse model arise for sufficiently fast synaptic processing. This is congruent with our result that short refractory times in (8) elicit coherent bump states, since both refractory times and synaptic processing timescales affect the average firing rate of the neuron. However, bumps cease to exist in our model if the refractory times are too long, whereas the lighthouse model supports stationary bumps for slow synapses, which highlights the subtle differences between the roles of refractoriness and synaptic processing in neural networks. It should also be noted that the meandering observed, for instance, in Fig. 2 is due to noise, and that all bumps will tend to wander; on the other hand, the meandering described by Chow and Coombes arises from the deterministic dynamics of the lighthouse model, and it is triggered by a sufficently fast synaptic process. We also remark that, without modification, the lighthouse model does not support travelling wave solutions, and so we cannot make comparisons regarding these solutions.
Travelling waves and bumps have almost identical meso-and macroscopic profiles: if microscopic data were removed from Figs. 2a and 4a, the profiles and activity sets of these two patterns would be indistinguishable. We have shown that a disambiguation is however possible if the meso-and macroscopic descriptions take into account microscopic traits of the patterns: in the deterministic limit of the system, where mathematical analysis is possible, the microscopic structure is used in the partition sets of Propositions 1 and 2; in the stochastic setting with Heaviside firing rates and an infinite number of neurons, the microscopic structure is reflected in the approximate probability mass functions appearing in Sect. 7; in the full stochastic finite-size setting, where an analytical description is unavailable, the microscopic structure is hardwired in the lifting operators of the coarse time-steppers (Sect. 8).
An essential ingredient in our analysis is the dependence of the Markov chain transition probability matrix upon the global activity of the network, via the firing rate function f . Since this hypothesis is used to construct rate models as Markov processes (Bressloff 2010), our lifting strategy could be used in equation-free schemes for more general large-dimensional neural networks. An apparent limitation of the procedure presented here is its inability to lift strongly unstable patterns with low activity, as pointed out in Sect. 10. This limitation, however, seems to be specific to the model studied here: when Δ → 0, bumps destabilise with transients that are too short to be captured by the coarse time-stepper.
A possible remedy would be to represent the pattern via a low-dimensional, spatially-extended, spectral discretisation of the mesoscopic profile (see Laing (2006)), which would allow us to represent the synaptic activity below the threshold h. This would lead to a larger-dimensional coarse system, in which noise would pollute the Jacobian-vector evaluation and the convergence of the Newton method. Variancereduction techniques (Rousset and Samaey 2013) have been recently proposed for equation-free methods in the context of agent-based models , and we aim to adapt them to large neural networks in subsequent publications.