Dynamics, Stability and Bifurcation Phenomena in the Nonlocal Model of Cortical Activity

Dynamics, Stability and Bifurcation Phenomena in the Nonlocal Model of Cortical Activity

Konstantin Doubrovinski

U.U.D.M. Project Report 2005:8

Examensarbete i matematik, 20 poäng Handledare och examinator: Sten Kaijser

Juni 2005

Department of Mathematics

Uppsala University

Dynamics, Stability and Bifurcation Phenomena in the Nonlocal Model of Cortical Activity

1 Introduction

1.1 Background

Computational neuroscience is a novel rapidly developing area of physical research, aiming to uncover the mechanisms behind the functionality of the brain. In recent years numerous detailed models of the dynamics of individual neurons were put forward. An arbitrarily accurate description thereof is readily obtained by incorporating sufficient number of structural features [Dayan 2003, Koch 1989]. However, when trying to understand the behavior of the nervous system, the requirement of simulating huge neuronal populations constitutes the major obstacle. In order to remedy this problem the concept of neural field was introduced in the early seventies: instead of concentrating on the dynamics of individual neurons, one considers the cell population as an isotropic homogeneous medium. Two major models, utilizing this approach, have been put forward. The first one is Wilson-Cowan model [Wilson 1973]. It is based on few thoroughly justified physiological assumptions, comprehensively recapitulated in the original paper. The pattern-forming properties of the Wilson-Cowan model received substantial attention [Enculescu 2003, Tass 1997, Ermentrout 1998]. Those of the second, so called, Amari model [Amari 1977] remain poorly understood. Here we aim to continue the study of the properties of the solutions of the Amari equations.

Recent developments in the area of experimental neuroscience, in particular those, utilizing voltage sensitive dyes allow to asses the validity of the model predictions. A typical experimental setup is described in [Grinvald 1994]. Craniotomy was performed on anesthetized macaque monkeys. Exposed region of the cortex was stained with a fluorescent dye which changes the intensity of fluorescence in response to varying membrane potentials of the cortical neurons. (According to [Grinvald 1994] the signal from the dye effectively measures “fast changes in membrane potential”). Drifting black and white gratings were presented to an animal on a computer screen, simultaneously monitoring the cortical activity with a macroscope through an optical chamber. Such imaging procedures demonstrated the presence of reproducible spatially confined bump- like structures, emerging in response to spatially localized visual stimuli, as depicted in Figure 1. However, the interpretation of these experiments remains highly arguable and we believe that mathematical modeling is the appropriate tool for understanding the underlying dynamics.

Figure 1

Spatial distribution of cortical activity, revealed by real-time optical imaging of macaque monkey visual cortex. Adopted from [Grinvald 1994].

We start with a brief derivation of the Amari equations.

The quantity characterizing the state of a neuron at a certain time instant is the membrane potential which is due to the ion concentration difference between the interior and the exterior of the cell. When modeling a neuronal population consisting of m cells one aims to deduce the membrane potential ui(t) of each cell i at time instant t. It is physiologically plausible to assume that in the absence of external input (either sensory or stemming from the rest of the neuronal population) it is to approach some resting potential h:

( )

du

i =− i +

1.1 Being connected to the other members of the population each neuron can be assumed to interact with the surrounding cells according to:

( ) _∑ [ ]

=

+ +

−

= ^m

j

j ij

i ui t h f u

dt du

1

ω

1.2 This is nothing but the most conventional neural net, frequently encountered in the literature, termed Hopfield network [Haykin 1994]. ωij are the so called synaptic weights, f is the nonlinear output function. Assume that the neuronal population under consideration constitutes a flat, spatially extended sheet, consisting of a single cell layer.

Now every neuron i can be assigned a spatial coordinate ri. Such a population can be regarded as a continuous medium, each cell being represented as a point of the Euclidean plane. (1.2) reads:

( )

( )

∫ ( ) ( ) [ ]

∂

∂

5

r' r' r' r r r

d t u f ω h t t u

t

u , , , ,

1.3 This is the two-dimensional one-layer lateral-inhibition type Amari model.

Viewing the cortex as a two-dimensional sheet is biologically motivated: the increase in the number of cortical neurons of vertebrates in the course of evolution is almost entirely due to the increase in the surface area of the brain. Furthermore, as the membrane potential is known to be approximately constant across any sufficiently small cortical region, comprising a very large number of neurons (the cortex consists of as many as 10¹⁰ cells) it is not implausible to view it as continuous pattern forming medium.

Suppose that the population is isotropic and homogeneous. That is, cell i influences cell j in the same way as cell j influences cell i, and two cells, situated equally far from a third one, exert the same effect on its dynamics. In other words, ω(r,r’)=w(||r- r’||). Following [Amari 1977] we will be using w(r)=Kexp(-kr²)-Mexp(-mr²), K>M, k>m unless otherwise stated. A typical profile of such a kernel is depicted in Figure 2. The neurons in the immediate proximity of a given cell enhance their activity, while those, situated further away counteract the increase of its potential.

Figure 2

Typical profile of an interaction kernel

Why is equation (1.3) of interest?

First of all (1.3) can be thought of as a natural continuous extension of a classical discrete Hopfield model.

Furthermore, (1.3) is one of the two major neural field models of the cortical activity, (the second one being the above mentioned Wilson-Cowan model). The major difference between the latter and (1.3) is that the nonlinear output function is “outside of the integral” rather then “inside”. The Wilson-Cowan model can be viewed as continuous limit of another discrete neural net, frequently encountered in the literature, sometimes termed activity-based model. For further discussion on the assumptions behind the two models see [Ermentrout 1998].

It was suggested that spatially periodic patterns, exhibited by non-local models such as (1.3) might explain the emergence of visual hallucinations after mescaline or LSD ingestion [Tass 1997]. A direct test of this hypothesis could be obtained by recording cortical activity after the ingestion of a hallucinogenic drug, in order to see whether any periodic patterns emerge. However, to our best knowledge, such experiments have not been conducted.

Supplementing (1.3) with an extra term s, modeling the external sensory input, it proved possible to rationalize the results of psychophysical experients.

( )

( )

( ) ( )

[

]

( )

t

u , , , , ,

r r' r' r' r r r

+ +

+

−

∂ =

∂

∫

5

1.4 s is interpreted as the intensity of visual stimulus and u as a quantity which is proportional to the perceived signal. Let us briefly recapitulate some of the examples described in [Wilson 1973].

Imagine that a grating with consecutive white and black stripes is presented to a subject. Decreasing the contrast, monitoring its level at which the pattern is perceived as homogeneous, a measure of sensitivity is obtained. This was seen to decrease upon increasing spatial frequency of the grating. This effect is readily captured by models such as (1.4). Applying s(r,t) of the form s(x,y,t)=k(1+cos(2nπx/L))/2, and plotting the amplitude of u(r,t) (which in this case becomes spatially periodic with the same spatial period as that of s) against n yields a curve very similar to the one, measured psychophysically.

A visual stimulus is known to be perceived as being weaker if yet another masking signal is present in the visual field. This is termed metacontrast phenomenon.

Such behavior was demonstrated in the Wilson-Cowan model: choosing the term, corresponding to s(r,t) as a spatially confined peak during some initial time interval and subsequently stimulating a neighboring zone, one observes that u(r,t) experiences a decrease in the initially stimulated region.

Yet another example is the so called spatial hysteresis. Assume that a singe point in the visual field of a subject splits in two. At some critical distance they will be perceived as separate. If they are let to approach shortly thereafter, they will be perceived as one whole at a distance well below that at which they originally seemed to split. This behavior is readily captured in neural field models such as (1.4). Representing the two

spots by two peaks of s(r,t), u(r,t) evolves in a way which would be expected from the results of the psychophysical experiments.

There have been attempts to apply (1.4) to autonomous robots [Giese 1999].

There are a number of possible implementations. Usually one two-dimensional field maps the explored environment of the robots, modeling rudimentary memory. The activated regions simply represent positions of the obstacles. The maximum of another one-dimensional field specifies the angle of the movement relative to the lateral axis of the robot. These algorithms were successfully applied to the problem of finding a way through the maze [Giese 1999].

Furthermore, (1.4) is interesting in its own right as a simple (infinite-dimensional) dynamical pattern-forming system.

Finally, we would like to make a brief subjective comment on the plausibility of (1.4). To us it seems that models such as (1.4) are incapable of explaining how human brain is able to implement pattern recognition, discerning among the different objects, encountered in the environment or to perform numerical calculations. However, when the excitability of the neurons is increased (for example under the influence of psychoactive drug), it could become perfectly plausible to regard the cortex as a continuous pattern- forming medium. In this case it is natural to represent the long-range axodendritic interactions with a spatial convolution. This way an interesting and experimentally verifiable description of cortical activity is obtained.

Considering numerous applications, mentioned above, it is of interest to solve (1.3). Particularly interesting are the stationary solutions, satisfying:

( )

∫ (

) [ ( ) ]

5

r' r' r'

r

r t h w f u t d

u , ,

1.5 Suppose f is chosen to be the Heaviside step, as is done in the classical Hopfield model. In this case it would be possible to rewrite (1.5) as:

( ) ( )

[ ]

,

R u

u ^r t = +h

∫

1.6 R[u] designates the region of activation: the subset of the plane where u is positive. It is readily seen that a stationary solution is uniquely determined by its region of activation.

Exploiting this fact it is possible to give a fairly exhaustive classification of the possible solutions of the one-dimensional variant of (1.5). The later is obtained by replacing the integral over ϒ² with an integral over ϒ. A similar approach was adopted for constructing circular solutions of (1.5). Assuming that the region of activation is a disc of some given radius, it proved possible to construct a fairly exhaustive classification of such rotationally invariant solutions [Werner 2001]. Moreover Werner et al. propose that (1.5) admits only one type of bifurcation, schematically depicted in Figure 3. This proposition implies that solutions arise in pairs, one being a stable and the other one being unstable.

The stability properties of the solutions persist for all values of the bifurcation parameter.

Figure 3

Earlier conjectured bifurcation diagram.

The aim of this work is to continue the search for new types of solutions of (1.5) and to further examine their properties. Section (2.1) states the main result: it proved possible to solve the eigenvalue problem governing the stability of circular solutions constructed in [Werner 2001]. This result allowed to disprove the earlier stability conjecture: the branch, corresponding to the stable circular solution was shown to loose stability as the bifurcation parameter increased past some critical value, as is explained in Figure 4.

Figure 4

Actual appearance of the bifurcation diagram.

Sections (2.2) and (2.3) demonstrate how a similar stability analysis can be applied to some other types of solutions: namely the one with an annular region of activation and the one with a stripe-shaped region of activation. These are readily constructed in the same fashion as circular solutions. Their existence was previously suggested in [Werner 2001], but it was never thoroughly investigated.

There is a long standing general question associated with (1.5). Assume that we define a (nonlinear) functional T in the following way:

( ) ( )

( ) ⁰

f

f w d

>

=

∫

r'

T r r r' r'

1.7 for some particular choice of w. Consider some generic f. Assume it is positive on a circle. Clearly, Tf will be constant on its boundary, due to the symmetry properties of the convolution kernel. At this point it is natural to wonder whether the opposite implication holds. That is, for any f define S={r∈ϒ²|f(r)>0}. Assume that S is bounded and simply connected. Suppose Tf in constant on the boundary of S. Is the boundary of S necessarily a circle? This longstanding question remained unanswered during the last thirty years.

Clearly, if (1.5) admits non-rotationally invariant, bounded solutions with simply connected region of activation the answer would be no. Section (4) can be seen as an attempt to tackle this problem. Indeed, as the eigenfunction, corresponding to the largest eigenvalue of the operator, governing stability of solutions of (1.5) is not rotationally invariant, it is natural to expect that there will be a non rotationally invariant solution in the vicinity of the one, undergoing stability loss.

Figure 5

If present, the bifurcating solution can be approximated by a power series in the vicinity of the bifurcation point.

If so, it should be possible to derive an approximation of this solution in terms of power series by exploiting customary techniques of the bifurcation theory. Yet, the expansion is valid only provided that there is an additional branch in the vicinity of the one, turning unstable.

Finally, the last section describes some results on traveling waves and spatiotemporal chaos, encountered in the two-layer modification of (1.5).

1.2 Recapitulation of former results 1.2.1 One-dimensional model

A particularly simple one layer variant of (1.4) was introduced in [Amari 1977]:

( )

(

) ( )

[ ]

( )

t x

u , ,

1

+ +

− +

−

∂ =

∂

∫

5

1.8

The non-linear output function f(u) is modeled by the Heaviside step. s denotes external input which could stem from sensory organs. h is the resting potential: clearly, u is to approach this value if s=f=0. w specifies synaptic weights and is taken to be a difference of two Gaussians: w(r)=Kexp(-kr²)-Mexp(-mr²), K>M>0, k>m>0 (short range excitation, long range inhibition).

A fairly exhaustive classification of possible stationary solutions of (1.8) is given in the original paper. Clearly, they are to satisfy:

( ) ( )

u=

∫

1.9 Assume that the region of excitation (that is, the subset of the real line where u is positive) is an interval (0,a). As (1.9) is translationally invariant, the choice of origin is arbitrary. At this point it is convenient to introduce:

( ) ( ) ( ) ( )

0

0

max lim

x

m x

x

W x w y dy

W W x

W W x

>

∞ →∞

=

=

=

∫

1.10 Substituting the assumption into (1.9) one arrives at an algebraic equation for a:

( )

W a + =h

1.11 This necessary condition for existence of one-bump solution can indeed be shown to be sufficient. Choosing the kernel according to the above assumptions it turns out that W(x) has two intervals of monotony. Being zero at x=0 it increases monotonically, attaining a global maximum Wm. It then decreases monotonically, approaching the value W∞ from above.

By studying the properties of (1.11), using h as an order parameter, the author arrives at a natural classification of the different stationary solutions.

Stability of a stationary one-bump (further denoted by ū) is determined by the dynamics at the boundary. Assume u(xi(t),t)=0, where xi(t) consequently denotes the coordinates of the boundary points of an excited region. Taking partial derivative with respect to t one arrives at:

( )

( )

1,2

2 1

0

1

x

dx W x x h

dt c

c u

x ₌

= −

= ∂

∂

m +

1.12 Using (1.12) together with translational invariance and the fact that the length of an excited region is given by a=x2-x1 it is readily seen that one-bump is stable precisely when

( )

x a

dW w a

dx ₌ = <

1.13 These results together with the existence condition recapitulated above yield three nontrivial types of dynamics:

I: {ϕ, ǎ, ∞} II: {ϕ, ǎ1, a2, ∞} III: {prd}

ǎ1 denotes an unstable one-bump, a denotes a stable bump, ϕ and ∞ are the trivial cases of having no activated region and having all of the domain activated respectively, while prd denotes the presence of infinite families of periodic solutions in the absence of localized bumps.

By considering the dynamics at the boundary of a one-bump solution, the response to stationary input can be predicted. Assuming that s is sufficiently small, (1.12) becomes:

( )

( )

1,2

2 1 1,2

0

1 ( )

x

dx W x x h s x

dt c

c u

x ₌

= − + +

= ∂

∂ m

1.14 The motion of the center of the excited region is determined by:

( ) ( )

( )

2 1

2 1

( ) 1

2

d x x

s x s x

dt c

+ = −

1.15 This implies that the bump is to shift in the direction of increasing input.

However, this optimization process is nonlocal, meaning that the bump will not, in

general, get trapped at the nearest local maximum, as was confirmed by numerical simulations. This suggests that the concept of neural field is potentially applicable to optimization problems.

It is of interest to consider the interactions between two proximal excited regions.

Adopting the same reasoning as when discussing the response to a stationary input it can be shown that two bumps are to approach each other at small distance while they are to repel, provided that the distance between them is sufficiently large.

1.2.2 Two dimensional extension

Cortex, being essentially a two dimensional sheet should be modeled as such in order to obtain a more plausible and physiologically applicable description. The first attempt to extend the above results to two dimensions was undertaken by Taylor. The dynamics is given by:

( )

( ) ( ( )

) ( )

u t

u h w f u t d s t

t

∂ = − + + − +

∂^r

∫

1.16 Now u is a function of two spatial coordinates. The idea is the same as above:

every neuron located at r is assumed to respond to nonlinear input from the rest of the population, convolved with an appropriate kernel.

The starting point is to assume the existence of rotationally invariant stationary solution with a region of activation being a circle of radius R centered at origin. Clearly, such a solution is to satisfy:

( ) ( ) ( ( ) )

( ( ) )

∫∫

∫

− + +

=

=

− +

=

π θ θ

2

0 0

2

2 ' 2 'cos '

'

,

R

d dr rr

r r w r h

d t u f w

h u

5

r' r' r'

r r

1.17 Now, the integral on the right hand side of (1.17) is a function of r=||r|| only. Let it be denoted by CR(r). The necessary condition for existence of a rotationally invariant one-bump solution is given by:

( )

C RR + =h

1.18 which is a two-dimensional counterpart of (1.11). However it is no longer possible to show the sufficiency of this condition.

In [Taylor 1999] the author claims that the global behavior of CR(R) in (1.18) is the same as that of W(x) in (1.11), implying that the classification of localized one-bump solutions carries over to the two-dimensional case essentially without modifications. This was however disproved in the later work by Werner and Richter [Werner 1997].

Although being necessary, condition (1.18) does not exclude the possibility of u(r) given by (1.17) turning negative for r<R. This eliminates certain classes of possible solutions.

Moreover, in the two-dimensional case CR(R) can have as many as three (and, indeed, even more) monotony intervals. In [Werner 1997], Werner and coworkers construct more accurate classification of the possible stationary one-bump solutions. It turns out to be convenient to subdivide parameter range in five different cases, three of which are essentially identical to their one-dimensional counterparts, while others could exhibit new types of dynamics. In those new cases CR(R), having three monotony intervals, allowed for coexistence of three rotationally invariant one-bump solutions.

It is immediate from the above discussion that the stability of such one-bump solution requires ∂CR(r)/∂r to be negative at R. This basically implies stability with respect to rotationally invariant perturbations. Yet, as there is a possibility of noncircular perturbation, the issue of stability has to be considered further.

Finally [Taylor 1999] includes a brief discussion on the dynamics of bump formation and disappearances in the presence of time varying input.

1.2.3 Two bump solutions

As mentioned above, two activated regions are to attract at small distances and to repel when initially situated sufficiently far apart. Therefore the existence of stable two-bump solutions is not to be expected. This issue is further addressed in a recent work by Laing and Troy [Laing 2003]. Assuming that the region of activation comprises two intervals of length a, algebraic equations for a and the inter-bump separation are readily derived. By considering dynamics of the boundary points, stability of the two-bump solution is settled. It was demonstrated that a two-bump solution remains unstable for large regions of parameter space. However, the possibility of having a stable two-bump for arbitrary parameters K, M, k and m, specifying the shape of the kernel, has not been rigorously disproved. A different interaction function such as

( )

(

)

w r = e⁻ −d x +d x −d x⁶

1.19 attaining three zeros for positive r does admit a stable two-bump solution.

Note that the idea behind this analysis is essentially identical to the one, adopted by Amari when considering existence and stability of one-bump solutions.

The existence of multi-bump solutions in the two dimensions has not been treated analytically. Yet, the possibility of having stable two-bumps in the two dimensional Amari model is suggested by numerical experiments, recapitulated in [Laing 2002].

1.2.4 PDE methods

There have been some recent attempts to approach the problem by transforming the original integro-differential equation into partial differential equation which was proved possible by Laing ad coworkers in [Laing 2003]. The authors considered a modified version of (1.16), given by

( ) ( ) ( ( ) )

( )

( )

, ,

2

r u th

u x t

u w x y f u y t d t

f u e H u th

∞

−∞

− −

∂ = − + −

∂

= −

∫

1.20 where H denotes the Heaviside step function. r controls the “steepness” of f at th. The exponential prefactor ensures continuity and differentiability of f at u=th. The resting potential has been incorporated into f. This, however, does not reduce the physical plausibility of the model. Finally an oscillatory kernel was adopted:

( )

(

)

w x =e⁻ b x + x

1.21 On one hand this choice of the interaction function will simplify the analysis below, on the other hand it allows for a rich set of possible stable multi-bump solutions as indicated in the above discussion. Finally, the authors argue that this form of interaction function might provide more faithful representations of “patchy” structure of interneuronal connections in certain regions of the cortex.

Setting time derivative to zero, calculating Fourier transform of (1.20) and making use of the special form of the connectivity function:

( ) ( ) ( )

( ) ( ) ^{( ( ) )}

2

4 2 2 2 2

4 1

2 1 1

i x b b

F u e u x dx F f u

b b

α

α α

∞

−

−∞

= = +

+ − + +

∫

1.22 Multiplying both sides by the denominator of the prefactor to the right and taking inverse transform gives:

( ) ( ) ( ) ^{( )}

( ) ( )

2 2 2 2

'''' 2 1 '' 1 4 1

lim_x , ', '', ''' 0, 0, 0, 0

u b u b u b b f

u u u u

→±∞

⎧ − − + + = +

⎪⎨

⎪⎩ =

u

1.23 The original integral equation (1.20) has been transformed into an equivalent ODE. Stationary solution is clearly given by a homoclinic orbit of (1.23). (1.23) is reversible and Hamiltonian as can be demonstrated by integrating both sides. Such equations gained considerable attention in recent years due to their applicability in fluid dynamics, continuum mechanics and quantum optics [Champneys 1998].

Using (1.23) Taylor and coauthors were able to construct bifurcation diagrams of stationary multi-bump solutions of (1.20).

Moreover, by carefully examining the properties of (1.23) it is possible to deduce an upper bound for b above which no non-constant stationary solutions of (1.20) exist:

4 16 2

0 th

b th

+ −

< ≤

1.24 The procedure of constructing solutions of neural field equations by considering the corresponding PDE are readily generalized to two dimensions by utilizing the properties of two-dimensional Fourier transforms. The starting point is the two- dimensional counterpart of (1.20):

( ) ( ^{( )}

)

u u+ =t

∫∫

1.25 Assuming differentiability of f as in (1.20), applying two-dimensional Fourier- transform to both sides, using F(∇²f)=-(α²+β²)F(f):

( )

( ) ( ) ( ( ) )

(

)

^{( ( ) )} (

) ^{( ( ) )}

1 ,

2

i x y

F u ut e u x y dxdy F w F f u

A F f u G F f u

B M

α β

π

α β α β

∞ ∞ − +

−∞ −∞

+ = =

= ≡ +

+ + −

∫ ∫

s

1.26 The last expression is obtained by replacing F(w) by its Padé approximant (a rational function that matches F(w) and as many of its derivatives as possible). (1.26) is no longer equivalent to (1.25). Yet, it can either be thought of as an approximation of (1.25) or it can be considered as an exact solution of the original problem where w is replaced by the inverse transform of the Padé approximant determined by:

( ) ( )

wˆ sG s J rs d

∞

=

∫

1.27 J0 denotes the Bessel function of the first kind of order zero.

Proceeding as in the one-dimensional case one arrives at:

( ) ( ) ( ) ^{( )} ( ^{( )} )

4 2 2

2 ,

t t t

u u M u u B M u u Af u x y t th

∇ + + ∇ + + + + = , −

1.28 Assuming rotational invariance and resorting to polar coordinates:

( )

( )

4 3 2 2

2

4 3 2 2 3 2

2 1 1 1

2M B M u u_t

r r r r r r r r r r

Af u th

⎡∂ + ∂ − ∂ + ∂ + ⎛ ∂ + ∂ ⎞+ + ⎤ +

⎢∂ ∂ ∂ ∂ ⎜⎝∂ ∂ ⎟⎠ ⎥

⎣ ⎦

= −

1.29 Most importantly, the authors derive stationary rotationally invariant solutions of (1.25) by solving (1.29). Being unstable, such solutions would not have been approached by numerically integrating the original equation.

Numerical experiments suggested that when perturbed such rotationally invariant solutions are to develop into multi-bumps. It is of interest to be able to predict the number of bumps emerging. In order to do so the authors assume a solution of the form:

u(r,θ,t)=ũ(r)+μυ(r,t)cos(mθ). (ũ(r) is rotationally invariant stationary solution.) Substituting into (1.29) and linearizing in μ yields:

( ) _{( )}

4 3 2 2 2 2 2

4 3 2 2 3

4 2 2 4 2 2

4

2 2 2 1 2 1 2

4 2

'

Mr m m Mr

r r r r r r r

m m B M r Mm r

Af u th

r t

υ υ υ

⎡∂ + ∂ +⎛ − − ⎞ ∂ +⎛ + + ∂

⎢∂ ∂ ⎜⎝ ⎟⎠∂ ⎜⎝ ∂

⎣

− + + − ⎤⎛ ∂ ⎞

+ ⎥ ⎜⎥⎦⎝ + ∂ ⎟⎠ %

⎞⎟

⎠

= −

1.30 Being linear (1.30) admits solutions of the form υ(r,t)=ζ(r)e^λt. By numerically integrating (1.30) for each given m and plotting the rate of growth of υ (that is,

∂sup(υ)/∂t≈λ ) versus m an eigenvaue spectrum is obtained. Integer m, maximizing this quantity, is assumed to specify the number of bumps in the multi-bump solution which is being approached. This assumption proved consistent with the results of numerical experiments.

1.2.5 Multilayer extension

The neural field models can be extended to incorporate the dynamics of several neuronal populations. The most immediate extension of the one-layer model would read:

( ) ( ( ) ) ( ) ( ( ) )

( ) ( ( ) ) ( ) ( ( ) )

1

1 11 1 12 2 1

2

2 21 1 22 2 2

' ', ' ' ', '

' ', ' ' ', '

u u w x x f u x t dx w x x f u x t dx h s

t

u u w x x f u x t dx w x x f u x t dx h

t

∂ = − + − − − + +

∂

∂ = − + − − − + +

∂

∫ ∫

∫ ∫

1

s2

1.31 Now the neural tissue is assumed to consist of two distinct neuronal populations.

The activity of excitatory population (u1) enhances that of the neighboring cells, while inhibitory neurons are to suppress the firing rate of their neighbors. Importantly, the one- layer model is nothing but one limiting case of (1.31). In order to see this, let w denote

w11-w12, assume further w11=w21, w12=w22. It is easily seen that u1-u2 is to approach zero exponentially. Consequently, for sufficiently large t substituting u2=u1 into the first line of (1.31), effectively yields (1.8).

In [Amari 1977] Amari undertakes an attempt to consider dynamics of (1.31). It is convenient to set w21 to Dirac delta function, which would imply that the excitatory cells are only capable of exciting the members of the inhibitory population in their immediate proximity. Furthermore, omitting the w22 term significantly simplifies analytical treatment. Under these assumptions the author demonstrates existence of traveling waves.

1.2.6 Lyapunov functional

Finally, it is of interest to point out that Amari model under some mild assumptions admits a Lyapunov functional [Giese 1999]. Namely, requiring that the domain Ω of u is bounded, that the nonlinear output function θ(u) is twice differentiable and strictly monotonically increasing, that w(z) and s(z) are bounded and that w(z) is symmetric (meaning that w(z)= w(-z)):

( )

( ) ( ( )

) ^{( )}

u u t w u t d s h

t θ

Ω

∂ = − + − + +

∂ ^z

∫

1.32 admits a Lyapunov functional of the form:

[ ] ( ) ( ( ) ) ( ^{( )} )

( ) ( ( ) ) ( ( ) )

( ), 0

1 , ,

2

' ,

u t

E u w u t u t d d

d s h u t d

θ θ

ηθ η η θ

Ω Ω

Ω

= − − +

⎡ − + ⎤

⎢ ⎥

⎣ ⎦

∫ ∫

∫ ∫

z z' z z' z z'

z z z

1.33 which is readily seen for example by varying (1.33) with respect to u and obtaining (1.5), the stationary Amari equation.

1.2.7 Wilson-Cowan Model: One-dimensional case

In a classical paper from 1973 Wilson and Cowan [Wilson 1973] put forward a two-layer neural-field model, similar to the one developed by Amari:

[ ]

[ ]

E EE EI

I IE II

E E S w E w I P

t

I I S w E w I

t

∂ = − + ∗ − ∗ +

∂

∂ = − + ∗ − ∗

∂

1.34 E(r,t) and I(r,t) are now interpreted as proportions of excitatory and inhibitory neurons respectively, becoming active at a certain time instant. The “star-product” denotes spatial

convolution. Kernels are taken as exponentially decaying functions of the radial distance.

SI and SE are the nonlinear output functions assumed to be nondecaying, bounded in ±∞

and having unique inflection points (frequently replaced with Heaviside steps in the many later modifications of the model). P denotes external input. The model is based on statistical considerations and a few neurophysiologically motivated assumptions, thoroughly outlined in [Wilson 1973].

The most significant difference between this model and (1.31) is that now the nonlinearity is “outside” of the convolution.

Both models are nothing but limiting cases of the two types of conventional neural nets with finite number of interacting elements. In the customary neural network theory both cases are frequently considered. For a more extensive discussion of physical interpretations of both types of models see [Ermentrout 1998].

Wilson and Cowan conducted extensive numerical studies of (1.34), demonstrating the presence of the following types of behavior:

(1)

Active transients: A brief localized stimulus of the form:

( )

2

0 ,

0

x L P x t p L x L x L

⎧ <

=⎪⎨ ≤ ≤

⎪ >

⎩

1.35 applied during a short time interval continues to grow even after the cessations of input.

The maximum value of E, attained after the onset of stimulations depends on the width of the stimulated region and the amplitude of the stimulus.

(2)

Edge enhancement. The spatial profile of the evoked active transient develops a local minimum in the middle of the stimulated regions, while two local maxima flank the borders of the stimulated zone.

(3)

Spatially localized limit cycles: A stimulated region undergoes a spatially confined oscillation. The frequency of the limit cycle increases monotonically with the stimulus intensity and it is effectively independent of the stimulus width. This mechanism might be utilized by the nervous system in order to encode the amplitude of the stimulus in terms of the limit cycle frequency.

(4)

Frequency demultiplication. When input P is taken to be a square wave, the frequency of the output E is a multiple fraction of that of the stimulus. Interestingly, this phenomenon was shown to occur in the thalamic tissue.

(5)

Traveling waves. Under certain conditions a stimulated region undergoes edge enhancement gradually developing into two distinct peaks which start to move apart in opposite directions.

(6)

Spatially inhomogeneous steady states.

1.2.8 Oscillatory dynamics in two-dimensional case

Spontaneous oscillatory dynamics in the two-dimensional Wilson-Cowan model is considered in [Tass 1997].

Substituting

( ) ( ) ( ) ( )

2 2

ˆ ˆ

, , ⁱ , , ⁱ

E t E t e I t I t e

∈ ∈

=

∑

∑

k Z k Z

r k ^g r k ^g

1.36 into (1.36) and linearizing, the dynamics of a single Fourier mode is given by:

( )

( ) ( ) ( )

( )

ˆ , ˆ ,

ˆ , ˆ ,

E t E t

d B

dt I t I

⎡ ⎤ ⎡ ⎤

⎢ ⎥= ⎢ ⎥

⎢ ⎥ ⎢ ⎥

⎣ ⎦ ⎣ ⎦

k k

k k k t

1.37 where B(k) is a 2*2-matrix. Suppose that an imaginary eigenvalue of B(k) crosses into the positive complex half-plane as one of the model parameters is varied. In this case an oscillatory state is expected to emerge as a consequence of Hopf bifurcation. By utilizing symmetry properties of (1.34) (namely, the rotational invariance of the convolution kernels), it can be shown that three types of stability loss are possible: a single complex pair of eigenvalue crosses into the positive complex half-plane; four complex pairs of eigenvalues cross the imaginary axis; as many as eight complex pairs of eigenvalues cross imaginary axis.

Expressing [Ê(kj,t),Î(kj,t)] (kj denotes the Fourier mode corresponding to eigenvalue pair j) as a linear combination of the eigenvalues of B(k) at the critical point:

( )

( )

ˆ ,

ˆ ,

j

j j

j

E t

I t ξ μ

⎡ ⎤

⎢ ⎥ = +

⎢ ⎥

⎣ ⎦

k Φ Φ

k

1.38

and by utilizing customary projection techniques describe in [Iooss 1989], the author obtains ordinary differential equations, governing the time evolution of ξj and μj for each of the three different cases. This way the full system of two itegro-differential equations can be reduced to an ODE, by effectively assuming that only linearly unstable modes contribute to the resulting oscillatory solution. This assumption can be rigorously justified, provided that the bifurcation parameter is in the vicinity of the critical point.

1.2.9 Traveling waves

In an article by Enculescu et al. [Enculescu 2003] the authors consider the emergence of traveling waves in the two dimensional Wilson-Cowan model. Now, numerical as well as analytical treatment is complicated by introducing an extra dimension. Let us briefly recapitulate the approach, undertaken by Enculescu and coauthors.

Numerical simulations suggested that (1.34) has a solution in the form of a radial traveling wave, emanating from the origin of excitation. That is, assume that the input is given by:

( )

0, 0

P t ⎧P =

= ⎨⎩ ≠ r r

r

1.39 Then a rotationally invariant wave pulse, originating at r=0 will propagate away from the origin of excitation. Consequently, the solution can be approximated as:

( ) (

( ) ( )

1

2

, ,

)

E t g r vt

I t g r vt

= −

= −

r r

1.40 Where r=||r|| ant v denotes the speed of the wave propagation. g1 and g2 are assumed to have the form of a localized peaks.

Formally, we can not substitute (1.40) into (1.34) in order to obtain a particular solution as (1.34) is not invariant under the transformation r→r-vt. However, if one assumes that the wave has propagated far enough from the center of excitation so that E(r,t) does not effectively depend on the y coordinate, the convolution terms can be approximated with:

( )

( ) (

)

0 0

* , ' ' cos '

w E t w r g r vt r r d dr

π

φ φ '

∞

=

∫ ∫

r

1.41 Now, substitution of (1.40) becomes valid. In terms of g1 and g2 (1.34) reads:

[ ]

[ ]

1

1 1

2

2 1

E EE EI

E IE II

vdg g S w g w

dz

vdg g S w g w

dz

− = − + −

− = − + −

2

2

g g

1.42 One assumes further:

( ) ( )

( ) ( )

E E

I I

S z H z V

S z H z V

= −

= −

1.43 where H denotes the Haeviside step function, VE and VI are positive constants, modeling activation thresholds.

Clearly, the left hand side of (1.42) becomes piecewise linear and one can further assume:

( ) ( )

( )

( ) ( )

( )

1 1 2

1 1 2

2 1 2

2 1 2

* * 0,

* * 0, otherw

* * 0, ,

* * 0, otherwise

EE EI E

EE EI E

IE II I

IE II I

G z w g w g V z r vt a

G z w g w g V

G z w g w g V z b c

G z w g w g V

= − − > ≡ − ∈

⎧⎨ = − − <

⎩

= − − > ∈

⎧⎨ = − − <

⎩

, 0 ise

1.44 Using the fact that the fields E and I are unperturbed for z≡r-vt>0 (1.42) is easily integrated, yielding the shape of the traveling wave front. However, parameters v, a, b and c remain undetermined. Their values are fixed by substituting the obtained solution into (1.44) and requiring:

( ) ( ) ( ) ( )

1 1

2 2

0 0

0

G G a

G b G c

= =

= =

1.45

2 Bifurcation and stability

2.1 Bifurcation from a one-bump solution Consider the one-layer Amari-type model:

( )

(

( ) )

t u

u =− + − +

∂

∂

∫

5

r' r' r'

r

2.1 S denotes the Heaviside step function. Below we will frequently be using w(||r- r’||)=w(r,θ,r’,θ’)=w(x,y,x’,y’) to signify the change of variables ((r,θ) and (r’,θ’) are polar coordinates of r and r’ respectively, while (x,y) and (x’,y’) are Cartesian coordinates thereof). This somewhat informal notion is frequently adopted in physics literature.

Consider a stationary solution of (2.1). Obviously, it is to satisfy:

( )

( )

(

( ) )

u =

∫

5

r' r' r'

r r

2.2 Assuming rotational invariance (that is u(r)=u(r)), a one-bump solution ū with a region of activation determined by u(r)>0 ⇔ ||r||≡r<R is readily constructed utilizing the methods recapitulated in [Werner 2001].

Consider a small perturbation of such a solution:

( )

( ) ( )

u r t =u r +εη r t

2.3 (ε - small). Substituting (2.3) into (2.1), keeping terms of order one in ε:

( )

( )

(

) ( ( )

) ( )

u w

t t

t , , ,

, η δ η

η _{=− + −}

∂

∂

∫

5

2.4 δ denotes Diric’s delta function. Recall:

( ( ) ) ⁽ _{( )} ⁾

i '

i

t i

f t t t

f t

δ =

∑

2.5

where ti are zeros of f and f’(ti) are derivatives at the corresponding points. Resorting to polar coordinates, using (2.5), the integral on the right hand side of (2.4) can be rewritten as:

( ) ( ( ) ) ( ) ( ) ( ) ( ) ( )

( )

( ) ( )

(

) (

)

' , ' 1 ,

' , , ,

' ' ' , ' , ' ' ' , ' , ,

,

2

0 2

0 2

0 0

θ θ η θ θ

θ θ

η ρ

θ ρ θ

θ θ

η δ

θ θ

η δ

π π

ρ π

d t R R

r w

Rd t u R

R r w

d dr r t r r u r

r w

d t u

w

R

∫

∫

∫∫

∫

Γ

≡

∂

∂

=

=

−

=

∞ 5

r' r' r' r'

r

2.6 Using (2.6), substituting η(r,t)=e^λtξ(r,θ) into (2.4) yields an eigenvalue problem, governing the stability of ū:

( ) ( )

( ) ( )

0

, , , , , ' , '

r r w r R R

π

d '

λξ θ = −ξ θ + Γ

∫

2.7 Utilizing that ξ(r,θ) enters the right hand side of (2.7) exclusively as ξ(R,θ), denoting ξ(R,θ) by ξ(θ) it is possible to rewrite (2.7) as:

( ) ( ) ( ) ( )

( ) ( ) ( )

2

0 2

0

, , , ' ' '

' ' '

w R R d

g d

π

π

λξ θ ξ θ θ θ ξ θ θ

ξ θ θ θ ξ θ θ

= − + Γ =

− + Γ −

∫

∫

2.8 where Γ=R/|∂ū(R)/∂r| and g is given by:

( )

( )

( ) ( ( ))

( ) ( )

( )

2 ' 2 '

2 2

2 2

2 2

2 1 cos ' 2 1 cos '

2 cos ' 2 cos '

2 2

, , , '

'

i i i i

kR e e mR e e

kR mR

kR mR

kR mR

w R R Ke Me

Ke Me

Ke e Me e g

θ θ θ θ

θ θ θ θ

θ θ θ θ

θ θ

θ θ

− − − −

− − − − − −

− −

− −

= − =

− =

− ≡ −

2.9 Clearly, g(θ) is 2π-periodic. Fourier-expanding the integrand in (2.8):

( ) ( ) ( ) ( )

( ) ( ) ( ) ( )

( ) ( ) ( )

2

0 0 2

0 0

2

0 0

cos ' ' '

cos cos ' ' '

sin sin ' ' '

n n

n n

n n

a n n d

a n n d

a n n d

π

π

π

λξ θ ξ θ θ θ ξ θ θ

ξ θ θ θ ξ θ θ

θ θ ξ θ θ

∞

=

∞

=

∞

=

⎡ ⎤

= − + Γ ⎢⎣ − ⎥⎦ =

− + Γ +

Γ

∫ ∑

∑ ∫

∑ ∫

2.10 Sine-terms are omitted from the expansion of g as g is even. From (2.10) it is readily seen that the eigenvalues and the corresponding eigenfunctions are:

( ) ( ) ( )

2

0

1 cos

cos

n

n

g n

n

π

λ θ θ d

ξ θ

= − + Γ

=

∫

2.11 Clearly, every eigenvalue is two-fold degenerate as sines could have been taken instead. Figure 6 depicts such eigenvalue spectra for a stable one-bump and an unstable one-bump. In the former case all eigenvalues are readily seen to be negative (except the one which corresponds to spatial translation as will be explained below). The eigenvalue spectrum corresponding to a linearization around an unstable solution contains positive eigenvalues, as would be expected. These stability properties are in full accordance with the results of numerical experiments.

(a) (b)

(c) (d)

Figure 6

Stable (a) and unstable (b) one-bump solutions of the Amari-type neural field equation together with their corresponding spectra. Parameters chosen according to: K=2.5, k=5, M=0.5, m=1.5, h=-0.281 in the case of the stable solution and K=2.5, k=5, M=0.5, m=1.5, h=-0.294 for the unstable one.

(a) (b) (c)

(d) (e) (f)

(g) (h) (i)

(j) (k) (l)

Figure 8

Bifurcation from the one-bump solution ((a)-(k)), initially perturbed by a small 4π-periodic deviation (l). K=1.5, k=5, M=0.5, m=1.5, h=-0.146347

Now that ξ(R,θ) is known, it can be substituted into the right hand side of (2.7) in order to obtain the ξ(r,θ) for arbitrary r. We have:

( )

( ) ( )

0

, , , ', ' cos

n r An w r r n d

π

' '

ξ θ =

∫

2.12 Where An is a normalization constant which can be explicitly calculated after defining an appropriate scalar product.

It is important to note that one of the eigenvalues must inevitably be zero, reflecting translational invariance of (2.2). It is readily seen that in the case of a one- bump solution this holds for n=1. In order to see this, imagine that we shift ū(r) by dr where dr is an infinitesimal vector. Restriction of ū(r-dr)- ū(r) to the circle of radius R centered at origin is clearly a multiple of ξ1(R,θ). Consequently:

( ) ( )

2 1

0

0 1 g cos 1 d

π

λ = = − + Γ

∫

2.13 This allows us to reexpress Γ as:

( ) ( )

2

0

1

g cos d

π θ θ θ

Γ =

∫

2.14 Therefore:

( ) ( ) ( ) ( )

2

0 2

0

cos 1

cos

n

g n

g d

π

π

θ θ θd λ

θ θ θ

= − +

∫

∫

2.15 This expression specifies each eigenvalue as a function of model parameters K, k, M, m and the radius of the region of activation of the hypothetical one-bump solution. It will prove extremely useful when deducing the stability properties of one-bump solutions.

Equation (2.15) has a number of interesting immediate implications.

First of all it settles the issue of stability of rotationally invariant one-bump solution of the Amari-type neural field model with respect to arbitrary perturbations. The solution is linearly stable as long as λn given by (2.15) is strictly negative for all n (except n=1).

In [Laing 2002] the authors conduct a numerical study of stability properties of rotationally invariant solutions of neural fields. Their investigation involves numerical derivation of the eigenvalue spectrum as was briefly recapitulated in the introduction.

Using (2.15) it is possible to derive this very spectrum analytically. We did not intend to do so for the case of the modified version of Amari-model, given in [Laing 2002]. Yet, we would like to point out that this is indeed possible. Such a derivation is outlined in the treatment of the stability properties of annular solutions of the original model given below.

Finally, the right hand side of (2.2) can be thought of as a nonlinear functional, defined on a set of some sufficiently well-behaved real-valued functions, having ϒ² as their domain. A stationary solution of Amari-type neural field is therefore a fixed point of this very functional, further denoted by T. Let us equip the domain of T with a scalar product defined by:

( )

( ) ( ) ( )

∫

=

5

r r r r u v d u

v

u, δ

2.16 This is not an appropriate definition as any function v(r) such that ū(r)=0 ⇒ v(r)=0 will have norm zero. However, if we restrict the domains of u and v to the subset of the real plane where ū is zero, this problem is overcome. Since w(||r-r’||)=w(||r’-r||) it is easy to show that the linear part of T becomes selfadjoint if regarded as a functional with the same domain as that of the scalar product (2.16). Consequently, the spectrum of the linear part of T is real. This excludes the possibility of Hopf bifurcations from a stationary solution of (2.1).

Moreover, consider the following proposition: the integral on the right hand side of (2.8) is positive. An immediate corollary would be: any stable circular solution of (2.2) is a stable fixed point of T. The proof is fairly trivial. Stability of the stationary solution implies that the largest eigenvalue of the linear part of T is smaller then one, as is readily seen from (2.11). Assuming that the above proposition holds, we know that all eigenvalues are strictly larger then zero. Consequently, the absolute value of each eigenvalue is smaller then one, implying the above corollary. We did not succeed in proving the proposition, yet, numerical investigations strongly suggest that it holds.

Moreover, numerically calculating Tⁿ(u) for some finite n, the stationary one bump is indeed frequently approached.

2.2 Bifurcation from an annulus-solution

Proceeding in the same way as when constructing a one-bump it can be shown that (2.2) admits stationary solutions with annular regions of activation (u(r)>0 ⇔ ||r||>R1 and

||r||<R2 for some R1, R2). Defining:

( ) ( )

( ) ( )

2

1 2

1

2

1 1 2 1

0 2

2 1 2 2

0

, , 0, ,

, , 0, ,

R

R R

R

W R R w R r rdrd

W R R w R r rdrd

π

π

θ θ

θ θ

=

=

∫ ∫

∫ ∫

2.17 It follows that in order to have an annular solution, W1(R1,R2)=W2(R1,R2)=-h must hold for some R1, R2. Figure 7c depicts a plot of W1-W2 for appropriately chosen model parameters, specified in the figure caption. However, although necessary, W1-W2=0 is not sufficient for the existence. Indeed, it must be checked that u(r), given by

( )

( ( ) )

1

2

2 2

0

' ' 2 ' cos

R

R

u h r w r r rr dr d

π

θ ' θ

= +

∫ ∫

r

2.18 is positive exclusively on the interior of the specified annulus. Although we do not intend to construct a general proof for arbitrarily chosen model parameters, this can be checked explicitly for each set of parameter values of interest. In this way an annular solution can be constructed.

Due to its rotational invariance, the stability of an annular solution (further denoted by ū) can be examined in the same way as that of the one-bump solution.

Consider a small perturbation εη of ū (ε small). Substituting into (2.2) and linearizing in ε:

( ) ( ) ( ) ( ( ) ) ( ) ( ) ∫∫ ( ) ( ( ) ) ( )

∫

∞

+

−

=

− +

−

∂ =

∂

π θ θ δ η θ

η

η δ

η η

2

0 0

' ' , ' ' ' ' , ' , , ,

, ,

,

d dr t r r u r r r w t

t u

w t t t

r

r' r' r'

r r

r

5

2.19 In the last line we resorted to polar coordinates. Using ū(r)=0 ⇔ r=R1,R2 together with

( ( ) ) ^{( )}

( )

i

i

t

i

f t t t

df t dt

δ =

∑