Critical phenomena and noiseinduced phase transitions in neuronal networks
Abstract
We study numerically and analytically first and secondorder phase transitions in neuronal networks stimulated by shot noise (a flow of random spikes bombarding neurons). Using an exactly solvable cortical model of neuronal networks on classical random networks, we find critical phenomena accompanying the transitions and their dependence on the shot noise intensity. We show that a pattern of spontaneous neuronal activity near a critical point of a phase transition is a characteristic property that can be used to identify the bifurcation mechanism of the transition. We demonstrate that bursts and avalanches are precursors of a firstorder phase transition, paroxysmallike spikes of activity precede a secondorder phase transition caused by a saddlenode bifurcation, while irregular spindle oscillations represent spontaneous activity near a secondorder phase transition caused by a supercritical Hopf bifurcation. Our most interesting result is the observation of the paroxysmallike spikes. We show that a paroxysmallike spike is a single nonlinear event that appears instantly from a low background activity with a rapid onset, reaches a large amplitude, and ends up with an abrupt return to lower activity. These spikes are similar to single paroxysmal spikes and sharp waves observed in EEG measurements. Our analysis shows that above the saddlenode bifurcation, sustained network oscillations appear with a large amplitude but a small frequency in contrast to network oscillations near the Hopf bifurcation that have a small amplitude but a large frequency. We discuss an amazing similarity between excitability of the cortical model stimulated by shot noise and excitability of the MorrisLecar neuron stimulated by an applied current.
pacs:
87.19.lj, 87.19.ln, 87.19.lc, 05.70.FhI Introduction
In the brain, interactions among neurons lead to diverse collective phenomena such as, for example, selforganization, phase transitions, brain rhythms, and avalanches Kelso_1995 (); Chialvo_2006 (); Buzsaki_2006 (). Among phase transitions, one can mention a nonequilibrium secondorder phase transition observed in human bimanual coordination Kelso_1984 (); hkb1985 (); Kelso_1986 (). Brain rhythms, epileptic seizures, and the ultraslow oscillations of BOLD fMRI patterns may also emerge as a result of nonequilibrium secondorder phase transitions SteynRoss_2010 (). Living neural networks stimulated by an electric field undergo a firstorder phase transition that can be seen as a jump of neuronal activity at a certain applied voltage Breskin_2006 (). Taking into account the role played by the collective phenomena in brain dynamics, it is very important to understand their nature and mechanisms. It is well known that bifurcations are responsible for the emergence of oscillations in nonlinear dynamic models Strogatz_book1994 (); Kuznetsov_1998 (), for example, the HodgkinHuxley model of a biological neuron Rinzel1989 (); Izhikevich2000 (). In the context of brain rhythms, the Hopf bifurcation was discussed in the framework of meanfield cortical models SteynRoss_2010 (), models of randomly connected integrateandfire neurons amit97 (); bh1999 (); b2000 (); obh2009 (); lb2011 (); ldl2005 (); bh2008 (); mmkn10 (), and networks of stochastic spiking neurons Benayoun_2010 (); Wallace_2011 (). However, when studying a phase transition, it is not enough to identify the type of bifurcation. It is also necessary to reveal and study critical phenomena accompanying the transition Stanley_book (). In the brain, various patterns of spontaneous activity representing collective phenomena were observed, such as neuronal avalanches Beggs_2003 (); Chialvo_2006 (); Plenz_2007 (), paroxysmal activity Steriade_1995 (); Timofeev_2004 (), spindle oscillations cdss1996 (), and many others. Despite a significant progress, understanding of collective phenomena in the brain and bifurcation mechanisms of phase transitions is elusive.
A neuronal network undergoes a phase transition from one to another state when a control parameter, such as an applied voltage or a flow of spikes bombarding neurons, reaches a critical value. In some cases (for example, for epileptic seizures), it is necessary to foresee that a neuronal network approaches to the critical point. An analysis of patterns and statistics of spontaneous neuronal activity and critical phenomena near the critical point may be a useful method for solving the problem. Nowadays, a comprehensive analysis of the critical phenomena in neuronal networks is far from to be complete.
In statistical physics, exactly solved models largely help us to understand phase transitions and critical phenomena Baxter_1982 (). Unfortunately, even simple versions of neuronal networks composed of integrateandfire neurons are very complex for an analytical consideration amit97 (); bh1999 (); b2000 (); obh2009 (); lb2011 (); ldl2005 (); bh2008 (); mmkn10 (). In the present paper, we study analytically and numerically an exactly solvable cortical model with stochastic excitatory and inhibitory neurons on complex networks. In the framework of this model, we consider first and secondorder phase transitions stimulated by shot noise (a flow of random spikes bombarding neurons). We also study critical phenomena accompanying the transitions and patterns of spontaneous activity signaling the transitions. First, we study a noiseinduced firstorder phase transition from low to high neuronal activity. The transition occurs if inhibitory neurons respond faster on stimuli than excitatory neurons. We demonstrate that bursts and avalanches of neuronal activity precede the transition. Second, we study two noiseinduced secondorder phase transitions that occur if inhibitory neurons respond slower on stimuli than excitatory neurons. The transitions represent two scenarios of appearance and disappearance of sustained network oscillations. We show that, when increasing the shot noise intensity, at first, sustained network oscillations appear due to a saddlenode bifurcation and then, at a higher shot noise intensity, the oscillations disappear due to a supercritical Hopf bifurcation. We study patterns of spontaneous neuronal activity near the bifurcations. We show that sharp paroxysmallike spikes of activity precede the secondorder phase transition caused by the saddlenode bifurcation. Above the Hopf bifurcation, spontaneous activity appears in a form of irregular spindles formed by damped oscillations. We also study analytically and numerically sustained network oscillations near the critical points of the bifurcations. Furthermore, we analyze the power spectral density (PSD) of spontaneous neuronal activity and its dependence on the noise intensity. We show that the PSD depends strongly on the bifurcation mechanism and the closeness to the critical point. We compare our results with experimental data and previous theoretical investigations. Finally, we discuss an amazing similarity between excitability of the considered cortical model stimulated by shot noise and excitability of the MorrisLecar neuron ML81 () stimulated by an applied current.
Ii 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 treelike and has the smallworld properties Albert_2002 (); dg2002 (); newman2003 () similar to ones found in brain networks sporns04 (). Our model also takes into account noise playing an important role in the brain dynamics white00 (); lgns04 (); faisal08 (); ermentrout08 (). We assume that neurons are bombarded by random spikes represented by Dirac delta functions,
(1) 
where are arrival times of spikes and is their amplitude. This kind of random input is socalled shot noise. The flow of random spikes bombarding neurons represents a combined effect of synaptic noise (spontaneous release of neurotransmitters), stimuli from other brain areas or sensory stimuli. 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, . In the present paper, we assume that the probability to receive random spikes during the integration time is Gaussian,
(2) 
where is the variance, is the mean number of spikes arriving during the time interval , and is the normalization constant, . We use as the control parameter characterizing the shot noise intensity.
Neurons also receive deltalike spikes from active neighbors. The spikes mediate interaction among neurons. We assume that efficacies of synaptic connections with excitatory and inhibitory neurons are uniform and equal to () and (), respectively. The total input includes spikes from shot noise and excitatory and inhibitory presynaptic neurons. We define the input to a neuron with index , , as the integral of over the time interval . It gives
(3) 
where , , and are the numbers of spikes arriving during the time interval from shot noise, active presynaptic excitatory and inhibitory neurons, 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
(4) 
where is a threshold. The frequency 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 the response of a neuron 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:

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.

An active excitatory (inhibitory) neuron is inactivated with probability () if .
We introduce a dimensionless activation threshold . is of the order of 1530 in living neuronal networks Eckmann07 (); breskin06 (); Soriano08 () and about in the brain. In our model, and are of the order of the firstspike latencies of excitatory and inhibitory neurons (from 6 to 100 ms in the cortex Heil04 (); Swadlow03 (); fmi04 (); mfmt05 ()). We introduce a parameter,
(5) 
If inhibitory neurons respond faster on stimuli than inhibitory neurons, i.e., the response time of an inhibitory neuron is smaller than the response time of an excitatory neuron, then . If excitatory neurons respond faster, i.e., , then . In the cortex, may be both larger and smaller than 1 Heil04 (); Swadlow03 (); fmi04 (); mfmt05 ().
ii.1 Rate equations
The behavior of the cortical 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 particular, the method of generating functions goltsev13 (), in the limit , we find explicit rate equations,
(6) 
where . is the probability that, at given activities and , input to a randomly chosen excitatory (inhibitory) neuron is at least . and represent fields acting on excitatory and inhibitory neurons. Note that the rate equations (6) are similar to the WilsonCowan equations Wilson_1972 () and rate equations derived for a stochastic rate model in Benayoun_2010 (); Wallace_2011 (). In the case of the classical random graph, we find
(7) 
where . , , and are the probabilities that, during the time interval , a randomly chosen neuron receives random spikes from shot noise, spikes from excitatory neurons, and spikes from inhibitory neurons, respectively. Note that the Poisson function is the probability that a randomly chosen neuron has presynaptic connections. Below we will study analytically and numerically Eqs. (6) and compare with simulations of the cortical model.
Our cortical model based on goltsev10 () is similar to the stochastic model of spiking neurons proposed by Benayoun et al. Benayoun_2010 (). Both models consider networks of stochastic neurons (‘inputdependent stochastic switches’ by Benayoun_2010 ()). The difference between the models is in some details about how to describe activation and deactivation processes and external input. Benayoun et al. Benayoun_2010 () assume that each neuron spikes with a rate dependent on its total synaptic input, while the resulting spiking activity decays at a constant rate independent on the input. In our model, we use a similar activation rule, while spiking activity decays with a certain rate only if the input becomes smaller than a threshold. The rates for activation and decay are different in Benayoun_2010 (), in contrast to our model where they are the same. Benayoun et al. assume that external input to each neuron is fixed in contrast to our model where external input is represented by shot noise. It is not surprising that, despite these differences, these models demonstrate similar dynamics. The advantage of the models with stochastic neurons is that they can be solved explicitly. Benayoun et al. Benayoun_2010 () and Wallace et al. Wallace_2011 () derived explicit rate equations for networks with alltoall connections while sparse randomly connected networks (classical random networks) were studied numerically. Methods of complex network theory dorogovtsev08 () allowed us to find explicit rate equations for neuronal networks on classical random graphs goltsev13 () and scalefree networks goltsev13 () and apply the model to study stochastic resonance Lopes2012 () and the role of synaptic plasticity Lee2012 ().
ii.2 Algorithm of simulations and parameters
In simulations, we built 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. (3) 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 our paper, we present numerical calculations for parameters , , , , , and . We analyze dynamical behavior of the cortical model in dependence on only two parameters: the parameter and the shot noise intensity. The latter parameter is the control parameter. Throughout this paper we use as time unit and as input unit. Following amit97 (), we choose . We use and for the amplitude and variance of shot noise.
ii.3 Steady states
The shot noise intensity determines activities and of excitatory and inhibitory populations at given model parameters. At zero fields , from Eqs. (6) we obtain in a steady state (). is a solution of the steady state equation,
(8) 
A graphical solution of the equation is displayed in Fig. 1. If the shot noise intensity is either sufficiently small or large, then there is only one solution, either point 1 or point 3. These fixed points correspond to steady states with low and high neuronal activity, respectively. In an intermediate range there are three fixed points (1,2, and 3). The critical point is the point where fixed points 2 and 3 coalesce. Fixed points 1 and 2 coalesce at . From Fig. 1 one sees that the coalescence occurs when
(9) 
Together with the steady state equation (8), the condition (9) determines the critical points and .
While the fixed points depend on , but not on , their local stability with respect to small perturbations depends on both and . It is determined by eigenvalues of the Jacobian of Eqs. (6),
(10) 
calculated at the fixed points. The eigenvalues are
(11) 
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 fixed points and their stability determine a phase portrait of Eqs. (6).
If the neuronal network is weakly perturbed from an equilibrium state corresponding to a stable fixed point , then the real and imaginary parts of at this point determines a relaxation rate to the state,
(12) 
and the angular frequency of damped oscillations about the fixed point,
(13) 
Ia  Ib  Ic  Id  Ie  IIa  IIb  IIIa  IIIb  
1  s  s  s  s  s  –  –  –  – 
2  –  sd  sd  sd  sd  –  –  –  – 
3  –  s  s sp  u sp  u  s  s sp  u sp & lc  u & lc 
ii.4 Phase diagram
Analyzing the local stability of the fixed points 1, 2, and 3 in the plane (see Table 1), we find the phase diagram of the cortical model displayed in Fig. 2. According to Table 1, in regions IaIe, the network relaxes exponentially to the stable fixed point 1 (of course, if a perturbation is 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 about the fixed point 3. 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. Nonlinear Eqs. (6) have different phase portraits in phase regions Ia–IIIb in Fig. 2. The phase portraits in the phase can be found by use of the standard methods Strogatz_book1994 (); Kuznetsov_1998 (). They determine the patterns of collective neuronal activity and response of the network on stimuli.
In Fig. 2, the phase boundaries are represented by the dashed and solid lines. The vertical lines and are determined by the selfconsistent solutions of Eqs. (8) and (9) discussed in Sec. II.3. The boundary between regions IIa and IIb and regions IIIa and IIIb is determined by the condition,
(14) 
(see the dashed lines in Fig. 2). The phase boundary between regions Ic and Id and regions IIb and IIIa is determined by the condition,
(15) 
(see solid line in Fig. 2). According to Eq. (15), on the boundary between regions IIb and IIIa, the relaxation rate is zero, i.e., critical slowing down occurs. The point in Fig. 2 is a tricritical point of coexistence of three phases: the low activity state (regions Ic and Id), the high activity state (region IIb), and the state with sustained network oscillations (region IIIa). At the point , the line of the firstorder phase transition meets the lines of two continuous phase transitions. The point is the common point of regions Ia, Ic, and Id. For the parameters used in our paper, we find , , , and .
Iii Firstorder phase transition
In this section, we study critical phenomena accompanying the firstorder phase transition stimulated by shot noise. In particular, we study neuronal bursts and avalanches as precursors of the transition. Though bursts and avalanches have been broadly studied both experimentally and theoretically, understanding of their mechanism in the brain is elusive Beggs_2003 (); Chialvo_2006 (); Plenz_2007 (); mmkn10 (); Orlandi_2013 (). Here, apart the standard measurements of the distribution function of avalanches over size, we also study critical behavior of the relaxation rate, a dependence of the power spectral density (PSD) of activity fluctuations on the shot noise intensity, and discuss finitesize effects. We find a dramatic increase of the zero frequency peak of the PSD when the shot noise intensity tends to a critical point while above the point the relaxation rate is nonzero and there are no critical fluctuations.
The firstorder phase transition occurs if , i.e., when the response time of an inhibitory neuron to stimulus is small enough in comparison with the response time of an excitatory neuron. In simulations and numerical solution of Eqs. (6), 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). When increasing the noise intensity , the neuronal activity undergoes a jump at ( in Fig. 3(b)). Therefore, the critical point is the limiting point of the firstorder phase transition. This phase transition is caused by a saddlenode bifurcation that corresponds to coalescence of the stable point 1 and the saddle point 2. Simultaneously, at , the eigenvalue becomes zero while remains to be negative. The firstorder phase transition was also found in Benayoun_2010 (). The line of the firstorder phase transition ends up at the point on the phase diagram. If , the neuronal networks undergoes a secondorder phase transition at that will be discussed in Sec. IV.
iii.1 Avalanches
In simulations, at , we observe bursts of neuronal activity (see Fig. 3(a)). When the mean interburst interval decreases while the mean burst duration increases. The bursts are caused by avalanches (activation of a single neuron triggers activation of a cluster of neurons). These activation processes are stochastic. In our model, in networks of finite size, bursts are generated by finitesize fluctuations. 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 pvalue is (the closeness of to 1 shows that the fit is good). Our estimation is close to the value 1.62 found in a stochastic rate model Benayoun_2010 (). Avalanches with the exponent about 1.5 were also found near a saddlenode bifurcation in networks of leaky integrateandfire neurons with shortterm synaptic depression mmkn10 (). Our estimation also agrees with experimental data Beggs_2003 (); Friedman_2012 () and the standard meanfield exponent obtained for other exactly solved models goltsev10 (); Sethna_1993 (); Sethna_2001 (); Goltsev_2006 (); Dorogovtsev_2006 (); Baxter_2012 ().
iii.2 Hysteresis
At a given , if decreases from a value above to a value below , the network activity remains as high as it was above (see Fig. 3(b)). The activity falls to a low value only at a critical intensity , which, in the general case, depends on and (see Fig. 2). If , where is the coordinate of the point on Fig. 2, hysteresis occurs in the range . If , hysteresis occurs in a smaller range of shot noise intensity, where is coordinate of the interception point of the line 1 with the phase boundary between region Ic and Id ending up at points S and T on the phase diagram in Fig. 2. The width of the hysteresis region, e.i., , tends to zero when . At , hysteresis is absent because, in regions Id and Ie, the fixed point 3 is unstable and there is only one stable fixed point 1. One notes that critical slowing down occurs at both limiting points of the firstorder phase transition, i.e., at in the low activity state and at (or ) in the high activity state. Hysteresis was observed, for example, in living neural networks Soriano () and in simulations of thalamocortical systems izhikevich08 ().
iii.3 Critical slowing down of neuronal dynamics
For deeper understanding of the firstorder phase transition, we now find analytically the relaxation rate to the low activity state. Writing Eq. (9) in the form and substituting it into Eq. (11), we find that at the eigenvalue is zero at the fixed point 1. Therefore, the relaxation rate, Eq. (12), to the low activity state is also zero,
(16) 
This phenomenon is socalled critical slowing down. Note that it takes place on the line at all , both above and below (see Fig. 2).
We now find dependence of the relaxation rate on at . We use the Taylor expansion of over small in Eq. (11) and obtain
(17) 
The first term is zero. Using Eq. (29) for from the Appendix A, in the leading order, we obtain
(18) 
This behavior occurs both at and .
If a neuronal network has a finite but large size , then according to the scaling law hypothesis, near a critical point of a continuous phase transition, the relaxation rate is described by the general scaling law,
(19) 
with exponents and which can be found by use of renormalization group techniques StaufferAharony (); Stanley_book (); Binder_1987 (). One assumes that the scaling law also is valid near the limiting point of the firstorder phase transition Sethna_2001 ():
(20)  
where . Thus, at a finite but large size , the relaxation rate is nonzero at any due to finite size effects that smear the critical singularity. This agrees with results of our simulations.
iii.4 Power spectral density of fluctuations near the firstorder phase transition
We now find the power spectral density (PSD) of activity fluctuations in the low activity state when is close to . In simulations, we measured the PSD of excitatory and inhibitory activities. We also solved analytically Eqs. (6) with weak whitenoise forces and , where . The forces mimic forces caused by finitesize effects (this method was also used in bh1999 ()). Our calculations are represented in Appendix B. We find that, in the low activity state, the PSD defined as
(21) 
where , has a sharp zero frequency peak described by the following shape function (see Eq. (40)):
(22) 
The peak maximum is . Fig. 4(a) displays the PSD measured in our+ simulations in the low activity state in region Ic. In Fig. 4(b) we compare simulations with the theoretical prediction. One sees that Eq. (22) describes well the measured frequency dependence of the PSD. According to Eq. (18), at , the peak maximum increases as
(23) 
Our simulations support the predicted increase of the zero frequency peak when . This behavior occurs if is not very close to . Due to finitesize effects, remains nonzero, though very small, even at the critical point (see Eq. (20)).
The Lorentzian behavior of the PSD of synaptic currents has been observed in cat cortex during wakefulness Bedard_2006a (). In Ref. Bedard_2006a (), it was suggested that this behavior may be driven by a white noise process. During slowwave sleep, the PSD deviates from the Lorentzian Bedard_2006a (). This deviation suggests that, in a general case, stochastic forces may be statistically different from white noise.
Thus, the cortical model shows that bursts and avalanches appear near the limiting point of metastable states of the firstorder phase transition caused by a saddlenode bifurcation in agreement with other network models Sethna_1993 (); Sethna_2001 (); mmkn10 (); Benayoun_2010 (). Critical phenomena (powerlaw statistics for avalanches and sharp zerofrequency peak of the PSD) due to critical slowing down in the low activity state (when approaching the critical point from below), the absence of the critical phenomena above the point (because, in the high activity state, the relaxation rate is nonzero at the critical point), and hysteresis are the characteristic properties of the first order phase transition, which can be experimentally tested. Another mechanism of avalanches based on ideas of selforganized criticality by Per Bak Bak () is discussed in Chialvo_2006 (). From our point of view, at the present time, there is no direct experimental evidence that supports one approach over the other. Further experimental and theoretical investigations of these two approaches are necessary for understanding avalanches in the brain.
Iv Secondorder nonequilibrium phase transitions
We now consider the case , i.e., when excitatory neurons respond faster on stimuli compared to inhibitory neurons. We show that, when increasing the shot noise intensity, the cortical model undergoes successively two secondorder phase transitions. We find that sustained network oscillations emerge at a saddlenode bifurcation and disappear at a Hopf bifurcation. We study properties of the phase transitions, critical phenomena, patterns of spontaneous activity, and sustained network oscillations near the critical intensities of shot noise.
iv.1 Saddlenode bifurcation
At a given , we increase shot noise intensity from (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: the stable point 1, the saddle point 2, and the unstable point 3 (see Table 1). At the points 1 and 2 coalesce and the network undergoes a secondorder phase transition due to a saddlenode bifurcation from a state with a low activity and shortrange temporal correlations between neurons into a state with regular sustained network oscillations (regions IIIa or IIIb). In regions IIIa and IIIb, dynamics of neuronal networks is determined by the unstable fixed point 3 surrounded by a limit cycle. At , the network oscillations emerge with a large amplitude (see Fig. 5(a)) and their frequency increases from zero as (see Fig. 5(c)). This frequency dependence is typical for the saddlenode bifurcation in nonlinear dynamic equations Strogatz_book1994 (); Rinzel1989 (); Izhikevich2000 (). Note however, that in our model, we deal with a phase transition, i.e., a collective phenomenon in neuronal networks. We suggest that for this kind of continuous phase transition the frequency is the order parameter.
In simulations, at below , we observed irregular almost identical sharp spikes of neuronal activity (see Fig. 6(a)). The mean frequency of the spikes is very small and increases when the shot noise intensity tends to the critical point while the spike duration is almost constant and much larger than the period () of oscillations generated by a single neuron. This kind of activity differs sharply from bursts found near the firstorder phase transition (compare Figs. 3 and 6). The sharp spikes emerge from a low background activity with a rapid onset (Fig. 6(b)). They reach a large amplitude, involve in synchronized activity about 90 of neurons, and end up with an abrupt return to lower activity. In Fig. 6, the spike duration is about 0.2 s and the mean interspike interval is about 34 s at ms.
In order to understand the mechanism of generation of the sharp spikes, we performed numerical integration of Eqs. (6) with nonzero stochastic forces and representing finitesize effects at the same parameters as in simulations. The numerical integration also reveals sharp spikes that are identical to ones observed in simulations (Fig. 6(b)). Our analysis of the phase portrait of Eqs. (6) in regions Id and Ie shows that the sharp spikes are strongly nonlinear events in neuronal activity generated by fluctuations. In the ( phase plane, their trajectories are topologically equivalent to the heteroclinic orbits found in the MorrisLecar model (see Fig. 7.4 in Ref. Rinzel1989 ()).
At near , the relaxation rate is given by Eq. (18). This result is in contrast to the standard mean field theory (the Landau theory) that predicts for a secondorder phase transition. The nonstandard scaling behavior and emergence of paroxysmallike spikes near the saddlenode bifurcation show an unusual character of the phase transition.
Analyzing properties of the sharp spikes, such as emergence conditions, course of the events, their shape, amplitude, duration, and low frequency oscillations, we find that this kind of spontaneous neuronal activity is similar to such epileptiform activity as the paroxysmal spikes observed in EEG activity Steriade_1995 (); Timofeev_2004 (). Based on this similarity we suggest that the paroxysmal spikes and other seizurelike events, such as slowwave oscillations Timofeev_2004 () or sharp waves in hippocampus Buzsaki_1985 (); Buzsaki_2006 (), are possible strongly nonlinear waves appearing in neuronal networks near a saddlenode bifurcation. Of course, in order to describe in detail the events, a realistic network structure and realistic singleneuron dynamics must be taken into account. As far as we know, paroxysmallike spikes as collective nonlinear objects were not studied analytically within a neuronal network model. A detailed investigation of the nature and mechanism of generation of the paroxysmallike spikes will be published elsewhere.
iv.2 Supercritical Hopf bifurcation
We now study the secondorder phase transition due to the supercritical Hopf bifurcation. For this purpose we perform simulations of the cortical model, numerical integration, and analytical analysis of Eqs. (6). We find critical behavior and demonstrate a difference in critical properties between the saddlenode and supercritical Hopf bifurcations.
When increasing the shot noise intensity above , the frequency of sustained oscillations increases while their amplitude decreases (see Fig. 5(c)). The oscillations disappear at a critical noise intensity which depends on (see the line 2 in Fig. 2). At , the network undergoes a phase transition from a state with the unstable point 3 surrounded by a limit cycle (region IIIa) into a state in which the fixed point 3 is a stable spiral (region IIb). From the stability analysis in Sec. II.3, it follows that this transition is due to the supercritical Hopf bifurcation. Above the network enters region IIb with damped network oscillations about the fixed point 3. Note also that network oscillations taking place near the saddlenode and supercritical Hopf bifurcations have different patterns (compare Figs. 5(a) and 5(b)). Oscillations emerging due to a Hopf bifurcation were also found in a stochastic rate model Wallace_2011 ().
First, we study sustained network oscillations at below . We expand in Eqs. (6) in a series in around the fixed point 3 and hold terms up to inclusively. Then, we solve Eqs. (6) in region IIIa, using the averaging theory Strogatz_book1994 (). Details of our calculations are in Appendix C. When , we find a decrease of the oscillation amplitude ,
(24) 
and a decrease of the relaxation rate ,
(25) 
that signals the supercritical Hopf bifurcation (see Fig. 5(c)). This behavior is typical for this kind of bifurcation Strogatz_book1994 (). The amplitude is the order parameter for the transition. In Appendix C, we show that a phase lag between synchronized activities of excitatory and inhibitory populations is proportional to , . is zero at . Comparing Eq. (25) with Eq. (18), we conclude that the continuous phase transitions corresponding to the saddlenode and supercritical Hopf bifurcations have different critical behaviors and, therefore, belong to different classes of universality.
We now analyze analytically the critical behavior of the cortical model at above . Our simulations show that, above , spontaneous activity has a form of spindle oscillations (see the inset in Fig. 7). The spindle oscillations are similar to spindles observed, for example, in thalamus cdss1996 (). Damped oscillations were observed in an instance of epilepsy (see, for example, Babloyantz_1986 ()). If tends to from above, then the amplitude of spindle oscillations increases. This results in an increase of the peak of the power spectral density of activity fluctuations at the frequency of damped oscillations (see Fig. 7). This critical phenomenon signals an approach to the Hopf bifurcation. In order to understand the phenomenon, we use simulations and analytical calculations. According to Appendix B.2, the PSD, , has a resonance peak,
(26) 
where and . The parameter ,
(27) 
is the damping ratio of damped oscillations. The peak maximum is . This behavior of is due to the fact that the amplitude of damped oscillations increases as when (see Eq. (44)). In turn, the amplitude decreases when increases and the network goes away from the supercritical Hopf bifurcation. Note that the relaxation rate determines the time decay of the damped oscillations and can be found from data analysis of a time dependence of the autocorrelation function, Eqs. (34) and (44). From this analysis one finds the dimensionless parameter that is an important characteristic of the closeness of the network to the critical point . The smaller is the close is the network to the critical point. In the infinite size limit, is zero at . A similar resonance peak of the PSD was also found within the integrateandfire model in b2000 (); obh2009 (); lb2011 ().
In Fig. 8 we represent the PSD of activity fluctuations measured in simulations. In agreement with the theoretical prediction, the measured PSD, , reveals a sharp maximum at the frequency of damped oscillations. Fig. 7 shows that, when , the maximum value first strongly increases and then saturates at a certain value due to the finitesize effects, Eq. (20). Fig. 8 shows that the shape of this maximum is well described by the shape function Eq. (26).
The critical behavior of the cortical model near the supercritical Hopf bifurcation helps to understand the attenuation of alpha rhythms by visual or auditory stimuli (the Berger effect) Hari_1997 (); Lehtel_1997 (). Recall that the Berger effect manifests itself in activation of alpha waves on the electroencephalogram when the eyes are closed and diminution of alpha waves when they are opened (see, for example, the review of Hari_1997 ()). Based on the cortical model, we suggest that opening eyes may result in an increase of the flow of spikes bombarding neurons in the area of the cortex that is responsible for the alpha waves. As a result, the neuronal network goes away from the Hopf bifurcation and the amplitude of damped oscillations decreases. A similar phenomenon was also observed in the auditory cortex where the tau rhythm (the tau rhythm belongs to the family of alpha rhythms) was transiently suppressed by auditory stimuli Lehtel_1997 ().
V Discussion
In Sec. IV, we discussed local stability of fixed points and bifurcations of nonlinear Eqs. (6) in the cortical model in dependence on the shot noise intensity and the parameter . Based on these results, one builds phase portraits of Eqs. (6). In the case , we revealed that the phase portraits in regions Id, Ie, IIIa, and IIIb are topologically equivalent (in other words, homeomorphic) to the phase portraits found in the MorrisLecar model stimulated by an applied current in the case when the relation is Nshaped Rinzel1989 (). Recall that the MorrisLecar model is a simplified version of the fourdimensional HodgkinHuxley model. Within the MorrisLecar model, a system of two nonlinear equations describes a relationship between the membrane potential and the activation of ion channels within the membrane. It is wellknown that the topological equivalence of phase portraits of two dynamical systems results in similar dynamics and similar responses on stimuli Kuznetsov_1998 (). Therefore, the dynamic behavior of the cortical model stimulated by shot noise (a flow of random spikes bombarding neurons) must be similar in some respects to the dynamic behavior of the MorrisLecar model stimulated by an applied current. In this case, we can apply results obtained for the wellstudied MorrisLecar model to the cortical model. Izhikevich Izhikevich2000 () showed that the MorrisLecar neuron acts as an ‘integrator’, when it is close to the saddlenode bifurcation, and as a ‘resonator’, when it is close to the Hopf bifurcation. Based on the topological equivalence, we can conclude that the cortical model acts in a similar way near the bifurcations. Indeed, in Sec. IV.1, we showed that if the mean frequency of incoming random spikes is a little bit larger than the critical frequency corresponding to the saddlenode bifurcation, then a neuronal network oscillates with an arbitrary low frequency. The higher the mean frequency of incoming random spikes, the higher the frequency of sustained network oscillations. Thus, we can say that the network acts as an ‘integrator’. In contrast, when the network is in the rest state near the supercritical Hopf bifurcation, it acts as a ‘resonator’ because it responds preferentially to a certain (resonant) frequency of input (see Sec. IV.2). Furthermore, in Sec. IV.1, the topological equivalence helped us to understand the nature of paroxysmallike spikes observed near the saddlenode bifurcation because similar nonlinear spikes were found in the MorrisLecar model by Rinzel and Ermentrout Rinzel1989 ().
Vi Conclusion
In conclusion, within an exactly solvable cortical model of neuronal networks with stochastic excitatory and inhibitory neurons, we studied first and secondorder phase transitions stimulated by shot noise (a flow of random spikes bombarding neurons). We performed simulations, numerical integration, and analytical analysis of nonlinear dynamical equations. These methods gave results in complete agreement with each other. The advantage of our model is that it gives a possibility to study both noiseinduced first and secondorder phase transitions in neuronal networks by use of a unified approach and standard analytical physical and mathematical methods. This unified approach allowed us to compare qualitatively and quantitatively critical phenomena accompanying the phase transitions, patterns of spontaneous neuronal activity, and their dependence on the shot noise intensity. Furthermore, the rate equations derived for the model allowed us to study strongly nonlinear events, such as paroxysmallike spikes and slow waves observed in neuronal activity, that cannot be described by a linear theory. Our results support the idea that collective behavior of neuronal networks may have universal properties that do not depend on details of single neuron dynamics. The universal collective phenomena are determined by general properties of neuronal networks, such as the network structure, balance between excitatory and inhibitory neurons, the presence of noise, and interaction between neurons.
We showed that if inhibitory neurons respond faster on stimuli than excitatory neurons, then a firstorder phase transition manifests itself as a jump from low to high neuronal activity at a critical noise intensity. We found the mechanism of the transition and showed that it occurs due to a saddlenode bifurcation. We studied in detail critical phenomena that accompany the transition and patterns of spontaneous activity near the critical point. In particular, we showed that bursts and avalanches are precursors of the firstorder phase transition. When the shot noise intensity tends to the limiting point of the metastable states then the size distribution of neuronal avalanches becomes a power law with the exponent about 1.5. Moreover, at the critical point, critical slowing down occurs in the infinite network, i.e., the relaxation rate is zero at the critical noise intensity. Our simulations revealed that finitesize effects smear the phase transition and make the relaxation rate to be nonzero at the critical point. Critical phenomena (powerlaw statistics for avalanches and sharp zerofrequency peak of the PSD) due to critical slowing down in the low activity state (when approaching the critical point from below), the absence of the critical phenomena above the point (because, in the high activity state, the relaxation rate is nonzero at the critical point), and hysteresis are the characteristic properties of the first order phase transition, which can be experimentally tested.
We also studied two noiseinduced secondorder phase transitions that occur if inhibitory neurons respond slower on stimuli than excitatory neurons. These transitions represent two scenarios of appearance and disappearance of network oscillations. When increasing the shot noise intensity, at first, sustained network oscillations appear due to a saddlenode bifurcation, and then, at a higher shot noise intensity, the oscillations disappear due to a supercritical Hopf bifurcation. Our analysis showed that the continuous phase transitions caused by the saddlenode and supercritical Hopf bifurcations are accompanied by different critical phenomena and different patterns of spontaneous neuronal activity. The transitions are characterized by different order parameters and belong to different classes of universality.
We analyzed patterns of spontaneous neuronal activity near the saddlenode and Hopf bifurcations. Our most interesting result is the observation of paroxysmallike spikes that precede the secondorder phase transition caused by the saddlenode bifurcation. We found that the spikes are strongly nonlinear objects that appear instantly from a low background activity with a rapid onset, reaches a large amplitude, involve in synchronized activity about 90 of neurons, and end up with an abrupt return to lower activity. These spikes are similar to single paroxysmal spikes and sharp waves observed in EEG measurements. With increasing the shot noise intensity above the critical point of the saddlenode bifurcation, low frequency network oscillations follow the irregular spikes. They appear with a large amplitude but a small frequency (at the critical point, the frequency is zero). The pattern of the oscillations resembles sharpslow waves Timofeev_2004 () or sharp waves in hippocampus Buzsaki_2006 (); Buzsaki_1985 (); Rex2009 (). In contrast to the saddlenode bifurcation, spontaneous activity above the Hopf bifurcation is represented by irregular spindles formed by damped oscillations. Sustained network oscillations below the supercritical Hopf bifurcation have a small amplitude (at the critical point, the amplitude is zero in the infinite size limit) and a finite frequency. These oscillations are also nonlinear and have properties like ones of the Van der Pol oscillator.
We also analyzed the power spectral density (PSD) of spontaneous neuronal activity near the critical points of the phase transitions. We showed that the frequency dependence of the PSD and its dependence on the shot noise intensity give a rich information about the kind of bifurcation and the closeness of the network to the critical point. In particular, the PSD has a zero frequency peak near the firstorder phase transition while above the supercritical Hopf bifurcation the PSD has a peak at the frequency of damped oscillations. The peaks are strongly enhanced when the noise intensity tends to the critical points of the phase transitions. These results may be applied to an analysis of spectral properties of EEG recording in order to predict an approach to a critical point in neuronal activity.
Finally, we discussed an amazing similarity between excitability of the considered cortical model stimulated by shot noise and excitability of the MorrisLecar neuron stimulated by an applied current ML81 (); Rinzel1989 (); Izhikevich2000 () . This similarity results from the fact that the cortical model of neuronal networks and the MorrisLecar model have topologically equivalent phase portraits. This similarity allowed us to conclude that a neuronal network acts as ‘integrator’, when it is close to the saddlenode bifurcation, and as a ‘resonator’, when it is close to the supercritical Hopf bifurcation. We believe this similarity may be useful for understanding many nonlinear phenomena in dynamics of neuronal networks.
In our model, a flow of random spikes bombarding neurons represents a combined effect of synaptic noise (spontaneous release of neurotransmitters), stimuli from other brain areas, and sensory stimuli. At given model parameters, the flow controls dynamics of the neuronal network. If the flow intensity is close to a critical value, then even a small change in the flow intensity can switch the network from one to another state. In other words, a small change of activity of neuronal networks to which a considered network is connected may strongly impact on a dynamic state of the network under consideration. This represents one of important mechanisms of interaction between neuronal networks Abbott_2006 ().
Vii Acknowledgements
We thank S. N. Dorogovtsev for stimulating discussions. This work was supported by the PTDC Projects No. SAUNEU/ 103904/2008, No. FIS/ 108476 /2008, No. MAT/ 114515 /2009, the project PEstC / 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.
Appendix A Neuronal activity near the critical points
Let us find the activity in the low activity state near the critical point of the saddlenode bifurcation, i.e., at . In simulations, can be found by measuring neuronal activity and averaging it over a sufficiently large observation time. In Eq. (8), we use the Taylor expansion of the function over and up to the second order in . Then Eq. (8) takes a form,
(28) 
where the function and its derivatives are calculated at . Using Eqs. (8) and (9), we find a solution,
(29) 
where
(30) 
The singular behavior Eq. (29) is a general attribute of hybrid and firstorder phase transitions (Dorogovtsev_2006, ; Goltsev_2006, ; Baxter_2012, ). Note that at the fixed point 2 is
(31) 
because, at , the point 1 and 2 coalesce and .
Neuronal activity near the Hopf bifurcation can be found at by use of the Taylor expansion Eq. (28) with and , where the function and its derivatives are calculated at . In this case, the linear terms give the leading contribution to a solution,
(32) 
in contrast to the square root dependence in Eq. (29).
Appendix B Power spectral density of activity fluctuations
The power spectral density (PSD) of fluctuations of neuronal activity encodes rich information about critical phenomena. According to the WienerKhintchine theorem, the power spectral density of activity fluctuations of the excitatory population is the Fourier transform of the autocorrelation function ,
(33) 
The autocorrelation function,
(34) 
where , describes fluctuations of activity of population , , around the averaged value . is a measure of correlations between values of and at two different instants separated by a lag and averaged over an arbitrary large time window (see, for example, in Gardiner_2002 ()).
In order to calculate the PSD, we assume that activity fluctuations are driven by weak whitenoise forces that mimic forces caused by finitesize effects,
(35) 
where . In this case, one can use the linear response theory and find from the linearized Eq. (6),
(36) 
where is a steady state solution of Eq. (8), , and . The Jacobian is given by Eq. (10). Making the Fourier transformation,
(37) 
we find the linear response,
(38) 
Substituting this result into Eq. (34), we find the PSD for excitatory neurons,
(39) 
The PSD of inhibitory neurons is obtained from this equation after replacements: and
b.1 PSD near the firstorder phase transition
At first, let us consider the power spectral density (PSD) of activity fluctuations in the low activity state (the fixed point 1) in regions Ib and Ic in Fig. 2. In these regions, eigenvalues and are real. When the noise intensity tends to the critical point of the firstorder phase transition, the eigenvalue tends to zero according to Eq. (18) while the eigenvalue remains finite. Therefore, at small , equation (39) takes a form
(40) 
b.2 PSD near the Hopf bifurcation
Now we consider the PSD of activity fluctuations in the high activity state (the fixed point 3) at (region IIb in Fig. 2). In this region the eigenvalues are complex. Their real and imaginary parts determines the relaxation rate and the frequency of damped oscillations, respectively (see Eqs. (12) and (13)). In this case, equation (39) can be written in a form,
(41) 
where , , and . is the damping ratio of the damped oscillations.
In the case when the shot noise intensity tends from above to the critical point , the relaxation rate tends to zero (see Eq. (25)). If , then the PSD has a sharp peak at the resonance frequency . The peak maximum is
(42) 
Near the resonance frequency , is described by a shape function ,
(43) 
Appendix C Oscillations near the supercritical Hopf bifurcation. nonlinear analysis.
In this appendix we study analytically oscillations around the fixed point 3 near the Hopf bifurcation in a noise intensity range (a range around the boundary between regions IIIa and IIb in Fig. 2). In this range, the oscillations have a small amplitude that allows us to use the Taylor expansion over in Eqs. (6). Assuming and taking into account terms up to the third order in , we obtain two coupled nonlinear equations,
(45) 
where and
(46) 
In Fig. 9(a) we compare results of numerical integration of the reduced equations, Eqs. (45), with the exact equations (6). In the numerical integration, we studied relaxation of the system to a state with sustained oscillations (see Fig. 9(a)) from an initial point . One sees that the frequency of the oscillations described by the reduced equations (45) is very close to the frequency of oscillations from exact Eqs. (6) though the amplitude of the sustained oscillations from Eqs. (45) is a little bit larger. These results evidence that the reduced equations (45) are a good approximation to the exact Eqs. (6). A similar analysis based on a reduced equation was used in bh1999 (); b2000 () to study analytically oscillations near the Hopf bifurcation in networks of integrateandfire neurons. Below we use the reduced equations to study a critical behavior of the amplitude of sustained oscillations, a relaxation rate to the state with the oscillations, and the phase lag between activities of excitatory and inhibitory populations.
It is convenient to rewrite Eqs. (45) in a vector form,
(47) 
where
is the Jacobian Eq. (10) and is a matrix which introduces nonlinear terms,
(48)  
The Jacobian , Eq. (10), can be represented in a form,
(49) 
where is the identity matrix. The parameter is determined by Eqs. (12) and (11) at the fixed point 3, i.e., . In regions IIb and IIIa, Eq. (11) gives . Furthermore, is a complex vector