Chaos in neural systems

Volume 114A, number 7

PHYSICS LETTERS

10 March 1986

C H A O S IN N E U R A L S Y S T E M S K.E. K O R T E N lnstitut fflr Theoreti~'che Physik, Uni~ersiAit zu KT~ln, 5000 Cologne. FRG

and J.W. C L A R K McDonnell Center for the Space Sciences and Department of Physics, Washington Unit)ersiO,, St. Louis, MO 63130, USA Received 13 August 1985; accepted for publication 12 December 1985

The occurrence of chaos in continuous-time neural network models is demonstrated through two examples, (i) a single neuron with an unusual kind of periodic input and (ii) a randomly connected network of 26 neurons with 7 incoming synapses per neuron, half the neurons being subject to a steady external stimulus. The former case is susceptible to exhaustive analysis, while the latter, discovered by computer simulation, is more attuned to neurobiological reality.

We wish here to present a rather general model of networks of neuron-like elements operating in continuous time and to demonstrate certain novel properties of this model. The existence of chaotic or unusually complex solutions in this model is of considerable intrinsic interest within the theory of nonlinear dynamical systems; there may also be important biological implications. The microscopic dynamics of a system of N interconnected neurons may be described with some realism by a set of NK coupled first-order nonlinear differential equations, the primary dynamical variables being the individual neuronal firing rates. This description invokes the common assumption that the mean interval between action-potential pulses emitted by any neuron is short compared to the time scale for appreciable changes in the imposed stimulus. We shall work in terms of a generalization of equations originally constructed by Stein et al. [1,2] to reproduce observed properties and response characteristics of individual neurons. The generalized equations read

a~-1 dX li/dt = -X li N K + rTla(fi(t) -- fOi + ]~ l Ci]Xl/ + k~=2bkiXki ) 0.375-9601/86/$ 03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

d X k i / d t = g k i [ X l i ( t); t] -- P k i X k i ,

(1)

with i = 1, ..., N and k = 2 ..... K. The function o(u) is supposed to be of sigmoid shape, being monotonic over _oo <( u < oo, and approaching the limiting values 0 and 1, respectively, as the argument u goes to its indicated extremes. To be definite, we take o(u) = [1 + exp(-u)] -1 , although the form 1/2 + 7r-1 t a n - l ( u ) would serve equaly well. The ingredients of these equations may be interpreted as follows. The principle variable x 1i - x i is taken as the firing rate of neuron i, while the auxiliary variables x/a , k f> 2, trace the time-course of adaptation of neuron i to its inputs. The total input F i to cell i is given by the sum of the first, third and fourth terms in the argument of the sigma function. The first term, f/, is the external stimulus or input; the third is the total synaptic input from all cells j of the net (including 0; and the fourth brings in a kind of latent self-inhibition and/or self-excitation. So far as singleneuron behavior is concerned, the proposed dynamical model may account for (a) the restricted range of firing rates due to an absolute refractory period (one parameter, ri), (b) decreased sensitivity of neuronal response at large excitatory or inhibitory inputs 413

Volume l14A, number 7

PHYSICS LETTERS

(through the sigmoid shape of the function o which summarizes the accumulation and processing of stimuli), (c) spontaneous firing, i.e, firing in the absence of input (one parameter,f0i), (d) response time for the increase of firing rate following a step increase of input (one parameter, a/-1 ,which lumps together synaptic delays, rise times of postsynaptic membrane potentials, time for passive spread of excitation to the axon hillock, etc.) and (e) the phenomenon of adaptation (accommodation) of neuronal firing rate to maintained inputs (one or more variables X k i , o n e or more level parameters bki and one or more rate constants Pki)" We note that explicit reference to the refractory periods r i may be eliminated by working with normalized, dimensionless impulse frequencies 2 i = rixi, lying between 0 and 1 and designated again simply as x i. We note also that the spontaneity parameter foi, a threshold value of input at which the steady firing rate is 1/2, may be absorbed into F i in the interests of economy. A third technical remark, which bears on property (b), is that one could adjust the width of the domain of linear response around F i - foi = 0 by scaling the argument of o with an additional parameter ui; however, this scaling constant may be gathered into Fi - foi without sacrificing any. mathematical generality. The topology of the neural network and the strengths of the interactions among the neurons are specified by the N 2 coupling coefficients cq. More specifically, ci] measures the efficacy of the activity x! of neuron ] in producing input to neuron i, represented as Giciix/. Patently, a positive (negative) coupling coefficient c i / m e a n s that the postsynaptic effect on i, of active ], is excitatory (inhibitory), and of course ell = 0 ifj has no synapses on i. In general, eli =¢=Cii; action and reaction are not equal and opposite in neural systems. It is conceptually helpful to think of the input to neuron i from a neuron] of the net, i.e. the contribution o f ] to Fi, as corresponding to the actionpotential-induced emission of transmitter substance into relevant afferent synapses of i, although in that case the use of a single warm-up constant a i to describe response both to external and internal inputs might be questioned. At any rate, other interpretations of the proposed equations are tenable at the level of cellular neurophysiology. For example, the quantity F i defined above might be considered as a mean soma potential of neuron i, measured relative to a threshold value (cf. ref. [3]). Another degree of latitude in inter414

10 March 1986

pretation is associated with the quantities gki [xt(t); t], which are meant to be arbitrary functionals o f x t i involving some arbitrary additional explicit time dependence. (We assume in any case that g is a smooth function of t.) Normally gki is just taken to be X l i , corres p o n d i n g to a simple mechanism for adaptation [1,2]. However, more elaborate modeling of neuronal behavior might exploit the element of arbitrary explicit time dependence to introduce a second or different channel for input to the cell, a channel which introduces response times longer (say) than a~ 1 . More broadly we may stress that the above system of equations is extremely general and flexible, with applications not only to the modeling of neural networks but also to condensed matter systems such as the current filaments in a semiconductor [4-6], to ecological systems and to economics. Within the context of continuous-time neural-net theories, the models of Cowan [7] and Hopfield [3] are more widely known than that sketched above. It can be argued that the present model is more flexible and therefore potentially more realistic than either of these. We point out that the steady states of all three models are formally the same, provided that (a) the same choice is made for the respective sigmoid functions and (b) all the adaptation channels are turned off(all b/k set zero) in the present formulation. A particularly elegant formal structure which embraces many of the existing continuous-time models has been presented by an der Heiden in ref. [8]. In refs. [ 1,2 ] the dynamical behavior of neural-network models based on eqs. (1), with gki = X l i , has been investigated in detail up to fourth order in N K . In particular, conditions for the existence of stable versus unstable steady-state solutions and for the occurrence of stable, sustained periodic oscillations of the neuronal f'tring frequencies have been derived analytically, based on linearization of the equations of motion [1,2,9]. Additionally, a significant numerical effort has been focused on special models of direct neurobiological interest. Among these are models designed to mimic invertebrate nervous systems [10] and simple cortical structures. (The latter include the mammalian olfactory bulb [11 ] and, in an extension involving time lags, the vertebrate retina [12] .) Numerical results have been obtained for the eigenvalues of the linearized system, for input-output relations and for typical steady-state and periodic solutions of

Volume 114A, number 7

PHYSICS LETTERS

the full nonlinear equations. However, as yet there has been no explicit consideration of the possibility that the equations admit chaotic solutions under appropriate conditions of input and/or connectivity. In view of the variety of nonlinear systems in which choas is known to occur (hydrodynamical problems, chemical reaction-diffusion systems, Josephson junctions, semiconductors, nonlinear optical systems, lasers ...), it would be surprising if neural networks were immune from this phenomenon. We shall establish here that certain neural systems described by the eqs. (1) of Stein et al. do indeed have chaotic solutions. Our first example exploits recent work of Aoki et al. [6] on the bifurcation routes to chaos in the firing wave instability of the current filaments of a semiconductor, based on a model involving four nonlinear differential equations. With a slight change of notation, their equations read x1 = - X l + °(2elXl - 2 y l ) , .~1 = (e2Xl + e3Y2) - 71Yl , x2 = -72x2 + ° ( 2 e 4 x 2 - 2y2), Y2 = (eSx2 + 3'3.vO) -- "Y3Y2 ,

(2)

wherein the x i and the Y i are considered, respectively, as the fftring activities and ignition thresholds of two subsystems i = 1,2; the 7l, as damping rates; and the e t, as temporal weights associated with memory of earlier firing. We note that the last two equations (for the motion of subsystem 2) decouple from the f'trst two (for subsystem 1). Thus they may be solved separately, the solution for Y2(t) serving as a modulating input to subsystem 1. In the circumstances leading to the demonstrated chaotic behavior of element 1, the second element or subsystem acts as a driven oscillator (with frequency adjusted by es), the background threshold levely 0 and the parameters 3'1- 3, e l - 4 being chosen such that there are three equilibrium points (two stable, the third an unstable saddle point) for element 1. The nature of the chaotic motion of the variable x 1(t) was then mapped out (with bifurcation diagram, Lyapunov characteristic exponents, phase portraits and power spectra) as a function of the control parameter e 5 . It is transparent that the set of equations (2) studied by Aoki et al. is in fact a special case of eqs. (1)

10 March 1986

w i t h N = 1 a n d K = 2. One setsxll = x I and x21 = Y l , takes g21 = e2Xl + e3Y2 and makes the identifications a 1 = 1,Cll = 2 e l , b 1 = - 2 , P l = 71,f1( t ) = 0 . T h u s , chaotic activity can be exhibited already at the level of a single "neuronal" element with periodically oscillating input. (The involvement of a time-varying input accords with the known fact that at least three differential equations are needed in this sort of nonlinear problem for a chaotic solution to exist; another option would be to insert time lags into the dynamics of the one-neuron system. An interesting example of single-neuron chaos which exploits the latter option has been presented by Mackey and an der Heiden [13], with the further important distinction that a nonmonotonic nonlinearity is assumed rather than a monotonic (sigmoid) nonlinearity as in system (1).) As established in ref. [6], the example we have constructed possesses the universality exhibited by the forced Lorenz model [14] : period-doubling cascade,inverse cascade and the usual tangent bifurcation occur in one of the asymmetric basins of attraction, while chaoschaos transitions and doubly tangent bifurcation result from the random switchings from one basin of attraction to another. It must be understood, however, that the neuronal element required for this example has rather unusual properties compared to the model neurons studied by Stein et al. and other authors. In particular it makes use of the rather complicated sort of second-order input alluded to earlier in discussing the flexibility allowed by the g functionals. The relevant input function, namely Y2(t), is to be fed into g21(t) from the solution of the second (decoupled) pair of equations in (2); the ordinary input f l ( t ) is zero. Moreover, the specific neuronal parameters for this case, a 1 = 1 (s), Cll = 12 and b 1 = - 2 , are certainly atypical (see below), and would imply rather sluggish behavior (activity variations on a scale of seconds rather than tens of milliseconds). The neurobiological interpretation given this system by Aoki et al. may, accordingly, be more appropriate: one envisions two neuron pools described by statistical, macroscopic activities Xl, x 2 and ignition thresholds y 1, Y2, with a special one-way interaction between the pools. Let us turn now to examples of neural systems which have the advantage of greater claim to realism but which, involving much higher dimensionality, have the disadvantage of being less amenable to compre415

Volume 114A, number 7

PHYSICS LETTERS

hensive and incisive analysis. In proceeding to larger systems one could of course try to connect together a number of individually chaotic elements like that examined above; however the composite system need not in general manifest chaotic behavior. In fact, we first encountered behavior resembling chaos h~ the neural context during a systematic investigation of the stability of steady states of randomly connected nets of simple, identical model neurons with parameters quite different from those of the first example. The class of randomly connected networks chosen for study is specified as follows. (1) Single-neuron properties. The dynamical behavior of neuron i is modeled by eqs. (1) with K = 2 and g[xi] = x i. As indicated above, r i is effectively replaced by unity. All neurons are taken to have the same intrinsic parameters, namely foi = f 0 = 0, ai = a = 100 s 1, bl i = b = - 2 0 0 , P l i = p = 10 s -1. These values produce a reasonable match of observed neuronal responses; they define a standard model which has received extensive numerical documentation. The individual neuronal elements are inherently stable [1 ]. (2) Connectivity o f network. The pattern of synaptic connections, their signs and their absolute strengths are chosen with the aid of some randomnumber generators. (i) Each neuron i is assigned exactly m nonzero ci/'s, where 0 <~m <.N, the distinct neurons/ feeding information to i being picked at random and democratically from the N neurons of the pool. (ii) A prescribed fraction h of the totality of effective n e u r o n - n e u r o n linksj -+ i in the net is taken inhibitory, each nonzero ci! having an equal chance of being negative. (iii) The magnitudes of the nonzero c#'s are decided by sampling a uniform distribution on the interval (0, L]. Typical choices for the gross network parameters (guided by studies of the basic circuit of the olfactory bulb [11,15] and used in the cases cited below) are: h = 0.5, L = 90. We have considered N up to 80, and examples with m ranging from 1 (very sparsely connected) to N (fully connected net). (3) Input. The network is stimulated by steady external inputsfi to a randomly selected set of neurons i, N s in number, with zero external input to the remainder. In the cases to be cited, Ns/N = 1/2. We take all nonzero inputs to have the same value f and consider a wide range in this control parameter. So far we have found at least two particular quasi416

10 March 1986

random nets which display genuine chaotic behavior. The example which has been scrutinized most intensively h a s N = 26 with m = 7 incoming connections per neuron. A thorough numerical search (using a damped N e w t o n - R a p h s o n method [16]) has revealed only one fixed point (steady state), which is unstable, while in some of the other nets considered we are able to find as many as 3 0 4 0 fixed points. Upon varying the external input ]; we observe instances of putatively chaotic output of a subset of neurons (some 9 in number) in the range 36 < f < 43, where the limiting values stated are not meant to be precise. These neurons. nontrivially active, are dubbed eligible, the remainder of the neurons being locked into a condition of either total inactivity (zero impulse frequency) or saturated activity (with firing rate unity) across the regime of interest. Increasing the stimulusfpast the lower critical boundary, the indicated subset of neurons undergoes a transition from simple periodic motion of the x 1i, with two oscillations per period and a period of roughly 40 ms, to a condition resembling intermittency; increasing f through the upper boundary, there is a transition from apparent chaos to periodic motion. By f = 42.8 the solutions are already periodic, with exceptionally long transients in the existing runs. A t f = 42.8, there are three oscillations per period; at f = 44.8, seven; and at f = 46.8, eleven (with a period of some 800 ms). By f = 54.8 we fred once more a simple mode with two oscillations per period (and a period of about 40 ms). It seems that different branches can be reached depending on the intial conditions, and different neurons generally have solutions of quite different structure. Solutions for the individual-neuron firing rates have been computed as far out as 25 s. An obvious general remark which is especially pertinent to systems with such high dimensionality (52 and 160, in the cases at hand) is that in many situations where the numerical solution looks chaotic, one may in reality just be seeing a very long transient, or else a periodic mode o f extremely long period [ 17]. Ideally, one would like to have available the Lyapunov exponents to remove any doubt, but these are extremely difficult to determine. Here we must be content with examination of power spectra and (more definitively) estimation of fractal dimensionality. Fig. 1 provides a selection of numerical solutions, in both ordered and chaotic domains, for the firing

Volume 114A, number 7

PHYSICS LETTERS

10 March 1986 I

I

I

-2

2.0

5.0

Z5

10.0

~

12.5 ~-8 0

!~ 'l Iii'

--

~ 40 I

I 80 I

t

120

160

I

~0

20.0

22.5

2510

-9

1 0

0

m

18.0

20.0

22.5

25.0

40

l 80 w (sec-I )

120

160

Fig. 3. Power spectra at two values of external input: (a) f = 41.8 (chaotic case) and (b) f = 46.8 (periodic). Here, P(to) = J'~3c(t) exp(itot) dr, where k(t) = x(t) - T -1 f l x ( u ) du and T is the upper limit of the computed time series.

t (sec)

Fig. 1. Representative solutions for the firing rate x(t) versus time t, of an eligible neuron at three values of the external input (control parameter): (a)f= 41.8 (chaotic case), (b)f = 46.8 (periodic) and (c)f = 36.3 (chaotic). rate of a representative eligible neuron (called neuron 1). Two sample phase projections, in which the output of neuron 1 is plotted against the simultaneous output of another chosen eligible neuron (neuron 2), are presented in fig. 2. Illustrative power spectra, based on 8000 data points, are juxtaposed in fig. 3 (note the difference in vertical scales). A t f = 36.3 it appears that we are just within the boundary o f the chaotic regime, while the solution at f = 46.8 is definitely

0.5

0.1

O.

0.5

t.0

0.2

0.5

1,0

Fig. 2. Phase portraits at two values of external input: (a) f = 41.8 (chaotic case) and (b)f= 46.8 (periodic). Firing rate xl (t) of cell 1 is plotted against firing rate x2(t) of cell 2. The initial transient is not included (b).

periodic. For f = 40.8, in a run very similar to that at input 41.8, chaos prevails. This has been verified convincingly by implementing the algorithm proposed in ref. [18] for the characterization of an attractor using a single-variable time series. The correlation exponent v (which provides a (generally tight) lower bound of the fractal dimension D) is found to be 3.1, whereas the expected value of unity is obtained in the periodic case at f = 46.8. Further (though less compelling) evidence of the chaotic nature of the f = 40.8 and f = 41.8 solutions is furnished by the associated phase portraits (see, e.g., fig. 2a) and especially by the power spectra we have derived. Referring to fig. 3a, we see the characteristic broad continuum at low frequency. There is an obvious contrast with the periodic example ( f = 46.8) of fig. 3b, which is distinguished by the occurrence of sharp peaks. The power spectrum for the run at f = 36.3 - on the edge of the chaotic regime - has an intermediate appearance. The high dimensionality of this example precludes a full understanding of its mathematical nature, in terms of bifurcation routes to chaos. Nevertheless, it will be interesting to accumulate information - derived numerically or otherwise - on the fine structure of the chaotic domain and of the adjacent regions, on arty windows o f limit cycles which might be embedded 417

Volume 114A, number 7

PHYSICS LETTERS

in chaotic bands, on the occurrence of tangent bifurcation, on the number of changes o f sign o f Lyapunov characteristic exponents (cf. the many changes of sign occurring in the first example), etc. At the same time it is important to discover other examples o f quasirandom nets displaying chaotic behavior. To date we have identified one other case where chaos is present, involving 80 neurons with 40 incoming connections per neuron, the f'mal certification again being made on the basis of the estimated fractal dimension. A further excellent candidate with N = 80 and m = 40 is under investigation. F r o m our experience we expect that chaos is not an uncommon occurrence for systems o f sizeable dimension, although it may be difficult to authenticate. There remians the important question o f the biological implications of chaos in neural networks. One's immediate judgement is that in general a tendency toward chaos would be harmful, an impediment to reliable response and to organized behavior conducive to survival. Regarding performance o f the mature animal, there would appear to be a need to suppress chaotic modes if the synaptic interactions (the cij ) and other conditions o f structure and stimulus permit it, by modifying the network parameters in some way. On the other hand, during development of the animal there would appear to be a need for some mechanism for avoiding, insofar as is possible, the patterns of synaptic and neuroanatomical organization which are favorable to chaos. The latter end is presumably achieved, at least in part, by genetic programming. However, it is not possible to fulfill either need by a mechanism which depends only on local monitoring o f activity - local in space and time - at the neuronal or synaptic level. Chaos is an emergent global, collective phenomenon; hence some global external controller, attuned to collective behavior, would appear to be required. (We may note in this connection that in some instances chaos can be eradicated by connecting together two or more systems undergoing chaos [ 19] .) In conclusion, one should not rule out the possibility that chaos could, in moderation, be beneficial, particularly in the context of higher mental processing where it is advantageous to have the potential of breaking out o f rigid, stereotyped behavior. Thus chaos, like neuronal spontaneity, opens access to a wider variety of responses and behavior, and may, if mechanisms exist for keeping it in check, be a valuable attribute o f complex biological systems. For extensive and per418

10 March 1986

ceptive commentary on the occurrence of chaos in neural systems and its role in neurobiology, the reader should consult refs. [20,21 ]. Funding for this work was provided in part by the National Science Foundation under Grant No. DMR83-04213 and in part by BRSG S07 RR0705420 awarded by the Biomedical Research Support Grant Program, Division of Research Resources, National Institutes o f Health. We thank J.W. Chen, U. an der Heiden, H. Moraal, K. Nildaus and H.G. Schuster for informative discussions. KEK acknowledges the kind hospitality o f the McDonnell Center for the Space Sciences and Department of Physics o f Washington University during several visits.

References [1 ] R.B. Stein, K.V. Leung, D. Mangeron and M.N. Oguztor~li, Kybernetik 15 (1974) 1. [2] R.B. Stein, K.V. Leung, M.N. Oguztor6li and D.W. Williams, Kybernetik 14 (1974) 223, and references therein. [3] J.J. Hopfield, Proc. US Natl. Acad. Sci. 81 (1984) 3088. [4] K. Aoki and K. Yamamoto, Phys. Lett. 98A (1983) 72. [5] K. Aoki, O. Ikezawa and K. Yamamoto, J. Phys. Soc. Japan 53 (1984) 5. [6] K. Aoki, O. Ikezawa and K. Yamamoto, Phys. Lett. 106A (1984) 343. [7] J.D. Cowan, in: Lectures on mathematics in the life sciences, Vol. 2, ed. M. Gerstenhaber (Am. Math. Soc., Providence, 1970). [8] U. an der Heiden, Analysis of neural networks (Springer, Berlin, 1980). [9] U. an der Heiden, Biol. Cybernetics 21 (1976) 37. [10] G.H. Paine, Ph.D. thesis, Washington University (1982), unpublished. [11] J.W. Chen, J.W. Clark and K.E. Kiirten, to be published. [12] M.N. Oguztor61i, Biol. Cybernetics 44 (1982) 1, and references therein. [13J M.C. Mackey and U. an der Heiden, J. Math. Biol. 19 (1984) 211. [14] E.N. Lorenz, J. Atm. Sci. 20 (1963) 130; Y. Aizawa and T. Uezu, Prog. Theor. Phys. 68 (1982) 64, 1864. [15 ] G.H. Shepherd, The synaptic organization of the brain, 2nd Ed. (Oxford Univ. Press, Oxford, 1979). [16] P Deuflhard, Num. Math. 22 (1974) 289. [17] S. Newhouse, Publ. Math. IHES 50 (1980) 101. [18] P. Grassberger and I. Procaccia, Physica 9D (1983) 189. [ 19] H.G. Schuster, Deterministic chaos: an introduction (Physik-Verlag, Weinheim, 1984); private communication. [20] E. Harth, IEEE Trans. Systems Man Cybernetics SMC-13 (1983) 782. [21] M.R. Guevara, L. Glass, M.C. Mackey and A. Shrier, IEEE Trans. Systems Man Cybernetics SMC-13 (1983) 790.