A derivation of a macroscopic field theory of the brain from the quasi-microscopic neural dynamics

ELSEVIER

Physica D 99 (1997) 503-526

A derivation of a macroscopic field theory of the brain from the quasi-microscopic neural dynamics V.K. Jirsa *, H. Haken Institute for Theoretical Physics and Synergetics, University of Stuttgart, Pfaffenwaldring 57/4, 70550 Stuttgart, Germany Received 23 February 1996; revised 12 April 1996; accepted 15 May 1996 Communicated by Y. Kuramoto

Abstract

The purpose of this paper is twofold: First, we present a semi-quantitative nonlinear field theory of the brain under realistic anatomical connectivity conditions describing the interaction between functional units within the brain. This macroscopic field theory is derived from the quasi-microscopic conversion properties of neural populations occurring at synapses and somas. The quasi-microscopic models by Wilson-Cowan (1972,1973) and Nunez (1974) can be derived from these. Functional units are treated as inhomogeneities within a nonlinear one-dimensional neural tissue. Second, for the case of the Kelso experiment the field equation is treated analytically and numerically and can be reduced to a set of ordinary differential equations which corresponds to a model by Jirsa et al. (1994,1995). This phenomenological model reproduces the spatio-temporal phenomena experimentally observed. Here the most prominent property of the neural tissue is the parametric excitation. The macroscopic field parameters can be expressed by quasi-microscopic neural parameters.

1. I n t r o d u c t i o n The brain operates as an open complex system and exhibits spatio-temporal behavior. A necessary condition for this pattern forming character of the brain is a nonlinear dynamics and a spatial interconnection. For a recent review of pattern formation in nonequlibrium systems with many further references see [4]. The functional behavior of the brain is encoded in these spatio-temporal structures and can be extracted from the dynamics of the macroscopic quantities measured by the EEG and MEG [8-10,20]. According to synergetics [14-16] this extraction contains all the relevant information about the spatio-temporal behavior of the brain and has, in general, a small number of degrees of freedom. This idea has been formalized to the order parameter concept based on circular causality: the order parameters are determined and created by the cooperation of microscopic quantities, but at the same time the order parameters govern the qualitative behavior of the whole system. Based on this approach phenomenological models were set up in the past for different experiments in order to find evolution equations that describe the experimentally observed macroscopic dynamics [8,21]. At the quasi-microscopic level ensembles of neurons are gathered together * Corresponding author. Tel.: 0049/711/685-4986; e-mail: viktor@theo, physik, uni-Stuttgart.de. 0167-2789/97/$17.00 Copyright © 1997 Elsevier Science B.V. All rights reserved PII S0167-27 89(96)00166-2

504

V.K. Jirsa, H. Haken/Physica D 99 (1997) 503-526

in units, which tend to behave as single entities. The dynamics of such neural ensembles distinguishing excitatory and inhibitory populations has been modeled by Wilson and Cowan [33,34] by means of action potential densities and by Nunez [29] by means of synaptic activities. The purpose of this paper is twofold: First, the present paper aims at bridging the gap between a quasi-microscopic and a macroscopic description. We define a mathematical description of the input-output behavior, the conversion operations, at the synapses and the somas of neurons in an ensemble. Taking physiological and anatomical conditions into account the quasi-microscopic models of Wilson-Cowan and Nunez can be derived. Here the idea of separation of time scales and spatial scales is utilized such that a spatio-temporal description is obtained relevant for the macroscopic level obtained from EEG and MEG measurements. This description is reduced to a one-variable field equation given as a nonlinear partial differential equation. Here the brain is considered as a nonlinear medium with specific dispersive properties. Functional units like the auditory cortex areas are assumed to interact as inhomogeneities within this medium. The coupling between these inhomogeneities and the neural sheet is nonlinear and given by the field equation. The geometry of the brain, given by the dimension and the boundary conditions of the brain, is an open question and nontrivial. Fuchs et al. [ 10] successfully used the geometry of a closed sphere in the case of EEG measurements of a-waves. Nunez [29] proposes two geometries: a closed sphere and a closed one-dimensional loop. Though we think that the closed sphere might be the best geometrical assumption, here we use the closed one-dimensional loop, since we are mainly interested in extracting the properties of interaction of functional units within the brain. In order to do so the chosen simple geometry suffices. Second, on the basis of the proposed field equation we tackle the Kelso experiment [23]. In this experiment, when hearing an acoustic stimulus, a subject had to press a button in between two consecutive tones, i.e., to perform a syncopated coordination pattern. When the stimulus frequency was increased the subject was no longer able to syncopate and spontaneously switched to the phase of the stimulus, i.e., to a synchronized motion. The sudden and completely involuntary transition from one sensori-motor coordination pattern to another is intimately connected with a change of the brain signals as detected by an array of 37 SQUIDs (Superconducting Quantum Interference Device). The patterns of brain activity are spatially coherent, but their temporal evolution is complex. The spatio-temporal dynamics observed in this experiment was successfully modeled by a set of phenomenological ordinary differential equations [21 ] governing the temporal behavior of the spatial modes. From our proposed field equation we can derive these mode equations under simplest assumptions of the localizations of the auditory and motor areas in the neural sheet. It turns out that the model in [21 ] represents a reduced normal form of that obtained here, just keeping the terms relevant for the phenomena observed in the Kelso experiment. A parametric excitation of the neural tissue is the most fundamental result of this analysis. Since the field equation is obtained from a quasi-microscopic description the macroscopic constant field parameters can be expressed by quasi-microscopic neural parameters. The paper is organized as follows: First, we give a brief survey of physiological and anatomical parameter ranges used here. We then briefly describe the quasi-microscopic Wilson-Cowan model and the Nunez model. After that the idea of conversion operations is introduced from which the macroscopic field equation is derived. We introduce the observed phenomena of the Kelso experiment and discuss analytically and numerically the field equation in this context.

2. Physiological parameter ranges The elementary unit of which the nervous system is composed is the neuron which is divided into three basic components [7]: the dendrites, the cell body and the axon. The dendrites act as the receptive side of the neuron and provide a large surface area for the synapses. The synapses act as small batteries and convert the inputs from other neurons by initiating electric currents on the dendrites which are spatially integrated at the cell body. There are

505

V.K. Jirsa, H. Haken/Physica D 99 (1997) 503-526

Table I Time scales Synaptic delay Neuronal refractory time Inhibitory delay due to propagation

Membrane time constant Excitatory delay due to propagation

1 ms I ms 1

ms

10ms 100 ms

Spatial scales

Spatial range of intracortical fibers Spatial range of corticocortical fibers

0.1 cm 10cm

mainly two kinds of synapses [1,3]: An excitatory synapse causes current to flow into the dendrite at the synapse, along the dendritic cable to the cell body. An inhibitory synapse causes a current to flow in the opposite direction. By Kirchoff's law these currents are linearly summed and cause a resulting current. If the resulting current at the cell body exceeds a certain threshold, it will be converted into a pulse train along the axon without attenuation. Pools of neurons in local neighborhoods tend to share activity and thus can be regarded as the quasi-microscopic entities that perform spatially coherent behavior. The EEG measures macroscopic quantities which mainly correspond to summed dendritic potentials [7] and the MEG measures macroscopic quantities mainly corresponding to summed dendritic currents [30]. Note that the precise relationship between the electric potentials of single dendrites and the resulting summed dendritic potential of an ensemble of neurons remains a topic of future investigation. The cortical neurons are connected via intracortical fibers over short distances which can be excitatory or inhibitory, but also via corticocortical fibers through the white matter over long distances [ 1]. The latter are exclusively excitatory. The spatial range of the dendritic trees is comparable to the intracortical connection range [ 1]. The ranges of physiological parameters vary a lot depending on the cortical areas and types of fibers considered, We emphasize that here we only want to give an impression of the relation of the parameter scales to each other. Detailed discussions of parameter ranges can be found in [ 1,3,25,30]. Synaptic delays and refractory times are of the order 1 ms, the neuronal membrane constant is in the range of several ms [3]. Corticocortical propagation velocities have a wide range from 0.2 rn/s [25] to 6-9 m/s [29]. Here we will use 1 m/s as an average propagation velocity. The lengths of corticocortical fibers range from 1 cm to 20 cm [29] which yields a broad range of delays from 10 to 200 ms. In Table 1 we give a survey of the parameter scales on which our subsequent discussion will be based. Single neurons have two main state variables [7]: first, the dendritic currents which are generated by active synapses and cause the waves of extracellular fields; second, the axon pulse frequency. The conversion of pulse frequency to current amplitude occurs at synapses, the dendritic wave amplitude is converted to a pulse frequency at the somas. In ensembles of cortical neurons activity densities are defined over spatial distributions of neurons. Here the pulse-to-wave conversion at the synapses is constrained to a linear small-signal range, whereas the wave-to-pulse conversion shows a sigmoidal behavior [7]. These two conversion operations for ensembles are shown in Fig. 1.

3. Models of brain activity There has been a vast variety of attempts in mathematical modeling of the dynamical phenomena observed on various spatial scales in the brain. On the microscopic level the nonlinear Hodgkin-Huxley equations serve as a model for the spatially localized generation of action potentials at the neuronal membrane [ 19], whereas the core conductor model, a second-order partial differential equation in space and first-order in time, describes the transfer properties of axons and dendrites [7]. The Hodgkin-Huxley equations were reduced by FitzHugh and Nagumo (FHN) [5,27] into a two-variable model which captures the key phenomena and is analytically tractable. By introducing spatial

506

V.K. Jirsa, H. Haken/Physica D 99 (1997) 503-526

Conversion o p e r a t i o n s

PO - - -~

Wave-to-pulse

i I ....

~.

,

.

., . .

Pulse-to-wave

Fig. 1. On the left the wave-to-pulse conversion operation performed at the somas of neurons in an ensemble is shown. This conversion shows a sigmoidal behavior. On the right the same situation for the pulse-to-wave conversion performed at the synapses in a neuronal ensemble is shown. Here the conversion is constrained to a small signal range. The functional forms of these conversions presented here are only valid for neuronal ensembles and differ from one of the single neurons.

diffusion into these equations the spatio-temporal FHN system is obtained which demonstrates the existence of traveling pulses that only propagate if a certain threshold perturbation is exceeded [26]. Mesoscopic theories aim at describing the spatio-temporal development of the overall mean level of electric or magnetic activity in neuronal populations. First models of this type were published by Beurle [2] in 1956 and Griffith [12,13] in 1963. Both used partial differential equations to describe the propagation of a field that describes the overall excitation of the nervous network. In the early 1970s Freeman [6] set up a model of the olfactory system (KIII model) mainly consisting of three oscillators. These oscillators are described by ordinary second-order differential equations that are coupled via excitatory and inhibitory connections by long pathways with distributed delays. In the following we first present in more detail two established mathematical models describing the spatio-temporal dynamics of localized neuronal populations. In the first model by Wilson and Cowan [34] the relevant variable is represented by the proportion of neurons becoming active per time unit which corresponds to the generation of action potentials at the somas of neurons in an ensemble. The second model by Nunez [28] uses the number of active synapses per unit volume at each instant in time. Both models distinguish between excitatory and inhibitory activities of neuron ensembles. In Section 4 we present our more general approach.

3.1. The Wilson-Cowan model In 1972 Wilson and Cowan [33] set up equations governing the dynamics of the generation, the so-called firing, of action potentials in localized populations of neurons distinguishing excitatory and inhibitory subpopulations. The proportion of excitatory cells in the corresponding subpopulation firing per unit time at the instant time t is denoted by E(t), the proportion of inhibitory cells in the corresponding subpopulation is denoted by l(t). The proposed evolution equations of the variables E(t) and l(t) obtain the following form:

VK. Jirsa, H. Haken/Physica D 99 (1997) 503-526

507

Iz[~(t) = - E + Se(Ne),

(1)

/zl(t) = - I + Si(Ni)

(2)

Nj : a ( c j l E - cj21 q- Pj).

(3)

with

Here the dot denotes the first derivative with respect to the time. The indices e and i denote excitatory and inhibitory subpopulations, Sj (Nj) the sigmoid function, /z the neural membrane time constant, a is the coefficient of the membrane impulse response, cjk are constant parameters, Pj is some external input to a population of class j neurons and Nj the resulting value of excitation of a class j population. Eqs. (1)-(3) were extended by Wilson and Cowan [34] for spatially interacting excitatory and inhibitory neural populations in a homogeneous one-dimensional neural sheet and read as follows: Iz

Iz

OE(x, t + r) Ot Ol(x, t + r) Ot

--

E ( x , t) + Se(Ne),

(4)

--

l ( x , t) + Si(Ni),

(5)

where

Ne(x,t)=otlz

peflee(x-X)

E X,t

-- f Pi~ie(x-X) Ni(x,t)=ot#

IOeflei(x-X) E X,t

Ix-X__ I

Ue

dX

[ x - g ' ) dX"l-Pe(x,t)l

Ix-Xl_ 13e

(6)

dX

-f piflii(x-X)l(X,t Ix-XI)dx+Pi(x,t)

(7)

--OG

Here the following symbols and expressions were used: x, X denote the location in space, r the synaptic delay time, pj is the space constant for connectivity and vj the propagation velocity, fljk denotes the connection density of the probability that neurons of class k are connected with those of class j at a distance I x - X I. 3.2. The Nunez model In 1974 Nunez [28] published a model which links synaptic action to action potential firings. Here synaptic action is to be understood as the proportion of active synapses at the time t in either an excitatory or inhibitory subpopulation. Since single synapses cause dendritic currents which are linearly summed along a dendrite (see Section 2), synaptic action of a neural ensemble can here be identified with the corresponding summed dendritic

V.K. Jirsa, H. Haken/Physica D 99 (1997) 503-526

508

currents. The excitatory synaptic action ~Pe(X,t) and the inhibitory synaptic action ~i(X, t) are governed by the following set of equations: Oo

~Pe(X,t)=u(x,t)+fdvfdeXRe(x,X,v)g(X,t o F

~i(X, t) = u0 + f dxX gi(x, X) g(X, t),

IX-Xv I),

(8) (9)

F where F represents the unfolded surface area of the brain, g(x, t) the action potential firings, u (x, t) the excitatory input to the neural sheet and u0 the corresponding inhibitory input which is assumed to be constant during a fixed physiological state. These equations express the idea that action potential firings g(X, t) at cortical location X cause synaptic activity at location x at later times depending on the separation distance I x - X I and action potential velocity v. The number density of corticocortical fibers connecting locations X with x with the propagation velocity v is described by the distribution function Re(x, X, v). The length of intracortical fibers which lead to inhibitory neurons is short, thus the delay due to propagation is negligible and the corresponding distribution function Ri (x, X) is independent of the velocity v. The action potential firings g(x, t) are dependent on the synaptic inputs ~Pe(X,t) and ~i (X, t). The relationship between action potentials and synaptic inputs is known to be sigmoidal. Here Nunez makes a linear expansion of this sigmoid function around an assumed fixed physiological state go and obtains the following relationship:

g(x, t) = go + Oe~e(X, t) - Oi~i(x, t),

(lO)

where 0g ae = - ~ e g=go'

0g ai = --~-~ g=go"

(11)

4. Derivation of a field theory of the brain

In this section we first want to discuss the possible conversion operations at the synapses and the somas of neural ensembles and derive from these the Wilson-Cowan model for the dynamics of the action potentials and an extended Nunez model for the dynamics of the synaptic activities. We will consider a one-dimensional neural tissue with periodic boundaries. We simplify the derived extended Nunez model and obtain a one-variable neural tissue equation. Second, we derive a nonlinear partial differential equation from the simplified neural tissue equation.

4.1. Conversion operations In Section 2 we pointed out that there are two conversion operations in the neural tissue: the pulse-wave conversion at the synapses and the wave-pulse conversion at the somas. The term pulse corresponds to the action potential densities E(x, t) and I (x, t) of the Wilson-Cowan model and the term wave to the synaptic activities ~Pe(x, t) and ~Pi(x, t) of the Nunez model. Here we interpret these four quantities as deviations from a fixed physiological state and rewrite them in the following form:

V.K.Jirsa,H.Haken/PhysicaD 99(1997)503-526 = f ox e(X.

t).

~i(X,

f

t) = dX fi(x, X)Hi(x, X, t),

F E(x. t) = f dX

5O9 (12)

F

I(x, t) = f dX fl(x, X)Hl(X,X, t).

X. t).

F

(13)

F

fk (x, X) the

The functions Hk (x, X, t) represent the output of a conversion operation and corresponding distribution function depending on the spatial connectivity. The considered surface area of the brain is denoted by F. In order to set up equations for the conversion output we make the following considerations (see Section 2): - Excitatory neurons only have excitatory synapses, inhibitory neurons only inhibitory synapses [ 1]. - In ensembles of neurons the pulse-wave conversion at the synapses is linear, the wave-pulse conversion at the somas nonlinear and sigmoidal. - The spatial distribution range of the dendrites and intracortical fibers is very short. Only the corticocortial connections cause a significant delay via propagation. - External input is realized such that afferent fibers make synaptic connections. These items lead to the following relations between conversion output and summed action potentials:

He(x,X,t)=S(E(X,t

Ix-Xve I ) ) ~ a e E ( X , t

ni(x,X,t)=S(l(X,t

Ix-Xl))'mail(

Ix-X__ve ') '

(14)

Ix-X[)vi

(15)

and between conversion output and summed synaptic activities

HE(x,X,t)=Se(~re(X,t

[x-Xve

I)-Oi(

Hl(x,X,t)=Si(~e(X, t

I x-X_ 1)e

l) - ~ i ( X , /

X'/

Ix-X_vi I)-k-Pe( X't

Ix-Xv

I x - X I)-~-pi(X,t Pi

I x - X ')) 13 '

pj (X, t)

')) ,

(16) (17)

Sj

where ae and ai are constant parameters, external input to the neural sheet and S and the sigmoid functions of a class j ensemble. The propagation velocities re, vi and v are assumed to be fixed with a small variance. The distribution functions fi(x, X), fE(x, X) and fi(x, X) are of short range and can be assumed to be g-like. Inserting these relations into (12) and (13) we obtain for the dynamics of the excitatory and inhibitory action potential densities of a neuronal ensemble the following set of equations:

E(x,t)=Se(f dXfe(x,X)aeE(X,

Ix--XI)-f r dXfi(x'X)ail(X't)+pe(x't))

'(x,,)=s,(f

'x-Xi)-f " r dXA(x' e X)a"(X't)+pi(x")) r

(18)

(19)

These equations correspond to the Wilson-Cowan model for the case that the time scale # of the neural membrane is much smaller than the one given by the delay via corticocortical fibers and that the spatial range of intracortical fibers is very short. Comparing (18) and (19) with (4) and (5) we can identify the distribution functions

fe(X, X) : flee(X -- X) ~---flei(X -- X),

(20)

j~(x, X) : flii(x - X) : ~ie(X - X)

(21)

V.K.Jirsa, H. Haken/PhysicaD 99 (1997) 503-526

510 and the parameters ae = ot/zpe,

ai = o~/zpi.

(22)

For the dynamics of the synaptic activities we obtain the following equations:

ape(x,t)=ae f d X f e ( x , X ) S e ( a p e ( X , t -

Ix-Xve l ) - a p i ( X , t

Ix-X__ve ])

F

api(X, t) = a i

/ dX fi(x, X)Si(ape(X, t) - api(X, t) + pi(X, t) ).

(24)

F The linearized versions of these equations correspond to the original Nunez model (8-10), which has recently been modified to include some nonlinear effects [30]. Here we provide an alternate approach to the nonlinearity. Note that in contrast to (8) and (9) here the external input pj (X, t) is included into the nonlinear sigmoid function.

4.2. The neural field equation As pointed out in Section 2 the EEG and MEG mainly measures macroscopic quantities generated by dendritic potentials and currents which correspond to synaptic activities. Thus we focus in the following on the evolution equations (23) and (24). We specify the sigmoid function Sj (nj) as the logistic curve 1

Sj = 1 + e x p ( - v j n j + vjOj)

1

1 +exp(vjOj)'

(25)

where Oj denotes a fixed physiological state, usually the excitation threshold, and vj denotes the sensitivity coefficient of response of the corresponding neural subset. Since we interpret apj as deviations from Oj, we can set Oj = 0. We expand (25) into a Taylor series up to third-order in nj and obtain

Sj(nj) ~ ctjnj - ~ot)nj43 3 ,

(26)

where otj = ¼vj. Here the sigmoid function is approximated by odd orders in nj due to the choice Oj = 0. With increasing distance to this inflexion point second-order terms of nj will turn up which can be eliminated by a linear transformation. Next, we look at (24) in more detail. The time scale of the intrinsic dynamics of api(x, t) is given by the synaptic decay time and the delay via propagation along intracortical fibers. But this system is also coupled to the excitatory synaptic activity ape(x, t - I x - X I ~re) which operates on a much slower time scale. We assume the connectivity functions to be of the following form:

fi(x, X) = 1- - e -Ix-XI/ai ~ 8(x - X),

(27)

fe(x, X) = -1 - e -ix-xa/'°,

(28)

2tri

2ae

V.K. Jirsa, H. Haken / Physica D 99 (1997) 503-526

511

where (27) takes into account that intracortical connections are mainly local. Taking only linear contributions of (24) into account we obtain the following behavior of the inhibitory synaptic activity: aioti ~ri(x, t) ~ - (~re(X, t) q- pi(x, t)). 1 + aioti

(29)

Here the dynamics of lPi(x, t) is expressed in terms of the leading order of the slowly varying quantities ~Pe(X,t) and Pi (x, t), which means that on this time scale the intrinsic dynamics of ~Pi(x, t) is negligible. The higher-order contributions of these quantities cause small modifications of the corresponding parameters and are neglected here. Inserting (29) into (23), we readily obtain for the dynamics of the excitatory synaptic activity:

qZe(X,t)=aef dX f e ( x , X )

Se(/5~re(X,t

Ix-XI)+P(

Ix-X__oe I ) ) ,

(30)

F where

p(X, T) = pe(X, t)

aioti - pi(X, t) 1 + aioti

(31)

and /5 = 1

aioti 1 + ait~i

(32)

Eq. (30) represents a one-variable neural tissue which will serve as the starting point of the following analysis. We will follow the idea that the input signals to the neural sheet can be understood as inhomogeneities which are embedded into the neural sheet. These inhomogeneities represent the localized areas in the brain that perform different tasks. The coupling of these inhomogeneities to the neural sheet is described by (30). This yields the following form of the stimulus p(x, t):

p(x, t) = iS(x) s(t),

(33)

where/~(x) defines the spatial properties and s(t) the temporal behavior of the input signal. In the case of m input signals at different locations with a different temporal behavior we obtain the following formulation: m

m

p(x, t) = ~-'~pi(x, t) = ~_,fli(x) si(t). i=1

(34)

i=1

Expressing the time delay via propagation along the corticocortical fibers by a delta function ~ (t - T - I x - X I /re) we can rewrite (30) as follows: OO

(35)

qle(x,t)=ffG(x-X,t-T)p(X,T)dXdT F -o¢~

with the Greens function

G(x - X, t - T) = ~ ( t - T

I x - X I]

1 e_lX_Xi/~re

J

(36)

and 4 3 p(X, T) = ae[Cte(/5 ~e(X, T) -t- p(X, T)) - ~t~ e (/5 ~e(X, T) d- p(X, T))3].

(37)

V.K.Jirsa, H. Haken/PhysicaD 99 (1997)503-526

512

We perform the following Fourier transformations: OG

OO

lffeikx-i°~tlpe(k,09)dkdog, ~Pe(X,t) -- (2:n') 2

(38)

--OO --OO

p(x, t) --

lffeikx-i°'tp(k,

(2rr) 2

09)dkd09,

(39)

09)dkd09,

(40)

--OO --OO OO

(3O

lffeik~-i°~tog(k, G(se, to) - (2rr)2 --OO --OO

where

~=x-X

to=t-T

(41)

and obtain the relation

~e(k, o)) = g(k, 09)p(k, o9).

(42)

The Fourier transform g(k, 09) of the Greens function G(~, to) can be determined as

g(k, 09) =

f/

d~

--OO

dto G(~, to)e -ik~+i°Jt°

(43)

(v2k2~ (090- i09)2)

--O~

with the parameter I) e

too = - - .

(44)

O"e

Rewriting (42) in the space and time domain the following partial differential equation is obtained: ~e+(092--o2A) lPe+2090~te =

092 "1-090~- p(x,t)

(45)

with

p(x,t)=aeSe /5 ~e(X, t) -t-

pi(x,t) i=1

~ a e Ore(/5~Pe(X,t) + Z p i ( x , i=1

(~ ~e(X,

Zpi(x,

.

(46)

i=l

Here the Laplacian is denoted by .4. Our field equation governs the spatio-temporal behavior of the electromagnetic brain activity. Functional units Pi (x, t) convey diverse input to the neural sheet. The output signals ~j (t) from the neural sheet that govern the behavior of other parts of the body are defined by

~j(t) = f

dx flj(X)~e(X, t),

F where flj (x) defines the spatial localization of the jth functional output unit.

(47)

V.K. Jirsa, H. Haken/Physica D 99 (1997) 503-526

513

5. The neural field under periodic stimulation In this section we want to investigate the wavelike properties of the neural field equation (45) for the case of the MEG experiment of Kelso et al. [23]. We first summarize the experimental phenomena and the conducted analyses of this experiment in Section 5.1 and recall a phenomenological model [21 ] reproducing the observed phenomena. In Section 5.2 we specify our field equation (45) for the case of the Kelso experiment and derive the experimentally observed properties. 5.1. The motor-brain-behavior experiment by Kelso et al.

In this experiment a subject was exposed to acoustic stimuli and pressed a button in a syncopated motion. The stimulus frequency at the beginning was set to 1 Hz and was increased by 0.25 Hz after l0 stimulus repetitions up to 2.25 Hz. We shall refer to a domain of constant stimulus frequency as a plateau. Around the frequency of 1.75 Hz, at the border of plateaus III and IV, the subject switched spontaneously to a synchronized motion. During this experiment the magnetic field data were recorded over the left parieto-temporal cortex, mainly covering the motor and auditory areas. Detailed analyses conducted by Fuchs et al. [ l l] revealed that in the pretransition region of plateaus I-III the registered brain signals oscillate mainly with the stimulus frequency. Stimulus and motor response frequency are locked l: 1. Then a switch in behavior occurs at the border of plateaus III and IV. The brain signals oscillate with twice the stimulus frequency in the post-transition region of plateaus IV-VI. Applying a Fourier transformation to the signals and looking at the component of the same frequency as the stimulus one observes that the phase difference between stimulus signal and the time series is stable in the pretransition region, performs an abrupt change of Jr at the transition point and remains stable again in the post-transition region. Near the transition point, predicted features of nonequilibrium phase transitions like critical slowing down and fluctuation enhancement are observed in both the behavioral data and the brain signals. In order to obtain information about the spatio-temporal behavior of the brain signals from the entire SQUID array a Karhunen-Lobve-expansion (KL expansion) was applied [ll ]. This procedure is also known as Principal Component Analysis or Singular Value Decomposition. The values of the magnetic field from the different channels at a specific time form a 37-dimensional vector H(t), which may be linearly expanded into a set of orthonormal spatial modes v ~k) with corresponding amplitudes xk (t) that depend on time: 37

H(t) = Z

xk (t)v (k).

(48)

k=l

The modes v tk) are chosen as the eigenvectors of the covariance matrix. By construction these KL modes have the property that the eigenvalues Xk are a measure of the contribution of the spatial mode v tk) to the whole signal. These modes are ordered such that the mode with the largest eigenvalue is k = 1, the mode with the next smaller eigenvalue is k = 2, etc.; thus the main parts of the spatio-temporal dynamics are represented by a few modes that belong to the lower eigenvalues (the eigenvalues with large amplitudes), whereas the higher modes (the eigenvalues with small amplitudes) represent the distortions or noise. The neglect of higher modes and consideration of only the few modes that are relevant for the dynamics corresponds to a reduction o f dimension in space. The KL modes v ~k) of the experimental time series, taken over all plateaus, were determined. Then the brain signals of the single plateaus were projected onto the modes v (k) and the following observation was made: In the pretransition region the main dynamics is dominated by the first KL mode v (l) oscillating with the stimulus frequency, then at the transition point a switch occurs and the dynamics of the signals are dominated by the second

514

V.K. Jirsa, H. Haken/Physica D 99 (1997) 503-526

KL mode v ~2) oscillating with twice the stimulus frequency [23]. The relative phase between the stimulus signal and the first Fourier component of the first KL mode is locked in the pretransition region. Here two possible stable stationary states of the relative phase are observed, thus a bistable situation is present. At the transition point a transition by Jr in the relative phase from one stationary state to the other is observed and a monostable situation is present: Only one stationary state of the relative phase exists. It also appears that the first Fourier component of the second KL mode performs such a transition in the phase relative to the stimulus signal. But since here the amplitude is very small, this transition is not very clear. In 1994 Friedrich, Kelso and we set up a phenomenological model [21] consisting of two ordinary second-order differential equations. In this model we interpret the two first spatial KL modes as order parameters whose temporal dynamics is described by their time-dependent amplitudes, denoted here by ~0 and ~l. In 1995 we presented a more detailed spatio-temporal analysis of the experimental MEG signal [20] based on the variational method by Uhl et al. [31,32]. This approach showed that in the case of the Kelso experiment the KL modes serve as a first linear approximation of the spatial structure. It turned out that the first-order parameter state has to be represented by a superposition of two spatial modes, where the one with the larger contribution to the entire MEG signal is similar to the first KL mode and the second has a higher symmetry. Their temporal dynamics is given by the oscillator ~0 which is determined by the first model equation (49) that we will present below. The second-order parameter state is described by a spatial pattern similar to one of the second KL mode, but being distorted. Its temporal behavior is governed by the oscillator ~l determined by the second model equation (53). The finer spatial and temporal structures were also determined by means of enslaved variables. Here we will focus on the gross features of the experimentally observed phenomena which are satisfactorily described by the two first leading KL modes. The phenomenological model system [21] describing the temporal dynamics of ~0 and ~l has the following form:

~0 + 09~0 + Yl~O + bll~2~O + b12~2~0 + fc1 + [¢2 + [¢3 = O,

(49)

[q = el sin 2.f2t ~o,

(50)

fez = Cl sin I2t

~2,

(51)

f~3 = dl~l

(52)

~1 q- 092~1 + Y2~I -t- b21~2~1 q- b22~2~1 + k4 q- f~5 + [¢6 = O,

(53)

f¢4 = E2 sin 212t ~l,

(54)

f¢5 "~"C2 sin 12t ~02,

(55)

f¢6 = d2~0,

(56)

and

where ogj, yj, bjk, Ej, cj and dj represent constant parameters. These phenomenological parameters were chosen such that the phenomena observed in the Kelso experiment could be reproduced. This model will be discussed in comparison to the field equation in the following Section 5.2.

5. 2. The field equation under periodic stimulation In this section we want to treat the field equation (45) analytically with respect to the Kelso experiment. In this paper we introduced the idea to view the brain as a nonlinear medium where functional units are embedded as inhomogeneities. For the case of the Kelso experiment, a motor-brain-behavior experiment, the following situation

515

V.K. J irsa, H. Haken / Physica D 99 (1997) 503-526

External

[ Finger oscillator L

Auditory stimulus

Functional units

Cortexwith neuralfield

Internal

~e (x, t)

I I Corticocortical

lit

X

Fig. 2. Three functionalunits, the auditory,the motor and the sensori-motorareas, are embeddedas inhomogeneitiesin the neural sheet. is given and sketched in Fig. 2 The spatio-temporal behavior of the brain is described by the quantity ~Pe(x, t) which represents an internal variable and whose dynamics is governed by the field equation (45). Three external quantities are present: First, the auditory periodic input pa(x, t) to the neural sheet whose temporal dynamics is given by Sa(t) and the localization of the auditory area in the neural sheet by tia(X). Second, the motor output ~m(t) from the neural sheet to the finger oscillator which performs a periodic movement. This output signal is picked up from the motor areas in the brain whose localization is given by tim(X). Third, the finger performs an oscillation with the amplitude z(t) which is transferred to the brain as a sensori-motor input signal Psm(X, t) via the sensori-motor areas defined by tism(X) with the temporal behavior given by Ssm(t). These three functional units, one output and two input components, interact with the neural field via the nonlinear coupling (46). First we consider only the neural sheet under acoustic periodic stimulation. Then we also introduce the motor and sensori-motor behavior.

5.2.1. Input: acoustic periodic stimulus The input function Pa (x, t) of the external acoustic periodic stimulus depends on the stimulus localization tia (x) and the temporal behavior Sa(t). Here the temporal stimulus behavior is taken to be sinusoidal, Sa(t) = sin 12t, where ~2 denotes the stimulus frequency of the acoustic stimulus. Using (33) we write the auditory stimulus signal as

pa(X, t) = tia(X) sin 12t.

(57)

Inserting (57) into (45) we obtain the following field equation after enumerating the nonlinearities in the right-hand side of (45): ~e + (~2~ - vZA) ~e + Y0 ~e + a ~3 + B ~e2"~e -t'- g l q-" K2 -t'- K3 = 0,

(58)

516

V.K. Jirsa, H. H a k e n / P h y s i c a D 99 (1997) 5 0 3 - 5 2 6

Table 2 Parameters of the field equation

wo

Ve/ae

/22 Y0 A

092 (1 - Ote~ae) + 2fla(X)2aeOt3O92/5 too (2 -- Ue~ae) + 2fla(X)2aeOt3W0p 4/3ae~3W2,63

B E

4aet~e3W0t33 2~a (x)2aet/3a)0p

CO0

Ve / tTe

C D1

4/~a (x)aea3toOp 2 -ae(~a(X)Oteto0 - ~a(X)3a3to0)

D2

- fla (x ) 3ae ~ 3 O.~O

where KI = e ( 2 ~ sin 2$2taPe - cos 212t d e -- tO0 COS2Qt ~Pe),

(59)

K2 = C(~(2 cos I2t ~p2 + 2 sin I 2 t ~e~te + O90sin ~2t ~pe2),

(60)

K3 = DI (/2 cos I2t + o90 sin S2t) + D2 (I2 cos 312t + ½o90sin 312t).

(61)

Here KI represents a parametric excitation of the neural tissue with twice the stimulus frequency, K2 a nonlinear term with a parametric coefficient oscillating with the stimulus frequency and K3 a linear periodic driving of the neural tissue with the frequencies/2 and 312. Eq. (58) determines the dynamics of a field containing linear and nonlinear damping terms, if Y0, B > 0, and an amplitude dependent frequency. The parameters in (58) can be expressed in terms of quasi-microscopic anatomical and physiological parameters and are given in Table 2.

5.2.2. Analytical treatment of the field equation Here we want to treat the field equation (58) analytically, where we first consider only the auditorily excited neural sheet. Since we do not know the auditory stimulus localization we make the approximation that global stimulation is only considered fla(X) ~" fla ~---constant.

(62)

This assumption is certainly not realistic, but it already provides some insight into the dynamics of the neural sheet. Taking a finite distribution length of the stimulus into account, the impact of the stimulus will become dependent on the distribution length. Further, in a two-dimensional neural sheet the stimulus localization will not only show radial, but also angular dependencies. These properties will be an area for further study. We will consider only standing waves in the analytical discussion and perform a mode decomposition such that only the first two modes are taken into account assuming that all the higher modes are damped: 1

~e(X, t) = Z

~n(t)e inkx,

(63)

n=--I

where ~n (t) = ~-*n(t) and the asterisk indicates the complex conjugate. In our present one-dimensional description the first experimental KL mode corresponds to the mode with n = 0 in (63) and the second KL mode to n = 1. Standing waves will be obtained if the time-dependent amplitudes ~n (t) are real times a factor e ix which corresponds to a spatial constant phase shift for n # 0. For reasons of convenience we set X = 0. Inserting (63) into (45) we obtain the following system of ordinary differential equations:

V.K. Jirsa, H. Haken/Physica D 99 (1997) 503-526

517

3 ~0 + ~2se0 + ~'o~o + a ( se2 + 6~ 2) ~o + B(~2~o + 2~2~0 + 4~0~I~1) + Z ki = O, i=l

(64)

kl = e ( 2 ~ sin 252t s~0 -- cos 2 ~ t ~o - if)0 cos 212t ~0),

(65)

k2 = C(5"2 cos S2t(~ 0 + 2~?) + sin ~2t (2s%~0 + 4s~l~l + O9o (sed + 2s~2))),

(66)

k3 -----Dr (~2 cos 12t + coo sin Y2t) + D2(f2 cos 3f2t + ½wo sin 3~2t)

(67)

and

~1 + (I2~ + v2k 2) ~1 + Vo~I + 3A(~02 + ~ ) ~l + B(3~2~ + ~g~l + 2~0~J~0) + k4 + k5 = 0,

(68)

k4 = ~(2$'2 sin 2S2t sel - cos 2I2t ~1 - ~0 cos 212t ~ ),

(69)

k5 = C(2S2 cos I2t ~0~1 + 2 sin Y2t (~0~l + ~0~l + w0~0~l)).

(70)

Let us turn to a comparison of the model system (64), (68) derived from the quasi-microscopic dynamics and the phenomenological model system I (49),(53). The equations (64), (68) and (49), (53) qualitatively contain the same terms except of an amplitude-dependent frequency contribution in (64) and (68) given by the term with the parameter A. We will see later that this contribution becomes negligible in the case of the Kelso experiment. The terms kl, k4, kl and k4 each represent a parametric excitation with twice the stimulus frequency, the terms k2, ks, k2 and [c5 a second-order term multiplied by the periodic stimulus. The linear driving term k3 does not appear in phenomenological model. It can be shown that this term is of lower significance for the phenomena in the Ke]so experiment, since it causes purely monostable behavior in the relative phase between brain and stimulus signal. The terms/c3 and k6 in the phenomenological model cause the transition in the relative phase of the second KL mode. As pointed out in Section 5.2 this phenomenon is not clear and the terms k'3 and/:6 actually do not appear in Eqs. (64) and (68), but note that these terms would appear in (64) and (68), if the simplification of global stimulation (62) was not made and a localized stimulus was considered. The most prominent feature of the above equations is the parametric excitation of the neural sheet with twice the stimulus frequency. The parametric excitation has two main characteristics: First, it provides a frequency selection by enforcing stable and unstable solutions dependent on the relation between coi and Y2. The unstable solutions are obtained for o)i/2S2 = k,

(71 )

where k = ½, 1, 3 .... Second, the parametric excitation causes a bistable situation for the phase of the unstable Fourier component. In contrast a linear driving term like k3 causes a purely monostable situation which is not observed in the Kelso experiment. Thus the parametric excitation must be the dominant excitation. The dynamics of the model system (49),(53) is as follows: In the pretransition region the oscillator ~0 is in the unstable region of k = ½ and oscillates with the stimulus frequency I2. The oscillator ~l is in the stable region and hence is damped. When the stimulus frequency S2 is increased the oscillator sel reaches the unstable region k = 1 at a critical frequency and starts to increase. Due to this increase of~el via the coupling terms ~:2 and/:5 the relative phase between the first Fourier component of ~0 and the stimulus signal performs a sudden transition by zr. The amplitude of the oscillator ~1, mainly oscillating with twice the stimulus frequency, becomes large and the oscillator ~0 becomes damped via 1Note that the parameters d I and c2 were set equal 0 in [21], since the contribution of the corresponding terms does not qualitatively change the behavior of the model system. In [20] they were introduced for reasons of symmetry.

V.K. Jirsa, H. Haken/Physica D 99 (1997) 503-526

518

the cubic cross coupling to ~l. In [21] the model system (49), (53) was treated analytically and numerically. Here we will follow the same lines as in the analytical treatment of (64) and (68). We make the following ansatz for the amplitudes ~0 and ~1: 2

~j = Z

~(m)eimI2t +c.c.

with j = 0, 1,

(72)

m=0

where c.c. denotes the complex conjugate. The complex amplitude ~(m) is slowly varying in time and higher harmonics are neglected. We perform two approximations well known in the theory of nonlinear oscillators. The slowly varying amplitude approximation means that we neglect terms ~i) compared to 52~j(i). In the present context the rotating wave approximation means that we may neglect terms of higher order than e 2i~2t . Next, we use the method of adiabatic elimination to obtain a lower-dimensional description of the system. This approach allows us to express the stable variables in terms of the unstable ones, respectively, setting the nonlinear terms containing stable variables to zero. We thus obtain equations solely in terms of the unstable variables ~0(l) and ~2). In order to examine the behavior of the phases and amplitudes separately we further make the hypothesis ~0(±1) = XI0 e ±i4~

(73)

~+2) = Y20 e±iO,

(74)

and

where xl0, y20 and 4,, 6) may be time-dependent. These procedures lead to the following set of equations for the amplitudes: )C10 =

('-~Yo+~Ecos24' '

~

.Y20 =

( - 7}'0+ l ~ee0(1 l +cos

20) -

)

sinZq~ Xlo--B

xZo+Z~y

220)

XlO,

(75)

~Eo~ ) 8--~ sin20 Y2o

-B((I+eo)x210+(~+6e2)y20)Y20

(76)

and for the phases q~ __ 092 2-~'-2-522

- COS24, + C~ ( cos 4, - ~o90 y20 16 sin 2q~ -- ~a.~Z - sin~) --,xlo

0 = o92 - 4522 -- 1/2690EE0 1 690eE0 452 + ~Ee0 sin20 -- 8--~ cos20,

(77)

(78)

where we used

o92 = 522,

o92 = 522 + v2k 2,

~ = 1 + 2Eo + 2%2,

Go = E69o/2w2.

(79)

Contributions of the additive driving terms with Dl, D2 were neglected. We want to discuss the qualitative behavior of (75)-(78) at the threshold where the phase transition is performed. This region is defined by the stationary solution 0st of the phase 0 which is given by ~b0 s i n 20st -- ~ c o s 20st

~

]sin 20st 1=

o92 -- 452 2 -- 1/2ff-,'OEE0 52~0

< --

1.

(80)

V.K. Jirsa, H. Haken / Physica D 99 (1997) 503-526

510

This region corresponds to the instability area of (71) for k = 1. The phase 4, has two stationary values at ~ Taking these into account we determine the time-independent stationary solutions Xst of xlo and yst of 3'20 and perform a linear stability analysis. Here it turns out that the only linearly stable stationary solutions are 2 ~5-- 2}'0 Xst -2B '

y2 = 0, st

(81)

which will also be shown numerically in Section 5.2.6. The nontrivial stationary solution Yst is unstable, thus the oscillator sel will remain damped even beyond the critical stimulus frequency and the phase ~ will remain bistable with two stationary values at 0 and ;r.

5.2.3. Output: motor signal In 1985 the periodic movement of the index finger has successfully been modeled by a Van der Pol-Rayleighoscillator [18]. In order to keep our analysis as simple as possible here we assume a linearly damped oscillator model for the amplitude of the displacement z(t) of the finger tip. The motion of the finger tip is driven by the motor areas in the neural tissue. These areas convey the signal ~m(t) to the finger oscillator z(t). According to (47) this motor signal is defined by dx tim(X)lPe(X, t),

~m(t) ~- I

(82)

F where tim (x) represents the spatial distribution of the motor areas. Taking only linear terms into account the dynamics of the finger displacement z(t) is given by

z(t) 4- y~(t) 4- o92z(t) ~- C~m(t - tf),

(83)

where c is a constant parameter and tf represents the time delay caused by the propagation from the motor areas to the finger tip. Here we assumed the finger oscillation to be linearly driven by the motor areas in the brain. Note that it might be interesting to investigate whether the coupling between motor signal ~m (t) and finger displacement z(t) could also be described to be sigmoidal as the coupling of the nervous tissue in (46). Peripheral nerve axons carry action potentials in the 100 m / s range, which is about 100 times faster than corticocortical axons. Thus the time delay tf is in the 10 ms range which is much smaller than the time scale 2zr/Y2 given by the periodic stimulus and we neglect tf in the following. In the brain the transition of the relative phase between brain signal and stimulus mainly occurs for the first KL mode which performs an oscillation with the stimulus frequency. These phenomena are also observed in the motor behavior which gives the motivation to choose the spatial distribution tim (x) such that ~m(t) represents the first KL mode which can be realized in the simplest case via a global coupling tim(X) ~ tim = constant,

(84)

which causes ~m (t) ~ se0(t). In the case of linear damping, y > 0 in (83), the oscillator z (t) will perform an enforced movement with a fixed parameter-dependent phase shift 3 between z(t) and ~m(t). Neglecting the homogeneous solution of (83) the inhomogeneous solution is OO

C

z(t) = ~

f

~m (O~)

j

o92 --co---2+ iyo9

ei~Ot &o.

(85)

--OO

Assuming that ~m (t) mainly oscillates with the stimulus frequency I2, i.e., l~m (t) -~ 1~0 COS a'?t, then the inhomogeneous solution in (85) can be written as z ( t ) = C 1~0 ((O92 -- $22) 2 4- (yS'2)2)-1/2 COS($-2t -- 3)

with S = }'.(2(W 2 -- ,Q2).

(86)

520

V.K. Jirsa, H. Haken/ Physica D 99 (1997) 503-526

Assuming ~ ~ 0 we obtain z(t) -~ CO~m(t) with co --

constant.

(87)

~,/(O02 _ ,.Q2)2 _1._(ya"2)2 5.2.4. Input: sensori-motor signal The information about the finger motion z(t) is conveyed to the sensori-motor areas in the brain and serves according to (34) as an additional input Psm(X, t) given by Psm(X, t) = J~sm(X)Ssm(t)

with Ssm(t) = z(t - tf),

(88)

where/~sm (x) represents the spatial distribution of the sensori-motor areas and Ssm(t) the temporal behavior of the sensori-motor signal. Taking into account the sensori-motor input Psm(X, t) as an additional input to the neural field ~e(X, t) we obtain according to (45) the following field equation ~e + (oo2 - v2A) 0e + 2o00@e = ae ( o92 + m o ~0 ) Se[p_ 0e(X, t) + pa(X, t) + Psm(X, t)],

(89)

which governs the spatio-temporal dynamics of the neural field under periodic acoustic stimulation and periodic finger movement. In order to simplify the subsequent calculations we assume that the acoustic input to the neural field is dominant and thus take only linear contributions of the sensori-motor input Psm(x, t) into account. The time delay tf caused by the propagation from the finger tip to the sensori-motor areas will again be neglected. Using 0902<< woO/Ot we can also drop the small linear frequency contribution of the sensori-motor signal. With (87) and (88) we obtain the following neural field equation: 7)e + (o002- v~z3) Oe + 2 o ~ e = ae

(

°)

o)02+ o0o-~ &[~ Oe(X, t) + pa(x, t)] - V~(X)~m(t),

(90)

where gl (x) = aeOleCOWO~sm(x). Under the simplifying assumption of global coupling ¢lsm(X) = ¢lsm -= 1

(91)

the term Vl (x) becomes a constant parameter F1 (x) = Vl = -aeaeCOo00. 5.2.5. Further analytical treatment of the field equation In Section 5.2.2 we determined the stationary solutions of the amplitudes xl0 and Y20 and the phases ¢ and 0 at the threshold, as well as conditions for their existence under purely auditory stimulation. Here we analytically discuss the field equation (90) which takes auditory and sensori-motor input into account. We consider an increased contribution of the parametric driving with ff,~ >mo and neglect the contributions of the additive driving terms with D1, D2. We obtain the following evolution equations of the amplitudes: - ~ - sin2~ Y20 =

(

- 1Vo + ¼eeo(1 + cos20)

)

~-~ sin20

(12 ) ((2)

xlo- B

Y2o - B

~Xlo + 2 ( ~ - %)Y2o xlo, 1+

x20 + 3(1 + 2'2)Y

(92) Y2o

(93)

and the phases ~ = 092-a'22 2~

,ffgO ( ¼Esin2~b-- -~-~-cos2~b+C cos~b

) y2 wO sinq~ 2__.o xlo'

(94)

V.K.Jirsa, H. Haken/PhysicaD 99 (1997)503-526 --

EEOff)

0

0 = 0)2 - 4122412 1/2cbo~Eo + ¼E~0 sin 20 - - - 8 1 2 cos 20.

521 (95)

The increased contribution by the term ~ causes a stronger deviation from the stationary values 0, Jr of ~b. This deviation can be calculated analytically and it proves that it does not cause any changes in the qualitative behavior of q~. Below the threshold, 12 < 12c, the following stable stationary solutions are obtained: 2 Echo - 2~()'o + )'1) Xst = 212B ,

2 Yst = 0.

(96)

The phase 0 becomes locked if cos20

-

7-2Y2sin 20 0)0

~rcos201

=

L0) 2 --4122 -- 1/25)OEEO 2 " -_ EE00)0

< 1.

(97)

The trivial stationary solution of Y20 becomes unstable if ½Yo + YI + ¼ " o + ½~o()'o + )'1) > -

~Eo - 812

E~

+ -212

(98) '

where the equality sign holds at the critical frequency 12 = 12c. Here the stable amplitude xlo can be adiabatically eliminated and expressed in terms of the unstable amplitude Y20, Xl0 =

(Ed~O-212()'O+)'I))l/2[l(212B + )'0 +B(1 0 . 5-,+00) ). 5' 1' ~ 1 2 -1 y 2 0 ) ] "

(99)

Above the threshold, 12 > ~c, the following stable stationary solutions are obtained: 2

Xst

----

0,

y2 = E~0(212 + 5~0) - 4 1 2 ~ 1212B(1 + 2E2)

(100)

Here no coexistence of the two states ~0 and ~l is present, but rather a selection between these states occurs. The amplitude xm decreases to 0 above the threshold. In the order of approximation we used in the present analytical discussion the terms K2 in (60) do not appear in the amplitude equation (92) of xm. These terms cause a small, but nonzero amplitude xl0 for 12 > 12c- Note that for stationary values of 4~, 0 the amplitude equations (92) and (93) correspond to the selection equations of the synergetic computer [ 17] that Haken introduced as a paradigm for pattern recognition in the brain. The phase ~b will show bistability if the terms sin 2~b and cos 2q~ dominate the dynamics in (94) which corresponds to a strong parametric excitation of the neural sheet. Similar phase equations were obtained and studied by Kuramoto [24] for discrete coupled populations of integrate-and-fire neurons. In contrast here the phase q~represents a collective variable. The stationary solutions of q~ correspond to minima in a potential field V which is given by a gradient dynamics:

$ = -OVla4~.

(lO1)

The terms cos ~b and sin ~b cause distortions of this potential and deviations of the minima from 0, Jr. Here we have the following situation: Below the threshold the terms cos q~ and sin q~ are negligible and bistability is present. At the threshold the amplitude of Y20 increases and xl0 decreases. These changes cause an increase of the contributions of the terms cos q~ and sin ~ which enforce monostability of V. Neglecting all kinds of distorsions and deviations the potential V can qualitatively be written as V(~b) = - a -ys2t - cos~b - b cos 2q~ Xst

(102)

522

V.K. Jirsa, H. Haken/Physica D 99 (1997) 503-526 V

V

7

/2 o; ......i

V

AV

V

/7-_

Fig. 3. Starting in the top row the potential V is plotted from left to right for increasing values of y2t/Xst. and is plotted in Fig. 3 in dependence of varying values of the quotient y2/xst. This behavior of ~ will be numerically tested in Section 5.2.6. To summarize, under the assumption of global coupling of the inhomogeneities to the neural tissue the field equation (90) was treated analytically such that it was reduced to a set of ordinary first-order differential equations describing the dynamics of amplitude and phase of the first and second field mode. The stationary solutions of these amplitudes and phases were determined and their stability investigated. Parameter ranges can be determined in which the phase transition observed in the Kelso experiment occurs. The derived equations (64) and (68) qualitatively correspond to the phenomenological mode equations (49) and (53) in [21] which can be considered as a reduced normal form of the spatio-temporal dynamics of the brain in the case of the Kelso experiment. ! 2.6. Numerical treatment of the field equation Here we want to present a numerical discussion of the field equation (90) under periodic stimulation and periodic finger tapping analytically treated in Section 5.2.5. So far we treated the neural tissue as a one-dimensional medium with periodic boundaries. These assumptions are obviously unrealistic. Thus this numerical study is not concerned with a detailed analysis of the possible phenomena in the parameter space, but rather with the proof that the proposed field equation is capable of producing phenoma like the ones observed in the Kelso experiment. We will use one fixed set of quasi-microscopic neural parameters and determine the consequent field parameters used in (90). We will use this set of field parameters in the numerical analysis and only consider minimal modifications, e.g. an increased value of &0, in order to ensure the existence and stability of the desired nontrivial stationary solutions. These parameter modifications can be interpreted, for instance, as caused by different stimulus localizations, a higher dimension or different geometry of the neural tissue or nonlinear contributions of the sensori-motor signals. Under the assumption of global coupling as in (62) and (91) the field parameters can be expressed by quasimicroscopic parameters as given in Table 2. In Tables 3 and 4 the sets of quasi-microscopic neural and consequent field parameters are given where time and space are scaled such that the time unit corresponds to 20 ms and the space unit to 1 cm. In the brackets the values are given that were used for the numerical simulations.

523

V.K. Jirsa, H. Haken/Physica D 99 (1997) 503-526

Table 3 Neural parameters v

1.52

fl0

1

/~a JOe c~

3.6 1 0.7

tre R0 t5 # C~e

10 10 0.9 1 1.2

Table 4 Field parameters

A F0 D1

0.15 0.85 (0.5) 5 x 10-3 1 (0.03) 1.4 (0)

,(2o d~0 B C D2

0.35 0.15 (2.3) 0.1 (0.4) 0.42 (0.1) - 1.6 (0)

Here R0 represents the considered extension of the neural sheet. Note that the contributions of the linear driving forces given by the neural parameters are strong which are a consequence of the chosen stimulus localization, the global coupling in (62). As discussed in the preceding section these linear driving forces destroy the bistability in the relative phase between the field mode ~0 and the stimulus signal, thus we assume these forces to be small under different stimulus localization or geometric conditions of the neural sheet and set them equal 0 in the numerical simulations. Since these terms are not apparent or only weakly present in the phenomena of the Kelso experiment, it is also worth investigating whether such linear driving forces are not quantitatively cancelled, but qualitatively averaged out, for instance via spatial averaging within an ensemble of neurons. All the following numerical simulations are based on (90) where the sigmoid function was expanded according to (46). The boundary conditions are periodic and only standing waves are observed. We use a semi-implicit FTCS procedure (Forward-Time-Central-Space) for integration. First we present the numerical simulations of (90) in the case without the sensori-motor modifications, i.e., Yl = 0. In Fig. 4 the time series of the field modes ~0 (bold line) and ~l (slim line) are plotted together with the stimulus signal (dotted line). The first row shows the pretransition region with I2 = 0.31 corresponding to the pretransition region where ~0 dominates oscillating with the stimulus frequency in in-phase, the second row the same situation in anti-phase. Note that here the terms in-phase and anti-phase refer to the situation where the motor behavior is in in-phase or anti-phase with respect to the stimulus signal. Due to possible fixed phase shifts between ~0 and the motor signal the field mode ~0 may be shifted with respect to the auditory stimulus. Both stationary solutions are stable and all the higher modes are damped, thus bistability in the relative phase between ~0 and the stimulus signal is present. Increasing the stimulus frequency up to 12 = 0.45 corresponding to the post-transition region the dynamic state of the field remains the same, only the amplitude of the dominating mode ~0 decreases. Here without the sensori-motor feedback ~1 does not have a stable nontrivial stationary solution. This is a prediction of the spatio-temporal dynamics in the brain for the case of no motor movements and purely acoustic periodic stimulation, which can be experimentally tested. Fig. 5 shows the following time series with Yl = 0.9 : The first and second rows present the pretransition situation at I2 = 0.31 where ~0 dominates oscillating with ~ in in-phase (first row) and in anti-phase (second row). All higher modes are damped. The third row shows the transition regime at I-2 = 0.4 where the second mode ~l comes up oscillating with 212 and the first mode ~0 performs a transition by Jr in the relative phase to the stimulus signal. The bottom row shows the stationary situation in the post-transition region at I2 = 0.47 dominated by ~l. Here ~0 has a smaller amplitude and is monostable in the relative phase to the stimulus signal. Note that the temporal behavior of ~0 directly corresponds to the dynamics of the motor behavior in the case of the Kelso experiment, since the finger

V.K.Jirsa,H. Haken/PhysicaD 99 (1997)503-526

524

0.00.':"

.:"

.:::

.:::" 1.00

:/

.:" 2.00

.;"

:':' '":".~ ":""'V :"""~ 3.00

i..,.. A.,.,A.,-IA..,,. A.-, i,..,i..,., i...,i J Iv.. .j. v...... V....j, V..i-. v....., v.....-.V.". V..J. v./. V..., 1.00

2.00

3.00

Fig. 4. The numerical integration of the field equation under periodic stimulation without sensori-motor modifications is shown. Here the amplitudes are scaled in arbitrary units and the time in s. The presented time series show the temporal behavior of the field modes ~0 (bold line) and ~1 (slim line) and the external stimulus signal (dotted line). All the higher field modes are damped. In the first row in the pretransition region ~0 dominates oscillating with the stimulus frequency in the in-phase state, in the second row in the anti-phase state. This situation remains the same, even when the stimulus frequency is increased into the post-transition region.

oscillator has been modeled in Section 5.2.3 as linearly driven by the motor areas in the brain which are localized such that ~0 is selected as the driving force. This statement is only valid for the case of the Kelso experiment as can be seen in Fig. 4 where ~0 performs a dynamic behavior, even when the sensori-motor modification is switched off.

6. Discussion, summary and conclusions The present paper bridges the gap between a quasi-microscopic and a macroscopic description of the spatiotemporal dynamics of the brain. We introduced a mathematical description of the notion of conversion operations at the synapses and the somas of neurons in an ensemble which represents the quasi-microscopic level of description. Considering a time scale hierarchy as well as a spatial scale hierarchy we could derive equations governing the spatio-temporal behavior of synaptic activities and of action potential densities. The considered time scale of 100 ms is due to delays via propagation of action potentials along corticocortical fibers, the considered spatial scale of several cm is due to the spatial range of excitatory corticocortical fibers and short range inhibitory fibers. These scales are anatomically realistic. On these time and spatial scales the derived equations of synaptic activity correspond to an extension of the model by Nunez [28] and the ones of pulse activity to the model by WilsonCowan [34]. The macroscopic level of description is given by EEG and MEG measurements that mainly pick up the signals of dendritic potentials and currents which correspond to synaptic activities. We reduced the equations of the synaptic activities to a one-variable field equation given as a nonlinear partial differential equation. In this order of approximation around a fixed physiological state inhibition plays a minor role and could be eliminated. It contributes to the constant field parameters as inhibitory background effects. This approximation will not be valid any more, if higher order patterns, i.e. smaller spatial scales, in the EEG or MEG are considered. We followed the paradigm that the brain acts as a nonlinear medium with dispersive properties. Functional units like the auditory and motor cortex areas are assumed to interact as inhomogeneities with this medium. These inhomogeneities obviously have a strong impact onto the dispersive properties of the brain. Here we made the approximation of global coupling of the inhomogeneities to the medium brain. One of the major tasks in the future will be the handling of the localization of these inhomogeneities. We discussed the field equation for the case of the Kelso experiment and could derive from this the phenomenological model in [21] which governs the spatio-temporal brain dynamics in this case. Here we

525

VK. Jirsa, tl. Haken / Physica D 99 (1997) 503-526

, 1.00

A1

2.00

3.00

100 " '".

:"" :"" ' ""'" "' ' '"'' /' " :"'. ' :"" :"" A"... A"., A"... A,. A', A":. A',. A".. .100 -V -V/ V / V / -V/ -V/ V: 1 1.00

1.00L

,

'

2.00

,

'

0.000"50 '~ "-. ~ ~: "'- - : ' '" : ." ':. . :' . "" ' "°'5°b V

~

.

-

~

0.50

1.00L ,

0.s0[~ ~

'

:"

""

-

1.00

,

3.00

',

"

/ 1.50

i

'

L

"'

',

/:

~

~ :

{

"

'"

" -/

2.00

/

"

1 "'

- .-t

2.50

i

,

A

A~q

A 'A

4

A

-1.00 0.50

1.00

1.50

2.00

Fig. 5. The first two rows show the same situation as in Fig. 4, but here with the inclusion of the sensori-motor modification. The third row shows the transition regime at I2 = 0.4 where the second mode ~! comes up oscillating with 2.(2 and the first mode 40 performs a transition by zr in the relative phase to the stimulus signal. The bottom row shows the stationary situation in the post-transition region at .c2 = 0.47 dominated by ~1. c o u l d prove that the parametric excitation o f the neural tissue postulated in [21 ] is a fundamental property o f the brain. Our p r o p o s e d n o n l i n e a r field theory opens a vast variety o f possible future theoretical and experimental investigations on the basis o f the ideas presented here. Particularly periodic p h e n o m e n a like periodic driving or self-generated oscillations like in the case o f a - w a v e s are interesting to v i e w under this aspect.

Acknowledgements

We gratefully a c k n o w l e d g e useful discussions with M. Bestehorn, R. Friedrich and A. Fuchs. We wish to thank A. Schiiz and P. N u n e z for helpful c o m m e n t s about a n a t o m y and physiology.

References

[1] [2] [3] [4] [5] [6] [7]

M. Abeles, Corticonics (Cambridge University Press, Cambridge, 1991 ). R.L Beude, Properties of a mass of cells capable of regenerating pulses, Phil. Trans. Soc. London Set. A 240 (1956) 55-94. V. Braitenberg and A. Schiiz, Anatomy of the Cortex. Statistics and Geometry (Springer, Berlin, 1991). M.C. Cross and EC. Hohenberg, Pattern formation outside of equilibrium, Rev. MOd. Phys. 65 (1993) 851. R. FitzHugh, Impulses and physiological states in theoretical models of nerve membrane, Biophys. J. l (1961) 445-466. W.J. Freeman, Mass Action in the Nervous System (Academic Press, New York, 1975). W.J. Freeman, Tutorial on neurobiology: From single neurons to brain chaos, Int. J. Bifurc. and Chaos 2 (3) (1992) 451-482.

526

V.K. Jirsa, H. Haken/Physica D 99 (1997) 503-526

[8] R. Friedrich, A. Fuchs and H. Haken, Spatio-temporal EEG patterns, in: Rhythms in Physiological Systems, eds. H. Haken and H.P. Koepchen (Springer, Berlin). [9] R. Friedrich and C. Uhl, Synergetic analysis of human electroencephalograms: Petit-Mal epilepsy, in: Evolution of Dynamical Structures in Complex Systems, eds. R. Friedrich and A. Wunderlin (Springer, Berlin, 1992). [10] A. Fuchs, R. Friedrich, H. Haken and D. Lehmann, Spatio-temporal analysis of multichannel a-EEG map series, in: Computational Systems - Natural and Artificial, ed. H. Haken (Springer, Berlin). [11] A. Fuchs, J.A.S. Kelso and H. Haken, Phase transitions in the human brain: Spatial mode dynamics, Int. J. Bifurc. and Chaos 2 (4) (1992) 917-939. [12] J.S. Griffith, A field theory of neural nets: I: Derivation of field equations, Bull. Math. Biophys. 25 (1963) 111-120. [13] J.S. Griffith, A field theory of neural nets: II: Properties of field equations, Bull. Math. Biophys. 27 (1965) 187-195. [14] H. Haken, Synergetics. An Introduction, 3rd Ed. (Springer, Berlin, 1983). [ 15] H. Haken, Advanced Synergetics, 2nd Ed. (Springer, Berlin, 1987). [16] H. Haken, Information and Self-Organization (Springer, Berlin, 1988). [ 17] H. Haken, Synergetic Computers and Cognition (Springer, Berlin, 1991). [ 18] H. Haken, J.A.S. Kelso and H. Bunz, A theoretical model of phase transitions in human hand movements, Biol. Cybernet. 51 (1985) 347-356. [19] A.L. Hodgkin and A.E Huxley, A quantitative description of membrane current and its application to conduction and excitation in nerve, J. Physiol. 117 (1952) 500-544. [20] V.K. Jirsa, R. Friedrich and H. Haken, Reconstruction of the spatio-temporal dynamics of a human magnetoencephalogram, Physica D 89 (1995) 100-122. [21] V.K. Jirsa and R. Friedrich, H. Haken and J.A.S. Kelso, A theoretical model of phase transitions in the human brain, Biol. Cybernet. 71 (1994) 27-35. [22] J.A.S. Keiso, On the oscillatory basis of movement, Bull. Psychon. Soc. 18 (1981) 63. [23] J.A.S. Kelso, S.L. Bressler, S. Buchanan, G.C. DeGuzman, M. Ding, A. Fuchs and T. Holroyd, A phase transition in human brain and behavior, Phys. Lett. A 169 (1992) 134--144. [24] Y. Kuramoto, Collective synchronization of pulse-coupled oscillators and excitable units, Physica D 50 (1991) 15-30. [25] R. Miller, Representation of brief temporal patterns, Hebbian synapses, and the left-hemisphere dominance for phoneme recognition, Psychobiology 15 (3) (1987) 241-247. [26] J.D. Murray, Mathematical Biology, 2nd Ed. (Springer, Berlin, 1993). [27] J. Nagumo, S. Arimoto and S. Yoshizawa, An active pulse transmission line simulating nerve axon, Proc. IRE 50 (1962) 2061-2070. [28] P.L. Nunez, The brain wave equation: A model for the EEG, Math. Biosci. 21 (1974) 279-297. [29] P.L. Nunez, Electric Fields of the Brain (Oxford University Press, Oxford, 1981). [30] P.L. Nunez, Neocortical Dynamics and Human EEG Rhythms (Oxford University Press, Oxford, 1995). [31 ] C. Uhl, R. Friedrich and H. Haken, Reconstruction of spatio-temporal signals of complex systems, Z. Phys. B 92 (1993) 211-219. [32] C. Uhl, R. Friedrich and H. Haken, Analysis of spatio-temporal signals of complex systems, Phys. Rev. E 51 (5) (1995). [33] H.R. Wilson and J.D. Cowan, Excitatory and inhibitory interactions in localized populations of model neurons, Biophys. J. 12 (1972) 1-24. [34] H.R. Wilson and J.D. Cowan, A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue, Kybernetik 13 (1973) 55-80.