- Email: [email protected]
Mechanisms for activity spread in a neural field model
ARTICLE IN PRESS
Neurocomputing 70 (2007) 2996–3000 www.elsevier.com/locate/neucom
Mechanisms for activity spread in a neural field model Mihaela Enculescu, Michael Bestehorn Brandenburgische Technische Universitaet Cottbus, Lehrstuhl fuer Theorestische Physik II, Erich-Weiner-Strasse 1, 03046 Cottbus, Germany Received 9 March 2006; received in revised form 12 September 2006; accepted 13 September 2006 Communicated by V. Jirsa Available online 11 November 2006
Abstract Field models of continuous neural networks incorporate nonlocal connectivities as well as finite axonal propagation velocities and lead therefore to delayed integral equations. For special choices of the synaptic footprint it is possible to reduce the integral model to a system of partial differential equations. One example is that of the inhomogeneous damped wave equation in one space dimension derived by Jirsa and Haken for exponential synaptic footprint. We show that this equation can be put into the form of a conservation law with nonlinear source, and explore numerically this representation. We find two mechanisms for the spread of the activity from an initially excited region. r 2006 Elsevier B.V. All rights reserved. Keywords: Neural field models; Traveling waves
A mathematical model of a continuously distributed neural population was first proposed 30 years ago in Wilson and Cowan [17]. They described the activity of the cortical and thalamic nervous tissue by means of the local axonal pulse rates of excitatory and inhibitory neurons, for which they derived a system of coupled time evolution equations. The time relaxation of the firing rate as well as the nonlinear conversion of incoming pre-synaptic pulses to post-synaptic firing rates were considered. The coupling between neurons was assumed to be dependent on the distance, fact that gives their model the mathematical form of nonlinear integral equations. Since then, a lot of continuum models for the electrical activity of neural tissue have been proposed and studied intensively. Interesting results concerning several aspects like spatio–temporal pattern formation and traveling wave solutions were achieved. Localized solutions with simple and multiple bumps as mechanisms for working memory were analyzed in [1,2,13]. Bifurcation from trivial equilibria and pattern formation in neural field models were studied in the context of visual hallucinations, for example in [6,7,5,15,16,10]. Also the propagation of excitation through
the neural tissue can be modeled by traveling wave solutions of equations of this type (e.g. in [1,2,8,4,3]). For an autonomous, single-layer neural network, continuously distributed over the real line, the local depolarization uðx; tÞ can be written as integral over all post-synaptic currents in the form Z t Z 1 uðx; tÞ ¼ Zðt sÞ wðx yÞ 1 1 jx yj ds dy, ð1Þ f u y; s v where f ðuÞ denotes the firing rate function. The synaptic kernel Zðt sÞ models the post-synaptic current initiated at time t by a single action potential arriving at time s, while the synaptic footprint wðx yÞ stands for the synaptic connectivity strength between neurons at x and y. The delay term in (1) arises from the finite conduction velocity v of action potentials along axons. The firing rate function is monotonically increasing and saturates for large values of u. A reasonable fit to experimental data is considered to be f ðuÞ ¼ ð1 þ expðbðu hÞÞÞ1 ,
Corresponding author.
E-mail address: [email protected] (M. Enculescu). 0925-2312/$ - see front matter r 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.neucom.2006.09.012
were h is a positive threshold, and b a positive steepness parameter.
ARTICLE IN PRESS M. Enculescu, M. Bestehorn / Neurocomputing 70 (2007) 2996–3000
2997
The synaptic kernel ZðtÞ is chosen for simplicity to be the Green’s function of a linear differential operator L. A common choice is the one-sided exponential function:
The Fourier transform of (5) reads
ZðtÞ ¼ a expðatÞYðtÞ with L ¼ 1 þ a1 qt ,
For exponential synaptic footprint, the Fourier transform ~ oÞ is given by Gðk;
(2)
where YðtÞ denotes the Heaviside step function. The synaptic footprint is usually assumed to be a symmetric function wðjxjÞ. For one-dimensional systems with delays, the reformulation of (1) in terms of a PDE has been demonstrated in [2] for special functional forms of wðxÞ. With the notations Cðx; tÞ:¼Luðx; tÞ
and
rðx; tÞ:¼f ðuðx; tÞÞ,
(3)
and for exponential synaptic footprint wðxÞ ¼ expðjxj=sÞ=ð2sÞ,
(4)
we have
1
0
Last equation can be written in the form Z 1Z 1 Gðx x0 ; t t0 Þrðx0 ; t0 Þ dx0 dt0 , Cðx; tÞ ¼ 1
(5)
1
with the Green’s function x Gðx; tÞ ¼ d t wðxÞ. v
o 1þi o 0 ~ oÞ ¼ , Gðk; o 2 2 2 1þi þs k o0 with o0 ¼ v=s. This allow us to write (6) in the form " # o 2 o 2 2 ~ ~ oÞ, þ s k cðk; oÞ ¼ 1 þ i 1þi rðk; o0 o0 and after performing the inverse Fourier transform we get
Ct vjx ¼ o0 C þ o0 r;
u
u
1 t=0
0.5
u
u
t = 13.5
0.5
t = 27
0.5
u
u
t = 27
0.5 0 1
t = 40.5
0.5
u
u
t = 13.5
0.5 0 1
0 1 0 1
t = 40.5
0.5 0 1
t = 54
0.5
u
u
t=0
0.5 0 1
0 1
0 1
t = 54
0.5 0 1
t = 67.5
0.5
u
u
(8)
jt vCx ¼ o0 j;
0 1
0 1
t = 67.5
0.5 0 1
t = 81
0.5
u
u
(7)
This equation has first been derived in [11] for Zðt sÞ ¼ dðt sÞ, that is a special case of (2) when a ! 1. Their model has been successfully applied to give a theoretical description to the movement control experiments in [12]. With an auxiliary function jðx; tÞ, (7) can be transformed into the system
1
0 1
t = 81
0.5 0 1
t = 94.5
0.5
u
u
(6)
Ctt v2 Cxx þ 2o0 Ct þ o20 C ¼ o20 r þ o0 rt .
jx x j Cðx; tÞ ¼ wðx x0 Þr x0 ; t dx0 . v 1 Z
~ oÞrðk; ~ ~ oÞ. Cðk; oÞ ¼ Gðk;
0
t = 94.5
0.5 0
0
5
10
15 x
20
25
30
0
5
10
15 x
20
25
30
Fig. 1. Numerical solution of (8) and (3) for v ¼ 0:1, h ¼ 0:35, b ¼ 1 and Z ¼ a eat YðtÞ. The time unit is 1=a, x is measured in s units. On the left: the width of the initially excited interval is 0:672, two fronts are created. On the right: the width of the initially excited interval is 0:67, two pulses are produced.
ARTICLE IN PRESS M. Enculescu, M. Bestehorn / Neurocomputing 70 (2007) 2996–3000
u u
0 1 u
u u
t = 43.2
0.5
0 1 0 1
t = 75.6
0.5 0
5
10
15
20
25
0 1 0.5
t = 32.4
t = 43.2
t = 54
0.5 0 1
t = 64.8
0.5
0.5
0 1
t = 54
0.5
t = 21.6
0.5 0 1
t = 32.4
0.5 0 1
u
u
t = 21.6
0.5 0 1
u
0 1
u
u
0 1
t = 10.8
0.5
u
0.5
t=0
0.5 0 1
t = 10.8
u
u
0 1
0
1
t=0
u
u
1 0.5
0.5
u
2998
0 1 0.5
t = 64.8
t = 75.6
0
30
0
5
10
15
20
25
30
x Fig. 2. Numerical solution of (8) and (3) for v ¼ 0:1, h ¼ 0:35, b ¼ 20 and Z ¼ a eat YðtÞ. On the left: the width of the initially excited interval is 0:71, two fronts are created. On the right: the width of the initially excited interval is 0:708, two pulses are produced.
0.5
u
t = 40.5
0.5 0 1
t = 54
u
0.5 0 1
t = 67.5
t = 67.5
0.5
u
u u
0 1
t = 81
0.5
t = 81
0.5
u
u
0 1
t = 54
0 1
0 1
t = 94.5
0.5
t = 94.5
0.5
u
u
t = 27
0.5
t = 40.5
0 1
0
0 1
t = 27
0 1 0 1 0.5
t = 13.5
0.5
u
u u
0 1
0.5
0 1
t = 13.5
0.5
0 1 0.5
t=0
0.5
u
u
0 1
0.5
1
t=0
u
u
1
0
5
10
15 x
20
25
30
0
0
5
10
15
20
25
30
x
Fig. 3. Numerical solution of (8) and (3) for v ¼ 0:1, h ¼ 0:35, b ¼ 1 and Z ¼ a2 t eat YðtÞ. On the left: the width of the initially excited interval is 0:684, two fronts are created. On the right: the width of the initially excited interval is 0:682, two pulses are produced.
ARTICLE IN PRESS M. Enculescu, M. Bestehorn / Neurocomputing 70 (2007) 2996–3000 1 u
t =0 0.5 0 1 t = 2.25 u
that exhibits the form of a linear conservation law with a nonlinear source term. We have simulated (8) by splitting the source term and alternatively solving the homogeneous linear part using Godunov’s Method, and the nonlinear part by an Euler step. The source term has to be evaluated by solving (3). In order to achieve numerical stability, the discretization must satisfy the CFL condition (see [14])
2999
vDt p1. Dx
0.5 0 1
1
u
t = 4.5 0.5 0 1 u
t = 6.75 0.5 0 1 u
t=9 0.5 0 1 u
t = 11.25 0.5 0 1 u
t = 13.5 0.5 0 1 t = 15.75 u
We analyzed the dynamics of the field u starting with an excited bounded interval of variable length D. Depending on D, we found two mechanisms for the spread of activity from the initially excited interval. When the width of this interval is big enough, the excited area grows by the formation of two fronts, that travel without changing shape. If the initially excited area is smaller, two pulses are created, propagate and undergo attenuation. Fig. 1 shows the results of a numerical simulation in the special case b ! 1, when the firing rate function is a Heaviside step function. For finite values of b we obtain qualitatively the same picture, but the attenuation of the pulses increases (see Fig. 2). Also the time kernel influences the propagation of pulses. For example, the choice Z ¼ a2 t expðatÞYðtÞ also leads to stronger attenuation (see Fig. 3). Interestingly, formation of pulses occurs only for small values of the axonal propagation velocity v, so that no pulses are observed when the propagation delay is neglected. This agrees to the analytical result in [3]. For bigger values of v either two fronts are created, or the solution in the initially excited region monotonically falls down into the ground state. Thus, in this last case, the excitation remains localized and does not spread to the neighboring regions. Fig. 4 shows the simulation results in this case. We remark that this behavior can be energetically explained in the limit case v ! 1. As shown in [9], a free energy functional can be found in this limit. This is given by Z Z 1 1 1 F ½u ¼ wðx yÞff ½uðx; tÞ 4 1 1 Z 1 f ½uðy; tÞg2 dx dy þ W ½uðx; tÞ dx
0.5 0 0
5
10
15
20
25
30
x
Fig. 4. Numerical solution of (8) and (3) for v ¼ 0:6, h ¼ 0:35, b ¼ 1 and Z ¼ a eat YðtÞ. The system reaches the homogeneous state u 0 without the creation of traveling pulses.
leads to F 0 ½u ¼ sðuþ u Þ2 ð1 eD=s Þ þ 2RW ðu Þ þ DðW ðuþ Þ W ðu ÞÞ.
with Z
u 0
f ðsÞ½s f ðsÞ ds.
W ½u ¼ 0
Evaluation of the energy functional for our initial conditions ( uþ ; for x 2 ½0; D; uðx; 0Þ ¼ u ; for x 2 ½R; R ½0; D
The derivative dF 0 ½u 1 ¼ ðuþ u Þ2 eD=s þ W ðuþ Þ W ðu Þ dD 2 is a monotonic decreasing function changing sign at ðuþ u Þ2 Dcr ¼ s ln . 2½W ðuþ Þ W ðu Þ Thus, a growth of D is energetically favorable only for DXDcr . Therefore, activity spread from the initially excited
ARTICLE IN PRESS M. Enculescu, M. Bestehorn / Neurocomputing 70 (2007) 2996–3000
3000
u
1
necessary when the model is reformulated in terms of the PDE system (8), that is local in time and space.
t=0
0.5
References
u
0 1 t = 4.5
0.5
u
0 1 t=9
0.5
u
0 1 t = 13.5
0.5
u
0 1 t = 18
0.5
u
0 1 t = 22.5
0.5
u
0 1 t = 27
0.5
u
0 1 t = 31.5
0.5 0
0
5
10
15
20
25
30
x
Fig. 5. Numerical solution of (8) and (3) for v ¼ 0:6, h ¼ 0:5, b ¼ 1 and Z ¼ a eat YðtÞ. The system reaches a stationary inhomogeneous state. The length of the excited interval depends on the initial conditions.
interval is possible only if the width of the interval exceeds a threshold value. In all simulations, the field approaches either the homogeneous resting state u 0 or the homogeneous excited state uþ 1 for t ! 1. These states are the fixed points of (8). Thus, the width of the excited interval in the initial field decides about which fixed point is chosen for t ! 1. The only exception from this behavior is the case when f is symmetric in respect to the two fixed points (h ¼ 0:5). In this case, a stationary solution is reached (see Fig. 5). Our method has several advantages, compared to a direct integration scheme, e.g. like the one presented in [10]. For example, the costly evaluation of the integral is no more necessary. Furthermore, the direct simulation of (1) requires the storage of the solution for previous times, due to the delay term under the integral. This is no more
[1] S. Amari, Dynamics of pattern formation in lateral-inhibition type neural fields, Biol. Cybern. 27 (1977) 77–87. [2] S. Coombes, G.L. Lord, M.R. Owen, Waves and bumps in neuronal networks with axo-dendritic synaptic interactions, Physica D 178 (2003) 219–241. [3] M. Enculescu, A note on traveling fronts and pulses in a firing rate model of a neuronal network, Physica D 196 (2004) 362–386. [4] M. Enculescu, M. Bestehorn, Activity dynamics in nonlocal interacting neural fields, Phys. Rev. E 67 (2003) 041904. [5] G.B. Ermentrout, Neural networks as spatio–temporal pattern formatting systems, Rep. Prog. Phys. 61 (1998) 353–430. [6] G.B. Ermentrout, J.D. Cowan, A mathematical theory of visual hallucination patterns, Biol. Cybern. 34 (1979) 137–150. [7] G.B. Ermentrout, M. Lewis, Pattern formation in systems with one spatially distributed species, Bull. Math. Biol. 59 (1997) 533–549. [8] G.B. Ermentrout, J.B. McLeod, Existence and uniqueness of traveling waves for a neural network, Proc. R. Soc. Edinburgh 123A (1993) 461–478. [9] D.A. French, Identification of a free energy functional in an integrodifferential equation model for neuronal network activity, Appl. Math. Lett. 17 (2004) 1047–1051. [10] A. Hutt, M. Bestehorn, T. Wennekers, Pattern formation in intracortical neuronal fields, Network 14 (2003) 351–368. [11] V.K. Jirsa, H. Haken, Field theory of electromagnetic brain activity, Phys. Rev. Lett. 77 (1996) 960–963. [12] A.S. Kelso, S.L. Bressler, G.C. DeGuzman, M. Ding, A. Fuchs, T. Holroyd, A phase transition in human brain and behavior, Phys. Lett. A 168 (1992) 134–144. [13] C.R. Laing, W.C. Gutkin, G.B. Ermentrout, Multiple bumps in a neuronal model of working memory, SIAM J. Appl. Math. 63 (2002) 62–97. [14] R. LeVeque, Numerical Methods for Conservation Laws, Birkha¨user-Verlag, Basel, 1992. [15] P. Tass, Cortical pattern formation during visual hallucinations, J. Biol. Phys. 21 (1995) 177–210. [16] P. Tass, Oscillatory cortical activity during visual hallucinations, J. Biol. Phys. 23 (1997) 21–66. [17] H.R. Wilson, J.D. Cowan, A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue, Kybernetik 13 (1973) 55–58. Mihaela Enculescu was born in Bucharest, Romania, in June 1979. She graduated in Physics from the Technical University of Cottbus in 2002. She is currently pursuing her Ph.D. degree and her research fields are pattern formation, neural networks and modeling of complex nonlinear systems.
Michael Bestehorn obtained his Diploma in Physics in 1983 at the University of Stuttgart. After his Ph.D. in Physics 1988 he spend one year at the University of Pamplona, Spain. He received his Habilitation 1995 at the Institute for Theoretical Physics and Synergetics in Stuttgart. Since 1999 he is Professor for Theoretical Physics at the University of Cottbus. His research interests are pattern formation, fluid mechanics, reduced models, and nonlinear systems.
Neurocomputing 70 (2007) 2996–3000 www.elsevier.com/locate/neucom
Mechanisms for activity spread in a neural field model Mihaela Enculescu, Michael Bestehorn Brandenburgische Technische Universitaet Cottbus, Lehrstuhl fuer Theorestische Physik II, Erich-Weiner-Strasse 1, 03046 Cottbus, Germany Received 9 March 2006; received in revised form 12 September 2006; accepted 13 September 2006 Communicated by V. Jirsa Available online 11 November 2006
Abstract Field models of continuous neural networks incorporate nonlocal connectivities as well as finite axonal propagation velocities and lead therefore to delayed integral equations. For special choices of the synaptic footprint it is possible to reduce the integral model to a system of partial differential equations. One example is that of the inhomogeneous damped wave equation in one space dimension derived by Jirsa and Haken for exponential synaptic footprint. We show that this equation can be put into the form of a conservation law with nonlinear source, and explore numerically this representation. We find two mechanisms for the spread of the activity from an initially excited region. r 2006 Elsevier B.V. All rights reserved. Keywords: Neural field models; Traveling waves
A mathematical model of a continuously distributed neural population was first proposed 30 years ago in Wilson and Cowan [17]. They described the activity of the cortical and thalamic nervous tissue by means of the local axonal pulse rates of excitatory and inhibitory neurons, for which they derived a system of coupled time evolution equations. The time relaxation of the firing rate as well as the nonlinear conversion of incoming pre-synaptic pulses to post-synaptic firing rates were considered. The coupling between neurons was assumed to be dependent on the distance, fact that gives their model the mathematical form of nonlinear integral equations. Since then, a lot of continuum models for the electrical activity of neural tissue have been proposed and studied intensively. Interesting results concerning several aspects like spatio–temporal pattern formation and traveling wave solutions were achieved. Localized solutions with simple and multiple bumps as mechanisms for working memory were analyzed in [1,2,13]. Bifurcation from trivial equilibria and pattern formation in neural field models were studied in the context of visual hallucinations, for example in [6,7,5,15,16,10]. Also the propagation of excitation through
the neural tissue can be modeled by traveling wave solutions of equations of this type (e.g. in [1,2,8,4,3]). For an autonomous, single-layer neural network, continuously distributed over the real line, the local depolarization uðx; tÞ can be written as integral over all post-synaptic currents in the form Z t Z 1 uðx; tÞ ¼ Zðt sÞ wðx yÞ 1 1 jx yj ds dy, ð1Þ f u y; s v where f ðuÞ denotes the firing rate function. The synaptic kernel Zðt sÞ models the post-synaptic current initiated at time t by a single action potential arriving at time s, while the synaptic footprint wðx yÞ stands for the synaptic connectivity strength between neurons at x and y. The delay term in (1) arises from the finite conduction velocity v of action potentials along axons. The firing rate function is monotonically increasing and saturates for large values of u. A reasonable fit to experimental data is considered to be f ðuÞ ¼ ð1 þ expðbðu hÞÞÞ1 ,
Corresponding author.
E-mail address: [email protected] (M. Enculescu). 0925-2312/$ - see front matter r 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.neucom.2006.09.012
were h is a positive threshold, and b a positive steepness parameter.
ARTICLE IN PRESS M. Enculescu, M. Bestehorn / Neurocomputing 70 (2007) 2996–3000
2997
The synaptic kernel ZðtÞ is chosen for simplicity to be the Green’s function of a linear differential operator L. A common choice is the one-sided exponential function:
The Fourier transform of (5) reads
ZðtÞ ¼ a expðatÞYðtÞ with L ¼ 1 þ a1 qt ,
For exponential synaptic footprint, the Fourier transform ~ oÞ is given by Gðk;
(2)
where YðtÞ denotes the Heaviside step function. The synaptic footprint is usually assumed to be a symmetric function wðjxjÞ. For one-dimensional systems with delays, the reformulation of (1) in terms of a PDE has been demonstrated in [2] for special functional forms of wðxÞ. With the notations Cðx; tÞ:¼Luðx; tÞ
and
rðx; tÞ:¼f ðuðx; tÞÞ,
(3)
and for exponential synaptic footprint wðxÞ ¼ expðjxj=sÞ=ð2sÞ,
(4)
we have
1
0
Last equation can be written in the form Z 1Z 1 Gðx x0 ; t t0 Þrðx0 ; t0 Þ dx0 dt0 , Cðx; tÞ ¼ 1
(5)
1
with the Green’s function x Gðx; tÞ ¼ d t wðxÞ. v
o 1þi o 0 ~ oÞ ¼ , Gðk; o 2 2 2 1þi þs k o0 with o0 ¼ v=s. This allow us to write (6) in the form " # o 2 o 2 2 ~ ~ oÞ, þ s k cðk; oÞ ¼ 1 þ i 1þi rðk; o0 o0 and after performing the inverse Fourier transform we get
Ct vjx ¼ o0 C þ o0 r;
u
u
1 t=0
0.5
u
u
t = 13.5
0.5
t = 27
0.5
u
u
t = 27
0.5 0 1
t = 40.5
0.5
u
u
t = 13.5
0.5 0 1
0 1 0 1
t = 40.5
0.5 0 1
t = 54
0.5
u
u
t=0
0.5 0 1
0 1
0 1
t = 54
0.5 0 1
t = 67.5
0.5
u
u
(8)
jt vCx ¼ o0 j;
0 1
0 1
t = 67.5
0.5 0 1
t = 81
0.5
u
u
(7)
This equation has first been derived in [11] for Zðt sÞ ¼ dðt sÞ, that is a special case of (2) when a ! 1. Their model has been successfully applied to give a theoretical description to the movement control experiments in [12]. With an auxiliary function jðx; tÞ, (7) can be transformed into the system
1
0 1
t = 81
0.5 0 1
t = 94.5
0.5
u
u
(6)
Ctt v2 Cxx þ 2o0 Ct þ o20 C ¼ o20 r þ o0 rt .
jx x j Cðx; tÞ ¼ wðx x0 Þr x0 ; t dx0 . v 1 Z
~ oÞrðk; ~ ~ oÞ. Cðk; oÞ ¼ Gðk;
0
t = 94.5
0.5 0
0
5
10
15 x
20
25
30
0
5
10
15 x
20
25
30
Fig. 1. Numerical solution of (8) and (3) for v ¼ 0:1, h ¼ 0:35, b ¼ 1 and Z ¼ a eat YðtÞ. The time unit is 1=a, x is measured in s units. On the left: the width of the initially excited interval is 0:672, two fronts are created. On the right: the width of the initially excited interval is 0:67, two pulses are produced.
ARTICLE IN PRESS M. Enculescu, M. Bestehorn / Neurocomputing 70 (2007) 2996–3000
u u
0 1 u
u u
t = 43.2
0.5
0 1 0 1
t = 75.6
0.5 0
5
10
15
20
25
0 1 0.5
t = 32.4
t = 43.2
t = 54
0.5 0 1
t = 64.8
0.5
0.5
0 1
t = 54
0.5
t = 21.6
0.5 0 1
t = 32.4
0.5 0 1
u
u
t = 21.6
0.5 0 1
u
0 1
u
u
0 1
t = 10.8
0.5
u
0.5
t=0
0.5 0 1
t = 10.8
u
u
0 1
0
1
t=0
u
u
1 0.5
0.5
u
2998
0 1 0.5
t = 64.8
t = 75.6
0
30
0
5
10
15
20
25
30
x Fig. 2. Numerical solution of (8) and (3) for v ¼ 0:1, h ¼ 0:35, b ¼ 20 and Z ¼ a eat YðtÞ. On the left: the width of the initially excited interval is 0:71, two fronts are created. On the right: the width of the initially excited interval is 0:708, two pulses are produced.
0.5
u
t = 40.5
0.5 0 1
t = 54
u
0.5 0 1
t = 67.5
t = 67.5
0.5
u
u u
0 1
t = 81
0.5
t = 81
0.5
u
u
0 1
t = 54
0 1
0 1
t = 94.5
0.5
t = 94.5
0.5
u
u
t = 27
0.5
t = 40.5
0 1
0
0 1
t = 27
0 1 0 1 0.5
t = 13.5
0.5
u
u u
0 1
0.5
0 1
t = 13.5
0.5
0 1 0.5
t=0
0.5
u
u
0 1
0.5
1
t=0
u
u
1
0
5
10
15 x
20
25
30
0
0
5
10
15
20
25
30
x
Fig. 3. Numerical solution of (8) and (3) for v ¼ 0:1, h ¼ 0:35, b ¼ 1 and Z ¼ a2 t eat YðtÞ. On the left: the width of the initially excited interval is 0:684, two fronts are created. On the right: the width of the initially excited interval is 0:682, two pulses are produced.
ARTICLE IN PRESS M. Enculescu, M. Bestehorn / Neurocomputing 70 (2007) 2996–3000 1 u
t =0 0.5 0 1 t = 2.25 u
that exhibits the form of a linear conservation law with a nonlinear source term. We have simulated (8) by splitting the source term and alternatively solving the homogeneous linear part using Godunov’s Method, and the nonlinear part by an Euler step. The source term has to be evaluated by solving (3). In order to achieve numerical stability, the discretization must satisfy the CFL condition (see [14])
2999
vDt p1. Dx
0.5 0 1
1
u
t = 4.5 0.5 0 1 u
t = 6.75 0.5 0 1 u
t=9 0.5 0 1 u
t = 11.25 0.5 0 1 u
t = 13.5 0.5 0 1 t = 15.75 u
We analyzed the dynamics of the field u starting with an excited bounded interval of variable length D. Depending on D, we found two mechanisms for the spread of activity from the initially excited interval. When the width of this interval is big enough, the excited area grows by the formation of two fronts, that travel without changing shape. If the initially excited area is smaller, two pulses are created, propagate and undergo attenuation. Fig. 1 shows the results of a numerical simulation in the special case b ! 1, when the firing rate function is a Heaviside step function. For finite values of b we obtain qualitatively the same picture, but the attenuation of the pulses increases (see Fig. 2). Also the time kernel influences the propagation of pulses. For example, the choice Z ¼ a2 t expðatÞYðtÞ also leads to stronger attenuation (see Fig. 3). Interestingly, formation of pulses occurs only for small values of the axonal propagation velocity v, so that no pulses are observed when the propagation delay is neglected. This agrees to the analytical result in [3]. For bigger values of v either two fronts are created, or the solution in the initially excited region monotonically falls down into the ground state. Thus, in this last case, the excitation remains localized and does not spread to the neighboring regions. Fig. 4 shows the simulation results in this case. We remark that this behavior can be energetically explained in the limit case v ! 1. As shown in [9], a free energy functional can be found in this limit. This is given by Z Z 1 1 1 F ½u ¼ wðx yÞff ½uðx; tÞ 4 1 1 Z 1 f ½uðy; tÞg2 dx dy þ W ½uðx; tÞ dx
0.5 0 0
5
10
15
20
25
30
x
Fig. 4. Numerical solution of (8) and (3) for v ¼ 0:6, h ¼ 0:35, b ¼ 1 and Z ¼ a eat YðtÞ. The system reaches the homogeneous state u 0 without the creation of traveling pulses.
leads to F 0 ½u ¼ sðuþ u Þ2 ð1 eD=s Þ þ 2RW ðu Þ þ DðW ðuþ Þ W ðu ÞÞ.
with Z
u 0
f ðsÞ½s f ðsÞ ds.
W ½u ¼ 0
Evaluation of the energy functional for our initial conditions ( uþ ; for x 2 ½0; D; uðx; 0Þ ¼ u ; for x 2 ½R; R ½0; D
The derivative dF 0 ½u 1 ¼ ðuþ u Þ2 eD=s þ W ðuþ Þ W ðu Þ dD 2 is a monotonic decreasing function changing sign at ðuþ u Þ2 Dcr ¼ s ln . 2½W ðuþ Þ W ðu Þ Thus, a growth of D is energetically favorable only for DXDcr . Therefore, activity spread from the initially excited
ARTICLE IN PRESS M. Enculescu, M. Bestehorn / Neurocomputing 70 (2007) 2996–3000
3000
u
1
necessary when the model is reformulated in terms of the PDE system (8), that is local in time and space.
t=0
0.5
References
u
0 1 t = 4.5
0.5
u
0 1 t=9
0.5
u
0 1 t = 13.5
0.5
u
0 1 t = 18
0.5
u
0 1 t = 22.5
0.5
u
0 1 t = 27
0.5
u
0 1 t = 31.5
0.5 0
0
5
10
15
20
25
30
x
Fig. 5. Numerical solution of (8) and (3) for v ¼ 0:6, h ¼ 0:5, b ¼ 1 and Z ¼ a eat YðtÞ. The system reaches a stationary inhomogeneous state. The length of the excited interval depends on the initial conditions.
interval is possible only if the width of the interval exceeds a threshold value. In all simulations, the field approaches either the homogeneous resting state u 0 or the homogeneous excited state uþ 1 for t ! 1. These states are the fixed points of (8). Thus, the width of the excited interval in the initial field decides about which fixed point is chosen for t ! 1. The only exception from this behavior is the case when f is symmetric in respect to the two fixed points (h ¼ 0:5). In this case, a stationary solution is reached (see Fig. 5). Our method has several advantages, compared to a direct integration scheme, e.g. like the one presented in [10]. For example, the costly evaluation of the integral is no more necessary. Furthermore, the direct simulation of (1) requires the storage of the solution for previous times, due to the delay term under the integral. This is no more
[1] S. Amari, Dynamics of pattern formation in lateral-inhibition type neural fields, Biol. Cybern. 27 (1977) 77–87. [2] S. Coombes, G.L. Lord, M.R. Owen, Waves and bumps in neuronal networks with axo-dendritic synaptic interactions, Physica D 178 (2003) 219–241. [3] M. Enculescu, A note on traveling fronts and pulses in a firing rate model of a neuronal network, Physica D 196 (2004) 362–386. [4] M. Enculescu, M. Bestehorn, Activity dynamics in nonlocal interacting neural fields, Phys. Rev. E 67 (2003) 041904. [5] G.B. Ermentrout, Neural networks as spatio–temporal pattern formatting systems, Rep. Prog. Phys. 61 (1998) 353–430. [6] G.B. Ermentrout, J.D. Cowan, A mathematical theory of visual hallucination patterns, Biol. Cybern. 34 (1979) 137–150. [7] G.B. Ermentrout, M. Lewis, Pattern formation in systems with one spatially distributed species, Bull. Math. Biol. 59 (1997) 533–549. [8] G.B. Ermentrout, J.B. McLeod, Existence and uniqueness of traveling waves for a neural network, Proc. R. Soc. Edinburgh 123A (1993) 461–478. [9] D.A. French, Identification of a free energy functional in an integrodifferential equation model for neuronal network activity, Appl. Math. Lett. 17 (2004) 1047–1051. [10] A. Hutt, M. Bestehorn, T. Wennekers, Pattern formation in intracortical neuronal fields, Network 14 (2003) 351–368. [11] V.K. Jirsa, H. Haken, Field theory of electromagnetic brain activity, Phys. Rev. Lett. 77 (1996) 960–963. [12] A.S. Kelso, S.L. Bressler, G.C. DeGuzman, M. Ding, A. Fuchs, T. Holroyd, A phase transition in human brain and behavior, Phys. Lett. A 168 (1992) 134–144. [13] C.R. Laing, W.C. Gutkin, G.B. Ermentrout, Multiple bumps in a neuronal model of working memory, SIAM J. Appl. Math. 63 (2002) 62–97. [14] R. LeVeque, Numerical Methods for Conservation Laws, Birkha¨user-Verlag, Basel, 1992. [15] P. Tass, Cortical pattern formation during visual hallucinations, J. Biol. Phys. 21 (1995) 177–210. [16] P. Tass, Oscillatory cortical activity during visual hallucinations, J. Biol. Phys. 23 (1997) 21–66. [17] H.R. Wilson, J.D. Cowan, A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue, Kybernetik 13 (1973) 55–58. Mihaela Enculescu was born in Bucharest, Romania, in June 1979. She graduated in Physics from the Technical University of Cottbus in 2002. She is currently pursuing her Ph.D. degree and her research fields are pattern formation, neural networks and modeling of complex nonlinear systems.
Michael Bestehorn obtained his Diploma in Physics in 1983 at the University of Stuttgart. After his Ph.D. in Physics 1988 he spend one year at the University of Pamplona, Spain. He received his Habilitation 1995 at the Institute for Theoretical Physics and Synergetics in Stuttgart. Since 1999 he is Professor for Theoretical Physics at the University of Cottbus. His research interests are pattern formation, fluid mechanics, reduced models, and nonlinear systems.