Noise-induced phase transitions in neuronal networks
Using an exactly solvable cortical model of a neuronal network, we show that, by increasing the intensity of shot noise (flow of random spikes bombarding neurons), the network undergoes first- and second-order non-equilibrium phase transitions. We study the nature of the transitions, bursts and avalanches of neuronal activity. Saddle-node and supercritical Hopf bifurcations are the mechanisms of emergence of sustained network oscillations. We show that the network stimulated by shot noise behaves similar to the Morris-Lecar model of a biological neuron stimulated by an applied current.
pacs:87.19.lj, 87.19.ln, 87.19.lc, 05.70.Fh
In the brain, interactions among neurons lead to diverse collective phenomena such as phase transitions, self-organization, avalanches, and brain rhythms Kelso_1995 (); Chialvo_2006 (). A non-equilibrium second-order phase transition was observed in human bimanual coordination Kelso_1984 (); hkb1985 (); Kelso_1986 (). Stimulation of living neural networks by electric fields causes a first-order phase transition Breskin_2006 (). There are evidences that brain rhythms, epileptic seizures, and the ultraslow oscillations of BOLD fMRI patterns also emerge as a result of non-equilibrium phase transitions Steyn-Ross_2010 (). The Hopf and saddle-node bifurcations are generic mechanisms for the emergence of oscillations in nonlinear differential equation models Strogatz_book1994 (). These mechanisms were found in the Morris-Lecar model of a biological neuron Rinzel1989 (); Izhikevich2000 (). In the context of brain rhythms, the Hopf bifurcations were discussed within mean-field cortical models Steyn-Ross_2010 () and for networks of randomly connected integrate-and-fire neurons bh1999 (); b2000 (); obh2009 (); lb2011 (). At the present time, understanding of nature of collective phenomena in the brain is elusive. For a complete description of a phase transition it is not enough to identify the bifurcation. It is also necessary to find critical phenomena that accompany it. In statistical physics, exactly solved models largely help us to understand phase transitions and critical phenomena Baxter_1982 ().
In the present paper, we propose an exactly solvable cortical model with stochastic excitatory and inhibitory neurons stimulated by shot noise (a flow of random spikes bombarding neurons). We show that shot noise stimulates first- and second-order non-equilibrium phase transitions in collective dynamics of neuronal networks. The first-order phase transition occurs as a discontinuous transition from low to high neuronal activity. Avalanches precede the transition. The saddle-node and supercritical Hopf bifurcations are the mechanism of emergence of sustained network oscillations. We find the order parameter for the continuous phase transitions and critical phenomena that accompany them in activity fluctuations. We show that the model exhibits collective excitability similar to excitability of the Morris-Lecar neuron stimulated by an applied current.
Model. We study a cortical model composed of excitatory and inhibitory neurons. is the network size, is the fraction of excitatory (inhibitory) neurons. Neurons are randomly connected with probability and form a classical random graph with Poisson degree distribution and the mean degree . The network is locally tree-like and has the small-world properties Albert_2002 (); dg2002 (); newman2003 () similar to ones found in brain networks sporns04 (). To take into account noise playing an important role in the brain dynamics lgns04 (); faisal08 (); ermentrout08 (), we assume that neurons are bombarded by random spikes represented by Dirac delta functions, , where are arrival times of spikes and is their amplitude. This kind of random input is a so-called shot noise. Random spikes represent spontaneous release of neurotransmitters in synapses (synaptic noise) and random spikes arriving from the remote part of the cortex. According to Schottky’s result, in the case of the Poisson distribution of interspike intervals, the power spectral density is proportional to the mean frequency of spikes, . We assume that the probability to receive random spikes during the integration time is Gaussian, where is the variance, is the mean number of spikes arriving during the time interval , and is the normalization constant, . We will use as the control parameter characterizing the shot noise intensity.
Neurons also receive delta-like spikes from active neighbors. The spikes mediate interaction among neurons. The total input includes spikes from shot noise and excitatory and inhibitory presynaptic neurons. We define the input to a neuron with index , , as integral of over the time interval . This gives
where , , and are the numbers of spikes arriving during the time interval from shot noise, active presynaptic excitatory and inhibitory neurons, respectively. We assume that efficacies of synaptic connections with excitatory and inhibitory neurons are uniform and equal to () and (), respectively. The numbers and are random and are determined by activity of presynaptic neurons during the interval . The network structure is encoded in the adjacency matrix.
In our model, neurons are tonic and the firing frequency versus input is the Heaviside function where is a threshold. is the same for both excitatory and inhibitory neurons. If and spike emission times of neurons are uncorrelated, then during the time interval , each active presynaptic neuron contributes to either one spike with probability or none with probability .
We consider stochastic neurons like those of Benayoun_2010 (); Wallace_2011 (); goltsev10 (); goltsev13 (). It means that response of neurons to an input is a stochastic process. Such a stochastic behavior might be caused by cellular noise and intensive bombardment by random spikes. Two rules determines dynamics of the cortical model: (1) If the input at an inactive excitatory (inhibitory) neuron at time is at least a certain threshold , then this neuron is activated with probability () and fires spikes. (2) An active excitatory (inhibitory) neuron is inactivated with probability () if . In our model, and are of the order of the first-spike latencies of excitatory and inhibitory neurons (from 6 to 100 ms, in the cortex). We introduce a parameter . If , inhibitory neurons have a larger response time than excitatory neurons. may be both larger and smaller than 1. We also introduce a dimensionless activation threshold . is of the order of 15-30 in living neuronal networks and about in the brain.
Rate equations. Behavior of the model is described by the fractions and of active excitatory and inhibitory neurons, respectively, at time . We will call them ‘activities’. We assume that activities are changed slightly during the integration time . Using the rules formulated above and the methods developed in goltsev10 (); goltsev13 (); dorogovtsev08 (), in the limit , we find explicit rate equations,
where . is the probability that, at given activities and , the input to a randomly chosen excitatory (inhibitory) neuron is at least . and represent fields acting on excitatory and inhibitory neurons. In the case of the classical random graph, we find
where , and are the probabilities that, during the time interval , a randomly chosen neuron receives spikes from excitatory and spikes from inhibitory neurons, respectively. Here .
Algorithm. In simulations, we builded a directed network, linking neurons with the probability . We divided time into intervals of width . At each time step, for each neuron we calculated input Eq. (1) given that each active presynaptic neuron contributes a spike with probability . The number of random spikes (shot noise) in this input was generated by the Gaussian process . Then, we updated states of neurons, using the rules formulated above. In the paper, numerical calculations are presented for parameters , , , , , and . We use as time unit and as input unit. Following amit97 (), we choose . We use and for the amplitude and variance of shot noise.
Steady states. If , then in a steady state, we have and is a solution of the steady state equation, . A graphical solution of the equation is displayed in Fig. 1. If shot noise intensity is either sufficiently small or sufficiently large, then there is only one solution, either point 1 or point 3. These fixed points correspond to the steady states with low and high neuronal activity, respectively. In an intermediate range there are three fixed points (1,2, and 3). At the critical point , points 2 and 3 coalesce. Points 1 and 2 coalesce at . From Fig. 1 one sees that the coalescence occurs when . Together with the steady state equation, this condition determines and . While the fixed points depend on , but not on , their local stability depends on both and . It is determined by eigenvalues of the Jacobian of Eqs. (2),
calculated at the fixed points. The eigenvalues are
where are the entries of the Jacobian. If at a fixed point, then this point is stable (attractor). If , then the point is unstable. If one of the eigenvalues is positive and the other is negative, then the point is saddle. If Re and Im, the point is a stable spiral. If Re and Im, the fixed point is an unstable spiral. The real and imaginary parts of determines the relaxation rate, , and the frequency, , of damped oscillations.
|3||–||s||s sp||u sp||u||s||s sp||u sp & lc||u & lc|
Phase diagram. Analyzing the local stability of the fixed points in the plane (see Table 1), we found the phase diagram displayed in Fig. 2. In regions Ia-Ie, the network relaxes exponentially to the stable fixed point 1 (if perturbations are small). In regions Ib and IIa, relaxation to the stable fixed point 3 is exponential while, in regions Ic and IIb, the relaxation occurs in the form of damped oscillations. In regions IIIa and IIIb, the fixed point 3 is an unstable point surrounded by a limit cycle. These are the regions with sustained network oscillations about the point 3. In Fig. 2, the phase boundaries represented by the dashed and solid lines are determined by the conditions and , respectively. Vertical lines represent lines and .
First-order phase transition. The pattern of collective behavior depends on and . First we consider the case ( is the critical value below which sustained oscillations appear). In simulations and numerical solution of Eqs. (2), we increased the noise level from zero (region Ia) to a value in region IIa (or IIb) above the critical point and afterwards decreased it again to a value below (see line 1 in Fig. 2). Neuronal activity undergoes a jump at ( in Fig. 3(b)).
Avalanches. Near , we observed bursts of neuronal activity caused by avalanches (see the inset in Fig. 3(a)). We studied avalanches, analyzing spike time series by use of the standard method (see Beggs_2003 () or the recent work Friedman_2012 ()). The avalanche size distribution is represented in Fig. 3(a). When is close to , is powerlaw, , in a broad range of . Using the maximum likelihood estimate (Touboul_2010, ), we found and the corresponding p-value is (the closeness of to 1 shows that the fit is good). This estimation agrees with the experimental data Beggs_2003 (); Friedman_2012 () and the standard mean-field exponent obtained for other exactly solved models Sethna_1993 (); Sethna_2001 (); Goltsev_2006 (). Thus, in the cortical model, bursts and avalanches are precursors of the first-order phase transition. Another mechanism of avalanches based on self-organized criticality is discussed in Chialvo_2006 ().
Hysteresis. If decreases from a value above to a value below , the network activity remaines as high as it was above (see Fig. 3(b)). The activity falls to a low value only at . Hysteresis occurs for and in regions Ib and Ic and is absent in regions Id and Ie where the fixed point 3 is unstable. Hysteresis was observed, for example, in living neural networks Soriano () and in simulations of thalamocortical systems izhikevich08 ().
An analytical analysis of the cortical model gives us deeper understanding of the first-order phase transition. Solving the equation at near and substituting into Eq. (5), we find that, in regions Ib and Ic the relaxation rate to the low activity state is , i.e., when . This phenomenon is known as critical slowing down. Solving Eqs. (2) with weak white-noise forces (), which mimic forces caused by finite-size effects, we found, using the linear perturbation analysis Thomas_1982 (), that the power spectral density (PSD) of the activity fluctuations, , has a sharp zero frequency peak, . If , the peak maximum increases as (note that in the high activity state, there is no singularity around ). If size is finite, then remains nonzero, though very small, even at the critical point due to finite-size effects that smooth the transition observed in simulations. This leads to a finite value of at .
Non-equilibrium phase transitions. An interesting scenario of collective dynamics takes place if inhibitory neurons respond slower than excitatory neurons and . In our analysis, we used simulations of the model, analytical and numerical solution of Eqs. (2). At a given , we start from and increase the shot noise intensity (see line 2 in Fig. 2). The neuronal network goes from region Ia with the single fixed point 1 into region Id or Ie where its dynamics is determined by three fixed points (see Table 1). At the points 1 and 2 coalesce and the network undergoes a non-equilibrium phase transition due to the saddle-node bifurcation. Above , the system is in region IIIa (or IIIb) with sustained oscillations (an unstable spiral or unstable fixed point surrounded by a limit cycle). Then, at on the boundary between regions IIIa and IIb, the network undergoes a phase transition due to the supercritical Hopf bifurcation. Above the network enters region IIb with damped network oscillations (a stable spiral in the point 3). The bifurcations were observed in the Morris-Lecar model stimulated by an applied current when the relation is N-shaped Rinzel1989 (). Thus, collective behavior of neuronal networks stimulated by a flow of random spikes bombarding neurons is similar to behavior of a single neuron stimulated by an applied current. From this similarity it follows that the network acts as an ‘integrator’, when it is close to the saddle-node bifurcation, and as a ‘resonator’, when it is close to the Hopf bifurcation Izhikevich2000 (). The important difference between the models is that, in our case, we deal with collective dynamics and phase transitions and it is necessary to find the order parameter and describe critical phenomena.
Saddle-node bifurcation. Above the saddle-node bifurcation, , the network oscillations emerge with a large amplitude (see Fig. 4(a)) and their frequency increases from zero as (see Fig. 4(c)). This behavior is generic for the bifurcation Strogatz_book1994 (); Rinzel1989 (); Izhikevich2000 (). We suggest that the frequency is the order parameter of the transition. In simulations, at below , we observed random almost identical bursts of activity that are the manifestation of critical fluctuations (inset of Fig. 5(a)). The bursts differ from ones found near the first-order phase transition in Fig. 3(a). We assume that the bursts are non linear excitations generated by finite-size effects Explanation (). The variance of the reciprocal of the interburst interval , , where is the mean burst frequency, has a peak at (see Fig. 5(a)) and behaves as about . Burst properties and generation mechanisms need further investigations.
Supercritical Hopf bifurcation. Near the supercritical Hopf bifurcation at , network oscillations differ from oscillations near the saddle-node bifurcation (see Figs. 4(a) and 4(b)). A decrease of the oscillation amplitude, , and the relaxation rate, , signals the supercritical Hopf bifurcation (see Fig. 4(c)). This behavior is generic for the bifurcation Strogatz_book1994 (). The amplitude is the order parameter for the transition. We found these results, expanding in Eqs. (2) in a series in around the fixed point 3 and holding terms up to inclusively. Then, we solved Eqs. (2) in region IIIa, using the averaging theory Strogatz_book1994 (). If tends to from the ‘paramagnetic’ region IIb with damped oscillations, then fluctuations of neuronal activity are increased, signaling the bifurcation. The PSD of the fluctuations has a resonance peak at the frequency of damped oscillations: , where , is the damping ratio, and the peak maximum (see also b2000 (); obh2009 (); lb2011 ()). These results agree well with our simulations in Fig. 5(b). They may help to understand the damping of alpha rhythms by visual or auditory stimuli Hari_1997 (); Lehtel_1997 ().
In conclusion, using the exactly solvable cortical model of neuronal networks with stochastic neurons, we showed that a flow of random spikes bombarding neurons controls behavior of the network. By increasing the flow, network oscillations appear at the saddle-node bifurcations and disappear at the supercritical Hopf bifurcation. Burst, avalanches, and critical fluctuations of neuronal activity are precursors of the non-equilibrium transitions. Within the model, neuronal networks are excitable systems having dynamical behavior similar to the Morris-Lecar model of a biological neuron.
Acknowledgements.We thank S. N. Dorogovtsev and J. F. F. Mendes for stimulating discussions. This work was supported by the PTDC Projects No. SAU-NEU/ 103904/2008, No. FIS/ 108476 /2008, No. MAT/ 114515 /2009, the project PEst-C / CTM / LA0025 / 2011, and the project ”New Strategies Applied to Neuropathological Disorders,” cofunded by QREN and EU. K. E. L. and M. A. L. were supported by the FCT Grants No. SFRH/ BPD/ 71883/2010 and No. SFRH/ BD/ 68743 /2010.
- (1) J. A. S. Kelso, Dynamic Patterns: The Self-Organization of Brain and Behavior (The MIT Press, Cambridge, 1995).
- (2) D. R. Chialvo, Nat. Phys. 2, 301 (2006); 6, 744 (2010).
- (3) J. A. Kelso, American Journal of Physiology - Regulatory, Integrative and Comparative Physiology 246, R1000 (1984).
- (4) H. Haken, J. A. S. Kelso, and H. Bunz, Biol. Cybern. 51, 347 (1985).
- (5) J. A. Kelso, J. Scholz, and G. Schöner, Physics Letters A, 118, 279 (1986).
- (6) I. Breskin, J. Soriano, E. Moses, and T. Tlusty, Phys. Rev. Lett. 97, 188102 (2006).
- (7) D. Steyn-Ross and M. Steyn-Ross, Modeling Phase Transitions in the Brain. Springer Series in Computational Neuroscience, v.4 (Springer, New York, 2010).
- (8) S. H. Strogatz, Nonlinear Dynamics And Chaos: With Applications To Physics, Biology, Chemistry, And Engineering ( Perseus Books Group, N.Y., 1994).
- (9) J. Rinzel, J. and G. B. Ermentrout, in Methods in Neuronal Modeling, edited by C. Koch and I. Segev, (The MIT Press, Cambridge, 1989), p. 251.
- (10) E. M. Izhikevich, Int. J. Bif. Chaos 10, 1171 (2000).
- (11) N. Brunel and V. Hakim, Neural Comput. 11, 1621 (1999).
- (12) N. Brunel, J. Comput. Neurosci. 8, 183 (2000).
- (13) S. Ostojic, N. Brunel, and V. Hakim, J. Comput. Neurosci. 26, 369 (2009).
- (14) E. Ledoux and N. Brunel, Front. Comput. Neurosci. 5, 25 (2011).
- (15) R. J. Baxter, Exactly solved models in statistical mechanics (Academic Press Inc., San Diego, 1982).
- (16) D. J. Amit and N. Brunel, Cereb. Cortex 7, 237 (1997).
- (17) R. Albert and A. L. Barabasi, Rev. Mod. Phys. 74, 47 (2002).
- (18) S. N. Dorogovtsev and J. F. F. Mendes, Adv. Phys. 51, 1079 (2002).
- (19) M. E. J. Newman, SIAM Review 45, 167 (2003).
- (20) O. Sporns, D. R. Chialvo, M. Kaiser, and C. C. Hilgetag, Trends Cogn. Science 8, 418 (2004).
- (21) B. Lindner, J. Garcia-Ojalvo, A. Neiman, and L. Schimansky-Geier, Phys. Rep. 392, 321 (2004).
- (22) A. A. Faisal, L. P. J. Selen, and D. M. Wolpert, Nat. Rev. 9, 292 (2008).
- (23) G. B. Ermentrout, R. F. Galan, and N. N. Urban, Trends Neurosci. 31, 428 (2008).
- (24) M. Benayoun, J. D. Cowan, W. van Drongelen, and E. Wallace, PLoS Comput. Biol. 6, e1000846 (2010).
- (25) E. Wallace, M. Benayoun, W. van Drongelen, and J. D. Cowan, PLoS ONE 6, e14804 (2011).
- (26) A. V. Goltsev, F. V. de Abreu, S. N. Dorogovtsev, and J. F. F. Mendes, Phys. Rev. E 81, 061921 (2010).
- (27) D. Holstein, A. V. Goltsev, and J. F. F. Mendes, Phys. Rev. E 87, 032717 (2013).
- (28) S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, Rev. Mod. Phys. 80, 1275 (2008).
- (29) J. M. Beggs and D. Plenz, J. Neurosci. 23, 11167 (2003).
- (30) N. Friedman, S. Ito, B. A. W. Brinkman, M. Shimono, R. E. L. DeVille, K. A. Dahmen, J. M. Beggs, and T. C. Butler, Phys. Rev. Lett. 108, 208102 (2012).
- (31) J. Touboul and A. Destexhe, PLoS ONE 5, e8982 (2010).
- (32) J. P. Sethna, K. Dahmen, S. Kartha, J. A. Krumhans, B. W. Roberts, and J. D. Shore, Phys. Rev. Lett. 70, 3347 (1993).
- (33) J. P. Sethna, K. Dahmen, C. R. Myers, Nature 410, 242 (2001).
- (34) A. V. Goltsev, S. N. Dorogovtsev, and J. F. F. Mendes, Phys. Rev. E 73, 056101 (2006).
- (35) J. Soriano and E. Moses (private communication).
- (36) E. M. Izhikevich and G. M. Edelman, Proc. Natl. Acad. Sci. U.S.A. 105, 3593 (2008).
- (37) P. Hanggi and H. Thomas, Phys. Rep. 88, 207 (1982).
- (38) The bursts are solutions of Eqs. (2). Their trajectories are heteroclinic orbits joining the points 1 and 2 around the unstable point 3 on the phase plane (see Ref. Rinzel1989 ()).
- (39) R. Hari and R. Salmelin, Trends Neurosci. 20, 44 (1997).
- (40) L. Lehtel, R. Salmelin, and R. Hari, Neurosci. Lett. 222, 111 (1997).