Finite-Size Effects on Traveling Wave Solutions to Neural Field Equations

Neural field equations are used to describe the spatio-temporal evolution of the activity in a network of synaptically coupled populations of neurons in the continuum limit. Their heuristic derivation involves two approximation steps. Under the assumption that each population in the network is large, the activity is described in terms of a population average. The discrete network is then approximated by a continuum. In this article we make the two approximation steps explicit. Extending a model by Bressloff and Newby, we describe the evolution of the activity in a discrete network of finite populations by a Markov chain. In order to determine finite-size effects—deviations from the mean-field limit due to the finite size of the populations in the network—we analyze the fluctuations of this Markov chain and set up an approximating system of diffusion processes. We show that a well-posed stochastic neural field equation with a noise term accounting for finite-size effects on traveling wave solutions is obtained as the strong continuum limit.


Introduction
The analysis of networks of neurons of growing size quickly becomes involved from a computational as well as from an analytic perspective when one tracks the spiking activity of every neuron in the network.It can therefore be useful to zoom out from the microscopic perspective and identify a population activity as an average over a certain group of neurons.
In the heuristic derivation of such population models it is usually assumed that each of the populations in the network is infinite, such that, in the spirit of the law of large numbers, the description of the activity in each population reduces to a description of the mean.By considering a spatially extended network and letting the density of populations go to infinity, neural field equations are obtained as the continuum limit of these models.Here we consider the voltage-based neural field equation, which is a nonlocal evolution equation of the form where u(x, t) describes the average membrane potential in the population at x at time t, w : R → [0, ∞) is a kernel describing the strengths of the synaptic connections between the populations, and the gain function F : R → [0, 1] relates the potential to the activity in the population.
Neural field equations were first introduced by Amari [1] and Wilson and Cowan [20,21] and have since been used extensively to study the spatio-temporal dynamics of the activity in coupled populations of neurons.While they are of a relatively simple form, they exhibit a variety of interesting spatio-temporal patterns.For an overview see for example [13,9,11,3,4].In this article we will concentrate on traveling wave solutions, modeling the propagation of activity, that were proven to exist in [12].
The communication of neurons is subject to noise.It is therefore crucial to study stochastic versions of (1).While several sources of noise have been identified on the single neuron level, it is not clear how noise translates to the level of populations.Since neural field equations are derived as mean field limits, the usual effects of noise should have averaged out on this level.However, the actual finite size of the populations leads to deviations from the mean field behavior, suggesting finite size effects as an intrinsic source of noise.
The (heuristic) derivation of neural field equations involves two approximation steps.First, the local dynamics in each population is reduced to a description of the mean activity.Second, the discrete network is approximated by a continuum.In this article we make these two approximation steps explicit.In order to describe deviations from the mean field behavior for finite population sizes, we set up a Markov chain to describe the evolution of the activity in the finite network, extending a model by Bressloff and Newby [5].The transition rates are chosen in such a way that we obtain the voltage-based neural network equation in the infinite population limit.We analyze the fluctuations of the Markov chain in order to determine a stochastic correction term describing finite-size effects.In the case of fluctuations around traveling wave solutions, we set up an approximating system of diffusion processes and prove that a well-posed stochastic neural field equation is obtained in the continuum limit.
In order to derive corrections to the neural field equation accounting for finite-size effects, in [2], Bressloff (following Buice and Cowan [6]) sets up a continuous time Markov chain describing the evolution of the activity in a finite network of populations of finite size N .The rates are chosen such that in the limit as N → ∞ one obtains the usual activity-based network equation.He then carries out a van Kampen system size expansion of the associated master equation in the small parameter 1/N to derive deterministic corrections of the neural field equation in the form of coupled differential equations for the moments.To first order, the finite-size effects can be characterized as Gaussian fluctuations around the mean field limit.
The model is considered from a mathematically rigorous perspective by Riedler and Buckwar in [18].They make use of limit theorems for Hilbert-space valued piecewise deterministic Markov processes recently obtained in [19] as an extension of Kurtz's convergence theorems for jump Markov processes to the infinite-dimensional setting.They derive a law of large numbers and a central limit theorem for the Markov chain, realizing the double limit (number of neurons per population to infinity and continuum limit) at the same time.They formally set up a stochastic neural field equation, but the question of well-posedness is left open.
In [5], Bressloff and Newby extend the original approach of [2] by including synaptic dynamics and consider a Markov chain modeling the activity coupled to a piecewise deterministic process describing the synaptic current (see also section 6.4 in [4] for a summary).In two different regimes, the model covers the case of Gaussian-like fluctuations around the mean-field limit as derived in [2], as well as a situation in which the activity has Poisson statistics as considered in [7].
Here we consider the question how finite-size effects can be included in the voltagebased neural field equation.We take up the approach of describing the dynamics of the activity in a finite-size network by a continuous-time Markov chain and motivate a choice of jump rates that will lead to the voltage-based network equation in the infinite-population limit.We derive a law of large numbers and a central limit theorem for the Markov chain.Instead of realizing the double limit as in [18], we split up the limiting procedure, which in particular allows us to insert further approximation steps.We follow the original approach by Kurtz to determine the limit of the fluctuations of the Markov chain.By linearizing the noise term around the traveling wave solution, we obtain an approximating system of diffusion processes.After introducing correlations between populations lying close together (cf.section 5.1) we obtain a well-posed L 2 (R)-valued stochastic evolution equation, with a noise term approximating finite-size effects on traveling waves, which we prove to be the strong continuum limit of the associated network.
The article is structured as follows.We recall how population models can be derived heuristically in section 2 and summarize the work on the description of finite-size effects that can be found in the literature so far.In section 3 we introduce our Markov chain model for determining finite-size effects in the voltage-based neural field equation and prove a law of large numbers and a central limit theorem for our choice of jump rates.We use it to set up a diffusion approximation with a noise term accounting for finite-size effects on traveling wave solutions in section 4. Finally, in section 5, we prove that a well-posed stochastic neural field equation is obtained in the continuum limit.

Assumptions on the Parameters
As usual, we take the gain function F : R → [0, 1] to be a sigmoid function, for example Our assumptions on the synaptic kernel w are the following Assumption (iv) on F implies that a 1 and a 2 are stable fixed points of (1), while a is an unstable fixed point.It has been shown in [12] that under these assumptions there exists a unique monotone traveling wave solution to (1) connecting the stable fixed points (and in [8], that traveling wave solutions are necessarily monotone).That is, there exists a unique wave profile û : R → [0, 1] and a unique wave speed and lim As also pointed out in [12], we can without loss of generality assume that c ≥ 0. Note that ûx ∈ L 2 (R) since in the case c > 0 and in the case c = 0, 2 Finite-Size Effects in Population Models

Population Models
In population models, or firing rate models, instead of tracking the spiking activity of every neuron in the network, neurons are grouped together and the activity is identified as a population average.We start by giving a heuristic derivation of population models, distinguishing as usual between an activity-based and a voltage-based regime.
We consider a population of N neurons.We say that a neuron is 'active' if it is in the process of firing an action potential such that its membrane potential is larger than some threshold value κ.If ∆ is the width of an action potential, then a neuron is active at time t if it fired a spike in the time interval (t − ∆, t].We define the population activity at a given time t as the proportion of active neurons, a N (t) = # neurons that are active at time t N ∈ {0, 1 N , 2 N , . . ., 1}.
We assume that all neurons in the population are identical and receive the same input.If the neurons fire independently from each other, then for a constant input current I, where F (I) is the probability that a neuron receiving constant stimulation I is active.In the infinite population limit, the population activity is thus related to the input current via the function F , called the gain function.Sometimes one also defines F as a function of the potential u, assuming that the potential is proportional to the current as in Ohm's law.F is typically a nonlinear function.It is usually modeled as a sigmoid, for example for some γ > 0 and some threshold κ > 0, imitating the threshold-like nature of spiking activity.
Sometimes a firing rate is considered instead of a probability.We define the population firing rate λ N as λ δ,N (t) = # spikes in the time interval (t − δ, t] δN . If δ = ∆, then ∆λ δ,N (t) = a N (t).At constant potential u, lim δ→0 lim N →∞ λ δ,N (t) = λ(u), where λ(u) is the single neuron firing rate.Note that λ ≤ 1 ∆ .The firing rate is related to the probability F (u) via If the stimulus varies in time, then the activity may track this stimulus with some delay such that a(t + τ a ) = F (u(t)) for some time constant τ a .Taylor expansion of the left-hand side gives an approximate description of the (infinite population) activity in terms of the differential equation to which we refer as the rate equation.
We now consider a network of P populations, each consisting of N neurons.We assume that each presynaptic spike in population j at time s causes a postsynaptic potential in population i at time t.Here the (w ij ) are weights characterizing the strength of the synaptic connections between populations i and j, and τ m is the membrane time constant, describing how fast the membrane potential relaxes back to its resting value.Under the assumption that all inputs add up linearly, the potential in population i at time t is given as In the infinite population limit we obtain where The behavior of the coupled system (u i , a i ) depends on the two time constants, τ m and τ a .We consider two different regimes in which the model can be reduced to just one of the two variables, u or a.
Case 1: In this regime we can assume that the activity reacts to changes in input immediately such that a j (t) = F (u j (t)).Then (3) can be closed in the variables u i and we obtain Differentiation yields the system of ordinary differential equations which we will call the voltage-based neural network equation.
Case 2: τ a ≫ τ m → 0 By (3), Letting τ m → 0 we obtain Using again that a i (t+τ a ) ≈ a i (t)+τ a a ′ i (t), we end up with the system of ordinary differential equations which we will call the activity-based neural network equation.

Finite-Size Effects in the Literature
In [5], Bressloff and Newby set up a model for the evolution of the activity in a network of finite populations.They define the activity in population j as where δ is a time window of variable size.If δ is chosen as the width of an action potential ∆, then we obtain our original notion of the activity, ∆a ∆,N = a N .Here the activity is modeled as a rate rather than a probability.Note that the number of spikes in the time interval (t − δ, t] is limited by n max := N ∨ δN ∆ .They describe the dynamics of a δ,N by a Markov chain with state space {0, 1 δN , ..., nmax δN } and jump rates where e i denotes the i-th unit vector, where λ(u) is the firing rate at potential u, related to the probability F (u) via ∆λ(u) = F (u), and where u N evolves according to (2.1), The idea is that the activation rate should be proportional to λ(u), while the inactivation rate should be proportional to the activity itself.The rates are chosen such that in the limit as N goes to infinity, we obtain the neural rate equation They consider two regimes.
In the first regime, the size of the time window δ is fixed, say δ = 1.If τ a ≫ τ m → 0, then as in section 2.1, u N i (t) = P j=1 w ij a δ,N j (t).The description of the Markov chain can thus be closed in the variables a δ,N i , leading to the model already considered in [2].In the limit N → ∞ one obtains the activity-based network equation By formally approximating to order 1 N in the associated master equation, they derive a stochastic correction to (9), leading to the diffusion approximation for independent Brownian motions B j .
In [18], Riedler and Buckwar rigorously derive a law of large numbers and a central limit theorem for the sequence of Markov chains as N tends to infinity.Note that the nature of the jump rates is such that the process has to be 'forced' to stay in its natural domain [0, nmax N ] by setting the jump rate to 0 at the boundary.As they point out, this discontinuous behavior is difficult to deal with mathematically.They therefore have to slightly modify the model and allow the activity to be larger than nmax N .They embed the Markov chain into L 2 (D) for a bounded domain D ⊂ R d and derive the LLN in L 2 (D) and the CLT in the Sobolev space H −α (D) for some α > d.
Case 2: δ = 1 N , τ m ≫ τ a In the second regime, the size of the time window δ goes to 0 as N goes to infinity such that δN = 1.In this case, They show that at fixed voltage u, the stationary distribution of the activity a δ,N evolving according to ( 8) is approximately Poisson with rate λ(u).This corresponds to the regime considered in [7].
In the limit N → ∞, a j (t) = λ(u j (t)) and the system reduces to the voltage-based network equation The third regime has not been considered explicitly in [5].It is the one which is relevant for us.
We go back to our original definition of the activity and fix the time window δ to be the length of an action potential ∆.We assume that the potential evolves slowly, τ m ≫ 0. Speeding up time, we define ũN For some large n, The potentials ũN i therefore only depend on the time-averaged activities given for We have If τ a ≪ τ m , the activity relaxes to its stationary distribution quickly on this time scale.At fixed voltage u, under the stationary distribution ν(u), If N < ∞, then the finite size of the populations causes deviations from (11).In order to determine these finite-size effects, in the next section we will set up a Markov chain X P,N to describe the evolution of the time-averaged activity ãN .

A Markov Chain Model for the Activity
We describe the evolution of the time-averaged activity by a Markov chain X P,N with state space E P,N = {0, 1 N , 2 N , . . ., 1} P .We define the jump rates as where for x ∈ R, x + := x ∨ 0, x − := −x ∨ 0, and where e i denotes the i-th unit vector.The idea behind this choice is the following: the time-averaged activity tends to jump up (down) if the potential in the population, which is approximately given by F −1 (x i ), is lower (higher) than the input from the other populations, which is given by P j=1 w ij x j .The probability that the activity jumps down (up) when the potential is lower (higher) than the input is assumed to be negligible.The jump rates are proportional to the difference between the two quantities, scaled by the factor F ′ (F −1 (x i )).They are therefore higher in the sensitive regime where F ′ ≫ 1, that is, where small changes in the potential have large effects on the activity.If a N i = F (u N i ) in all populations i, then the system is in balance.Note that the state space is naturally bounded since lim N , when N ≥ N 0 is large enough, and similarly at 0. We will see in Proposition 1 below that the Markov chain converges to the solution of (11) as the size of the populations N goes to infinity.
In [2] a different choice of jump rates was suggested in analogy to (8): Also this choice leads to (11) in the limit.In this picture, the jump rates are high in regions where the activity is high.Since, as explained above, one should think of the Markov chain as governing a slowly varying time-averaged activity, (12) seems like a more natural choice.
The generator of Q P,N of X P,N is given for bounded measurable f : E P,N → R by Let N 0 be such that the jump rates out of the interval 1 N0 , 1 − 1 N0 are 0.
Proposition 1.Let X P be the (deterministic) Feller process on [1/N 0 , 1 − 1/N 0 ] P with generator Proof.By a standard theorem on the convergence of Feller processes (cf.[15], Thm.19.25) it is enough to prove that for f ∈ C ∞ ([1/N 0 , 1 − 1/N 0 ] P ) there exist bounded measurable Then it is easy to see that uniformly in x.

Diffusion Approximation
We are now going to approximate X P,N by a diffusion process.To this end, we follow the standard approach due to Kurtz and derive a central limit theorem for the fluctuations of X P,N .This will give us a candidate for a stochastic correction term to (11).

A Central Limit Theorem
We write where π k : (0, 1) P → (0, 1), x → x k , is the projection onto the k-th coordinate, and is a martingale describing the fluctuations of the process.We start by determining the limit of these fluctuations.
w kj X P j (s) where B is a P -dimensional standard Brownian motion, and X P is the Feller process from Proposition 1 .
Proof.The bracket process of M P,N k is given in terms of the carré du champ operator as w kj X P,N j (s) w kj X P,N j (s) ds. Thus, since for y with q P,N (X P,N (s), y) > 0 at least one of y k − X P,N k (s) and y l − X P,N l (s) is always 0. Now and the statement follows by the martingale central limit theorem, see for example Theorem 1.4, Chapter 7 in [14] .
This suggests to approximate X P,N by the system of coupled diffusion processes Using Itô's formula, we formally obtain an approximation for u P,N k Since the square root function is not Lipschitz continuous near 0, we cannot apply standard existence theorems to obtain a solution to (13) with the full multiplicative noise term.Instead we will linearize around a deterministic solution to the neural field equation and approximate to a certain order of 1 √ N .

Fluctuations around the Traveling Wave
Let ū be a solution to the neural field equation (1).To determine the finite-size effects on ū, we consider a spatially extended network, that is, we look at populations distributed over an interval [−L, L] ⊂ R and use the stochastic integral derived in Proposition 2 to describe the local fluctuations on this interval.Let m ∈ N be the density of populations on [−L, L] and consider P = 2mL populations located at k m , k ∈ {−mL, −mL+1, ..., mL−1}.We choose the weights w kl as a discretization of the integral kernel w : R → [0, ∞), Since we think of the network as describing only a section of the actual domain R, we add to each population an input F (ū t (−L)) and F (ū t (L)), respectively, at the boundaries with corresponding weights Fix a population size We write and assume that v k is of order 1/ √ N .Linearizing (13) around (ū k ) we obtain the approximation to order 1/ √ N .Note that bm k (t, ū) ≈ ∂ t ūk (t) = 0 for a stationary solution ū, with equality if ū is constant.The finite-size effects are hence of smaller order.Since the square root function is not differentiable at 0 we cannot expand further.
However, the situation is different if we linearize around a moving pattern.We consider the traveling wave solution u T W t (x) = û(x − ct) to (1) and we assume without loss of generality that c > 0. Then bm k (t, ū) ≈ ∂ t ūk (t) = −c∂ x u T W t < 0. This monotonicity property allows us to approximate to order 1/N in (13).Indeed, note that since û and F are increasing, So for L large enough, − bm k (t, u T W t ) > 0 and we have, using Taylor's formula and ( 16), As a possible diffusion approximation in the case of traveling wave solutions we therefore obtain the system of stochastic differential equations for which there exists a unique solution as we will see in the next section.

The Continuum Limit
In this section we take the continuum limit of the network of diffusions (18), that is, we let the size of the domain and the density of populations go to infinity in order to obtain a stochastic neural field equation with a noise term describing the fluctuations around the deterministic traveling wave solution due to finite size effects.
We thus have to deal with functions that 'look almost like the wave' and choose to work in the space S := u : R → R : u − û ∈ L 2 .Note that since for u 1 , u 2 ∈ S, u 1 − u 2 < ∞, the L 2 -norm induces a topology on S.

A Word on Correlations
Recall the definition of the Markov chain introduced in section 3. Note that as long as we allow only single jumps in the evolution, meaning that there will not be any jumps in the activity in two populations at the same time, the martingales associated with any two populations will be uncorrelated, yielding independent driving Brownian motions in the diffusion limit (cf.Proposition 2).
This only makes sense for populations that are clearly distinguishable.In order to determine the fluctuations around traveling wave solutions, we consider spatially extended networks of populations.The population located at x ∈ R is to be understood as the ensemble of all neurons in the ǫ-neighborhood (x − ǫ, x + ǫ) of x for some ǫ > 0. If we consider two populations located at x, y ∈ R with |x − y| < 2ǫ, then they will overlap.Consequently, simultaneous jumps will occur, leading to correlations between the driving Brownian motions.
Thus the Markov chain model (and the associated diffusion approximation) is only appropriate as long as the distance between the individual populations is large enough.When taking the continuum limit, we therefore adapt the model by introducing correlations between the driving Brownian motions of populations lying close together.

The Stochastic Neural Field Equation
We start by defining the limiting object.For u ∈ S and t ∈ [0, T ] set , y)h(y)dy for some symmetric kernel q(x, y) with q(x, •) ∈ L 2 ∩ L 1 for all x ∈ R and sup x∈R ( q(x, •) + q(x, •) 1 ) < ∞. (Details on the theory of Q-Wiener processes can be found in [17,10].)We assume that the dispersion coefficient is given as the multiplication operator associated with σ : [0, T ] × S → L 2 (R), which we also denote by σ, where σ is Lipschitz continuous with respect to the second variable uniformly in t ≤ T , that is, we assume that there exists L σ > 0 such that for all u 1 , u 2 ∈ S and t ∈ [0, T ], The correlations are described by the kernel q.For f, g in L 2 (R), , y)", where we denote by q * q(x, y) the integral q(x, z)q(z, y)dz.We could for example take for some small ǫ > 0 (cf.section 5.1).
We have σ(t, v) ∈ L 0 2 since by Parseval's identity Note that for uncorrelated noise (i.e.Q=E), this is not the case.Therefore, in [18] Riedler and Buckwar derive the central limit theorem in the Sobolev space H −α .Splitting up the limiting procedures, N → ∞ and continuum limit, allows us to incorporate correlations and finally to work in the more natural function space L 2 .
Proposition 3.For any initial condition u 0 ∈ S, the stochastic evolution equation has a unique strong S-valued solution.u has a continuous modification.For any p ≥ 2, For a proof see for example Prop.6.5.1 in [16].

Embedding of the Diffusion Processes
As a next step we embed the systems of coupled diffusion processes (18) into L 2 (R).Let m ∈ N be the population density and L m ∈ N the length of the domain with ), and let Note that the Brownian motions are independent as long as m < 1 4ǫ .For m ∈ N let σm : [0, T ] × R P → R P and assume that there exists L σm > 0 such that for any t ∈ [0, T ] and Consider the system of coupled stochastic differential equations with weights as in ( 14) and (15).We identify u = (u k ) −mL m ≤k≤mL m −1 ∈ R P with its piecewise constant interpolation as an element of L 2 via the embedding Then where and where We assume joint continuity and Lipschitz continuity in the second variable uniformly in m and t ≤ T , that is, there exists Proposition 4. For any initial condition u 0 ∈ L 2 (R) there exists a unique strong L 2 -valued solution u m to (22).u m admits a continuous modification.For any p ≥ 2, Proof.Again we check that the drift and diffusion coefficients are Lipschitz continuous.Note that Therefore, for and for an orthonormal basis (e k ) of L 2 (R) we obtain, using Parseval's identity,

Convergence
We are now able to state the main convergence result.We will need the following assumption on the kernel w.
Assumption 5.There exists C w > 0 such that for x ≥ 0, That assumption is satisfied for classical choices of w such as w Theorem 6. Fix T > 0. Let u and u m be the solutions to (21) and ( 22), respectively.Assume that m→∞ − −−− → 0, and for all p ≥ 2, We postpone the proof to section 6.
Remark 7. Let ǫ > 0. The kernel q(x, y) = 1 2ǫ 1 (x−ǫ,x+ǫ) (y) satisfies assumption (i) of the theorem.Indeed, note that for x, z with |x − z| ≤ 1 m , |{y : m .Therefore we obtain that for all k and for any The theorem applies to the case of the fluctuations described in section 4.2.In order to ensure that the diffusion coefficients are in L 2 (R), we cut off the noise outside a compact set {∂ x u T W t ≥ δ}, δ > 0. Note that the neglected region moves with the wave such that we always retain the fluctuations in the relevant regime away from the fixed points.
Theorem 8. Assume that the wave speed is strictly positive, c > 0. Fix δ > 0. The diffusion coefficients as derived in Section 4.2, where are jointly continuous and Lipschitz continuous in the second variable with Lipschitz constant uniform in m and t ≤ T , and satisfy condition (ii) of Theorem 6.
6 Proof of Theorem 6 For the proof of the theorem it therefore suffices to show that Using the Cauchy-Schwarz inequality we get With (25) it follows that Another application of the Cauchy-Schwarz inequality yields Using integration by parts, (23), and assumption (24), we obtain Analogously, Last we observe that Finally we consider We have and, as in (25), The last summand satisfies, using (23) and (25),
Let (e k ) be an orthonormal basis of L 2 (R).Note that by Parseval's identity (

σ 2 . 2 .
last term using the Burkholder-Davis-Gundy inequality, (26), and Young's inequality:m (s, v m s + π m (u T W s )) • Φ m − σ(s, v s + u T W s s, v m s + π m (u T W s )) • Φ m − σ(s, v s + u T W s s, v m s + π m (u T W s )) • Φ m − σ(s, v s + u T W sBringing the first summand to the left-hand side of (27) this implies that t,v m t + π m (u T W t )) • Φ m − σ(t, v t + u T W tWe estimate the last term as before and obtainE T 0 σ m (s, v m s + π m (u T W s )) • Φ m − σ(s, v s + u T W s 26)Applying Itô's formula to the real-valued stochastic process v m t − v t 2 we obtain for p ≥ 2≤ K 1 p v m t − v t p dt + K 2 p v m t − v t p−2 R(t, v t , m)dt + p v m t − v t p−2 dM tEstimating the last term as above and using Young's inequality we obtain σ m (t, v m t + π m (u T W t )) • Φ m − σ(t, v t + u T W 2 p R(t, v t , m) p 2 dt + p v m t − v t p−2 dM t .Integrating, maximizing over t ≤ T , and taking expectations we get E sup