Pulsating fronts in periodically modulated neural field models.

We consider a coarse-grained neural field model for synaptic activity in spatially extended cortical tissue that possesses an underlying periodicity in its microstructure. The model is written as an integrodifferential equation with periodic modulation of a translationally invariant spatial kernel. This modulation can have a strong effect on wave propagation through the tissue, including the creation of pulsating fronts with widely varying speeds and wave-propagation failure. Here we develop a new analysis for the study of such phenomena, using two complementary techniques. The first uses linearized information from the leading edge of a traveling periodic wave to obtain wave speed estimates for pulsating fronts, and the second develops an interface description for waves in the full nonlinear model. For weak modulation and a Heaviside firing rate function the interface dynamics can be analyzed exactly and gives predictions that are in excellent agreement with direct numerical simulations. Importantly, the interface dynamics description improves on the standard homogenization calculation, which is restricted to modulation that is both fast and weak.


I. INTRODUCTION
The propagation of waves of neural activity across the surface of the brain is known to subserve both natural and pathological neurobiological phenomena. An example of the former is spreading excitation associated with sensory processing [1], whilst waves in epilepsy are a classic example of the latter [2]. There is now a long history of using integro-differential neural field models to understand the properties of such waves; for a review of this area we refer the reader to [3]. For mathematical convenience such models are often assumed to be spatially translationally-invariant, and the issue of treating heterogeneity in brain dynamics is avoided. This approach allows one to bring to bear powerful techniques from the study of other translationally-invariant physical and biological systems, typically modeled with partial differential equations (PDEs). A review of the analysis of traveling fronts in such local models can be found in [4]. However, it is hard even at a first approximation to view the brain as a homogeneous system and so there is a pressing need to develop a set of mathematical tools for the study of waves in heterogeneous media that can be used in brain modeling. One common tool used in the PDE community for understanding systems with rapidly oscillating coefficients is homogenization. A particularly nice application of this technique can be found in the paper by Keener [5], who studied wave propagation in models of cardiac dynamics with rapidly varying spatial structure. This multi-scale approach leads to the study of smoothed or homogenized systems that are averaged over some well identified micro-scale. Bressloff [6] was the first to extend this technique to neural field models and to show how to describe fronts that travel through a neural model with a periodically modulated micro-structure. However, it is important to emphasize that as a perturbation technique homogenization theory requires that modulation on the micro-scale be both small in amplitude and rapidly varying in space. Interestingly, there are other techniques in the applied mathematical sciences that have been developed for the study of waves in random media that might also translate well to heterogeneous neural field models. In particular asymptotic techniques, stochastic homogenization, variational methods and large deviation methods, as reviewed by Xin [7,8], are all worthy of consideration. Of all these techniques the one developed by Shigesada et al. [9] for treating periodic waves in a heterogeneous (periodically modulated) single species population model of reaction-diffusion type seems the most appealing. This method uses linearized information from the leading edge of a traveling periodic wave to obtain wave speed estimates, generalizing the approach commonly used for the study of homogeneous PDEs [10]. One of the results we present here is the generalization of this technique to cover non-local interactions and its subsequent application to neural field models. Another result is the development of a new technique for studying the evolution of interfaces in heterogeneous neural fields; a technique that may translate back to the PDE community.
In section II we introduce the mathematical model of a one dimensional neural field with a periodically modulated non-local interaction kernel. Numerical simulation is used to illustrate the possibility of a traveling pulsating front solution in an excitatory model with a sinusoidally modulated exponentially decaying connectivity. Next, in section III we use linear analysis to develop a minimum wave speed calculation of traveling pulsating fronts, for a specific piecewise-linear firing rate function. This analysis is related to mathematical approaches for studying Kolmogorov-Petrovskii-Piscuinov type equations [11,12], but is able to cope with the non-locality of neural field models. Even though the calculation uses information from only the leading edge of a wave it gives wave speed predictions in very good agreement with those from direct numerical simulations. In section IV we consider a Heaviside firing rate function and focus on time-independent spatially varying solutions and show how to determine existence and stability of pinned front solutions. A fully nonlinear analysis of traveling fronts is developed in section V based on tracking the evolution of level set contours of the neural activity. For a Heaviside firing rate this interface dynamics simplifies considerably and the dynamics for a traveling front is determined by a single scalar nonlinear ordinary differential equation (ODE). For weak modulation this ODE can be analyzed using standard perturbation arguments and the wave speed of a pulsating front is easily calculated. Once again the theory is in excellent agreement with numerical simulations and, in contrast to results obtained using homogenization theory, can describe wave behavior even when the period modulation of the spatial kernel is not rapid. Finally in section VI we summarize and discuss possible extensions of this work.

II. THE MODEL
Many current coarse grained models of neural tissue trace their origins back to seminal work in the 1970s by Wilson and Cowan [13,14] and Amari [15,16]. In one spatial dimension the simplest single population model that describes the evolution of neural activity u = u(x, t), with x ∈ R and t ∈ R + is The function f represents the firing rate of the tissue and is often chosen to have a sigmoidal form. The weight kernel W (x, y) represents anatomical connectivity between points x and y in the tissue and the presence of this function means that the model has a non-local structure. For a more comprehensive discussion of such models, their generalization and their use in neuroscience we refer the reader to [17]. In this paper we mostly focus on a very specific form of heterogeneity, namely one where the connectivity has the following product structure: where J is some σ-periodic function. The use of periodically modulated translationally-invariant kernels of this form is particularly relevant to modeling primary visual cortex, which is known to have a periodic micro-structure on the millimeter scale, reviewed in some detail in the paper by Bressloff [6]. For later convenience we shall use a Fourier representation for the real periodic function J and write where J n ∈ C and † denotes complex conjugation. To provide a sense of the type of solutions that the model defined by (1), (2) and (3) can support we first present a direct numerical simulation for the choice and where H is the Heaviside function. Figure 1 shows a space-time plot of u(x, t) and a wave that connects a low activity state to a high activity state. Although there is an easily identified front that separates these two states it does not travel with constant speed and in fact oscillates periodically in a comoving frame (whose speed is chosen in an appropriate fashion, described in section III). In the wake of the front one also sees a periodic spatial pattern. We call such a wave a pulsating front. Numerical simulations with such a firing rate function show that with increasing such fronts can fail to propagate, instead becoming pinned, in qualitative (though not quantitative) agreement with predictions from homogenization theory. We will provide an analysis of pinned fronts in section IV, but before that we consider their propagating counterparts, for a specific, piecewise-linear, firing rate function.

III. TRAVELING PERIODIC WAVE
We focus on solutions of the type shown in Fig. 1 where enough time has passed so that a periodic traveling wave can develop defined by the conditions u(x, t) = u(x + σ, t + T ), lim x→∞ u(x, t) = 0 and lim x→−∞ u(x, t) = q(x) where q(x) is a stationary (spatially periodic) function. Note that the solution repeats itself in time T if it is observed at two points a distance σ apart. We define a speed c of the wave according to c = σ/T . Waves of this type have previously been studied by Shigesada et al. [9] in a PDE model of population dynamics in a heterogeneous environment. These authors developed wave speed estimates after linearizing about the leading edge of the wave. Here we adapt their arguments and develop a method for studying non-local systems that can be applied to neural field equations. To apply their methods we will work with the piecewise-linear firing rate function given by Note that for a more general sigmoidal form of firing rate, estimates of the wave speed would rely on information from across the whole solution and not just the leading edge. Later in sections IV and V we shall work with a Heaviside firing rate, which is obtained from equation (5) in the limit γ → ∞.

A. Dispersion relationship and minimum speed formula
For an exponentially decaying W (x) and firing rate functions for which f (0) = 0 the leading edge of a pulsating front solution is given simply by u = 0. This is often called a pulled front (as it is pulled by the unstable state u = 0). Linearizing around this solution gives the following linear integro-differential equation: where γ = lim u 0 f (u). It is natural to look for solutions to (6) Here we recognize ξ as a co-moving co-ordinate (in a frame of unknown speed c), so that u(x, t) has a product structure that separates into a traveling wave part and a spatial part. Substitution of this ansatz into (6) gives For the choice w(ξ) = αe λξ , λ ∈ (−∞, 0) (so that solutions decay to u = 0 as ξ → ∞), and writing v(x) = n v n e 2πinx/σ we find where Note that W † (−k, λ) = W (k, λ). Projecting (8) onto e −2πilx/σ and integrating over x gives an infinite system of linear equations for the amplitudes v l as Introducing the infinite matrix W(c, λ) with components we see that the system defined by (10) has a non-trivial solution if zero is an eigenvalue of W(c, λ). This matrix represents an unbounded linear operator W : In order to carry out practical calculations we truncate the system of equations (10) and consider only indices l from −N, . . . , N . This reduces the infinite eigenvalue problem of operator W to the calculation of a finite determinant. This truncation technique has a sound mathematical foundation and has often been used to study Hill's infinite matrix, especially as it arises in the analysis of the Mathieu equation (see for example [18]) or its time-delayed extension [19]. (In practice we used N = 20 and checked that our results were insensitive to increases in N .) Following [11], we now introduce a minimum propagation speed, c * , as where A(c, λ) is interpreted as the truncated (finite) version of (11). From this we can also define a critical decay rate λ * via det A(c * , λ * ) = 0. If we consider a traveling wave that arises from initial data with lim x→∞ u(x, 0) = αe λx with λ < 0 then around u = 0 we expect solutions of the form where c(λ) is defined by det A(c, λ) = 0. If the initial data decays sufficiently fast, namely λ < λ * then the first term in equation (13) dominates and we expect waves to emerge with the minimum wave speed c * [20,21]. If the initial data has a slow decay then the resulting wave will have a minimum speed determined by c(λ).
Consider for example the choice (4) so that J 1 = /(2i) and , |λ| < 1. (14) In this case the (finite) matrix A is tridiagonal with entries and the fixed point at u = 0 is unstable for small when −1 + γJ 0 > 0 (as can be seen by considering spatially coherent perturbations of the form u(x, t) = β exp(λt) in equation (6) and demanding that Re λ < 0).

B. Numerical Results
A plot of c as a function of λ, defined by the relation det A(c, λ) = 0, is shown in Fig. 2 for various values of . Here it can be seen that c = c(λ) has a well defined minimum, c * , that increases with increasing . In Fig. 3 we trace c * as a function of for different values of γ and make a comparison of the wave speed estimate with direct numerical simulations. For small values of these agree very well, and even for ∼ O(1) we still find good agreement. In Fig. 4 we show some corresponding plots for c * as a function of σ with various values of , to emphasize that good agreement between analysis and simulations holds over a large window of parameter space.
In the next two sections we move away from the use of linearized information to estimate wave speeds, and consider both stationary and moving patterns, showing how to analyze wave-speeds in the full nonlinear model when the firing rate is a Heaviside function.

IV. TIME INDEPENDENT SOLUTIONS
A stationary (time-independent) solution of (1) satisfies the equation  The stability of this solution can be determined by considering perturbations to q(x) of the form u(x)e λt . In this case linearizing around q(x) gives the equation For bounded solutions u(x) of (17) we shall say that q(x) is stable if Re λ < 0.
Since pinned pulsating front solutions provide a connection between a periodic spatial pattern (in the wake of the front) and a steady state (ahead of the front) it is natural to first consider the construction of a pure spatially periodic solution. Since solutions of (16) are rarely avail- able in closed form for an arbitrary firing rate (though see [22] for a recent discussion of how to analyze a particular class of sigmoids) we focus on the special case of a Heaviside for which solutions are easily constructed.

A. Periodic solutions
For a Heaviside firing rate σ-periodic solutions with q(x) > h for all x are given simply by For example, with the choice (4) then To ensure q(x) > h we must have that J 0 − | |/[1 + (2π/σ) 2 ] > h. For this problem the right hand side of (17) vanishes (since the solution never crosses h) and we have simply that λ = −1 and the solution is stable.

B. Pinned front
For a Heaviside firing rate function f (u) = f 2 (u) connections between periodic solutions, described in section IV A, and the homogeneous state q(x) = 0 can easily be constructed, and describe pinned fronts. An example is shown in Fig. 5. If q(x) > h for x < η where 0 ≤ η < σ and q(x) ≤ h otherwise we have simply that Imposing the condition q(η) = h fixes the parameter conditions for the existence of such a solution. with the choice (4) then where and tan φ = 2π/σ subject to the constraint 2h = J 0 + 1 + (2π/σ) 2 sin(2πη/σ − φ).
As a function of η equation (25) has either two solutions or none. In the latter case traveling fronts arise, invading the region where u = 0 if h is too low and retreating if h is too high. A plot of q(η) and q = h = const. is shown in Fig. 6 to illustrate the possible solutions of q(η) = h (where the line and curve meet). To ensure q(x) > h for x < 0 we must have that J 0 > | |/ 1 + (2π/σ) 2 . To determine stability we note that f 2 (q(y)) = δ(q(y)−h) = δ(y − η)/|q (η)| and (17) reduces to Demanding non-trivial solutions (u(η) = 0) we must have that where Hence, solutions with dq(η)/dη < 0 are stable and those with dq(η)/dη > 0 are unstable. For the example of Fig. 5 we see from Fig. 6 that of the two possible solutions it is the one with largest η that is stable.

V. INTERFACE DYNAMICS
Motivated by the form of the pulsating front in Fig. 1 we seek to describe the properties of this solution solely in terms of the behavior at the front edge which separates high activity from low. If the front is not pulsating (which is the case in the absence of period modulation of the connectivity) then in a traveling wave frame (of the same speed as the wave) the rising edge of the front may be identified with a single (traveling wave) co-ordinate. For a pulsating front this point is no longer stationary in time and instead oscillates. We now show how to derive the dynamics for this interface between high and low activity states.
In a co-moving frame the model (1) takes the form u = u(ξ, t) where ξ = x − c 0 t for some fixed c 0 and where We define a moving interface (level set) according to for some constant h. Here we are assuming that there is only one point on the interface (though in principle we could consider a set of points). Differentiation of (31) gives an exact expression for the velocity of the interface in the formξ Focusing now on the case of a Heaviside firing rate with f (u) = f 2 (u) means that for a pulsating front solution with u > h for ξ < ξ 0 (30) takes the simple form

A. Perturbation analysis
Consider a perturbation around the unmodulated case and write W (x, y) = W (|x − y|)[1 + K(y)] (K(y) = K(y + σ)) for some small parameter . For = 0 there is a traveling front q(ξ) given by the solution of where the speed c 0 is determined by q(0) = h. For small we assume that the slope of the traveling front varies sufficiently slowly so that we may make the convenient approximation u ξ | ξ=ξ0(t) = q ξ | ξ=0 . In this case we have, using equations (29) and (34), that Substitution of equations (35) and (36) into equation (32) givesξ For the choice K(x) = sin(2πx/σ) and W (x) = e −|x| /2 the time-dependent speed of the front is then given by with A = 1 2h − 1 Pulsating fronts are T -periodic solutions of the nonautonomous ordinary differential equation (37) with ξ 0 (t) = ξ 0 (t + T ). Introducing x 0 = ξ 0 + c 0 t with x 0 ∈ [0, σ] we may solve for the trajectory using: Using a half angle substitution we may evaluate this to give where z 0 (t) = tan[(2πx 0 (t)/σ − φ)/2] and x 0 (0) = 0. A periodic pulsating front with speed c = σ/T can be found by demanding that σ = x 0 (T ). Substitution of this condition into (41) shows that the speed of the pulsating front is given by Hence, a propagating wave is only supported if | | < 1/|A|.

B. Numerical Results
A plot of the wave speed as a function of the spatial scale σ is shown in Fig. 7. We see excellent agreement between the theoretical prediction from the interface dynamics approach and direct numerical simulations. Although our analysis in this case is restricted to small we have not had to make any assumptions about the scale of periodic modulation as determined by the parameter σ. In contrast a homogenization analysis would require both small and a periodic modulation that occurs on a smaller length-scale than the correlation length of W (x) (which is set to unity here). Following the calculation by Bressloff in [6] a naive application of homogenization theory gives As expected the result from the interface approach equation (42) recovers that from homogenization theory, given by (43), in the limit σ → 0. However, as illustrated in Fig. 7 the homogenization result gives poor results for σ away from zero. This has already been pointed out by Schmidt et al. [23], who use different arguments to derive equation (42). In their perturbation analysis they assume that the front is approximately homogeneous in a very small time window. Several other applications of the ideas in this section are possible, and we explore these briefly below.

C. Other forms of spatial modulation
In the main part of this paper we only considered spatial modulation of the coupling function W . Here we consider several other possible forms of spatial modulation. Suppose we replace (1) by where K is periodic. Performing the same analysis as in section V A, one obtains the equatioṅ and using the same functions W and K as in that section results in the expression where A 1 = 2h/(2h − 1). Note that in this expression c is independent of σ, the spatial wavelength of the modulation. Comparison of equation (47) with results from numerical simulations of (44)-(45) are shown in Fig. 8. We see good agreement for small , as expected.
Another possible way of including spatial modulation is to study the system in which the time-scale for the dynamics is modulated.
Performing a similar analysis to that above one obtains In this case, the critical value of beyond which waves are predicted not to travel is = 1, independent of both the spatial wavelength of modulation and the value of h, the threshold in the firing rate function. Figure 9 shows a comparison of (50) with results from numerical simulations of (48)-(49). The agreement is very good, although appears to become worse as σ → 0.
A third alternative is the system modeling a spatially-dependent stimulus, as originally proposed by Amari [16]. For functions K and W as above, we have where A 2 = 2/(2h − 1). Note that in this expression c is independent of σ, the spatial wavelength of the modulation. Figure 10 shows a comparison of (53) with results from numerical simulations of (51) and (52).

VI. DISCUSSION
In this paper we have shown that the biologically interesting case of a periodically modulated neural field model is amenable to mathematical analysis using two complementary techniques. One of these techniques uses linearization to obtain wave speed estimates for a particular piecewise-linear firing rate function, and is a generalization of approaches used in the PDE community. The other is a novel interface dynamics, especially relevant to neural field models where the firing rate switches rapidly between low and high values about some threshold. Although we have focused on describing front propagation both techniques carry over to the study of traveling pulses. The linearization approach cannot distinguish between the speed of a front and the speed of a pulse, whilst the interface approach provides more information as it tracks both the leading edge and trailing back of the pulse (via a pair of ODEs). For a recent study of propagating pulses in inhomogeneous neural media using averaging and homogenization theory we refer the reader to Kilpatrick et al. [24].
The two firing rate functions we considered are the piecewise-linear f 1 (u) (eqn. (5)) and the shifted Heaviside f 2 (u) (eqn. (18)). Both are approximations to the type of biophysically more plausible sigmoidal function that is often used [13,25]. To obtain the results in section III on traveling fronts we needed f (0) = 0 and γ ≡ lim u 0 f (u) > 1, thus f = f 1 was an appropriate choice. Since the value of γ is the only information about f (u) for u positive that is used in the analysis in section III, the results there will equally apply to other firing rate functions such as (for which lim u→∞ f (u) = 1) or point at which q(y) = h. Comparing Figs. 4 and 7 it is interesting to note that the speed of a front increases as σ is increased when f = f 1 and decreases as σ is increased when f = f 2 . This demonstrates that the qualitative behavior of a model such as (1)-(2) can depend on the exact form of the nonlinear function f . The only reasonable way to compare results obtained for f = f 1 and f = f 2 would be to take γ → ∞ in f 1 and h → 0 in f 2 . However, the result would be f (u) = H(u), the unshifted Heaviside, for which the speed of a front in the unmodulated case is infinite [c 0 = (1 − 2h)/(2h)]. Thus, a meaningful comparison is not really possible and our methods and results can be regarded as complementary to one another. In general, of course, one would like to derive results for an arbitrary sigmoidal function f .
The shifted Heaviside function f 2 is more appropriate for use in studying spreading patterns in two spatial dimensions, where interfaces will be more complicated than simple points. For example, in the case of a spreading spot the interface would be a closed curve. Importantly, the formalism of the level set approach that we have used in one spatial dimension carries over and the evolution of interfaces can be solely prescribed in terms of a dynamics on the interface (at least when the firing rate is a Heaviside). Indeed an interface description is a very natural way to treat neural field dynamics with other forms of heterogeneity, such as patchy connections [25][26][27][28].
It is worth emphasising that fronts or interfaces in periodic media are deterministic problems. In contrast to this there are many instances in neuroscience where lack of data means that it is more natural to model the inhomogeneous environment as a random process. It is then an open challenge to understand responses of the corresponding stochastic tissue model, though ideas developed for studying PDE models with random parameters may be an appropriate starting point [8]. One recent advance in the analysis of random traveling waves in the Nagumo equation with multiplicative noise has been made by Lord and Thümmler [29], and it would be interesting to explore the extension of their work to neural field models.