A nonlocal eigenvalue problem for the stability of a traveling wave in a neuronal medium, Discrete Contin

Past work on stability analysis of traveling waves in neuronal media has mostly focused on 
linearization around perturbations of spike times and has been done in the context of 
a restricted class of models. 
In theory, stability of such solutions could be affected by more general forms of 
perturbations. 
In the main result of this paper, 
linearization about more general perturbations is used to derive 
an eigenvalue problem for the stability of a traveling wave 
solution in the biophysically derived theta model, for which stability of waves has 
not previously been considered. 
The resulting eigenvalue problem is a nonlocal equation. 
This can be integrated to yield a reduced integral equation relating eigenvalues and wave speed, which is itself related to the Evans function for the nonlocal eigenvalue problem. 
I show that one solution to the nonlocal equation is the derivative of the wave, corresponding to translation invariance. 
Further, I establish that there is no unstable essential spectrum for this problem, 
that the magnitude of eigenvalues is bounded, and that for a special but commonly assumed form of coupling, any possible eigenfunctions for real, positive 
eigenvalues are nonmonotone on $(-\infty,0)$.


1.
Introduction.A variety of mathematical models have been developed to describe neuronal dynamics in various levels of detail.Models derived from biophysical principles provide the most complete and accurate description of neuronal behavior; however, their complexity may preclude mathematical analysis.The theta model, which is more mathematically tractable, arises in certain limits from a class of biophysically based neural models for cells near an activity threshold [5,4,7].In particular, this model applies near a saddle-node bifurcation on a limit cycle, through which the model cell undergoes a transition from being excitable, or intrinsically at rest, to oscillatory, or spontaneously active.
For a single cell, the theta model takes the form dθ dt = 1 − cos θ + (1 + cos θ)Z(t). (1.1) Here, θ(t) evolves on a circle, and Z(t) represents all inputs to the neuron.When θ(t) crosses through π, this corresponds to the firing of an action potential, or spike, by the model neuron.If Z is a negative constant, then equation (1.1) has stable and unstable fixed points.These coalesce in a saddle-node bifurcation on a limit cycle at Z = 0, and for Z > 0 the θ-neuron fires with period π/ √ Z.
Next, consider a network of synaptically coupled neurons, each near a saddlenode bifurcation on a limit cycle.Consider the transformation leading to the theta 926 JONATHAN E. RUBIN model, applied in the continuum limit.If the synaptic decay in the network occurs on the same time scale as the θ-dynamics, then the network of neurons will be reduced to a scalar model ([9], Remark 2) ∂θ(x, t) ∂t = 1 − cos θ(x, t) + (1 + cos θ(x, t)) β + g syn Ω J(x − y)s(y, t) dy , (1.2) ∂s(x, t) ∂t = aδ(θ(x, t) − π) − bs , (1.3) where β ∈ (−1, 0) is a bias parameter that controls the excitability of the cell, Ω is the domain of the network, s(x, t) denotes synaptic input, and a, b are positive constants.(See [7], [9] for a discussion of the case that arises when the synaptic decay rate is faster than the θ-dynamics.)Note that if the coupling term g syn J * s in equation (1.2) equals zero at some time t, then equation (1.2) has fixed points θ rest := − cos −1 ((1 + β)/(1 − β)), which is stable, and θ T := cos −1 ((1 + β)/(1 − β)), which is unstable.
In [13], a theorem is proved about the existence of traveling wave solutions to equation (1.2) with (x, t) ∈ IR × IR, such that Ω = IR.These solutions take the form θ(x, t) = θ(ξ) for ξ = ct − x with c > 0, and they satisfy θ(ξ) → θ rest as ξ → −∞ and θ(0) = π. (1.4) The second condition in (1.4) specifies that spikes occur at ξ = 0, corresponding to spike times t * (x) = x/c.This is part of the traveling wave ansatz and selects a particular wave from a translationally invariant family.For the result in [13], the s-equation (1.3) is dropped and the s(x, t) in equation (1.2) is simply taken to be a function α(t − x/c) = α(ξ/c) that satisfies certain assumptions stated below, which are consistent with equation (1.3).With this substitution, equation (1.2) can be considered without equation (1.3).The corresponding traveling wave equation takes the form (1.5)In addition to integrability, it is assumed that J(x), α(t) satisfy the following assumptions: The theorem states that if (H1), (H2) hold, then for g syn sufficiently large, there exist (at least) two traveling wave solutions θ slow,f ast (ξ) satisfying (1.4), (1.5), and that in these solutions, each cell fires more than one spike, corresponding to lim ξ→∞ θ(ξ) > 3π.For this paper, I will also assume that J, α are differentiable for positive arguments, such that θ slow,f ast (ξ) are both at least twice differentiable.The function α is not assumed to be differentiable at 0; indeed, many commonly used forms of α are not differentiable, or even not continuous, at 0. Note, however, that even when α is not differentiable, the traveling wave solution is C 1 since J is continuous.
In continuum networks in which neurons are characterized as simpler integrateand-fire units rather than as theta cells, fast and slow traveling waves also exist, and stability analysis has shown that the fast wave is stable and the slow wave is unstable under various conditions [2,3,6].However, due to the simplified nature of the integrate-and-fire model, the set of neuronal firing times uniquely specifies the form of a traveling wave solution.Hence, this stability analysis proceeds by simply treating linearization about spike times, without requiring independent consideration of perturbations to the waveforms.It is of great interest to analyze how the added biophysical realism in the theta model, and in particular the possibility for independent spike time and waveform perturbations, affects the stability of traveling waves.The stability of traveling waves was not addressed in [13], however.Indeed, this issue is confounded by the fact that the exact number of spikes fired by each cell in each traveling wave solution found there is not known.
In some neurons, there exist biological mechanisms that tend to suppress the firing of multiple spikes.These include synaptic depression, in which inputs corresponding to non-leading spikes in spike trains are weaker than inputs due to leading spikes, and adaptation, which involves changes in the membrane currents of cells that prevent the firing of repeated spikes even when input persists.The theta model with adaptation is presented in [12].An approximate version of this model takes the form (1.6) where ψ(t) is differentiable and satisfies (H3) The existence proof for traveling waves from [13] immediately generalizes to (1.6), assuming (H3) holds.Moreover, when k is sufficiently large and ψ does not decay too quickly, adaptation is strong enough to ensure that each cell fires only a single spike in each wave.The primary goal of this paper is to present the derivation of the eigenvalue problem that will determine the linearized stability of the single spike traveling wave solutions to the theta network model with strong adaptation.This is done in Section 2 of the paper.The resulting equation is nonlocal; technically, it is a nonautonomous ordinary differential equation, but the eigenvalue parameter enters into this equation in a nonlocal way.Stability of a traveling wave solution to a certain continuum neuronal model, itself an integrodifferential equation of a somewhat similar form to (1.2)-(1.3),has been analyzed previously, through use of an Evans function [15]; however, in that case, the eigenvalue problem lacks the nonlocal eigenvalue term.It is hoped that the derivation given here will inspire future work to analyze the nonlocal eigenvalue equation that results.
After presenting the derivation, I confirm in Section 3 that the derivative of the wave is an eigenfunction.In Section 4, I prove several potentially useful facts about the eigenvalue problem, as well as possible corresponding spectrum and eigenfunctions.Note that although I refer to the equation derived via linearization about a wave as an eigenvalue problem, it will also admit essential spectrum.I use standard arguments to show that there are no essential instabilities in this problem and that there is an upper bound on the magnitude of eigenvalues.The eigenvalue equation can be integrated, and I use this to derive a reduced integral equation relating eigenvalues and wave speed.In fact, solutions to this reduced equation correspond precisely to zeros of the Evans function for the eigenvalue problem derived in Section 2. This Evans function can be defined within the abstract framework presented in [10]; although the nonlocal term in the eigenvalue problem here is not considered in [10], the theory there can be easily extended to handle this term.Finally, I check that, for a special form of J and α, any eigenfunction corresponding to an eigenvalue λ ∈ IR + cannot be monotone on ξ ∈ (−∞, 0), where the wave itself is monotone.This illustrates the potential importance of considering perturbations of waveforms, not just of spike times, in stability analysis for such solutions.
2. Derivation of the eigenvalue equation.Let θ 0 (ct − x) ≡ θ 0 (ξ) denote a traveling wave solution (fast or slow) of the theta model (1.6); recall that c > 0. Note that θ 0 also satisfies (1.5) on ξ ∈ (−∞, ), where ψ = 0 by (H3).Consider the strong adaptation limit, such that the synaptic input α(t) yields only a single spike from each cell and θ 0 (ξ) → θ rest + 2π as ξ → ∞.Recall that this spike occurs at ξ = 0, or t = x/c.To check stability, one must consider perturbations of the waveform, θ(x, t) = θ 0 (ct − x) + v(x, t), and of spike times, say t * (x) = x/c + δt(x).However, for linearization of (1.6), it is necessary to expand α(t − t * (x)) about t − x/c, and the sign of δt(x) is not known a priori.This appears to be an obstacle, since some commonly assumed forms of α(t) fail to be differentiable at 0.
To circumvent this problem, one can substitute the perturbations into the δfunction formulation of the underlying θ-equation, namely To reduce notation, I will henceforth absorb the positive constant g syn into the function which holds when g(x) has a unique zero at x 0 [11].For this rest of this derivation, I will omit the kψ term from (2.7), as it is easily treated.The corresponding terms will be restored at the end of the section.
An expansion of the term θ t (y, t * (y)) yields, to first order, From equation (1.5), and the fact that θ 0 (0) = π, it follows that θ 0 (0) = f (π)/c = 2/c.Hence, the expression in (2.9) is positive for small perturbations, and the absolute value in (2.7) can be dropped.Using (2.9), equation (2.7) can be rewritten as (2.10) Next, take the leading order expansions of terms in (2.10) and cancel the appropriate terms with the sum ∂θ0 ∂t + c ∂θ0 ∂ξ on the left hand side.When the δ-function in (2.10) is expanded, it must be differentiated formally; this will be resolved in the course of the derivation.Also, note that α(t − s) = 0 for s > t.Thus, the expansion and cancellation yield (2.11) The right hand side of (2.11) features three integrals with respect to y.To simplify these, the following Propositions are useful: Proposition 2.1.To leading order, δt(y) = −v(y, y/c) cθ 0 (0) .
Proof: Recall that t * (y) denotes the spike time for the cell at position y, such that θ(y, t * (y)) = π.This yields to leading order, where the relation θ 0 (0) = π contributes to the final equality.

Remarks
1.The fact that θ 0 (0) equals zero could have been applied earlier in the derivation.I chose to keep the θ 0 (0) term in the derivation to add generality, as other integrodifferential equations related to (1.6) may feature solutions that do not have this property.2. Common choices for α(t) on [0, ∞) include te −bt and e −bt for some positive constant b.Even if α(0 + ) = 0, the eigenvalue equation (2.18) clearly maintains its nonlocal integral term.3.In the Appendix, a different form of eigenvalue problem from equation (2.17) is attained through an alternative derivation that omits the spike time perturbation (as well as the adaptation terms, for simplicity) but allows for perturbations to the time course of synaptic input.The resulting equation is identical to (2.17) but lacks the term proportional to v(0)α(0 + )J(ξ).This comparison makes clear that it can be essential to explicitly allow for perturbations of spike times, and that such perturbations can contribute a term depending on v(0)α(0 + )J(ξ).
3. Eigenvalue equation for λ = 0.Here I prove the following result, which makes up the last piece of Theorem 2.3 and which is expected to hold due to translation invariance, Proposition 3.1.The C 1 function dθ 0 /dξ is an eigenfunction for equation (2.18), with eigenvalue λ = 0.

Analysis of eigenvalue equation (2.18).
4.1.Essential spectrum and large |λ|.Continue to assume that each cell fires only a single spike in the wave.It is not difficult to restrict the possible instabilities arising in equation (2.18).In particular, the following results show that there are no essential instabilities and that a bound exists on the magnitude of any eigenvalues.Note that by definition, the spectrum, including the essential spectrum, is a property of a family of linear operators.The family corresponding to equation (2.18) is given by T (λ) : Below, we simply refer to the spectrum and the essential spectrum of equation (2.18).In particular, (4.21) holds for Reλ ≥ 0. Suppose that λ is chosen such that (4.21) is true.Then the nonlocal part of (2.18) represents a compact perturbation of a local equation, and the essential spectrum for (2.18) is determined by the corresponding asymptotic systems (see e.g.[14]), which take the form where n = 0 as ξ → −∞ and n = 1 as ξ → ∞.Therefore, if λ is selected such that (4.21) holds, then for bias β < 0, λ is in the essential spectrum for (2.18) if and only if for n = 0 or 1.In particular, if Reλ ≥ 0, then λ cannot be in the essential spectrum for (2.18), so no essential instabilities arise.
Remark: As β ↓ 0, the essential spectrum of (2.18) approaches the imaginary axis.In this limit, the uncoupled cell approaches the oscillatory state; that is, its rest and threshold states coalesce in a saddle-node bifurcation on a limit cycle at β = 0. Note, however, that β is not necessarily small, but rather can take any value in (-1,0), as the transformation to the theta model ensures that the cell is in the neighborhood of a saddle-node on a limit cycle for all β in this interval.
Proof: Let y = ξ|λ| and transform the eigenvalue equation (2.18) so that the independent variable is y, rather than ξ (cf.[1]).In the limit as |λ| → ∞, the resulting form of equation (2.18) becomes c v = −e i arg λ v, where differentiation is with respect to y.This limiting equation has no nontrivial bounded solutions that decay asymptotically.Hence, there exists a constant M > 0 such that all eigenvalues satisfy |λ| < M .
Remark: In [10], the Evans function is defined for a class of nonlocal eigenvalue problems dv dξ where the nonlocal term [N v](ξ) may take the form v(ξ 0 )h(ξ) for fixed ξ 0 ∈ IR and continuous function h(ξ).The definition given does not explicitly allow for the e −λη/c term, present in equation ( 2 4.3.Nonmonotonicity of eigenfunctions.In this subsection, I consider a certain typical choice of J(ξ) and α(t).In this case, I show that for λ ∈ IR + , any eigenfunctions must be nonmonotone on (−∞, 0).This illustrates the potential importance of considering perturbations to waveforms, not just to spike times, in the stability analysis of traveling wave solutions to the theta network model.Let ζ = e ξ and, abusing notation with A and B, consider the autonomous eigenvalue system where B(ζ, λ) = v(0) B(ζ, λ).Since ψ(ξ/c) = ψ (ξ/c) = 0 on (−∞, 0), I will omit the ψ terms from A, B below.A nontrivial solution must have v(0) = 0; without loss of generality, one can scale v(0) = θ 0 (0) > 0. For consistency, every eigenfunction must be a nonzero scalar multiple of a solution to (4.26) that makes a connection from (0, 0) to (v(0), 1) for this fixed v(0).The same calculation that ruled out essential instabilities for (2.18) shows that the unstable manifold of (0,0) under (4.26) is 1-dimensional, and one can seek such connections by tracking this manifold as ζ increases from 0. The v-nullcline for system (4.26) is given by Consideration of this nullcline yields a result on possible eigenfunctions for a concrete choice of J, α.
Remark: As λ becomes large, the magnitude of the nullcline curve v n (ζ) in (4.28) becomes small for all ζ.Thus, for λ sufficiently large, the maximum v-value along v n (ζ) is less than v(0).Thus, for such λ, there cannot be a connection from (0, 1) to (v(0), 1) at all; see Figure 2.This is consistent with the fact that there exists an upper bound on |λ| for eigenvalues λ.

Conclusion.
Since essential instabilities have been ruled out, it follows that a single spike traveling wave solution θ 0 (ξ), ξ = ct − x, of (1.6) is linearly stable if and only if for every complex number λ = 0 with Reλ ≥ 0, there exists no nontrivial, bounded, C 1 solution v(ξ) of equation ( 2.18) such that v(ξ) → 0 as ξ → ±∞, and further dθ 0 /dξ is the unique such solution for λ = 0, up to constant multiplication.Such a solution exists exactly for those values of λ for which equation (4.22) is satisfied, where A(ξ, λ) is given by (4.23) and B(ξ, λ) = B(ξ, λ)/v(0) for B(ξ, λ) in (4.24).Because the existence proof for fast and slow traveling waves in [13] proceeded by shooting, the wave speeds c for these waves are not precisely known, which of course complicates the solution of equation (4.22).It is not surprising that speed plays a crucial role in the stability analysis here, since stability analysis based on spike time perturbations has shown that, in related models, fast waves are stable and slow waves are unstable [2,3,6].Indeed, numerical evidence suggests that the fast wave is stable and the slow wave is unstable in the θ-model as well [13].
In addition to deriving the eigenvalue equations (2.18) and (4.22) and ruling out essential instabilities for traveling waves, I have shown here that the derivative of a traveling wave is an eigenfunction for eigenvalue λ = 0, that there is an upper bound M > 0 on |λ| for eigenvalues λ, and that any eigenfunctions corresponding to λ ∈ IR + are nonmonotone on (−∞, 0), at least for a special choice of J and α.Complete characterization of the set of eigenvalues for (1.6) remains open.In addressing this task, it may be helpful to note that the eigenvalue equation (4.22) is equivalent to the equation E(λ) = 0 on the subset of C that constitutes the domain of the Evans function E(λ) for equation (2.18).This Evans function can be defined by an immediate extension of the construction given in [10].
In the traveling waves considered here, the function α(t) may be nondifferentiable, or even discontinuous, at t = 0.This corresponds to a synaptic input that turns on at a discrete spike time.While this feature may be relatively specialized, eigenvalue equations with the general form of equations (2.17), (2.18) also arise in other contexts.For example, nonlocal eigenvalue problems, featuring exponential dependence on the eigenvalue parameter in integral terms, appear in the study of certain partial function differential equations, such as nonlinear reaction-diffusion equations with delay [16], and in certain approaches to analyzing the stability of viscous shock waves [8].Thus, it is expected that the future study of equation (2.18) may yield insight into stability problems for a variety of dynamical systems.
Acknowledgments: I am very grateful to Remus Oşan, Bill Troy, and especially Bard Ermentrout for useful discussions.I also thank Bard for pointing me to reference [11], Björn Sandstede for directing me to references [10,16], and Peter Howard for his comments on viscous shock waves.where b > 0 denotes the synaptic decay rate, such that α(t) = e −bt for t ≥ 0. Here, I will linearize equation (6.29) about perturbations (v(ξ, t), σ(ξ, t)) to (θ 0 (ξ), α(ξ/c)), which constitutes a stationary solution to (6.29), for comparison to the linearization performed in Section 2. Note that such a linearization does not allow for perturbations to spike times, but does allow for perturbations to the time course of synaptic input.As previously, the constant g syn will henceforth be absorbed into J.

Theorem 4 . 3 .
A complex number λ ∈ C is an eigenvalue for equation (2.18) if and only if it satisfies the equation
.18), in the nonlocal part N , however.Nonetheless, if this term decays as ξ → ±∞, as stated in equation (4.21), then the nonlocal term still represents a compact perturbation of an appropriate local problem.This allows for the definition of the Evans function for equation (2.18) when (4.21) holds, which is true in particular for Reλ ≥ 0. This extension yields an Evans function E(λ) with an open domain D that includes {λ ∈ C : Reλ > 0}.The domain D can be extended to include {λ ∈ C : Reλ ≥ 0} for a choice of J with sufficiently rapid decay as ξ → ±∞.For λ ∈ D, E(λ) = 0 if and only if λ satisfies equation (4.22).