A geometric method for eigenvalue problems with low rank perturbations

TL;DR

This paper introduces a geometric approach leveraging differential geometry to analyze the spectrum of operators with low-rank non-normal perturbations, applicable to various models in applied mathematics.

Contribution

It presents a novel geometric method for explicitly determining spectra of low-rank perturbed operators, extending classical techniques to non-normal cases.

Findings

Successfully analyzed three applied models using the geometric method.

Provided explicit spectral characterizations for low-rank perturbations.

Demonstrated the method's effectiveness in non-normal operator contexts.

Abstract

We consider the problem of finding the spectrum of an operator taking the form of a low-rank (rank one or two) non-normal perturbation of a well-understood operator, motivated by a number of problems of applied interest which take this form. We use the fact that the system is a low rank perturbation of a solved problem, together with a simple idea of classical differential geometry (the envelope of a family of curves) to completely analyze the spectrum. We use these techniques to analyze three problems of this form: a model of the oculomotor integrator due to Anastasio and Gad (2007), a continuum integrator model, and a nonlocal model of phase separation due to Rubinstein and Sternberg (1992).

Equations

$$\tilde{\mathbf{M}}\vec{w}=\mathbf{M}\vec{w}+{\rho }_1{\vec{f}}_1⟨{\vec{g}}_1,\vec{w}⟩+{\rho }_2{\vec{f}}_2⟨{\vec{g}}_2,\vec{w}⟩=\lambda \vec{w}$$

$${\mathrm{disc}}_{\lambda }(\det (\tilde{\mathbf{M}}-\lambda \mathbf{I}))=P({\rho }_1,{\rho }_2),$$

$$\mathcal{V}=\{({\rho }_1,{\rho }_2)|P({\rho }_1,{\rho }_2)=0\}$$

$$\tilde{D}(\lambda )=\det (\widetilde{\text{mbox}\mathbf{M}}-\lambda \mathbf{I})$$

$$\det (\tilde{}\text{mbox}\mathbf{M}-\lambda \mathbf{I})=D(\lambda )+P_1(\lambda ){\rho }_1+P_2(\lambda ){\rho }_2+Q(\lambda ){\rho }_1{\rho }_2$$

$$P_1(\lambda )$$

$$P_2(\lambda )$$

$$Q(\lambda )$$

$$P_1(\lambda )$$

$$P_2(\lambda )$$

$$Q(\lambda )$$

$$0$$

$$0=1+{\rho }_1⟨{\vec{g}}_1,(\text{mbox}\mathbf{M}-\lambda \mathbf{I}{)}^{-1}{\vec{f}}_1⟩+{\rho }_2⟨{\vec{g}}_2,(\text{mbox}\mathbf{M}-\lambda \mathbf{I}{)}^{-1}{\vec{f}}_2⟩.$$

$$D(0)+{\rho }_1{P}_1(0)+{\rho }_2{P}_2(0)+{\rho }_1{\rho }_{2}Q(0)=0$$

$$\Rightarrow$$

$$D(\lambda )+{\rho }_1{P}_1(\lambda )+{\rho }_2{P}_2(\lambda )+{\rho }_1{\rho }_{2}Q(\lambda )=0$$

$$D^{\prime }(\lambda )+{\rho }_1{P}_1^{\prime }(\lambda )+{\rho }_2{P}_2^{\prime }(\lambda )+{\rho }_1{\rho }_2{Q}^{\prime }(\lambda )=0.$$

$${\rho }_1$$

$${\rho }_2$$

$${\rho }_1$$

$${\rho }_2$$

$$\operatorname{rank}\left(\begin{array}[]{ccc}P_{1}&P_{2}&Q\\ P_{1}^{\prime}&P_{2}^{\prime}&Q^{\prime}\end{array}\right)<\operatorname{rank}\left(\begin{array}[]{cccc}D&P_{1}&P_{2}&Q\\ D^{\prime}&P_{1}^{\prime}&P_{2}^{\prime}&Q^{\prime}\end{array}\right).$$

$$\displaystyle{\bf C1}~{}~{}~{}~{}~{}\operatorname{rank}\left(\begin{array}[]{cccc}D&P_{1}&P_{2}&Q\\ D^{\prime}&P_{1}^{\prime}&P_{2}^{\prime}&Q^{\prime}\end{array}\right)$$

$$\displaystyle{\bf C2}~{}~{}~{}~{}~{}\left\{\begin{array}[]{c}P_{1}\wedge Q(\lambda)\neq 0\\ P_{2}\wedge Q(\lambda)=0\\ D\wedge P_{1}(\lambda)=0\\ P_{1}\wedge P_{2}(\lambda)=D\wedge Q(\lambda)\end{array}\right.$$

$$\displaystyle{\bf C3}~{}~{}~{}~{}~{}\left\{\begin{array}[]{c}P_{2}\wedge Q(\lambda)\neq 0\\ P_{1}\wedge Q(\lambda)=0\\ D\wedge P_{2}(\lambda)=0\\ P_{2}\wedge P_{1}(\lambda)=D\wedge Q(\lambda)\end{array}\right.$$

$$F({\rho }_1,{\rho }_2,\lambda )=0$$

$$\frac{\partial F}{\partial \lambda }({\rho }_1,{\rho }_2,\lambda )=0.$$

$$\frac{1}{2}\frac{{\partial }^{2}F}{\partial {\lambda }^2}(\delta \lambda {)}^2+{\nabla }_{{\rho }_1,{\rho }_2}F\cdot (\delta {\rho }_1,\delta {\rho }_2)=O(\delta \rho \delta \lambda ,\delta {\rho }^2,\delta {\lambda }^3)$$

$$\frac{1}{2}\frac{{\partial }^{2}F}{\partial {\lambda }^2}(\delta \lambda {)}^2+{\nabla }_{{\rho }_1,{\rho }_2}F\cdot (\delta {\rho }_1,\delta {\rho }_2)$$

$$\left|\begin{array}[]{cc}P_{1}+\rho_{2}Q&P_{2}+\rho_{1}Q\\ P_{1}^{\prime}+\rho_{2}Q^{\prime}&P_{2}^{\prime}+\rho_{1}Q^{\prime}\end{array}\right|\neq 0.$$

$$(P_1\wedge P_2\wedge Q)(P_1\wedge P_2\wedge D)$$

$$P_1\wedge P_2\wedge Q$$

$$P_1\wedge P_2\wedge Q=P_1\wedge P_2\wedge D=D\wedge P_2\wedge Q=P_1\wedge D\wedge Q=0.$$

$${\rho }_1^{\prime }(P_1(\lambda )+{\rho }_{2}Q(\lambda ))+{\rho }_2^{\prime }(P_2(\lambda )+{\rho }_{1}Q(\lambda ))$$

$${\rho }_1^{\prime }(P_1^{\prime }(\lambda )+{\rho }_2{Q}^{\prime }(\lambda ))+{\rho }_2^{\prime }(P_2^{\prime }(\lambda )+{\rho }_1{Q}^{\prime }(\lambda ))$$

$$(D^{\prime \prime }(\lambda )+{\rho }_1{P}_1^{\prime \prime }(\lambda )+{\rho }_2{P}_2^{\prime \prime }(\lambda )+{\rho }_1{\rho }_2{Q}^{\prime \prime })=0.$$

$$\left|\begin{array}[]{cc}P_{1}+\rho_{2}Q&P_{2}+\rho_{1}Q\\ P_{1}^{\prime}+\rho_{2}Q^{\prime}&P_{2}^{\prime}+\rho_{1}Q^{\prime}\end{array}\right|\neq 0$$

Full text

Keywords: Bifurcation theory, Aronszajn-Krein formula, Rank one perturbations

Introduction

In this paper we analyze eigenvalue problems of the following form

[TABLE]

with $\rho_{1}$ and $\rho_{2}$ parameters, $\vec{f}_{i},\vec{g}_{i}$ fixed vectors and ${\bf M}$ is an operator with known spectrum. While Eqn. (1) might appear very specific, we are aware of a number of interesting eigenvalue problems which take this form. These include:

• •

A model due to Anastasio and Gad of the behavior of the oculomotor integrator[1].

• •

A non-local Allen-Cahn model due to Rubinstein and Sternberg for phase separation[2].

• •

The stability problem for spike solutions to activator-inhibitor models in the limit of slow activator diffusion[3, 4].

• •

Stability for models of runaway ohmic heating[5, 6, 7] and microwave heating [8].

• •

Stability for stationary solutions of model for phytoplankton growth.[9]

Four of the models take the form of a reaction-diffusion equation with a nonlocal term. The study of the stability of stationary solutions of such models naturally leads to an eigenvalue problem which takes the form of a self-adjoint Sturm-Liouville operator plus a finite rank perturbation coming from the nonlocal term. Freitas[3] has considered a similar problem and has some related results for a single perturbation (rather than a two parameter family); Bose and Kriegsmann[8] have some other related results, mainly in the case where $\vec{f}_{i}=\vec{g}_{i}$ where the problem is self-adjoint. We refer the interested reader to the review paper of Freitas[10], which details a number of models whose stability problems take the form of a “nice” operator with a low rank (typically rank one) perturbation.

In all of these problems the eigenvalue problem arises in the study of the stability of a particular steady state. In this situation one is typically interested in understanding qualitative properties of the spectrum as a function of the parameters $(\rho_{1},\rho_{2})$. In particular one might wish to understand, for any particular pair $(\rho_{1},\rho_{2})$,

• •

how many eigenvalues are in the right half-plane, and

• •

how many eigenvalues are real (vs. complex).

Here we give a direct way to construct a phase-diagram in the $(\rho_{1},\rho_{2})$ plane which answers these questions.

The main approach that we take here is to exploit the low rank nature of the perturbations, along with some geometric constructions for the quantities of interest. Since the eigenvalues will vary continuously as a function of the parameters $(\rho_{1},\rho_{2})$, the quantities noted above are constant on open sets, with the boundary of these sets being respectively

• •

the set of $(\rho_{1},\rho_{2})$ for which $\widetilde{\bf M}$ has a purely complex eigenvalue $\lambda=\imath\omega$ (including $\lambda=0$), and

• •

the set of $(\rho_{1},\rho_{2})$ for which $\widetilde{\bf M}$ has a double real eigenvalue.

Knowledge of these boundary sets, together with knowledge of the spectrum of the unperturbed operator ${\bf M}$, would therefore enable us to study stability in the entire plane.

In passing we note that, while the unperturbed operator ${\bf M}$ is self-adjoint in some of the examples presented here, self-adjointness is not strictly necessary. What is necessary is that the spectrum of the unperturbed operator ${\bf M}$ is known (at least qualitatively), so that we have a baseline with which to compare the perturbed operator (much as a constant of integration fixes a particular solution of a differential equation). In each example studied here, the unperturbed operator ${\bf M}$ will have a purely real spectrum.

Our main result in this paper will be to identify – and gave a recipe for computing – a set of geometric quantities associated with the spectrum of a rank-two perturbation of a well-known operator, which can be used to analyze the perturbed operator in the entire plane: this is done in §1. We will then apply our technique to three specific problems which can be written in the form of Eqn. (1). The first is a model of a coupled brainstem-cerebellum neuronal network called the oculomotor integrator (§2); the second is a continuum version of that model, in which the (relatively numerous) brainstem neurons are replaced by a neural “line” (§3). Finally, we analyze a stability problem that arises in a nonlocal reaction-diffusion equation [2] (§4). In this last problem, we also use an intermediate result from §1 (the Aronszajn-Krein formula and its consequences) to prove a new theorem about stability of stationary solutions.

1 Basic Calculations

We begin with some general dimension-counting arguments. The matrix111For purposes of exposition we will consider the case where $\widetilde{\bf M}$ is a matrix, but everything we say will apply equally to the case where $\widetilde{\bf M}$ is an operator of compact resolvent. $\widetilde{\bf M}$ is assumed to be a real $N\times N$ matrix. Real non-symmetric matrices will generically have a real eigenvalue of multiplicity higher than one on a set of codimension one. In a two-parameter model such as we are considering here this codimension one set divides the parameter space into open sets having a constant number of real eigenvalues. As one crosses this set the number of real eigenvalues changes by (generically) two. This motivates the following definition:

Definition 1**.**

We define the bifurcation curve to be the locus of points ${\mathcal{V}}=(\rho_{1},\rho_{2})$ for which $\widetilde{\bf M}$ has a real eigenvalue of multiplicity two or higher.

For matrix problems, of course, there exists an algebraic procedure for determining the values in the $(\rho_{1},\rho_{2})$ plane where the matrix has multiple eigenvalues. One can simply compute the discriminant (in $\lambda$) of the characteristic polynomial of the matrix $\widetilde{\bf M},$

[TABLE]

which gives a polynomial in the parameters $(\rho_{1},\rho_{2})$. The variety defined by the zero set of this polynomial

[TABLE]

determines the bifurcation curve. Unfortunately, this computation isn’t practical to carry out analytically for real problems: for a large matrix, $P(\rho_{1},\rho_{2})$ will be a polynomial of large degree and the zero set will be difficult to compute. For the case of operators, even very nice ones, it is not clear that the discriminant makes sense at all.

Instead we use the fact that the perturbations are of finite rank to give an explicit rational or algebraic parameterization of the bifurcation curve. We begin by stating a preliminary lemma, which is basically the Aronszajn-Krein formula for rank one perturbations:

Lemma 1**.**

Let $\widetilde{\,\mbox{\bf M}}$ be an $N\times N$ matrix defined as in Equation (1). The characteristic polynomial of $\widetilde{\,\mbox{\bf M}}$

[TABLE]

takes the following form:

[TABLE]

where $D(\lambda)=\det({\bf M}-\lambda{\bf I})$ is the determinant of the unperturbed problem, $P_{i}(\lambda)$ are polynomials of degree (at most) $(N-1)$ and $Q(\lambda)$ of degree (at most) $N-2$. In the case where $\vec{g}_{1}$ and $\vec{g}_{2}$ (or $\vec{f}_{1,2}$) are linearly dependent $Q(\lambda)=0$.

Proof.

Due to the rank two nature of the perturbation, the characteristic polynomial can contain no powers of $\rho_{1}$ or $\rho_{2}$ above the first. The easiest way to see this is via multilinear algebra. The determinant is clearly polynomial in $\lambda,\rho_{1},\rho_{2}$. A general term of the form $\rho_{1}^{j}\rho_{2}^{k}$ comes from the wedge product of $j$ factors of $\vec{g}_{1}$, $k$ factors of $\vec{g}_{2}$ and ($N-(j+k)$) columns from $\,\mbox{\bf M}-\lambda{\bf I}$. Any term with more than one factor of $\vec{g}_{1}$ and one factor of $\vec{g}_{2}$ must be zero. Hence the determinant is of the form given in Eqn. (2).

Finally, if $\vec{g}_{1}$ and $\vec{g}_{2}$ are linearly dependent, then any wedge product including both $\vec{g}_{1}$ and $\vec{g}_{2}$ will vanish, so that $Q(\lambda)\equiv 0$.

The explicit form of the polynomials $P_{i}(\lambda),Q(\lambda)$ is easy to compute from the above construction (a detailed derivation is provided in the Appendix). If we let $\operatorname{cof}^{t}\left(\,\mbox{\bf M}-\lambda{\bf I}\right)$ denote the transpose cofactor matrix of ${\bf M}-\lambda{\bf I}$ then one has the following formulae

[TABLE]

If in addition  M is self-adjoint, we have alternative formulae, in which the polynomial form of $Q(\lambda)$, etc., is even more evident:

[TABLE]

where $\lambda_{i}$ and $\vec{\phi}_{i}$ are the eigenvalues and eigenvectors of the unperturbed matrix ${\bf M}$. ∎

It is often more convenient to divide Eq. (2) by $\det(\,\mbox{\bf M}-\lambda{\bf I})$ to put the eigenvalue condition in the form

[TABLE]

or even more simply, when $Q(\lambda)=0$,

[TABLE]

This form has the advantage that it is expressed in terms of resolvents (i.e. $R_{\lambda}\equiv(\,\mbox{\bf M}-\lambda{\bf I})^{-1}$), which are defined very generally for operators, rather than determinants and cofactors, which are not.

One geometric way to interpret this characteristic polynomial is as defining a one parameter family of rational curves, the curves of constant eigenvalue. For each value of $\lambda$, Eqn. (2) defines a curve in the $(\rho_{1},\rho_{2})$ plane along which $\lambda$ is an eigenvalue. For instance the matrix has zero as an eigenvalue along the curve

[TABLE]

in the $(\rho_{1},\rho_{2})$ plane. In the special case where $Q(\lambda)\equiv 0$, Eqn. (2) defines a one-parameter family of lines.

Given a family of curves it is always fruitful to consider its envelope, which is the simultaneous solution to

[TABLE]

Whereas each curve in the family encodes information on the location of a particular eigenvalue, the envelope of the curve encodes information on eigenvalue coincidence: this is the content of the next lemma. To simplify notation, we use the wedge notation for Wronskians: $f\wedge g=fg^{\prime}-f^{\prime}g$.

Lemma 2**.**

The solutions to (9, 10) give the bifurcation curve ${\mathcal{V}}$. If $Q(\lambda)$ is not identically zero the bifurcation curve $\mathcal{V}$ is generically given by the union of the pair of parametric curves given by

[TABLE]

If $Q(\lambda)$ is identically zero then the bifurcation curve is generically given by the parametric curve

[TABLE]

Proof.

The equivalence of the envelope and the discriminant is standard – see for instance Bruce and Giblin[11] or Spivak[12]. Generically one just has to solve Equations (9, 10) for $\rho_{1},\rho_{2}$. This is equivalent to solving a linear and a quadratic equation in the general case and a pair of linear equations in the special case $Q(\lambda)=0.$

For certain values of $\lambda$ the system (9, 10) maybe be inconsistent, or consistent but underdetermined. Inconsistency indicates that this eigenvalue cannot be achieved by any choice of $\rho_{1,2}$; this will occur if

[TABLE]

If the system is consistent but underdetermined, then there may be a curve in the $(\rho_{1},\rho_{2})$ plane along which $\lambda$ is a multiple eigenvalue. The system (9, 10) is consistent and underdetermined for $\lambda$ if one of the following conditions ${\bf C1,C2,C3}$ hold:

[TABLE]

If any of these genericity conditions are satisfied for some value of $\lambda$ then (generically) there exists a curve in the $(\rho_{1},\rho_{2})$ plane along which that value of $\lambda$ is a multiple eigenvalue. If they are not satisfied for any $\lambda$, then the bifurcation curve is equal to the envelope curve. Note that one can always check whether or not a pair of polynomials have a common root by computing the resultant of the polynomials: one need not be able to explicitly factor the polynomials to test this condition. Thus the genericity condition is readily checkable. ∎

Lemma 3**.**

As $(\rho_{1},\rho_{2})$ are varied so as to cross the envelope the number of real eigenvalues generically changes by two.

Proof.

We assume that most readers are familiar with this phenomenon from the theory of first order quasilinear partial differential equations, where the characteristic curves form a one parameter family of curves and the envelope (caustic) marks the transition between regions which are (typically) singly and triply covered by characteristics, but we give a short proof.

The envelope curve is defined by the simultaneous solution to

[TABLE]

Assume that $(\rho_{1}^{*},\rho_{2}^{*})$ is a point along this curve with corresponding eigenvalue $\lambda^{*}$, and that the gradient $\nabla_{\rho_{1},\rho_{2}}F(\rho_{1}^{*},\rho_{2}^{*},\lambda^{*})$ and the second derivative $\frac{\partial^{2}F}{\partial\lambda^{2}}(\rho_{1}^{*},\rho_{2}^{*},\lambda^{*})$ are non-vanishing. Expanding in a neighborhood of this point, i.e. letting $\rho_{1}=\rho_{1}^{*}+\delta\rho_{1},\rho_{2}=\rho_{2}^{*}+\delta\rho_{2}$ and $\lambda=\lambda^{*}+\delta\lambda$, we find the normal form

[TABLE]

By the Weierstrauss preparation theorem, we can approximate solutions to $F(\rho_{1},\rho_{2},\lambda)=0$ in a neighborhood of the point $(\rho_{1}^{*},\rho_{2}^{*})$ by fixing $(\delta\rho_{1},\delta\rho_{2})$ and solving for $\delta\lambda$; i.e. by solving the quadratic equation

[TABLE]

for $\delta\lambda$. The tangent line to the envelope at $(\rho_{1}^{*},\rho_{2}^{*})$ is given by $\nabla_{\rho_{1},\rho_{2}}F\cdot(\delta\rho_{1},\delta\rho_{2})=0$; on one side of this line, Eqn. (28) has two distinct real roots $\pm\delta\lambda$, so each point lies on two constant eigenvalue curves, $\lambda^{*}\pm\delta\lambda$. On the other side of the tangent line there are a pair of complex conjugate roots. This general picture is illustrated in Figure (1), which shows a close-up of the envelope curve from an example to follow.

For more details on the envelope see the text of Bruce and Giblin[11]. ∎

**Figure 1 here

**

We would next like to consider the possibility of eigenvalues of higher multiplicity. The envelope curve is, in general, a well-behaved curve and admits a parametrization by arc length. However this may fail at an isolated set of points in $(\rho_{1},\rho_{2})$. The next lemma says that (modulo some genericity assumptions) the following are all equivalent, and occur on a codimension two set (isolated points in the $(\rho_{1},\rho_{2})$ plane):

• •

Points in the $(\rho_{1},\rho_{2})$ plane where the model has a real eigenvalue of multiplicity at least three.

• •

Points where the tangent vector to the envelope curve vanishes.

• •

Cusps in the envelope curve.

Lemma 4**.**

The vanishing of the tangent vector to the envelope curve at a point implies that $\widetilde{\bf M}$ has an eigenvalue of multiplicity (at least) three at that point. The converse holds as long as the following determinant is non-zero at the point in question:

[TABLE]

Alternatively, the model has a triple eigenvalue if and only if either

[TABLE]

or

[TABLE]

Proof.

The conditions for an eigenvalue of multiplicity (at least) two are given by (9,10). Differentiating each of these with respect to $\lambda$ gives the following equations for $\rho_{1}^{\prime},\rho_{2}^{\prime}$:

[TABLE]

The conditions for an eigenvalue of multiplicity at least three are given by (9,10) together with the condition

[TABLE]

Thus it is clear that $\rho_{1}^{\prime}=0,\rho_{2}^{\prime}=0$ implies that the eigenvalue is of multiplicity (at least) three. Further if

[TABLE]

then the above system can be solved uniquely for $\rho_{1}^{\prime},\rho_{2}^{\prime}$ and the existence of an eigenvalue of multiplicity three implies $\rho_{1}^{\prime}=0,\rho_{2}^{\prime}=0.$

A bit more algebra gives another characterization of points where the eigenvalue has multiplicity three or higher. Equations (9,10,29) form a system of three equations in three unknowns $\rho_{1},\rho_{2},$ and $\rho_{3}=\rho_{1}\rho_{2}$. Solving these three equations for $(\rho_{1},\rho_{2},\rho_{3})$ and imposing the consistency condition $\rho_{3}=\rho_{1}\rho_{2}$ shows that one has a root of multiplicity three if and only if either

[TABLE]

or all of the Wronskians

[TABLE]

vanish. The first possibility is typically of codimension two - it is expected to occur at isolated values of $\lambda$ corresponding to isolated values of $(\rho_{1},\rho_{2})$. The second does not typically happen at all, since it requires the simultaneous vanishing of several polynomials. However in Example (1) this case occurs because it is forced by a symmetry of the model.

Finally recall that a simple zero of the tangent vector represents a cusp, and generically an eigenvalue of multiplicity at least three will have multiplicity exactly three, so typically cusps in the envelope are equivalent to triple eigenvalues. ∎

The geometry of a bifurcation in the neighborhood of a triple eigenvalue is illustrated in Figure 4B. In a neighborhood of this point there are three dominant eigenvalues which participate in the bifurcation. The cusp of the envelope represents a transition between a bifurcation between the intermediate and the smallest eigenvalue in the trio, and a bifurcation between the intermediate and the largest eigenvalue in the trio. Emerging from the cusp-point is a curve that represents an exchange of dominance phenomenon, with a complex conjugate pair of eigenvalues crossing a single real one.

When considering questions of stability and the behavior of the dominant eigenvalue it is also important to understand the behavior of the complex eigenvalues. In particular one would like to understand the locus of points at which the matrix has purely imaginary eigenvalues, as this curve indicates where the model loses stability due to a Hopf bifurcation.

Definition 2**.**

The Hopf curve is the locus of points in the $(\rho_{1},\rho_{2})$ plane where $\widetilde{\bf M}$ has a pair of purely imaginary eigenvalues. Generically this curve is given parametrically by

[TABLE]

where ${\rm Re,Im}$ represent the real and imaginary parts respectively. The genericity conditions are the same as in Lemma (2) with the replacement of the Wronskians $f\wedge g$ by the quantities ${\rm Re}(f(i\omega)){\rm Im}(g(i\omega))-{\rm Re}(g(i\omega)){\rm Im}(f(i\omega)).$

The Hopf curve, the envelope, and the $\lambda=0$ eigenvalue curve will all intersect at a single point, the point at which there is a zero eigenvalue of higher multiplicity. We now present an illustrative example.

Example 1**.**

We consider the following model:

[TABLE]

It is straightforward to compute that the characteristic polynomial of this matrix is given by

[TABLE]

The zero eigenvalue curve is given by

[TABLE]

The bifurcation curve is given by the envelope

[TABLE]

together with the singular piece $\rho_{1}=-\frac{1}{2}\cup\rho_{2}=-\frac{1}{2}$, which is associated to the value $\lambda=-2+\frac{\sqrt{2}}{2}$, where the equations defining the envelope fail to have full rank. The envelope and the singular piece of the bifurcation curve meet tangentially at $(\rho_{2}=-\frac{1}{2},\rho_{1}=-\frac{3}{2})$ and $(\rho_{2}=-\frac{1}{2},\rho_{1}=-\frac{3}{2})$. Because of the symmetry we have $P_{1}=P_{2}$ and thus $P_{1}\wedge P_{2}\wedge Q\equiv 0$, so the condition for a triple eigenvalue reduces to simultaneous vanishing of $P_{1}\wedge D\wedge Q$ and $D\wedge P_{2}\wedge Q.$ Note that since $P_{1}=P_{2}$ these are not independent — $P_{1}\wedge D\wedge Q=-D\wedge P_{2}\wedge Q.$ Calculating we find that the triple eigenvalue condition becomes

[TABLE]

This cubic has three real roots: a double root at $\lambda=-2+\frac{\sqrt{2}}{2}$ and a simple root at $\lambda=-(\frac{1}{2}+\frac{\sqrt{2}}{4})\approx-2.35.$ The envelope is not defined for $\lambda\in(-\frac{5}{2},-\frac{3}{2}),$ so the root at $\lambda=-(\frac{1}{2}+\frac{\sqrt{2}}{4})$ does not correspond to a real multiple eigenvalue. Thus the only real eigenvalue of multiplicity higher than two is $\lambda=-2+\frac{\sqrt{2}}{2}.$ Since this eigenvalue is associated to the singular piece of the bifurcation curve we can potentially have many points where this is a triple eigenvalue. Along the curve $\rho_{1}=-\frac{1}{2}$ the eigenvalues are

[TABLE]

So the only triple eigenvalue is at $\rho_{2}=-\frac{3}{2}$, the point of intersection with the envelope curve. A similar calculation holds along $\rho_{2}=-\frac{1}{2}.$

The Hopf curve is given parametrically by

[TABLE]

where, as always, the signs are not independent. Note that the argument of the square root is strictly positive, so there exists purely imaginary eigenvalues corresponding to oscillations of any desired frequency. The Hopf curve, the envelope, and the zero eigenvalue line all meet at the points $(\rho_{1}=(4+\sqrt{2})(-1\pm\frac{\sqrt{30(18-8\sqrt{2})}}{14})=4+\sqrt{2}\pm\sqrt{30},\rho_{2}=(4+\sqrt{2})(-1\mp\frac{\sqrt{30(18-8\sqrt{2})}}{14})=-(4+\sqrt{2})\mp\sqrt{30}.$ The most interesting region of the stability diagram is depicted in Figure (2). The zero eigenvalue curve is depicted in dashed red, the envelope in blue (including a dot at the origin), the singular piece of the bifurcation curve in dot-dashed magenta, and the Hopf curve in solid dotted green.

From this information it is easy to derive the stability diagram. At the origin the eigenvalues are $\lambda=-3,\lambda=-2,\lambda=-2,\lambda=-1$. Since there is a degenerate eigenvalue one needs to do a local perturbation analysis near $\lambda=-2,\rho_{1}=0,\rho_{2}=0$ to determine if in the neighborhood of this point one has a real pair of eigenvalues or a complex conjugate pair. Letting $\lambda=-2+\delta$ shows that near this point one has

[TABLE]

where $O(3)$ denotes terms of order three or higher in $\delta,\rho_{i}.$ The discriminant of the above is $(\sqrt{2}\rho_{1}+\sqrt{2}\rho_{2})^{2}-4\rho_{1}\rho_{2}=2\rho_{1}^{2}+2\rho_{2}^{2}>0,$ indicating that in a neighborhood of the origin the double eigenvalue splits into a real (distinct) pair of eigenvalues. Thus in the region containing the origin and bounded by the singular pieces of the bifurcation curve and the upper branch of the envelope (labelled ${\bf A}$) there are four real eigenvalues in the left half-plane. As $\rho_{2}$ is decreased so as to cross the line $\rho_{2}=-\frac{1}{2}$ the first bifurcation occurs. Since this line corresponds to eigenvalue $\lambda=-2+\sqrt{2}$ and $-2<-2+\frac{\sqrt{2}}{2}<-1$ the bifurcation consists of the two dominant real eigenvalues bifurcating to a complex conjugate pair. Thus in region ${\bf B}$ we have two real eigenvalues and two complex eigenvalues, all in the left half-plane. As one leaves region B across the Hopf curve into the region labelled ${\bf E}$ the complex conjugate pair moves into the right half-plane, giving two complex eigenvalues in the right half-plane and two real eigenvalues in the left half-plane. Proceeding in this fashion the stability diagram can be labelled as follows:

• •

Region A: Four real eigenvalues in the left half-plane.

• •

Region B: Two real and two complex eigenvalues in the left half-plane.

• •

Region C: Four complex eigenvalues in the left half-plane.

• •

Region D: Two complex eigenvalues in right half-plane, two complex eigenvalues in the left half-plane.

• •

Region E: Two complex eigenvalues in the right half-plane, two real eigenvalues in the left half-plane.

• •

Region F: One real eigenvalue in the right half-plane, three real eigenvalues in the left half-plane.

• •

Region G: One real eigenvalue in the right half-plane, one real and two complex eigenvalues in the left half-plane.

Additionally there is a narrow region between the regions labelled E and F (to the left of $\rho_{2}=-(4+\sqrt{2})-\sqrt{30}\approx-10.9$ above the zero eigenvalue curve and below the envelope curve) where there are two real eigenvalues in the left half-plane and two real eigenvalues in the right half-plane. This region is not labelled since it is not visible on this scale.

One feature which we have not labelled are points where a real eigenvalue and a complex conjugate pair all have the same real part, corresponding to points where a real eigenvalue and a complex conjugate pair exchange dominance. Although it is easy to write down an implicit equation satisfied by these curves, it is generally not possible to find an explicit formula as for the Hopf curve, envelope, etc.

Since the characteristic polynomial of this problem is of order four it would in principle be possible to extract the above information directly from the solution formula for the quartic. In practice it would be exceedingly difficult to recover such detailed information. In the next section we will consider a model that arises from a differential equation of order eight; in this situation it is no longer possible even in principle to state a direct formula for the roots.

**Figure 2 here

**

In this section, we identified – and gave a recipe for computing – a set of geometric quantities associated with the spectrum of a rank-two perturbation of a well-known operator. Because some of these quantities (specifically the bifurcation curve, the $\lambda=0$ eigenvalue curve, and the Hopf curve) form curves which together partition the $(\rho_{1},\rho_{2})$ plane into open sets with qualitatively similar behavior, this information allows us to analyze the perturbed operator in the entire plane. We now apply this procedure to three specific problems which can be written in the form of Eqn. (1): a model of a coupled brainstem-cerebellum neuronal network called the oculomotor integrator (§2); a continuum version of that model, in which the (relatively numerous) brainstem neurons are replaced by a neural “line” (§3); and a stability problem that arises in a nonlocal reaction-diffusion equation [2] (§4).

2 Example: a model for the oculomotor integrator

The oculomotor integrator is a neural network that is essential for eye movement control. This network holds your gaze steady despite transient body motions, by using cues it receives from oculomotor subsystems such as the vestibular system (which processes sensory input from the semi-circular canals of the inner ear) [13, 14]. Here, we present a model for the oculomotor integrator which can either simulate normal integrator function, or a common eye movement disorder — congenital or infantile nystagmus (IN) — with a few small changes. We are able to make this determination by analysis of the envelope, constant eigenvalue curves and Hopf curve as described in §1. Some of this analysis was presented in [15], which focused on using the model to simulate various eye movements associated with IN; here, we focus on analyzing the full phase space generated by the low-rank perturbations.

We first summarize how the integrator works. The neurons that compose the integrator are located mainly in a brainstem region known as the vestibular nucleus. These neurons must integrate velocity signals into a desired position; thus, they perform the operation of integration with respect to time (or temporal integration). Single neurons produce a small amount of temporal integration, in that a firing rate increase due to a transient input will decay with a time constant of about 5 ms. The observed integrator time constant is closer to 20 seconds (s), so the integrator network must lengthen the single-neuron time constant by about 4000 times [14]; for a linear system, this corresponds to having an eigenvalue near zero [16]. Furthermore, the integrator should have the right gain in response to velocity signals; the ratio of its response to the input should be appropriate. It should also be plastic; i.e. it should be able to adjust, if injury or some other change occurs (for example, you adjust your oculomotor integrator gain when you get a new pair of glasses) [17, 18, 19, 20].

In order for the oculomotor integrator to function, the vestibular nucleus must be connected to a second major brain region, the cerebellum[21, 22, 23]. These connections are essential both for normal operation and for plasticity [24, 25, 26, 27]. Neuroanatomical research has shown that these connections are asymmetric in an important way: while the connections from the vestibular nucleus to the cerebellum are numerous and excitatory, the feedback connections from the cerebellum are sparse and are inhibitory [28, 29, 30, 31, 32, 33]. Anastasio and Gad [1] showed previously that this asymmetry permits the cerebellum to sensitively control both the time constant and the gain of the oculomotor integrator, despite (or perhaps because) the projections from the cerebellum back to the vestibular nuclei are so sparse. These sparse cerebellar-to-vestibular connections, which can plastically change their strength, are precisely the low-rank perturbations we analyze in detail here.

Anastasio and Gad [1] proposed a linear differential equation model that combined a vestibular network with sparse, asymmetric feedback connections from cerebellar Purkinje cells, i.e., the evolution of the system is given by:

[TABLE]

where $\vec{v}$ represents the response of the system (vestibular neurons and Purkinje cells together), $\widetilde{\,\mbox{\bf M}}$ represents a matrix of connections, $s(t)$ is a scalar input (velocity) signal to the integrator, and $\vec{b}$ a fixed vector representing the pattern in which the vestibular neurons receive the input signal.

The connection matrix $\widetilde{\,\mbox{\bf M}}$ proposed by [1] was:

[TABLE]

The vestibular neuron-to-Purkinje cell coupling vectors are given by the $\vec{w}_{i}$; the Purkinje cell-to vestibular neuron coupling vectors are given by $(\vec{u}_{1})_{j}=\delta_{j,k_{1}},(\vec{u}_{2})_{j}=\delta_{j,k_{2}}$, where $\delta_{j,k}$ is the Kronecker delta.

We state the remaining details for completeness; for a full explanation of the biophysical motivation, please see [1, 15]. The matrix ${\bf T}$ of effective connections between vestibular neurons is given by

[TABLE]

The parameter $\alpha$ is set to $200\,{\rm s}^{-1}$ (corresponding to a typical $5\,{\rm ms}$ membrane time constant); the parameter $\beta$ is fixed so that the vestibular sub-network has a time constant of $0.2\,{\rm s}$ in the absence of cerebellar interaction ($\rho_{1}=\rho_{2}=0$). The largest eigenvalue of ${\bf T}$ is given by $\lambda_{1}/\alpha=-1+\beta(1+2\cos(\frac{\pi}{N+1}))$, where $N$ is the number of vestibular cells (on each side of the network); therefore we choose $\beta=\frac{\lambda_{1}/\alpha+1}{1+2\cos\left(\frac{\pi}{N+1}\right)}$, where $\lambda_{1}=5\,{\rm s}^{-1}$.

As we have noted, we will treat the Purkinje-to-vestibular feedback connections as our low rank perturbations; to relate this model to Eq. (1), take

[TABLE]

where $N$ is the number of vestibular neurons (therefore  M is $(N+2)\times(N+2)$, and $\vec{f}_{i},\vec{g}_{i}$ are in $\mathbb{R}^{N+2}$). We will assume that the output of the integrator is a linear readout of the vestibular neuron responses, which for simplicity we take to be equal to $\vec{b}$: i.e. $\langle\vec{b},\vec{v}(t)\rangle$. Defining the eigenvectors $\vec{e}_{i}$ and adjoint eigenvectors $\vec{f}_{i}$ respectively by

[TABLE]

the linear readout at time $t$ (assuming $\widetilde{\,\mbox{\bf M}}$ diagonalizable) is given by

[TABLE]

When there is a separation of time scales (${\rm Re}(\lambda_{2})\ll{\rm Re}(\lambda_{1})$), the response is largely due to the dominant eigenvalue; we define the gain, $\gamma$, to be the ratio of this response, to the magnitude of the filtered input:

[TABLE]

Thus, $\gamma$ captures how the circuit amplifies — or suppresses — the incoming signal.

Setting this ratio correctly allows the organism to respond appropriately to its environment. However, injury or normal growth can alter the sensory systems that supply inputs to the integrator, and thus create a mismatch between input and desired output. To counteract this, the system must change this gain to compensate. We will allow our network to change, by adjusting the Purkinje-to-vestibular weights $\rho_{1}$ and $\rho_{2}$; this is biologically plausible, since one of the major functions of the cerebellum is to regulate motor plasticity and learning.

The asymmetry of the matrix $\widetilde{\,\mbox{\bf M}}$ is crucial, to allow gain to be adjusted freely; for a symmetric (more generally normal) matrix, $-1\leq\gamma\leq 1$. When can we get big changes in gain, with relatively small changes in $\rho_{1}$ and $\rho_{2}$? It would be ideal for the denominator of $\gamma$ to be near zero (assuming both $\vec{f}_{1}$, $\vec{e}_{1}$ are normalized to unit length): $\langle\vec{f}_{1},\vec{e}_{1}\rangle\ll 1$. However, this corresponds to near-orthogonality of the right and left eigenvectors, which occurs near a double eigenvalue (perfect orthogonality, $\langle\vec{f}_{1},\vec{e}_{1}\rangle=0$, can occur only when there is a degenerate double eigenvalue).

We now show an example of an integrator that can perform normal integration with arbitrary adjustment of gain: let ${\bf T}$ be $6\times 6$ and choose

[TABLE]

where $\vec{\bf e}_{j}$ is the $j$th identity vector (in this case $\vec{\bf e}_{j}\in\mathbb{R}^{6}$).

**Figure 3 here

**

With this choice, $\widetilde{\,\mbox{\bf M}}$ is the reduced matrix of a model with 6 vestibular neurons and 2 Purkinje cells on each side of the bilaterally symmetric network. The only parameters that may vary are $\rho_{1}$ and $\rho_{2}$, the strengths of the Purkinje-to-vestibular connections. In order to fix a certain time constant, $\rho_{2}$ and $\rho_{1}$ should be constrained to lie on the appropriate constant eigenvalue curve. The biologically appropriate time constant is in the neighborhood of $20\,{\rm s}$, so the appropriate eigenvalue is $\lambda=-\frac{1}{20}.$ From the results of §1 we know that this holds along the curve

[TABLE]

or (approximately)

[TABLE]

This constant eigenvalue curve is tangent to the envelope at the simultaneous solution of

[TABLE]

(Note: we have divided each equation through by a constant.) The biologically important root of the above pair of equations is the one in the first quadrant, $(\rho_{2},\rho_{1})=(1.22,2.23):$ negative values of $\rho$ would correspond to an excitatory Purkinje-to-vestibular connection, which is not known to occur.

The basic picture of integrator operation is as follows: let us suppose that $\rho_{2}$ is allowed to vary but that $\rho_{1}$ is given by Eqn. (54), so that $\lambda=-\frac{1}{20}$ is always an eigenvalue. As $\rho_{2}$ is increased from zero the dominant eigenvalue is fixed at $\lambda=-\frac{1}{20}$ and the (in this case real) subdominant eigenvalue increases. Because we are nearing the envelope curve (and thus a degenerate double eigenvalue $\lambda_{2}=\lambda_{1}=-\frac{1}{20}$), we expect the gain to increase. At $\rho_{2}=1.22$, where the constant eigenvalue curve is tangent to the envelope, there is a collision of eigenvalues and the dominant eigenvalue is degenerate. As $\rho_{2}$ is further increased the formerly subdominant eigenvalue is now dominant - it is real and larger than $-\frac{1}{20}$. Furthermore, we can compute the gain associated to the dominant mode along the $\lambda=-\frac{1}{20}$ curve and find

[TABLE]

Note that, as expected, the denominator of the gain diverges at $\rho_{2}=1.22,$ where the constant eigenvalue curve is tangent to the envelope and the $\lambda=-\frac{1}{20}$ eigenvalue is degenerate.

In Figure 3A we show the response of the network to an impulsive forcing of the form $f(t)=\delta(t)(1,1,1,1,1,1,0,0)^{t}$ as $\rho_{2}$, $\rho_{1}$ are varied to increase gain. (Note that the impulsive forcing is equivalent to free decay with a corresponding initial condition). Three responses are shown, with $\rho_{2}$ set to $0.65$,$0.955$, and $1.095$ respectively. For each value of $\rho_{2}$, $\rho_{1}$ is set so that the network lies on the constant time constant curve $\tau=20s$. The gains predicted based on a single dominant mode ($\gamma$ in Eqn. (55)) are $2.52,5.92$ and $12.88$. The corresponding measurements from the impulse response (given by dividing the maximum response by $\|\vec{b}\|^{2}$) yielded $2.54,5.87$ and $12.53$ respectively; so we see excellent agreement.

Figure 3B displays the interaction of the two dominant eigenvalues of the network in the vicinity of the current operating region. In order to increase gain, the network must climb up the $\lambda=-0.05$ curve in the vicinity of the double eigenvalue point. At the intersection of this curve with the envelope, the integrating eigenvalue exchanges dominance with another real eigenvalue, producing an unstable integrator. Note that the larger two gain cases straddle the point where the $\lambda=-\frac{1}{20}{\rm s}^{-1}$ curve crosses the envelope, indicating an eigenvalue bifurcation in a subdominant mode. In this case it is the mode(s) with the next largest real part. At about $\rho_{2}\approx 1.02$ the subdominant complex conjugate pair collides at the real axis, and for $\rho_{2}$ above this value the first three most dominant eigenvalues are all real. As $\rho_{2}$ increases the subdominant eigenvalue increases until $\rho_{2}\approx 1.22$ where there is an eigenvalue collision and exchange of dominance.

Figure 3B also uses letter labels to show the character of the dominant eigenvalue. Normal operation requires the network to remain in regions A or G. If the network is in error, it may wander into a region where the two eigenvalues are complex (F), or into the region where one or both eigenvalues are in the right half-plane (B,C, D or E). In neither case is normal operation possible.

Infantile nystagmus (IN) is a hereditary disorder characterized by involuntary, periodic eye movements. These movements (or waveforms) can be broadly classified into two forms, jerk and pendular. In jerk waveforms, the eye moves outward from the central position with increasing speed until interrupted by a sudden saccade; a pendular waveform resembles a sinusoidal oscillation. The presence of a jerk waveform suggests an unstable eigenvalue ($\lambda>0$); as noted in many integrator models, the need to maintain an eigenvalue very near zero implies that an unstable eigenvalue is a natural consequence of imprecision.

In our model, the ability to modulate gain is enhanced near the bifurcation curve; this suggests that the network is also near a point where the dominant eigenvalue is complex, which would generate the sinusoidal oscillations characteristic of pendular nystagmus. For example, consider the network specified in Eqn. (57). This is very similar to Eqn. (51); the only change is that the vestibular-to-Purkinje input from a handful of neurons has been altered. Here, as the cerebellum attempts to increase gain by adjusting $\rho_{1}$ and $\rho_{2}$, it enters an oscillating regime.

[TABLE]

**Figure 4 here

**

We see the impulse response of the network in Figure 4A. As in Fig. 3A, $\rho_{2}$ and $\rho_{1}$ are varied so as to remain along the $\lambda=-\frac{1}{20}{\rm s}^{-1}$ curve. As $\rho_{2}$ (and the gain) increase, the network enters a regime where an oscillation is superimposed on normal integration. Figure 4B illustrates why this behavior occurs. As we follow the $\lambda=-\frac{1}{20}{\rm s}^{-1}$ curve from left to right through Region A, gain will increase. However, the eigenvalue curve has an intersection with the Hopf curve far to the left of its intersection with the envelope. At this point the integrating eigenvalue exchanges dominance with a pair of imaginary eigenvalues. Beyond this point, the response of the network contains both the normal integrating mode and a superimposed oscillation (Region F).

3 Example: a continuum model of the oculomotor integrator

We next consider a continuum model that can be derived from the model of Anastasio and Gad by replacing the (relatively numerous) vestibular neurons with a continuum vestibular “line” denoted by $\psi(x)$, while the relatively few Purkinje cells remain discrete.

Suppose the network has $N$ vestibular cells on each side of the integrator, arranged in a row of length $L$. Recall that the vestibular interaction matrix ${\bf T}$ has nearest-neighbor structure; each row of ${\bf T}$ takes the form $\beta v_{j-1}+(-1+\beta)v_{j}+\beta v_{j+1}$ (Eqn. (39)). This will converge, in the continuum limit ($N\rightarrow\infty$), to a second derivative, because:

[TABLE]

The sum over vestibular cells in the equation for the Purkinje cells (two final rows of Eqn. (38)) can similarly be replaced by an integral:

[TABLE]

where $\psi(j\Delta x)=v_{j}$ and $\phi_{1}(j\Delta x)=(\vec{w}_{1}^{t})_{j}$.

This leads to the following integro-differential eigenvalue problem on $L_{2}[0,L]\times\mathbb{R}^{2}$

[TABLE]

where we have replaced $L/N$ with $\Delta x$. 222We have also scaled the entire problem by the common factor of $\alpha$. The delta function coupling reflects the sparseness of the Purkinje-to-vestibular connections, with $x_{1,2}$ denoting the points on the vestibular line innervated by the Purkinje cells, and the functions $\phi_{1,2}(x)$ represent the density of vestibular to Purkinje connections.

Eliminating $P_{1}$, $P_{2}$ from Eqn. (60) algebraically, we arrive at the single equation

[TABLE]

This maps to our original problem, Eqn. (1), as follows:

[TABLE]

This model can be solved in much the same way as the discrete model analyzed in §2. To illustrate we take $L=1$ and the vestibular-to-Purkinje connections to be $\phi_{1}(x)=\phi_{2}(x)=1.$ Note that when $\rho_{1}=0,\rho_{2}=0$, the vestibular neurons decouple from the Purkinje cells, and the eigenvectors $\psi$ are given by trigonometric functions with the appropriate boundary conditions:

[TABLE]

after which $P_{1}$, $P_{2}$ are obtained by using Eqns. (61):

[TABLE]

The corresponding eigenvalues are given by $\lambda_{n}=-1+3\beta-\beta\frac{n^{2}\pi^{2}}{N^{2}}$, together with $\lambda=-1$, an eigenvalue of multiplicity two corresponding to the Purkinje cells. Note that the even modes ($n=2k$) do not actually excite a Purkinje cell response (i.e. $P_{1}=P_{2}=0$), and thus the even modes are eigenfunctions of this problem for all values of $\rho_{1},\rho_{2}$: these modes do not change under perturbation by the Purkinje cells.

We now apply the techniques developed earlier in the paper to find the lines of constant eigenvalue. We will not reproduce the entire calculation, but merely note a few salient points. The first step is to act on Equation (60) with the appropriate resolvent (inverse) operator (here, $R_{\lambda}=(\beta(\Delta x)^{2}\partial_{xx}+3\beta-(1+\lambda))^{-1}$). The terms $R_{\lambda}\delta(x-x_{i})$ are simply the Green’s function for the operator $\beta(\Delta x)^{2}\partial_{xx}+3\beta-(1+\lambda)$ acting on $L_{2}[0,L]$ with Dirichlet boundary conditions, which can be easily calculated – see (for example) the text of Keener[34]; we also supply some details in the Appendix.

The main conclusion, is that the perturbed problem will have piecewise trigonometric eigenfunctions of the form $\sin(\omega x)$, etc. (with appropriate continuity conditions; see Eqn. (Derivation of Equation (8))), where

[TABLE]

We now fix the remaining constants in the model: $x_{1}=\frac{1}{3}$ and $x_{2}=\frac{1}{2}$, and perform the calculation thus described. We find that the resulting envelope curves have asymptotes when $\omega$ is a multiple of $4\pi$ and $6\pi$; this is a consequence of the location of the Purkinje cell innervation (i.e. $x_{1}$ and $x_{2}$).

One issue that becomes apparent, is that as the network becomes large we are unable to find an integrating eigenvalue using Eqn (66) in this way; observe that $\lambda$ is bounded above by $\lambda=-1+3\beta-\beta(\Delta x)^{2}\omega^{2}<=-1+3\beta\approx-\frac{1}{40}+\frac{39}{120}\left(\frac{\pi}{N+1}\right)^{2}$, which $\rightarrow-\frac{1}{40}$ as $N\rightarrow\infty$. However, the desired eigenvalue for a 20 s time constant (recall that we have scaled out the 5 ms single neuron time constant, $\alpha=200$) is $\lambda=-0.05/200$, which is much closer to zero than $-1/40$; as $N$ increases, the integrating eigenvalue will become out of reach. One way to view this observation is that as ratio of vestibular to cerebellar (Purkinje) cells increases, sparse cerebellar innervation becomes less able to alter the time constant of the vestibular network.

The alternative, is to consider piecewise exponential eigenfunctions; i.e.

[TABLE]

where now

[TABLE]

Although such functions cannot match the boundary conditions of the unperturbed problem, the discontinuities in the perturbed problem bring them back into consideration. We can plot the corresponding bifurcation curve; it begins near $(\rho_{2},\rho_{1})\approx(-0.09,0.09)$ (for $N=12$) and increases without bound into the second quadrant. In Fig. 5, we plot the phase plane for several values of $N$; $N=12$, 24, 50 and 100. In each panel, we show several pieces of the bifurcation curve in different colors: $(0,4\pi)$ (purple), $(4\pi,6\pi)$ (blue), $(6\pi,8\pi)$ (cyan), $(8\pi,12\pi)$ (green), and the curve for exponential eigenfunctions (yellow). We note that the structure is stable, for increasing $N$; however, the corresponding values of $\lambda$ are “compressed” (see Eqn. (66)). This is to be expected; in the finite-dimensional problem, the eigenvalues cluster together as $N$ increases; therefore we expect increasing density of constant eigenvalue curves as $N$ increases.

We now return to the question that originally motivated our analysis in §2; can this network behave as an integrator; i.e. can it achieve the correct time constant and gain, by adjusting its Purkinje-to-vestibular weights $\rho_{1}$, $\rho_{2}$?

Recall that we characterized the performance of an integrator by two quantities; the dominant eigenvalue $\lambda_{1}$ and the gain $\gamma$. In §2 we established that in order to achieve high gain, the parameters $\rho_{1},\rho_{2}$ must be set near the bifurcation point for the corresponding eigenvalue $\lambda_{1}$; this is where $\lambda_{1}$ has multiplicity two, and therefore (if degenerate) where the angle between left and right eigenvectors is near zero. Finally, recall that connections from the cerebellum are always inhibitory; by the sign convention in Eqn. (62), this means that both $\rho_{1},\rho_{2}$ must be positive. Therefore, a necessary condition for the continuum model to act as an integrator, is for the bifurcation point corresponding to $\lambda_{1}$ to lie in the first quadrant.

This clearly does not occur in Fig. 5, where in every panel the yellow curve lies entirely in the second quadrant $\rho_{2}<0$, $\rho_{1}>0$. Therefore this network — where $\phi_{1}=\phi_{2}=1$, $x_{1,2}=\{1/3,1/2\}$, and plasticity comes from manipulating the strength of the Purkinje-to-vestibular connections $\rho_{1}$, $\rho_{2}$ — cannot regulate both time constant and gain.

We now take one further step in generality, by asking how the network behaves as the Purkinje-to-vestibular connection locations — $x_{1}$ and $x_{2}$ — change. This would reflect a different source of network plasticity. We will show that for the vestibular-to-Purkinje projection patterns chosen here, the $\rho_{1},\rho_{2}$ coordinates of the bifurcation point will always be of opposite sign, and therefore will never occur in the first quadrant.

Lemma 5**.**

Consider the continuum integrator model defined in Eq. (62) with $L=1.$ Suppose $\phi_{1}=\phi_{2}=1$; then for sufficiently large $N$, the bifurcation point corresponding to the integrating eigenvalue $\lambda_{1}$ will not be found in the first quadrant ($\rho_{1},\rho_{2}>0$), for any $0<x_{1},x_{2}<1$.

Proof.

We assume that $N$ is large enough that the eigenvalue of interest requires an exponential eigenfunction; i.e. $\lambda_{1}>-1+3\beta-\beta(\Delta x)^{2}\omega^{2}$. We begin by directly computing the required terms in the characteristic polynomial:

[TABLE]

(refer to Eq. (8) for definitions and the Appendix for a more detailed calculation in terms of resolvent operators). Note that

[TABLE]

for the same function $f(\omega,x)$.

Then by Eqns. (13, 14) the coordinates of the bifurcation point corresponding to any desired $\omega=\omega_{0}$ will be:

[TABLE]

consequently,

[TABLE]

We will show that $\frac{\partial f}{\partial\omega}(\omega,x)>0$ for $\omega>0$ and $0<x<1$; therefore, the ratio of $\rho_{1}(\omega_{0})$ and $\rho_{2}(\omega_{0})$ must be negative.

Differentiating the function $f(\omega,x)$ given in Eqn. (72), we find that:

[TABLE]

It may not be obvious what the sign of this function is (the first line is positive and the second negative), but (re-)defining $F_{\omega}(x)\equiv\frac{\partial f}{\partial\omega}(\omega,x)$, we will show that $F_{\omega}(x)\geq 0$ for any $\omega>0$ and $x\in[0,1]$. We do this by confirming that $F_{\omega}(x)$ is concave in $(0,1)$, and that $F_{\omega}(0)=F_{\omega}(1)=0$. Therefore by the minimum principle for superharmonic functions, $F_{\omega}(x)>0$ on the interior $x\in(0,1)$.

To confirm concavity, we check that the second derivative in $x$ is negative:

[TABLE]

Checking that $F_{\omega}(0)=F_{\omega}(1)=0$ is a simple matter of substituting $x=0,1$ into Eqn. (73).

∎

In conclusion, we cannot make this network act as an integrator, by adjusting its Purkinje-to-vestibular connections. Instead, we would have to adjust the underlying vestibular-to-Purkinje projection pattern $\phi_{1}$.

**Figure 5 here

**

4 Example: The Rubinstein-Sternberg model

In this section we give a third, quite different, application for the technique based on low-rank perturbations. Rubinstein and Sternberg[2] introduced a nonlocal model for phase separation of the form

[TABLE]

We will compute the stability of a standing front type solution of this model, where the stability operator takes the form of a rank one perturbation to a standard Sturm-Liouville operator. This model has been analyzed by a number of authors, most notably Freitas[3, 35, 36, 10], and the closely related problem of the coarsening rate for the Cahn-Hilliard equation has been analyzed in the classic paper of Bates and Fife[37]. The main goal of this example is to illustrate the utility of treating the problem using the rank-one perturbation formula: while the final result appears to be new several of the intermediary results have analogues in the work of Freitas, and we will point these out where germane.

We begin with the equation

[TABLE]

This equation always admits a constant solution and for sufficiently large widths it admits front type solutions. At $L=\frac{n\pi\sqrt{f^{\prime}(0)}}{2}$ a bifurcation occurs giving rise to a solution containing $n$ fronts. In the absence of the non-local term these front solutions are (for $L$ finite) always unstable. The most unstable mode has a non-vanishing mean, so the instability is connected with non-conservation of mass. Rubinstein and Sternberg introduced the non-local term as a Lagrange multiplier to enforce mass conservation and remove this instability mechanism.

It is easy to see that, if $u$ represents a stationary solution to Eqn. (75) then the linearized evolution equation is given by

[TABLE]

Therefore the associated eigenvalue problem takes the form of a rank-one perturbation of a self-adjoint operator. To explicitly relate Eqn. (76) to Eqn. (1) we can define

[TABLE]

The problem with a bistable cubic nonlinearity, $f(u)=u-u^{3}$,

[TABLE]

is the simplest and most natural from a physical perspective, and can be analyzed rather explicitly. Assuming $L>\frac{\pi}{2}$ there is a front solution which can be expressed in terms of elliptic functions333For fixed period $L$ there is actually a one-parameter family of front solutions. The one given here has zero net mass and is the simplest. The rest will be discussed later.. After a simple rescaling ($u=(1+k^{2})^{-1}v$, $x=(1+k^{2})^{-1/2}y$, and $t=\frac{\sqrt{2}k}{\sqrt{1+k^{2}}}s$), the equation can be written in the form

[TABLE]

Here the quantity $k\in(0,1)$ denotes the elliptic modulus and $K(k)$ denotes the complete elliptic integral of the first kind

[TABLE]

The elliptic modulus $k$ is determined from $L$ by the relation

[TABLE]

The unperturbed problem now becomes

[TABLE]

a well-known identity for elliptic functions states that a stationary solution to this equation is given by

[TABLE]

where ${\rm sn}$ is the Jacobi elliptic sinus function.

Note that as the function ${\rm sn}(x,k)$ is odd the nonlocal term $\frac{1}{2K}\int_{-K(k)}^{K(k)}((1+k^{2})u-2k^{2}u^{3})\,dx$ vanishes; thus, this function solves both the perturbed and unperturbed problem (i.e. both Eqn. (78) and Eqn. (79)). However, the non-local term changes the stability problem; we will find that the stability properties of the two problems are not the same.

Linearizing around the elliptic function solution gives the following non-local evolution equation

[TABLE]

where we have identified $\widetilde{\bf H}$ as a perturbation from a “simpler” operator, ${\bf H}$. The unperturbed operator ${\bf H}$ is a two-gap Lamé operator, for which the spectral problem can be solved exactly [38, 39]. When subject to periodic boundary conditions on $[-2K(k),2K(k)]$ the largest five eigenvalues are simple and are given by

[TABLE]

with $a(k)=\sqrt{1-k^{2}+k^{4}}$. Here the superscript indicates whether the function satisfies a Neumann or a Dirichlet condition at $\pm K(k)$. Thus the unperturbed operator subject to Neumann boundary conditions has one positive eigenvalue $\lambda_{0}=-(1+k^{2}-2a(k))$ and the remainder of the eigenvalues on the negative real line.

To summarize, the stability problem (Eqn. (81), plus boundary conditions $v_{x}(\pm K(k))=0$) takes the form of a rank-one perturbation of a self-adjoint problem:

[TABLE]

with $G=-\frac{1}{2K(k)}{\bf 1}$ and $g=(1+k^{2}-6k^{2}{\rm sn}^{2}(x,k))$. While neither $G$ nor $g$ is an eigenvector of ${\bf H}$, both $G$ and $g$ lie in the span of the zeroth and second eigenfunctions $\phi_{0/2}^{(N)}(x)=k^{2}{\rm sn}^{2}(x,k)-\frac{1}{3}(1+k^{2}\mp a(k))$. Since the eigenfunctions of a self-adjoint operator are orthogonal this implies that the rank-one piece, $G\langle g,v\rangle$ vanishes on the span of all the remaining eigenfunctions. This implies that the perturbed operator decomposes as a direct sum of two operators, one that is self-adjoint and negative definite and one that is rank-two:

[TABLE]

We emphasize that the perturbation term here is very special, in that $G$ and $g$ are actually given by linear combinations of just two of the eigenfunctions of the unperturbed operator. In the generic case (i.e. for a different nonlinearity $f(u)$) one expects $G$ and $g$ to have non-trivial projections onto all eigenfunctions.

From the exact eigenvalues in Eqn. (82) the second term satisfies the coercivity estimate

[TABLE]

and the entire stability problem is reduced to understanding the two by two matrix eigenvalue problem defined by $\widetilde{\bf H}|_{{\rm span}(\phi_{0},\phi_{2})}.$ Furthermore, since the range of ${\widetilde{\bf H}}$ consists of mean zero functions it follows from the Fredholm alternative that $\widetilde{\bf H}|_{{\rm span}(\phi_{0},\phi_{2})}$ must have a zero eigenvalue.

This zero eigenvalue is connected with mass conservation. There is a one-parameter family of solutions to this equation of fixed spatial period $L$, which can be thought of as being related to the total mass of the stationary solution. Here we have only considered the simplest solutions – those with zero net mass – but there are analogous expressions in terms of elliptic functions for the general solution. Since we have a one parameter family of solutions, by Noethers theorem there must be an element in the kernel of linearized operator corresponding to the generator of this family.

It is straightforward to compute the restriction of the linearized operator $\widetilde{\bf H}|_{{\rm span}(\phi_{0},\phi_{2})}$ in the (independent but non-orthogonal) basis $\{1,{\rm sn}^{2}(x,k)\}$ in terms of complete elliptic integrals of the first and the second kind. By using the following identities:

[TABLE]

we derive the following expression for $\widetilde{\bf H}|_{{\rm span}{\phi_{0},\phi_{2}}}$, the restriction of the operator to the span of the zeroth and second eigenfunctions, in the basis $\{1,{\rm sn}^{2}(x,k)\}:$

[TABLE]

The eigenvalues of this restriction are given by

[TABLE]

It is easy to check directly from the definition of the elliptic integrals that the quantity $(3-3k^{2})K(k)-6E(k)$ is strictly negative for $k\in[0,1)$. This follows, for instance, from the Taylor series representations

[TABLE]

from which it is easy to see that the quantity $(1-k^{2})K(k)-2E(k)$ is even with only negative terms in its Taylor series, implying that it is strictly negative. Since the equation conserves mass it makes sense to consider only perturbations with zero net mass (those orthogonal to the kernel). In this case we have all eigenvalues strictly in the left half-plane and the stationary solution is non-linearly stable.

The problem for a general nonlinearity is slightly more complicated, since we do not have explicit formulae for the eigenvalues and eigenfunctions, but it can essentially be completely solved. Assuming that the equation

[TABLE]

has a stationary solution $u(x)$ we can compute the linearized operator as

[TABLE]

where ${\bf H}$ is the unperturbed operator: ${\bf H}v=v_{xx}+f^{\prime}(u)v$. It is convenient to introduce a coupling constant $\rho$ controlling the strength of the perturbation, although we are mainly interested in the special case $\rho=1$, and we thus consider a one-parameter family of eigenvalue problems

[TABLE]

Applying the Aronszjan-Krein formula gives the following eigenvalue condition:

[TABLE]

Thus the spectrum again decomposes into two pieces. Any eigenvectors of ${\bf H}$ which happen to be orthogonal to $f^{\prime}(u)$ ($\langle f^{\prime}(u),v\rangle=0$) remain eigenvectors of the perturbed problem. Eigenvectors which are not orthogonal to $f^{\prime}(u)$ must satisfy the Aronszajn-Krein eigenvalue condition

[TABLE]

(This condition was also identified by Freitas [3], who referred to the eigenvalues which satisfy Eqn. (86) as “moving eigenvalues”).

Now note that one has the identity ${\bf H}{\bf 1}=f^{\prime}(u)$ and thus we can write

[TABLE]

where we use the identity ${\bf 1}=\sum\langle{\bf 1},v_{i}\rangle v_{i}$ and therefore that $({\bf H}-\lambda{\bf I})^{-1}{\bf 1}=\sum\frac{\langle{\bf 1},v_{i}\rangle}{\lambda_{i}-\lambda}v_{i}$; here, $\lambda_{i}$ and $v_{i}$ are the eigenvalues and (orthogonal) eigenvectors of the unperturbed operator ${\bf H}$. Notice that in the last expression $h(\lambda)$ is a Herglotz function444Sometimes called a Nevanlinna or Nevanlinna-Pick function. — an analytic function that is real on the real axis and maps the open upper half-plane to itself — for $\rho\in[0,1].$ It is well known that the zeroes and poles of a Herglotz function are real, implying that the eigenvalues of the Rubinstein-Sternberg model for $\rho\in[0,1]$ are real (see Simon[40] page 920 for a more detailed discussion of Herglotz functions). This observation is analogous to Lemma 5.2 in the paper of Freitas[36], where reality of the eigenvalues for $\rho\in[0,1]$ is established using a combination of identities derived from the original equation. In later work Freitas [10] considers a general rank one perturbation of a self-adjoint operator: $\widetilde{\bf H}={\bf H}+|a\rangle\langle b|$, and associates to each eigenvalue $\lambda_{i}$ a signature given by ${\rm sgn}(\langle b,v_{i}\rangle\langle v_{i},a\rangle)$, where $v_{i}$ is the eigenfunction of the unperturbed problem — see, in particular, section 3 of [10]. The condition that all of these signatures are the same is equivalent to requiring that the Aronszajn-Krein function is Herglotz.

This calculation shows that, despite the fact that the linearized operator is not self-adjoint, the spectrum is purely real for $\rho\in[0,1]$. Now let us assume that we understand the spectrum of ${\bf H}_{0}={\bf H}$, the unperturbed operator, in particular that we know $n_{+}({\bf H})$, the number of positive eigenvalues of the unperturbed operator. Let us consider doing a homotopy in the parameter $\rho$. Since we have shown that the eigenvalues are real the only way that an eigenvalue can move from the right half-plane to the left half-plane (or vice-versa) is by passing through the origin. We can detect when this occurs by taking $\lambda=0$ in equation (86), which becomes

[TABLE]

Here we have used the fact that ${\bf H}{\bf 1}=f^{\prime}(u)$ so that ${\bf 1}={\bf H}^{-1}f^{\prime}(u)$. (We are also assuming that ${\bf H}$ is invertible. Minor changes are required in the case that ${\bf H}$ has a kernel; see Remark 1.) This calculation shows that the unique value of $\rho$ for which ${\bf H}_{\rho}$ has a zero eigenvalue is $\rho=1$.

We are now in a position to count the number of positive eigenvalues of $\widetilde{\bf H}$ using a continuation argument. The operator $\widetilde{\bf H}$ is bounded above and, by standard arguments, has a finite number of positive eigenvalues. We assume that the number of positive eigenvalues of the unperturbed operator is given by $n_{+}({\bf H})=k$: at $\rho=0$ there are $k$ positive eigenvalues and the remaining eigenvalues are negative. For $\rho\in(0,1)$ the kernel of ${\bf H}_{\rho}$ is empty, so no eigenvalues cross from the left half-line to the right. At $\rho=1$ there is an eigenvalue at $\lambda=0$, which either came from the left half-line or the right half-line. We can determine which by computing $\frac{d\lambda}{d\rho}$ and evaluating at $\lambda=0$. Doing so we find that

[TABLE]

If $\langle{\bf 1},{\bf H}^{-1}{\bf 1}\rangle>0$ the eigenvalue is moving from the positive half-line to the negative, while if $\langle{\bf 1},{\bf H}^{-1}{\bf 1}\rangle<0$ it is moving from the negative to the positive half-line.

Finally, the assumption that the unperturbed operator ${\bf H}$ is invertible implies that the perturbed operator $\widetilde{\bf H}$ has at most a one dimensional kernel, since if $\widetilde{\bf H}$ had a higher dimensional kernel one of the eigenfunctions could be chosen to be orthogonal to $f^{\prime}(u)$ and would thus lie in the kernel of the unperturbed operator ${\bf H}$. This completes the proof of the following theorem:

Theorem 1**.**

Suppose that the unperturbed operator ${\bf H}=\partial_{xx}+f^{\prime}(u)$ is non-singular and has $n_{+}({\bf H})=k$ positive eigenvalues. The perturbed operator $\widetilde{\bf H}$ has a simple kernel and $n_{+}({\bf{\widetilde{H}}})$, the number of positive eigenvalues of the linearized operator, is given by

[TABLE]

Thus a necessary and sufficient condition for spectral stability is that $n_{+}({\bf H})=0$ (from which it follows that $\langle{\bf 1},{\bf H}^{-1}{\bf 1}\rangle<0$) or $n_{+}({\bf H})=1$ and $\langle{\bf 1},{\bf H}^{-1}{\bf 1}\rangle>0.$

Remark 1**.**

The case where the unperturbed operator ${\bf H}$ has a kernel is somewhat more involved but can be addressed similarly. There one must do another perturbation calculation near $\rho=0$ to understand how the zero eigenvalue(s) move with $\rho$ in order to compute $n_{+}({\bf H}_{\rho})$ for $\rho$ small but non-zero. From there the calculation is the same: the number of positive eigenvalues can stay the same or decrease by one, and this is determined by the sign of $\langle{\bf 1},{\bf H}^{-1}{\bf 1}\rangle$. Here ${\bf H}$ is singular but ${\bf 1}$ is in the range of ${\bf H}$, so “${\bf H}^{-1}{\bf 1}$” may be interpreted in the sense of the Moore-Penrose pseudo-inverse.

As we remarked earlier there is actually a one-parameter family of stationary solutions to the nonlocal equation Eqn. (75). We now give an alternative characterization of the stability criterion by expressing it in terms of the dependence of the integrated reaction rate on the mass of the solution.

Corollary 1**.**

Two alternative characterizations of the stability criterion in Theorem (1) are as follows:

Let $M$ denote the total mass of the solution

[TABLE]

and $R$ denote the total reaction rate

[TABLE]

Assume that the family of stationary solutions can locally be parameterized by the total mass $M$, and that the number of positive eigenvalues of the unperturbed operator is given by $n_{+}({\bf H})=k$. Then the dimension of the unstable manifold is given by

[TABLE]

In particular a necessary condition for stability is that the total reaction rate must be an increasing function of the total mass. 2. 2.

Suppose that the stationary solution is given by the quadrature

[TABLE]

for appropriate constants $E,\kappa$. Let $\mu_{-}$ and $\mu_{+}$ be two turning points for the quadrature: i.e. $\mu_{\pm}$ are simple roots of $2E+2\kappa u-2F(u)=0$ such that $2E+2\kappa u-2F(u)>0$ for $u\in(\mu_{-},\mu_{+})$, with $F$ the antiderivative of the reaction rate, $F^{\prime}(u)=f(u)$. Define the period type integrals

[TABLE]

where $\Gamma$ is a simple closed contour containing the branch cut along the real axis from $\mu_{-}$ to $\mu_{+}$. Then the dimension of the unstable manifold is given by

[TABLE]

where $\tau$ is defined by

[TABLE]

and represents the rate of change of $M$ divided by the rate of change of $R$ along the one-parameter family of stationary solutions.

Proof.

The proof consists of computing the family of stationary solutions and observing that translation along the family of stationary solutions generates the appropriate element of the range of ${\bf H}$ needed to compute $\langle{\bf 1},{\bf H}^{-1}{\bf 1}\rangle$.

First, we show that we can find a one-parameter family of stationary solutions through quadrature: see, for instance, the text of Landau and Lifshitz[41]. A stationary solution $u$ must satisfy

[TABLE]

for some suitable $\kappa$; multiplying by $u_{x}$ we find that

[TABLE]

for some integration constant $E$, and $F^{\prime}(u)=f(u)$. Solving for $u_{x}$ yields

[TABLE]

The Neumann boundary condition $\frac{du}{dx}=0$ is satisfied at a turning point of the function $\sqrt{2E+2\kappa u-2F(u)}$; therefore we integrate along the real axis between two points $\mu_{-}=u(-L)$, $\mu_{+}=u(L)$, at which $\sqrt{2E+2\kappa u-2F(u)}=0$:

[TABLE]

or

[TABLE]

where $\Gamma$ is a simple closed contour that loops around the segment between $\mu_{-}$ and $\mu_{+}$. The factor of $\frac{1}{2}$ in the second equation arises in the usual way due to the fact that the contributions from the top and the bottom of the square root branch cut add. Condition (89) defines a curve in the $(E,\kappa)$ plane along which we have a stationary solution, defined locally by the vector field

[TABLE]

If we choose to parameterize the curve $(E(s),\kappa(s))$ in the $(E,\kappa)$ plane by arc length $s$ then we can take

[TABLE]

Note that along this curve we have the identities

[TABLE]

Also note that one has the identity $\kappa P(E,\kappa)=R(E,\kappa)$, so one could in principle eliminate one of these quantities although we have chosen not to do so here.

Having found a family of stationary solutions parameterized by the arc length, $u(x;s)$, we now proceed to compute the quantity $\langle{\bf 1},{\bf H}^{-1}{\bf 1}\rangle$ in terms of $M$ and $R$. First, take the equation for the stationary solution

[TABLE]

and differentiate with respect to the arc length parameter $s$, giving

[TABLE]

which we recognize as (recalling that ${\bf H}v=v_{xx}+f^{\prime}(u)v$)

[TABLE]

The right-hand side is a constant, thus we have

[TABLE]

Integrating this identity gives

[TABLE]

or

[TABLE]

Since $L>0$ is positive, the quantity $\langle{\bf 1},{\bf H}^{-1}{\bf 1}\rangle$ is positive if $R$ increases with increasing $M$ and is negative if $R$ decreases with increasing $M$. Applying the chain rule and Equations (90) and (91) we find that

[TABLE]

∎

While we are not aware of previous results of this type for dissipative equations there are a number of results for the stability of nonlinear dispersive waves that relate the index of some linearized operator to the sign of the derivative of some conserved quantity with respect to a parameter, analogous to the result presented here. The classical Vakhitov-Kolokolov criteria, relating the stability of solitary wave solutions $\psi(x,t)=e^{i\omega t}\phi(x;\omega)$ to the nonlinear Schrödinger equation

[TABLE]

to the sign of $\frac{d}{d\omega}\int|\phi|^{2}(x;\omega)dx$ falls into this category [42, 43, 44, 45]. The problem of the stability of periodic solutions to nonlinear dispersive waves is perhaps even more similar to the present case: in this situation the index is determined by the signs of certain determinants of derivatives of period integrals (see [46, 47, 48] for examples).

It is worth discussing the relationship of this result to the results of Freitas[36], who relates the dimension of the unstable manifold to the lap number of the underlying solution. Specifically Freitas shows (see Theorems 5.1 and 5.13 of [36]) that the dimension of the unstable manifold is related to $m$, the number of extrema of the underlying solution $u$ in $(-L,L)$ via

[TABLE]

Thus the dimension of the unstable manifold can differ from the number of extrema by up to two. To summarize the reasoning, the unstable manifold of the unperturbed operator has dimension equal to the lap number, $l(w)=m+1$; the corresponding dimension for the perturbed operator can then go up (to $m+2$) or down (to $m$) by at most one.

Here, we determine this direction (whether up to $m+2$, or down to $m$), by computing the index of the unperturbed operator; in terms of $m$, our result states that

[TABLE]

where the sign of $\langle{\bf 1},{\bf H}^{-1}{\bf 1}\rangle$ can be replaced by the sign of $\frac{dR}{dM}$ or any of the other equivalent forms discussed previously.

5 Discussion

In conclusion, we have presented a general method of analyzing low-rank perturbations of self-adjoint operators. We show how to use a simple idea of classical differential geometry (the envelope of a family of curves) to completely analyze the spectrum. When the rank of the perturbation is two, this allows us to view the system in a geometric way through a phase diagram in the perturbation strengths $(\rho_{1},\rho_{2})$. By locating constant eigenvalue and eigenvalue coincidence curves (both computable through simple formulas), we can determine where the perturbed operator is stable, and where double real eigenvalues bifurcate into complex pairs. This latter situation (bifurcation into a complex pair) coincides with a poorly conditioned eigenvalue, which in turn signals that small changes in the perturbation parameter will yield large changes in the operator behavior.

We used these techniques to analyze three problems of this form; a model of the oculomotor integrator due to Anastasio and Gad[1], a continuum version of the oculomotor integrator model, and a nonlocal model of phase separation due to Rubinstein and Sternberg[2]. In the first two problems, the physical interpretation of our model (a neural network that must maintain a steady eye position in the absence of input) required that we identify (a) where the perturbed system had a specific eigenvalue and (b) where this particular eigenvalue would be poorly conditioned. Our results in §1 show that both (a) and (b) can only occur in proximity to a specific point on the $(\rho_{1},\rho_{2})$ plane, which was then easy to visualize. In §2 and §3, some portions of the model are not completely specified by biology (such as the vestibular-to-Purkinje connections), but must be chosen arbitrarily (even randomly). In this paper, we analyzed a few carefully chosen examples. But, the geometric method we describe here also gives us rapid way to survey a large family of such models; using such a survey to draw conclusions about the vestibular-to-cerebellar pathway is an area for future work.

The problem analyzed in §4 involves a rank one (rather than rank two) perturbation, and so a phase plane approach is not applicable. Instead, we systematically exploit the rank one nature of the perturbation to characterize stability of stationary solutions in terms of the unperturbed operator. We further show how to construct a one-parameter family of stationary solutions, and relate the stability condition to the relative change in two integrated quantities (mass and reaction rate) as one travels along this family.

In analyzing these three problems, we have by no means exhausted the possible applications. For example, the eigenvalue problem in Eqn. (62) is very similar to the stability problem for spike solutions to activator-inhibitor models in the limit of slow activator diffusion[3, 4] (although the problem we study here differs because the eigenvalue enters in a non-linear way). Similar models of reaction-diffusion equations with non-local interactions have arisen in a number of other contexts including population dynamics[49], runaway ohmic heating[5, 6, 7], and microwave heating [8]. Therefore, we anticipate that the techniques presented here should be applicable to understanding these problems.

Competing Interests

The authors have no competing interests.

Data Accessibility

Not applicable: the paper contains sufficient detail to reproduce the results.

Authors’ Contributions

TJA, AKB, and JCB designed the project. TJA originated the integrator model described in §2; AKB and JCB originated the continuum integrator model (§3) and analyzed the Rubinstein-Sternberg model described in §4. AKB and JCB created figures; TJA, AKB and JCB wrote the paper.

Acknowledgements

We would like to thank two anonymous reviewers whose comments have helped us to improve the manuscript.

Funding

JCB and AKB received support from grant NSF-DMS0354462, and JCB from NSF-DMS0807584.

Appendix

Derivation of Equations (3) and (7)

Here we include details for the derivation of Eqn. (3), for matrices. We will make repeated use of the matrix determinant lemma:

[TABLE]

and the Sherman-Morrison formula:

[TABLE]

We first apply Eqn. (96) with ${\bf A}=\,\mbox{\bf M}-\lambda{\bf I}+\rho_{1}\vec{f}_{1}\vec{g}_{1}^{t}$, $\vec{u}=\rho_{2}\vec{f}_{2}$, and $\vec{v}=\vec{g}_{2}$, resulting in Eqn. (100). We then apply Sherman-Morrison to the inverse matrix, resulting in Eqn. (LABEL:eqn:applySherMorr). Finally we apply Eqn. (96) a second time, but with ${\bf A}=\,\mbox{\bf M}-\lambda{\bf I}$, $\vec{u}=\rho_{1}\vec{f}_{1}$, and $\vec{v}=\vec{g}_{1}$, and simplify to emphasize the polynomial form in $\rho_{1},\rho_{2}$:

[TABLE]

Using the cofactor formula $\operatorname{cof}^{t}{\bf A}=\det({\bf A}){\bf A}^{-1}$ to replace each instance of $\det({\bf M}-\lambda{\bf I})({\bf M}-\lambda{\bf I})^{-1}$ with $\operatorname{cof}^{t}({\bf M}-\lambda{\bf I})$ gives the formula in Eqn (3).

To derive Eqn. (7), use the spectral decomposition for self-adjoint operators:

[TABLE]

where $\lambda_{i}$ and $\vec{\phi}_{i}$ are the eigenvalue and eigenvectors of the unperturbed operator; substitute into Eqn. (Derivation of Equations (3) and (7)) and use the orthogonality of eigenfunctions ($\vec{\phi}_{i}^{t}\vec{\phi}_{j}=0$ for $i\not=j$).

Derivation of Equation (8)

We fill in a few more details for calculating eigenvalues of:

[TABLE]

where $\widetilde{\bf M}$ is an operator with compact resolvent.

First, act on this by the resolvent of the unperturbed operator, $R_{\lambda}=({\bf M}-\lambda{\bf I})^{-1}$:

[TABLE]

Now, act on the equation by both $\vec{g}_{1}$ and $\vec{g}_{2}$, to yield two consistency conditions:

[TABLE]

Substituting Eqn. (108) into Eqn. (106) and dividing by the now common factor $\langle\vec{g}_{1},\vec{w}\rangle$, will yield a single polynomial equation for $\rho_{1}$ and $\rho_{2}$:

[TABLE]

When $\vec{g}_{1}=\vec{g}_{2}$, the final $\rho_{1}\rho_{2}$ term is zero; this is what happens in the continuum example presented in §3.

We now explain how to apply this formalism to the eigenvalue problem in §3:

[TABLE]

Here, ${\bf M}=\beta\left(\Delta x\right)^{2}\psi_{xx}+(-1+3\beta)\psi$ and (as we expect negative eigenvalues) we will act on Eqn. (LABEL:reduce_appendix) with the resolvent operator $R_{\lambda}\equiv\left(\beta\left(\Delta x\right)^{2}\psi_{xx}+(-1+3\beta)\psi-\lambda\psi\right)^{-1}$. Here $\vec{f}_{1}\propto\delta(x-x_{1})$ and $\vec{f}_{2}\propto\delta(x-x_{2})$; therefore $R_{\lambda}\vec{f}_{1}$ is equivalent to solving the Green’s function problem

[TABLE]

Since $\psi$ solves the PDE $\psi_{xx}+\frac{-1-\lambda+3\beta}{\beta\Delta x^{2}}\psi=0$ for any $x\not=x_{1},x_{2}$, $\psi$ must have the following form:

[TABLE]

where $\omega^{2}=\frac{-1-\lambda+3\beta}{\beta(\Delta x)^{2}}$. By linearity of the resolvent, and the fact that $\vec{g}_{1}=\vec{g}_{2}$, we can solve them together: i.e. $\rho_{1}\langle\vec{g}_{1},R_{\lambda}\vec{f}_{1}\rangle+\rho_{2}\langle\vec{g}_{2},R_{\lambda}\vec{f}_{2}\rangle=\rho_{1}\langle\vec{g}_{1},R_{\lambda}\vec{f}_{1}\rangle+\rho_{2}\langle\vec{g}_{1},R_{\lambda}\vec{f}_{2}\rangle=\langle\vec{g}_{1},\rho_{1}R_{\lambda}\vec{f}_{1}+\rho_{2}R_{\lambda}\vec{f}_{2}\rangle=\langle\vec{g}_{1},R_{\lambda}(\rho_{1}\vec{f}_{1}+\rho_{2}\vec{f}_{2})\rangle$.

Therefore we solve the Green’s function problem only once; $\psi$ must satisfy the following conditions, which impose continuity of $\psi$ at $x_{1},x_{2}$ and the appropriate jump of the derivatives:

[TABLE]

This yields a system of 4 equations in the 4 unknowns $A$, $B$, $C$, and $D$. Similarly, we can express the action of $\vec{g}_{1}$ on any function of the form Eqn. (Derivation of Equation (8)) as a vector inner product: in this case,

[TABLE]

Thus our final polynomial is $0=1+\mathbf{g}^{T}\mathbf{w}$, where $\mathbf{w}$ is the vector of coefficients we obtained by solving Eqn. (113–116).

Figure captions