Irregular Dynamics in Up and Down Cortical States

Complex coherent dynamics is present in a wide variety of neural systems. A typical example is the voltage transitions between up and down states observed in cortical areas in the brain. In this work, we study this phenomenon via a biologically motivated stochastic model of up and down transitions. The model is constituted by a simple bistable rate dynamics, where the synaptic current is modulated by short-term synaptic processes which introduce stochasticity and temporal correlations. A complete analysis of our model, both with mean-field approaches and numerical simulations, shows the appearance of complex transitions between high (up) and low (down) neural activity states, driven by the synaptic noise, with permanence times in the up state distributed according to a power-law. We show that the experimentally observed large fluctuation in up and down permanence times can be explained as the result of sufficiently noisy dynamical synapses with sufficiently large recovery times. Static synapses cannot account for this behavior, nor can dynamical synapses in the absence of noise.


Introduction
Neural systems, even in the absence of external stimuli, can exhibit a wide variety of coherent collective behaviors, as in vivo and in vitro experiments show [1][2][3]. One of the most prominent examples is the spontaneous transition between two different voltage states, namely up and down states, observed in simultaneous individual single neuron recordings as well as in local field measures. Such behavior, which is generated within the cortex, may provide a framework for neural computations [4], and could also coordinate some sleep rhythms into a coherent rhythmic sequence of recurring cortical and thalamocortical activities [3,5,6]. The phenomenon of up and down transitions has been measured in a number of situations, such as in the primary visual cortex of anesthetized animals [7,8], during slow-wave sleep [1,5,6], in the somatosensory cortex of awake animals [9], or in slice preparation under different experimental protocols [3,10,11], to name a few. The origin of such structured neuronal activity is still unclear, although several studies have shown that both intrinsic cell properties [12][13][14] and the high level of recurrency present in actual neural circuits [3,15,16] may contribute to the generation of up and down transitions. In particular, the contribution that reverberations in recurrent neural networks may have in the appearance of these transitions could depend strongly on synaptic properties. It is known, for instance, that excitatory synapses with a slow dynamics (such as synapses mediated by NMDA receptors) may play a relevant role in the generation of persistent activity or up cortical states [17]. On the other hand, several modeling studies indicate that activitydependent synaptic mechanisms, such as short-term synaptic depression and facilitation, can induce voltage transitions between up and down neural states as well [16,[18][19][20].
Many crucial points about the understanding of up and down transitions are, however, still lacking. For instance, in vivo experiments in the cat visual cortex show that the permanence times in the depolarized or up state present a high variability, and can range from a scale of milliseconds to seconds [7]. A similar level of irregularity has also been recently found in in vivo recordings of updown transitions in the rat auditory cortex [21], as well as in sleepwake transitions [15,22,23], where power-law distributions in the duration of wake states have been measured. Such complexity in the time series of the neuron membrane potentials remains far to be explained, and could reflect scale invariance in permanence times, which could in turn be a (preliminary) indicative of criticality. In fact, there are many recent studies that have shown criticality in different contexts in the brain [24,25], as well as in neural network models which present self-organization and criticality properties [26][27][28], and even it has been reported to occur in sleep-wake transitions in in vivo conditions [22,23]. Although it is worth noting that the irregularity of the dynamics of up and down states is not a sufficient condition for criticality, a concrete characterization of such irregularity may be a convenient starting point for future works on this topic.
To study in detail the relevant issue of irregular up and down cortical dynamics, we propose in this work a minimal model for up and down transitions in neural media. We consider a simple bistable rate model whose stable solutions represent two possible voltage states of the mean membrane potential of the network. More precisely, such states correspond, respectively, to high and low levels of activity in the network (that is, the up and down cortical states). In addition, we consider that the synaptic connections between neurons of the network present short-term synaptic depression (STD) mechanisms, which introduce temporal correlations, as well as synaptic stochasticity, in the dynamics of the system [29][30][31][32]. A complete analysis of this simple mathematical model depicts, both numerically and within a theoretical probabilistic approach, the appearance of power-law dependences in the distribution of permanence times in the up state. Our results show that the appearance of such scale free distributions is due to the complex interplay between several factors including synaptic stochasticity and the temporal correlations introduced by STD. The emergence of power-law dependences could, indeed, explain the high variability in permanence times in the up state suggested by experiments [7,21].

Methods
Our starting point is a bistable rate model, which mimics the dynamics of the electrical activity of a population of interconnected excitatory neurons (although it can be easily extended to other situations) with two stable levels of activity. The model has the form [33] where n(t) is the mean firing rate of the (homogeneous) neural population, n m is the maximum level of activity which can be reached by the population (in absence of noise), J(w0) is the synaptic coupling strength in absence of STD, and h is the firing threshold of the neurons in the population. The variable f(t) is a Gaussian white noise of zero mean and standard deviation d, which takes into account the inner stochasticity of the neural population (caused by other sources of uncontrolled noise in the system). The parameter t n is the population time constant, which may be assumed to be around the duration of the synaptic current pulse [34,35]. For generality purposes, we set t n~1 , and therefore time and frequency are given in units of t n and t {1 n , respectively. The term S(z): 1 2 ½1z tanh (z) represents the transduction function, which gives the nonlinear effect that the mean postsynaptic current (coming from recurrent connections of the neural population) induces in the network mean firing rate. Employing this form for S(z), the up and down stable levels of activity correspond to n^n m and n^0, respectively.
On the other hand, the variable x(t) in equation (1) takes into account the dynamical modification of the strength of the synaptic connections during short time scales due to high network activity, and it is usually named short-term synaptic plasticity. Based on the model proposed in [29,36] for short-term depression, and following previous studies concerning the dynamics of neural populations [16], we assume that x(t) evolves according to where t r is the characteristic time scale of the STD mechanism, and u is a parameter related with the reliability of the synaptic transmission. According to experimental measurements for these parameters in the somatosensory cortex of the rat [36], we set t r~1 000 and u~0:6 unless specified otherwise. Assuming, for instance, a population time constant of t n~1 ms, which would approximately correspond to the duration of a fast synaptic current pulse mediated by AMPA receptors, we obtain t r~1 000 ms, which is within the physiological range measured in [36]. The last term on the right hand side of equation (2) is added to the original model in [36] to include some level of stochasticity in this, otherwise, deterministic description of synaptic transmission. The inclusion of such term constitutes a simple manner of considering the stochasticity due, for instance, to the unreliability of synaptic transmission [31,32], the stochastic properties of receptor-transmitter interactions [37], the sparse connectivity of cortical circuits [38,39], or other sources of noise not yet considered (see the Discussion Section for more details). The parameter D controls the strength of this fluctuating term, and j(t) is a Gaussian white noise with zero mean and variance one.
Equations (1) and (2) constitute our minimal model of an excitatory neural network with stochastic depressing synapses. The simplifications assumed by such model allows to obtain some analytical derivations for the quantities of interest, and concretely for the probability distributions of permanence times in the up state, denoted by P(T). Bistable systems in the presence of different sources of noise have been theoretically studied in detail in many works [40][41][42][43][44]. Here, however, we have employed a probabilistic approach which is very appropriate for the computation of the distribution of permanence times. In the following, we will derive an approximate expression for P(T) within this approach. First, we obtain the potential function and the conditions in which the dynamics of the system is driven by the variable x. After that, we compute the probability distribution of ruin times of x(t) which, as we will see, leads to the probability distribution of permanence times in the up state, namely P(T).

A. The potential function
In order to compute the potential function of the dynamics (1,2) (namely W(n,x)) one can see that, for realistic values of t r , the dynamics of x is very slow compared to that of n. We therefore can write equation (1) as where we have adiabatically eliminated x from the dynamics of n.
The extrema of W are given by the solutions of the equation In the following, we choose h~Jx 0 n 0 , with n 0 : 1 2 n m and x 0 :1=(1zut r n 0 ). With this choice, one can easily check from equation (3) that the potential becomes symmetric in n around n 0 when x^x 0 . Equation (4) may have one or three solutions, depending on the slope of the hyperbolic tangent and on the value of h. In order to obtain three solutions of (4) (that is, the bistable regime) the maximal slope of the hyperbolic tangent must be large enough, concretely the condition Jxn 0 w1 must be fulfilled. In addition, the threshold term must be not too small or too large so that f (n) has three crossing points with the straight line n rather than one. This last condition can be written, as a first approach, as f (n 1 )wn 1 and f (n 2 )vn 2 , where n 1,2 are the values where the curvature of the hyperbolic tangent is maximal and minimal, respectively. The points n 1,2 can be easily computed from the third derivative of f (n): f '''(n)~{n m J 3 x 3 1{3 tanh 2 (Jxn{Jx 0 n 0 ) cosh 2 (Jxn{Jx 0 n 0 ) : By setting f '''(n)~0 we obtain Using now these values for n 1,2 , the conditions f (n 1 )wn 1 and f (n 2 )vn 2 can be written as which implies that, in order to have one maxima and two minima in W(n,x), the variable x must be in the range x 1 vxvx 2 , where From equation (8), one can see that the range of x that allows to have three extrema in the potential is Dx: The condition Dxw0 implies Jx 0 n 0 * > 1:14 which is, therefore, a sufficient condition to obtain a double well potential for some value of x. One can find, however, a small discrepancy between this approximate prediction and the actual properties of W(n,x). The discrepancy appears because we have assumed that a sufficient condition for the existence of the three fixed point solutions of equation (4) is that f (n 1 )wn 1 and f (n 2 )vn 2 , and such assumption is only approximately correct. Plotting directly the potential as a function of n reveals that the condition to obtain a double well potential for x^x 0 is Jx 0 n 0 w1, rather than Jx 0 n 0 w1:14.
Assuming that the above condition (Dxw0) is satisfied, three different shapes for the potential function W(n,x) can be found, as the figure 1A illustrates. When xvx 1 the potential function presents only one minimum, located around n^0. Similarly, for x 2 vx the potential presents also a single minimum, but now located around n^n m . Finally, for x 1 vxvx 2 the potential will take a double well shape, with the maximum being located around n^n 0 and the minima located around n^0 and n^n m , respectively.
It is worth noting that x 1 vx 0 vx 2 , with x 0 being the mean value of x. Due to this, if the range Dx is small compared with the fluctuations of x, namely s x , the potential function will spend most of the time in the regimes xvx 1 and x 2 vx, with the double well regime appearing only when the system tries to jump from one of these regimes to the other (that is, when x^x 0 ). A direct consequence of this is that the mean firing rate will basically be switching between the up and down states (that is, n^0 and n^n m ), and that this switching will be driven by the dynamics of x, as the figure 2 illustrates. Therefore, one expects that the distribution of permanence times of n in the up (down) state, becomes approximately equal to the distribution of permanence times of x in the xwx 0 (xvx 0 ) regime, as long as Dx%s x is satisfied. Due to this equivalence, in order to compute P(T) we only need to compute the distribution of permanence times of the variable x in the xwx 0 regime, denoted as P x (T).
On the other hand, it should be noted that, since x is a fraction of available neurotransmitters, its value should be kept within the range ½0,1. In practice, this means that the value of s x must not be too large, so in order to make Dx%s x one has to restrict to Dx small. In the results presented here, x remain in its realistic range of values, and imposing ad hoc restrictions in such a way that x is always within the range ½0,1 does not affect the results obtained here.

B. Distribution of permanence times
In order to compute the distribution of permanence times of x in the xwx 0 (or xvx 0 ) regime, one can assume that the firing rate takes its mean value n^n 0 in equation (2). This is a reasonable approach since x is much slower than n for realistic values of the parameters. Considering this approach, and after the rescaling z:(1zut r n 0 )x{1, equation (2) can be written as which is the equation of the Ornstein-Uhlenbeck (OU) process (see [45] for details), with t:t r =(1zut r n 0 ) being the correlation time and z 0 :z(x 0 )~0. Therefore, computing the distribution of permanence times in the up state for our system is equivalent to obtain the distribution of the so called ruin times for the OU process [46,47], which may be defined as follows: if we consider a stochastic process y(t) starting at t~t 0 from y~y 0 , the ruin time is the interval t 1 {t 0 , where t 1 is the time at which y(t) returns to y 0 for the first time. Since y(t) is a stochastic process, the ruin times are stochastic quantities which follow a certain probability distribution. The strategy employed here to calculate the distribution of ruin times is based on the relation between the ruin time and the first passage time, which is the typical time that a stochastic process needs to arrive at a certain threshold value when starting from a certain initial condition [47]. Because of the symmetry of the OU process, the distribution of ruin times are equivalent when considering excursions of the variable z in the zv0 region or in the zw0 region. If we consider excursions in the zv0 region, we can set a small positive threshold e near zero (that is, 0vE%1), in such a way that the typical ruin time will be approximately equal to the corresponding first passage time, as the figure 1B illustrates. The excursions in the region zw0 typically lead to very short first passage times (since E is too small) which we will not take into account in our calculations by considering only large enough ruin times.
The first passage time for the OU process with a small threshold E can be performed by using the relation where P(a,t a Db,t b ) is the conditional probability distribution of the OU process, and r(t) is the first passage time distribution. This equation can be solved by taking into account the following property of the Laplace transformation wheref f i (s) is the Laplace transform of f i (t). By solving the Fokker-Planck equation associated with equation (10), one can obtain the conditional probability for the OU process where Dt:t 2 {t 1 w0, and s x :D= ffiffiffiffiffi 2t p being the standard deviation of x. From expression (13), and assuming that t is large enough (more precisely, assuming that t&Dt, which is a valid hypothesis since most of the permanence times in the up state are much lower than t), one arrives at We denote f 1 (T):P(E,TD0,0) and f 2 (T{t'):P(E,TDE,t'). Employing the Laplace transformation in f 1 (T) and f 2 (T{t') the following expressions are obtained Now, taking into account the property (12) in equation (11), the expression forr r(s) iŝ Finally, for small E one can approximate E 2 z4sts 2 x^4 sts 2 x . With this approximation, one can easily perform the inverse Laplace transformation to equation (16) and obtain the distribution of first passage times for the OU process A similar expression may be obtained if one considers more classical derivations of the first passage time of the OU process (see, for instance, [48]). In order to obtain the distribution of ruin times of the OU process, one has to consider a small (but positive) value of e, which leads to r(T)*T {3=2 . The distribution of ruin times of the variable x, namely P x (T), and therefore, the distribution of permanence times in the up state, namely P(T), for our system are also given by which corresponds to a power-law probability distribution for T. Summarizing, the three following conditions must be fulfilled to obtain a power-law dependence in P(T) with exponent {3=2: N Large enough values of t r . With this condition, we ensure that the dynamics of x(t) is much slower than that of n(t).
N Large enough values of D. In particular, we must have D&2tDx, according to the condition Dx%s x and the definitions s x :D= ffiffiffiffiffi 2t p and t:t r =(1zut r n 0 ). This condition can be achieved even with very small values of D, since Dx can be arbitrarily small (by increasing J, for instance).
N The condition Jx 0 n 0 w1 must hold to ensure the existence of two well defined (up-down) states.
All these conditions may be easily achieved (up to some point) with realistic values of the model parameters, indicating that power-law distributions of the permanence times in the up state are plausible to be found in actual cortical media.

Results
As we have stated in the previous sections, equations (1)(2) govern the dynamics of our simplified neural system. A typical time series of the dynamics of this model, for the case of deterministic synapses (that is, D~0), is depicted in figure 2A. In this case, the mean firing rate of the population is characterized by a periodic switching between up and down states. This type of periodic behavior was already found and analyzed in previous theoretical studies [14,16,18] and yields bimodal histograms for the mean firing rate of the neural population (see figure 2B), as the experiments indicate [3]. However, these approaches ignore the stochastic nature of synaptic transmission, and other forms of stochasticity at the synaptic level, which seem to be crucial for information processing in neural systems [31,32,49]. Considering a certain level of synaptic stochasticity in addition to STD in our model, one obtains a qualitatively different emergent behavior, as is shown in figure 2C for D~20. The mean firing rate presents then a complex switching between up and down states, and in particular involves a high variability in the permanence times in the up state. When deterministic synapses are considered (that is, D~0) the dynamics of the mean firing rate becomes quasi periodic, as it was reported in [16,18,19], for instance. This type of dynamics naturally leads to exponential distributions for the permanence times. More precisely, for D~0 our model is similar (except for the term f(t)) to the one analyzed in [16], which shows periodic oscillations of the network mean firing rate. In our case, however, the term f(t) introduces certain level of stochasticity which turns these periodic oscillations into quasi-periodic oscillations. This leads to the exponential distributions for the permanence times in the up state. When D is increased, on the other hand, the stochasticity of the synapses leads to the appearance of power-law distributions for the duration of the up states. This behavior is shown in figure 3A, where low values of D corresponds to exponential distributions for P(T), while larger values of D give P(T)*T {3=2 as predicted by our theoretical calculations. Such power-law distributions may explain the high variability of permanence times in the up state, which has been observed in a number of in vivo experiments, such as in the cat visual cortex [7] and rat auditory cortex [21], to name a few. Interestingly, similar power-law dependences have been observed during sleep-wake transitions in vivo when one measures the distribution of permanence times in the wake state [22,23]. On the other hand, exponential-like distributions, obtained for the case of having D~0, are not able to explain this variability of the duration of up states.
By looking at the data for D~20 in figure 3A, one can observe the existence of a small deviation of the numerical results (blue points) with respect to the theoretically predicted slope (solid line) for very large values of T. Such deviation is due to the fact that the separation of timescales between the dynamics of n(t) and x(t) (a necessary condition to obtain power-law dependences in P(T)) is only approximate when considering realistic values of the parameters (and in particular, realistic values for t r ). More precisely, the approximation fails when the activity of the system falls in the occasional periods of very long permanence times in the up state (that is, for large enough values of T, comparable with t r ). In order to study the effect of the separation of timescales between n(t) and x(t), we have computed P(T) for different (increasing) values of t r while keeping fixed values for Dx and D=t r (this can be done by properly modifying J and D with t r , respectively). As a consequence of this, the only effect of increasing t r will be a clearer separation of timescales between n(t) and x(t). The results are shown in figure 3B, where one can see that larger values of t r (that is, clearer separation of timescales) lead to a displacement of the effective cut-off towards higher values of T, as expected, and a clearer power-law distribution emerges.
It is worth noting that the appearance of an effective cut-off in T for realistic conditions does not represent an unrealistic feature of the model, but rather it constitutes a prediction about the effective range of permanence times which are expected to occur in actual neural systems. Indeed, for realistic values of the parameters, our results predict permanence times in the up state up to *1000 ms, which is the maximum permanence time observed in experimental realizations [7]. Larger permanence times in the up state (of about 10 seconds, for instance) should be expected to appear only as a consequence of input driven mechanisms (such as persistent activity associated with working memory tasks [50,51]), and not as a consequence of spontaneous transitions between different voltage levels, which are the matter of interest in this work.
For a better characterization of the dynamics of the system, one can use, for instance, other statistical magnitudes such as the autocorrelation function C(t') of n, which can be defined as Here, S Á Á ÁT indicates a temporal average. The autocorrelation function is depicted in figure 4A for the case of deterministic depressing synapses (D~0) and stochastic depressing synapses (D~20). C(t') presents, for D~0, two well located peaks at t'^+200, which indicates a strong periodicity of the time series (as can be seen in figure 2A). On the contrary, the inclusion of a certain level of intrinsic stochasticity in the dynamics of x introduces more pronounced temporal correlations in the dynamics of the system. This fact reflects the existence of long permanence stays in the up state, which occurs with more probability for high enough values of D, as we have already discussed.
The spectral properties of the dynamics can be analyzed as well, via the power spectrum defined as As one could expect, the power spectrum of the case D~0 presents a pronounced peak around a certain frequency, which in the particular case presented in the figure 4B is f *5 : 10 {3 . The power spectrum for higher values of D shows however different properties than the case D~0. For instance, the figure 4B (which considers D~20) indicates an approximated power-law behavior for the power spectrum, F (f )*f {b with b^1:7. This scale-free dependence can be understood by considering that, if P(T) is algebraic with exponent c, the corresponding power spectrum becomes also algebraic with exponent b, where the equation czb~3 relates both exponents [44]. In our particular case, since c^1:5, one obtains a theoretical prediction of b^1:5 for the exponent of the power spectrum. The theoretical relation between P(T) and F (f ) exposed above, however, is only valid under the so called single interval approximation, which implies that the integration variable t in equation (20) is smaller than the permanence time T (see [44] for details). This condition does not strictly hold for our system (where T ranges over several scales), and therefore it may introduce deviations in the theoretically predicted value of b (which is around b^1:5) with respect to the value found in simulations (of around b^1:7).
Besides the level of synaptic stochasticity, i.e. D, other parameters of the model could have an important effect on the dynamics as well. The parameter d, for instance, controls the level of stochasticity of the dynamics of n, and therefore one should expect that increasing its value could strongly influence the probability distribution P(T). This is shown in figure 5A, where an increase of d disrupts the appearance of power-law dependences, and exponential distributions appear instead. This change in P(T) is due to the fact that high levels of the additive noise d make the system to jump more frequently from one state to the other, and therefore long stays in the up state (and thus distributions with long power-law tails) rarely occur.
The parameters involving the dynamics of x also affect the probability distributions P(T). The parameter u, for instance, is responsible for the modulation of x via the mean firing rate n (see equation (2)), and therefore it can influence both the dynamics of x and n. As one may see in figure 5B, when u takes low values a bump in P(T) emerges for high T. Such deviation from the power-law dependence indicates that long stays in the up state occur more frequently than in the power-law case. Attending at equation (2), one can see that an increase of the mean firing rate n decreases the variable x via the parameter u. Therefore, if u takes lower values the decrement of x will be smaller. As a consequence, the stays of x in the x 0 %x regime (see Methods Section) will last longer, and the stays of the system in the up state will also last longer, causing the observed deviation from the power-law tendency. It should be noted, however, that the values of u which allow the appearance of power-law dependences in P(T) for our model agree with the values of u measured in actual cortical media where up and down transitions are observed [36]. We have also analyzed in detail the effect that varying t r has on the probability distribution of permanence times. Note that, contrarily to the previous study presented above, we have now varied the parameter t r while all the other parameters are kept fixed. This implies that the modification of t r will now have an effect on the separation of timescales between n(t) and x(t), but also on the concrete value of Dx and on the amplitude of the noisy term of equation (2) (namely D=t r ). The results are shown in figure 5C, where one can distinguish three different regimes as a function of the particular value of t r . For low t r (red region in the figure), the probability distributions show an exponential decay for large permanence times. The reason for this decay is that, for low t r , the variable x does not perform long excursions in the region x 0 %x (see Methods Section), and therefore the probability to have large values of T decreases and the power law behavior for P(T) is not obtained. As t r is increased, long excursions for x begin to occur, and we obtain a power law behavior P(T)*T {3=2 (green region in the figure). Finally, one can appreciate that, for even larger values of t r (blue region in the figure), the probability distribution of permanence times in the up state presents a power law dependence P(T)*T {c(tr) with c(t r )w3=2, being an increasing function of t r . Such dependence can not be explained by our previous theoretical predictions, based in the assumption that the system is in the bistable regime, and deserves a detailed analysis which will be exposed in the next section.

Further analysis
In the Methods Section, we established several conditions which had to be fulfilled in order to obtain power law dependences for P(T). In particular, our previous analysis indicates that the condition Jx 0 n 0 w1 must hold in order to have a potential function W(n,x) with three extrema (bistable regime). However, as we will see in the following, power law expressions for P(T) may appear even if the potential function has only one extremum in n (concretely, one minimum), although the origin of such power law distributions is different from the one considered in previous sections, as we will see.
When Jx 0 n 0 v1 (which occurs for J%1 or t r &1, for instance), the potential function W(n,x) has only one minimum in n, whose location strongly depends on x. An approximated expression for the location of this minimum as a function of x can be obtained by expanding the hyperbolic tangent of the fixed point expression of n(t) (see equation (4)) around its argument (which is small in this limit), yielding where n min is the value of n which corresponds to the minimum of the potential function. Therefore as x varies around x 0 , the location of the minimum of the potential n min also varies in the same way around n 0 . As an example, time series of both n and x are shown in figure 6A for a given set of parameters which satisfies Jx 0 n 0 v1. In this time series, the variable n fluctuates around the value n min , which is fully determined by x (that is, the variable n becomes a slave variable of x). The predictions of equation (21) agree approximately well with simulations and with the numerical evaluation of the fixed points of equation (1), as the figure 6B shows. Since n behaves now as a stochastic variable which does not present a clear bistable dynamics, the numerical computation of the distribution of the permanence times will depend on the exact value of n above which the system is considered to be in the up state. As we have seen before, this threshold value takes the form gn m (see caption of figure 3), where usually g may take a value between 0:6 and 0:9. While the results presented for Jx 0 n 0 w1 (that is, the bistable regime) are quite robust for different values of g, in the regime Jx 0 n 0 this parameter has indeed some effect on P(T), which indicates the difficulty to accurately analyze the up and down dynamics in this case.
In figure 7A, one observes that the distribution P(T) shows also a power law behavior P(T)*T {c for g~0:75 and different values of D, for a set of parameter values which satisfies Jx 0 n 0 v1 (that is the monostable regime). The concrete value of c depends strongly on D and it has also a weaker dependence with g, as the figure 7B illustrates. This type of power-law behavior appearing in the monostable regime corresponds to the blue region in figure 5C, as well.
It is worth noting that actual recordings of up and down transitions does not present a clear distinction between up and down states, and several nontrivial methods are commonly employed to discriminate between both states [52]. Therefore, the results found for the regime Jx 0 n 0 v1 could indeed reflect the behavior of actual cortical up-down transitions, showing power law dependences in P(T) with cw3=2 and indicating that the concrete nature of the transitions is a synaptic-driven monostable dynamics.
For a complete characterization of the model, one can summarize all the observed behaviors in a phase plot such as the one presented in figure 8A. A total of four different behaviors can be found in the (t r , D) space. The first one concerns the dynamics of n whose permanence times in the up state follows an exponential distribution (labeled as ''E'' in the figure). If the noise amplitude D is sufficiently high, one can increase the value of t r to reach the regime ''C'', in which the dependence P(T)*T {1:5 is obtained. By increasing t r even more, the probability distribution P(T) takes the form *T {c , with cw1:5 (regime denoted by ''S''), as we have already seen in figure 6. Finally, we also observe that when the depression time scale is not large enough (and D v * 3), a regime of quasi-periodic time series of n is obtained, with a welldefined duration of up states (regime denoted by ''P''). The lines  (4)). The inset shows the situation in which the system shows a bistable dynamics, analyzed in the previous section. (C) The potential function as a function of n for different values of x. One can appreciate the existence of only one minimum, whose location is controlled by x. (D) Histograms of the mean firing rate of the system for different values of J. For the cases showed in this panel, the condition Jx 0 n 0 v1 is only satisfied for the case J~0:55. For all panels, u~0:6, t r~1 000, D~20, d~0:3, n m~5 : 10 {3 , and J~0:55 V unless specifically specified. doi:10.1371/journal.pone.0013651.g006 between the different regimes have been obtained by visual inspection of P(T) for different values of t r and D. In particular, the regime ''P'' is characterized by the appearance of a bump in the probability distribution for some value of T (which reflects a preferred duration of the up state), and the existence of such bump has been used as a criterion to distinguish between regimes ''P'' and ''E''. Similarly, we assumed that the regimes ''C'' and ''S'' correspond to the situation in which a power-law behavior that extends for two decades or more is found for P(T). Such criterion, together with an estimation of the slope of the power-law via standard Levenberg-Marquardt fitting algorithms, allows to distinguish between regimes ''E'', ''C'' and ''S''.
It must be clarified, however, that actual up and down cortical transitions present most likely a richer repertoire of dynamical regimes than the one obtained with our simplified model. It is known, for instance, that attractor neural networks with dynamic synapses may exhibit different dynamics corresponding to memory, non-memory and switching regimes [18,19]. In this work, we have extensively explored different regimes of switching behavior, and its implications for the up and down dynamics observed in the cortex. The memory and non-memory regimes, however, can be also found in our simplified model by assuming that D, d?0. After taking these limits, the system will be in the memory regime if the potential function W(n,x) is bistable, or in the non-memory regime if W(n,x) is monostable.

Discussion
We have shown that the experimentally observed large fluctuations in up and down permanence times can be explained as the result of sufficiently noisy dynamical synapses with sufficiently large recovery times. Our study suggests that a power-law distribution for these permanence times may emerge as a consequence of these two ingredients. Static synapses cannot account for this behavior, nor can dynamical synapses in the absence of noise.
The origin of up and down cortical transitions is still unclear, although different factors that may influence their occurrence have been recently reported. It is known, for instance, that inhibitory GABAergic currents strongly contribute to the temporal coding and spike timing precision of cortical networks during up states of activity [3,53,54]. Several modeling studies also show the relevance of inhibitory interneurons in the generation of many types of oscillations in the brain (see for instance [55]). However, other studies indicate that most of the main features of up and down transitions depends strongly on synaptic plasticity mechanisms, both of long-term and short-term ones [16,56], and that the transitions appear even in the absence of inhibition [16]. In this work we have made the common assumption that the effects of inhibition can be treated as additive and can be incorporated in the threshold of the neuron. This is known to be a valid approximation in mean field neural network analysis, but may fail when precise timing and details of the dynamical aspects of the neuron affect the inhibition [57,58].
Regarding to synaptic characteristics, recent works show that synaptic fluctuations could have an important role in the generation of transitions between up and down states [14,59,60]. Since our model introduces stochasticity in the synaptic dynamics in a highly simplified manner, however, the last term in equation (2) should not be associated only with ureliability in synaptic transmission. Indeed, we have assumed that other sources of stochasticity may be contributing to this fluctuating term in the mean-field quantities n(t) and x(t). For instance, it is widely known that connectivity in actual cortical media is highly sparse. Such feature implies that, in order to obtain the mean-field quantity x(t), the average over synapses must be performed over a number of C synapses, with C ranging over 100*1000 connections per neuron [38]. In this situation, the fluctuations of x(t) would be of order 1= ffiffiffiffi C p , which leads to a range of 0:1*0:03 for the values given above. As we have seen, our results state that a value of D=t r~0 :02 is enough to obtain power-law distributions (see figure 3), which lies within this range. Therefore, topology-induced fluctuations constitute an important source of stochasticity which could be responsible of the appearance of power-law distributions in P(T). Other sources of stochasticity at synaptic level, such as the stochastic properties of receptor-transmitter interactions, may also contribute to the last term of equation (2). Moreover, the low activity rates typical from cortical media lead to a poor timeaveraging of the incoming input, and therefore the fluctuations at the postsynaptic level will be large at these short-time scales (of the order of the typical synaptic integration time constant).
On the other hand, the amplitude of the noisy term, D=t r , does not need to be very high to induce the appearance of power-law distributions in P(T). As we have stated above, a sparse connectivity already induces a level of stochasticity which is within the desired range, for instance. Furthermore, the noisy term could even be arbitrarily small: attending to our theoretical predictions, a necessary condition to have power-law distributions is that fluctuations of x(t) must be much larger than Dx (see Methods Section). Since Dx may be lowered to arbitrary levels (by increasing J, for instance), even a small noisy term in the dynamics of x(t) may induce power-law distributions.
It is also known that short-term synaptic mechanisms, such as short-term depression and facilitation, usually play a role in the efficient processing of information. In particular, they may be relevant in many tasks, such as in signal detection and coding [29,[61][62][63] or switching between different activity patterns previously stored [19,64]. However, their role on the transitions between cortical states has been pointed out only by a few studies [15,16,65], and their possible effects on the statistics of the transitions, which is the focus of our work, have been ignored. To the best of our knowledge, the present study is the first one which analyzes, even in a simplified manner, the strong effect of synaptic stochasticity-in a general sense-and dynamic synapses in the statistics of the up and down transitions. The possible role of other short-term synaptic mechanisms, such as STF, has not been addressed yet and constitutes a interesting issue still open.
In our analysis we assumed that the dynamics is symmetric in the up and down states. This is in contradiction with experimental evidences [66] which shows that power-law distributions are obtained for permanence times in the up state, while permanence times in the down state are exponentially distributed. However, this discrepancy disappears when one considers a more realistic transduction function which gives an asymmetric potential for the dynamics, and as a consequence the up-down symmetry is broken. More detailed studies considering, for instance, some of the biologically realistic aspects discussed above, should be performed to test our predictions. In particular, a more elaborated study considering realistic neuron models (such as Hodgkin-Huxley model [67]) and stochastic STD models (see [29,49], for instance) is necessary, as well as more detailed experimental studies which may confirm our predictions.
From a general point of view, evidences of criticality have been recently found in an increasing number of neural systems, such as in the functional connectivity of the living human brain [24], in critical avalanches of neuronal activity [25], or in sleep-wake transitions [23], to name a few. According to the results presented in this work, transitions between up and down cortical states could also present some relevant properties typical of systems at criticality. Some of these properties have been already measured in experiments, such as a high sensitivity of the system to external stimuli [8], or the presence of power-law dependences in the power spectra of the neural dynamics [53].
It is worth noting that other kind of probability distributions for P(T), such as a log-normal distribution, could also satisfactorily explain the irregularity in the up states found in experiments. Our study shows the importance of some biophysical factors, such as the neurotransmitter recovery time and the inherent synaptic stochasticity, and predicts a power-law dependence on P(T) as a consequence of such factors. However, further study is needed to investigate other mechanisms, not taken in account in this work, which could influence the permanence times in the up state. In a more general sense, our results may proportionate a new perspective of the phenomena of up and down transitions (and a theoretical framework) that could serve to conciliate the main experimental findings, and that could help for a deep understanding of this complex dynamics of the brain activity. Phase plot which shows the different behaviors found in our system. These behaviors corresponds to time series of n for which permanence times in the up state follow an exponential distribution (E), a power-law distribution P(T)*T {c with c~3=2 (C), or a power-law distribution with cw3=2 (S). In addition, a phase with a well-defined duration of the up state is found (P). In panel (B) some of these behaviors are depicted. From top to bottom one can see situations P, E and C. Other parameters are J~1:1 V , u~0:6, d~0:3 and n m~5 : 10 {3 . doi:10.1371/journal.pone.0013651.g008