Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://num-anal.srcc.msu.ru/list_wrk/ps1/ch4st3.ps Äàòà èçìåíåíèÿ: Tue Dec 17 13:00:58 2002 Äàòà èíäåêñèðîâàíèÿ: Mon Oct 1 20:46:36 2012 Êîäèðîâêà:

Mathematical Modeling of Complex Information Processing Systems 135
MATHEMATICAL AND NUMERICAL ANALYSIS OF PERIODIC SOLUTIONS NEAR
STATIONARY STATES IN A NONLINEAR MODEL OF BRAIN CORTEX ACTIVATION
A. Fraguela 1 , A. I. Greb'ennikov 2 , and J. F. Leyva 1
A nonlinear ordinary differential equation is considered as a macroscopic (global) model of brain
cortex activation. The nonlinearity has the form of a sigmoid function of two parameters a and b.
The zoning of the plane (b; a) with respect to the number and type of stationary states (critical
points) for the model under consideration is performed. The existence of the center­type solutions
in a particular domain of the plane (b; a) is proved analytically. A numerical algorithm for finding
periodic solutions is proposed; this algorithm was realized with the MATLAB system and was verified
using several model examples.
1. Introduction. The model of brain cortex activation presented in the form of a nonlinear ordinary
differential equation was first formulated by P. N~unez [1] as a result of analysis done on a system of two
integral equations that describe a space­time evolution of synaptic activity of the brain cortex. This model
was considered in [2] under some assumptions concerning anatomical and functional characteristics of the brain
cortex (such as its columnar organization, the distribution of synaptic connections, etc.) and was based on
the use of the activation variable u = u(y; Ü ). The authors of [2] proposed the following differential equation
describing the above model:
@ 2 u
@Ü 2
+ 2 @u
@Ü
+ u \Gamma
@ 2 u
@y 2
=
Ÿ
g(u) + @g(u)
@Ü
–
(1)
Here y and Ü are the space and time variables, respectively. The function g(u) (called a shot function) is
expressed by
g(u; a; b) = 1
1 + a exp(\Gammabu) (2)
where the parameters a and b are positive. This function represents a fraction of the total number of neurons
of a column that shoot the action potentials over a period of time t. Its behavior is sigmoidal [2], well known
in problems of this type. All variables and parameters of this model are dimensionless.
2. The equation of stationary state. Analysis of critical points. In order to obtain an equation
describing the stationary states, it is necessary to eliminate the time­dependent terms entering into (1). As a
result, we obtain the equation
d 2 u
dy 2
+ g(u; a; b) \Gamma u = 0 (3)
Let us consider its periodic solutions corresponding to the boundary conditions u (0) = u (l) and u 0 (0) = u 0 (l),
where l is a period.
2.1. Analysis of stationary states. Equation (3) can be transformed to the following Hamiltonian
system of differential equations:
du
dy
= z
dz
dy
= \GammaF (u)
(4)
Here the function F (u) = u \Gamma g(u; a; b) is auxiliary. The Hamiltonian (or the function of total energy) is of the
form [3]
H(u; z) = T (z) + U (u) (5)
where T (z) = z 2
2 is the kinetic energy and U (u) is the potential energy expressed by
U (u) = \Gamma
u
Z
0
F (s) ds = 1
b
ln exp(bu) + a
1 + a
\Gamma
u 2
2 (6)
1 Benem'erita Universidad Aut'onoma de Puebla, Facultad de Ciencias F'isico Matem'aticas, Apartado Postal
1152, Puebla, M'exico, e­mail: fraguela@fcfm.buap.mx, leyva@fcfm.buap.mx
2 Research Computing Center, Moscow State University, Moscow, 119899, Russian Federation, e­mail:
greben@iname.ru

136 Mathematical Modeling of Complex Information Processing Systems
The stationary states of the above Hamiltonian system are critical points that, in their turn, correspond
to the extreme points of the function U (u) [3]. Hence, the critical points can be determined from the equation
U 0 (u) = 0 that leads to g(u; a; b) = u, i.e., to the equation
1
1 + a exp(\Gammabu) = u (7)
The following cases of graphical representations of solutions to equation (7) are given in Figure 1, A for
different values of a and b: 1) the solid line corresponds to the case when a single critical point exists; 2) the
dashed line illustrates the case of two critical points, and 3) the broken line corresponds to the case of three
critical points.
Let us determine the domains on a plane (b; a) where the above cases are realized. The case of two critical
points follows from equation (7) and the condition of tangency g 0 (u) = 1:
1
1 + a exp(\Gammabu) = u
ab exp(\Gammabu)
\Gamma
1 + a exp(\Gammabu)
\Delta 2
= 1
(8)
Solving this system, we obtain the equation
bu 2 \Gamma bu + 1 = 0 (9)
whose solution yields the tangency points
u 0
\Sigma = 1
2 \Sigma
1
2
r
1 \Gamma
4
b
(10)
Since D = 1 \Gamma
4
b
should be positive, we get the inequality b ? 4. Substituting (10) into the first equation in (8),
we come to the following dependence between a and b:
a \Sigma =
/
b
2
/
1 \Upsilon
r
1 \Gamma
4
b
!
\Gamma 1
!
exp
/
b
2
/
1 \Sigma
r
1 \Gamma
4
b
!!
This dependence is represented in Figure 1, B. The curves corresponding to a+ and a \Gamma divide the plane
into three
domains\Omega k (k = 1; 2; 3) as follows: 1) the
domains\Omega 1
and\Omega 2 contain the values of a and b for
the case of a single critical point, 2) their boundaries correspond to the case of two critical points, and 3) the
domain\Omega 3 presents the case of three critical points.
2.2. Conditions for stability of the critical points. The center­type stable critical points correspond
to a minimum of the function U (u) [3]. We can see from Figure 1, B that these points are situated only in the
domain\Omega 3 (i.e., between the curves a+ and a \Gamma ). The function U (u) has three extrema when (b; a)
2\Omega 3 . Only
one of these extrema is a local minimum. Note that this minimum corresponds to the center­type case and
satisfies the inequality U 00 (u) = \GammaF 0 (u) ? 0, which, together with equation (7), yields the condition for the
existence of stable critical points:
1
1 + a exp(\Gammabu) = u
ab exp(\Gammabu)
\Gamma
1 + a exp(\Gammabu)
\Delta 2 ? 1
(11)
Solving (11), we obtain the inequality
bu 2 \Gamma bu + 1 ! 0 (12)
that is valid for b ? 4 and u 0
\Gamma ! u ! u 0
+ , where u 0
\Sigma are determined by (10). Expressing the above conditions
of stability in terms of the parameters a and b, we come to the conclusion
that\Omega 3 is the domain of stability.
Thus, we can formulate the following theorem.
Theorem. Ordinary differential equation (3) has a center­type critical point u 0 corresponding to given
values b 0 and a 0 of the parameters b and a if and only if the conditions (b 0 ; a 0 )
2\Omega 3 and u 0
\Gamma ! u 0 ! u 0
+ hold.
Here u 0
\Sigma are determined by expression (10).

Mathematical Modeling of Complex Information Processing Systems 137
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
A
u
g(u,a,b)
4 5 6 7
0
10
20
30
40
50
60
70
b
a(b)
B
W2
W1
W3
a -
a +
Fig. 1. Analysis of critical points. A: the number of critical points for different values of a and b.
B: the domains with the different number of critical points
3. Numerical results. Let us assume that u 0 is a center­type critical point of equation (3). Taken its
solution in the form u = u + u 0 , from (3) we obtain the equation
d 2 u
dy 2 +
`
1
1 + a exp(\Gammabu) exp(\Gammabu 0 ) \Gamma
1
1 + a exp(\Gammabu 0 )
'
\Gamma u = 0 (13)
Now we consider the boundary conditions u (0) = u (l), u 0 (0) = u 0 (l) (where l is an unknown period) and
construct the corresponding periodic solutions to equation (13). Note that near the center­type critical point
all such solutions can be represented by closed (periodic) trajectories on the phase plane. Hence, the boundary
value problem can be reduced to the Cauchy problem for equation (13) with the initial conditions
u(0) = u 0 ; u 0 (0) = u 0
0 (14)
where u 0 and u 0
0 are some given values. In order to solve problem (13), (14) numerically, we introduce a uniform
grid fy i g N
i=1 with a step length h and approximate the derivatives by the corresponding difference formulas.
As a result, we obtain a system of nonlinear algebraic equations, which can be solved by the least­squares
method. The above algorithm was realized with the MATLAB system, and some model examples were solved
numerically.
The periodic solutions obtained near the center u 0 = 0:472 are represented in Figure 2 for the following
cases: b = 6, a = 19 (b; a
2\Omega 3 ), u 0 = \Gamma0:01, u 0
0 = 0:01 (the solid line), and u 0 = \Gamma0:03, u 0
0 = 0:01 (the dashed

138 Mathematical Modeling of Complex Information Processing Systems
0.44 0.45 0.46 0.47 0.48 0.49 0.5 0.51
-0.025
-0.02
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
0.02
0.025
u
du/dy
Fig. 2. The behavior of periodic solutions near the center
line). The numerical experiments we performed confirm our theoretical conclusions. The above results can be
useful for further studies of basic rhythms of the brain cortex.
REFERENCES
1. P.L. N~unez, Neocortical Dynamics and Human EEG Rhythms, Oxford, 1995.
2. A.C. Fraguela and J.A. Escamilla, ``A mathematical model for studying the activation processes in the
cerebral cortex'', in: Numerical Analysis: Theory, Applications, Programs, pp. 75--83, Moscow, 1999.
3. L. Perko, Differential Equations and Dynamical Systems, New York, 1996.