UP-DOWN cortical dynamics reflect state transitions in a bistable balanced network

In the idling brain, neuronal circuits often exhibit transitions between periods of sustained firing (UP state) and quiescence (DOWN state). Although these dynamics occur across multiple areas and behavioral conditions, the underlying mechanisms remain unclear. Here we analyze spontaneous population activity from the somatosensory cortex of urethane-anesthetized rats. We find that UP and DOWN periods are variable (i.e. non-rhythmic) and that the population rate shows no significant decay during UP periods. We build a network model of excitatory (E) and inhibitory (I) neurons that exhibits a new bistability between a quiescent state and a balanced state of arbitrarily low rate. Fluctuating inputs trigger state transitions. Adaptation in E cells paradoxically causes a marginal decay of E-rate but a marked decay of I-rate, a signature of balanced bistability that we validate experimentally. Our findings provide evidence of a bistable balanced network that exhibits non-rhythmic state transitions when the brain rests.


Introduction
A ubiquitous pattern of spontaneous cortical activity during synchronized brain states consists in the 36 alternation between periods of tonic firing (UP states) and periods of quiescence (DOWN states) 37 (Luczak et al., 2007;Steriade et al., 1993a;Timofeev et al., 2001). Cortical UP and DOWN dynamics 38 take place during slow-wave-sleep (SWS) (Steriade et al., 1993a) and can also be induced by a 39 number of anesthetics (Steriade et al., 1993a). More recently however, similar synchronous cortical 40 dynamics have been observed not only in awake rodents during quiescence (Luczak et al., 2007;41 Petersen et al., 2003), but also in animals performing a perceptual task, both rodents 42 (Sachidhanandam et al., 2013;Vyazovskiy et al., 2011) and monkeys (Engel et al., 2013). 43 Spontaneous activity resembling UP and DOWN states has been found in cortical slices in vitro 44 (Cossart et al., 2003;Fanselow and Connors, 2010;Mann et al., 2009; Sanchez-Vives and 45 McCormick, 2000), in slabs (Timofeev et al., 2000) and in vivo under extensive thalamic lesions 46 (Steriade et al., 1993b). This has led to suggest that the underlying mechanism had an intracortical 47 origin. In such scenario, the standard hypothesis postulates that during UP periods a fatigue cellular 48 mechanism -e.g. spike frequency adaptation or synaptic short-term depression -decreases network 49 excitability until the state of tonic firing can no longer be sustained and the cortical network switches 50 into a DOWN state (Contreras et al., 1996;Sanchez-Vives and McCormick, 2000). During DOWN 51 periods, in the absence of firing, the fatigue variables recover until the circuit becomes self-excitable 52 and autonomously transitions into an UP state (Cunningham et   sequence by a lag k (Fig. 2C, right). The cross-correlation ( , + ) displayed consistently non-145 zero values for k=0 and k=1 (mean±SD: 0.21±0.09, 0.17±0.09, respectively; significant cross-146 correlation in 6/7 animals, P<0.05 permutation test), whereas remained close to zero for the rest of 147 lags, showing that period duration relationship is limited to adjacent periods (Fig 2C-D). The positive 148 correlation between adjacent periods was not due to slow changes in their duration, as we corrected 149 by the correlation obtained from surrogate D-U sequences obtained from shuffling the original 150 sequence within 30 second windows (see Methods). Positive correlations between consecutive 151 periods of activity and silence can be generated when fluctuation driven transitions are combined with an adaptive process such as activity-dependent adaptation currents (Lim and Rinzel, 2010;Tabak et al., 2000): if a fluctuation terminates a U period prematurely without much build-up in adaptation, the 154 consecutive D period also tends to be shorter as there is little adaptation to recover from. However, a 155 major role of adaptation currents in dictating UP-DOWN dynamics (Compte et al., 2003b)

Spiking activity during UP and DOWN states 174
We next searched for more direct evidence of an adaptive process by examining the time course of 175 the population firing rate R(t) during U and D periods (see Fig. 1C; see Methods). The mean firing rate 176 in U periods was low (mean±SD: 3.72±0.9 spikes/s, n=7). Moreover, D periods displayed occasional 177 spiking (mean±SD rate 0.018±0.007 spikes/s; see e.g. Fig. 3A-B and Supp. Fig. S2), in contrast with 178 the idea that DOWN periods do not display spiking activity (Chauvette et al., 2010), but see (Compte R(t) during Us and an increase during Ds, and this impact on R(t) dynamics should be more evident 181 during longer periods due to a larger accumulation (during Us) or recovery (during Ds) of the 182 adaptation. For each experiment, we aligned the rate R(t) at the DOWN-to-UP (DU) and UP-to-DOWN 183 (UD) transition times (Fig. 3A). We then computed the average rates ( ) and ( ) across all DU 184 and UD transitions, respectively, with = 0 representing the transition time (   Fig.1). (B-C) Each U period is aligned at the DU (B, left) and UD (B, right) transition times in order referred to as U-onset and U-offset windows, respectively. The windows were chosen 50 ms away 206 from τ = 0 to avoid the transient change due to the state transitions ( Fig. 3C-D). We found no 207 significant mean difference between population average rate at U-onset and U-offset windows across 208 our experiments (mean±SD onset minus offset population rate 0.04 ± 0.40 spikes/s, P=1, Wilcoxon 209 signed rank, n=7 animals). The equivalent analysis performed on D periods yielded a small but 210 significant mean increase in the population rate between the D-onset and D-offset windows (mean±SD 211 -0.014 ± 0.013 spikes/s, P=0.047, Wilcoxon signed rank test). To examine in more detail the lack of 212 population rate change during Us, we looked at the modulation of individual neuron rates normalized 213 by the overall temporal average of each unit (Fig. 3E). We found that the change between U-onset and 214 U-offset averaged across all our neurons (n=448 cells) was not significantly different from zero ( Fig.  215 3E right, mean±SD of the onset vs offset difference of normalized rates 0.057 ± 1.163, P=0.12, 216 Wilcoxon signed rank test) but that the recovery during D periods was significant (Fig. 3E Cowan, 1972). The excitatory-inhibitory (EI) network model described the dynamics of the mean 229 instantaneous rates r E and r I of each population in the presence of fluctuating external inputs. In 230 addition, the E population included an adaptation mechanism, an additive hyperpolarizing current a 231 that grew linearly with the rate r E (Fig. 4A; see Methods). We did not consider adaptation in the 232 inhibitory population for simplicity, and because inhibitory neurons show little or no spike-frequency 233 adaptation when depolarized with injected current (McCormick et al., 1985). Our aim was to search for 234 a regime in which, in the absence of adaptation and external input fluctuations, the network exhibited 235 bistability between a quiescent (D) and a low-rate state (U) fixed point. Although bistability in low-236 dimensional EI networks has been described since the seminal work of Wilson & Cowan (1972), 237 previous models primarily sought to explain bistability between a low-rate and a high-rate state, and 238 exploited the combination of expansive and contractive non-linearities produced by the transfer 239 function (Amit and Brunel, 1997;Renart et al., 2007; Wilson and Cowan, 1972), short-term synaptic robustly obtained solely using the expansive nonlinearity of the transfer function caused by the spiking 243 threshold. Given this, we choose the simplest possible transfer function with a threshold: a threshold-244 linear function (Fig.4B, see Methods). Our choice to only use an expansive threshold non-linearity 245 constrained strongly the way in which the network could exhibit bistability as can be deduced by 246 plotting the nullclines of the rates r E and r I (Fig. 4C): only when the I nullcline was shifted to the right 247 248 249 250 251 and had a larger slope than the E nullcline, the system exhibited two stable attractors (Eq. 20 in Methods). This configuration of the nullclines was readily obtained by setting the threshold and the 267 gain of the I transfer function larger than those of the E transfer function (Fig. 4B), a distinctive feature 268 previously reported when intracellularly characterizing the f -I curve of pyramidal and fast spiking 269 interneurons in the absence of background synaptic activity (Cruikshank et al., 2007;Schiff and 270 Reyes, 2012). This novel bistable regime yielded a quiescent D state, and arbitrarily low firing rates for 271 both E and I populations during U states, depending on the values of the thresholds and the synaptic 272 weights (Fig. 4C). This is remarkable as in most bistable network models the rate of the sustained 273 activity state is constrained to be above certain lower bound (see Discussion). Moreover, in this 274 bistable regime, the U state is an inhibition-stabilized state, a network dynamical condition in which the 275 excitatory feedback is so strong that would alone be unstable, but is balanced with fast and strong 276 inhibitory feedback to maintain the rates stable (Ozeki et al., 2009;Tsodyks et al., 1997) (see 277

Methods). 278
There are two ways in which transitions between U and D states can occur. On the one hand, 279 transitions could be driven by external input fluctuations, which were modeled as a stochastic process 280 with zero mean and short time constant (Fig. 4D). This fluctuating input reflected either afferents 281 coming from other brain areas whose neuronal activity was stochastic and uncorrelated with the 282 cortical circuit internal dynamics or the stochasticity of the spiking happening during U periods which 283 was not captured by the dynamics of the rates (Holcman and Tsodyks, 2006;Lim and Rinzel, 2010). 284 On the other hand, in the absence of fluctuations, state transitions could also occur solely driven by 285 adaptation currents (Fig. 4E). Because the adaptation time constant was much longer than the time 286 constants of the E and I rates, the dynamics of the rates r E (t) and r I (t) relaxing rapidly to their steady-287 state can be decoupled from the slow changes in a(t) (Latham et al., 2000;Rinzel and Lee, 1987).

288
The network dynamics can be described in the phase plane (r E (t), r I (t)) with variations in a(t) causing a 289 displacement of the E-nullcline. In particular, during U periods the build-up in adaptation produced a 290 downward displacement of the E-nullcline (Fig. 4E). If adaptation strength β was sufficiently large the 291 displacement increased until the U state was no longer a fixed point and the network transitioned to 292 the only stable fixed point D. Recovery of adaptation during D periods shifted the E-nullcline upwards 293 until the D state disappeared and there was a transition to the U state ( Fig 4E). In this limit cycle 294 regime the network exhibited an oscillatory behavior with a frequency close to the inverse of the 295 adaptation recovery time constant. When the two types of transitions are combined, two types of 296 stability in U and D states can be distinguished: (1) metastable, referred to a state that was stable to 297 the dynamics of both the rates and the adaptation but could transition away due to input fluctuations; 298 (2) quasi-stable, referred to a state that was stable for the fast rate dynamics but unstable for the slow To quantify the relative impact of activity fluctuations and adaptation in causing U-D transitions in the 303 data, we compared the dynamics of the model for different adaptation strengths β and different values 304 of the E effective threshold θ E . The (θ E , β) plane was divided into four regions with UD alternations, 305 corresponding to the four combinations of metastability and quasi-stability (Fig. 5A). Since only 306 metastable states tend to give exponentially distributed durations with CV~1, the large variability found 307 in both U and D durations ( show much difference whereas the average r I (t) showed a larger decrease over the U period (Fig. 6A). 358 Thus, although only the E and not the I population included intrinsic adaptation mechanisms, it was 359 r I (t) the one that exhibited the most pronounced decay during U periods. This was a direct 360 consequence of the specific conditions that gave rise to bistability in our model: the difference in 361 thresholds, i.e. θ I > θ E , and the fact that the I-nullcline has a higher slope than the E-nullcline (Eq. 21 362 in Methods). These features imposed that as adaptation built up during U periods, the downward 363 displacement of the E-nullcline caused a greater decrease in r I (t) than in r E (t) (compare "decay" colored 364 bands in Fig. 6B). With this arrangement the drop in r E (t) could be made arbitrarily small by increasing 365 the slope of the I-nullcline (Fig. 6B). During D periods the average r E (t) did show a substantial increase 366 due to the recovery of adaptation, whereas the r I (t) did not. This was because in the D state, the 367 quiescent network behaved as isolated neurons reflecting the dynamics of intrinsic adaptation which 368 was only present on E cells. In sum, if the majority of the neurons that we recorded experimentally 369 were excitatory, the model could explain why adaptation currents did not cause a significant decrease 370 in the average rate during U periods ( periods, E cells also showed a significant increase in rate (Wilcoxon signed rank test P=0.0092), just 384 like that observed in the whole cell population, whereas no rate change was found in I cells (not 385 shown). Although these changes observed during D periods were also predicted by the model, 386 properly testing the significance of this interaction would require a larger data set with more I cells. The 387 validation of the prediction on the counter-intuitive emergent dynamics of E and I rates during U 388 periods strongly suggests that the mechanism dissected by the model drives the spontaneous 389 dynamics of the cortical circuitry under urethane anesthesia. 390 392 393 (mean ~150 ms) and slightly more regular (CV~0.5) (unpublished observations). Such an asymmetry 441 in the duration and irregularity of U-D periods can be easily reproduced in our model by choosing 442 parameters in the mixed region where U is meta-stable and D is quasi-stable (Fig. 5A light orange). 443 In addition, we found non-zero correlations between consecutive D-U and U-D period durations, a 444 feature that was not observed previously in similar experimental conditions (Stern et al., 1997). An alternative mechanism to generate bistability between UP and DOWN states has been the 483 shunting or divisive effect of inhibitory synaptic conductances, a mechanism that can produce non-484 monotonic transfer functions and yield bistability between a zero rate state and a state of very low rate 485 Our model proposes a more parsimonious mechanism underlying UP-DOWN bistability: the 498 ubiquitous expansive threshold non-linearity of the transfer function plus the asymmetry in threshold 499 (θ I > θ E ) and gain (larger for I than E cells). We used a threshold-linear function for simplicity but other 500 more realistic choices (e.g. threshold-quadratic) produced the same qualitative results. The threshold 501 asymmetry is supported by in vitro patch clamp experiments revealing that firing threshold of inhibitory 502 fast-spiking neurons, measured as the lowest injected current causing spike firing, is higher than that circuits alternate between a quiescent state with no activity and a state of balanced irregular and that the rate of coordinated transitions to the DOWN state predicts the magnitude of correlations 514 across different brain states (Mochol et al., 2015). 515 A direct implication of the bistability obtained in our model was that intrinsic adaptation of excitatory 516 neurons (McCormick et al., 1985) did not cause a noticeable decrease in r E during the UP periods but 517 instead produced a significant decay in the inhibitory rate r I . We confirmed this prediction in our data 518 ( Fig. 6C-D). Interestingly, the same effect was also observed in ketamine anesthetized animals from 519 both extracellular (Luczak and Barthó, 2012) and intracellular recordings resolving synaptic 520 conductances (Haider et al., 2006). During DOWN periods in contrast, the network is not in a balanced 521 state and recovery from adaptation caused a significant increase in the rate of putative excitatory 522 neurons, as predicted by the model. In sum, our results present the first EI network model with linearly 523 summed inputs exhibiting bistability between a quiescent state and a balanced state with arbitrary low 524 rate. 525 526

The role and origin of fluctuations in UP-DOWN switching 527
Our findings stress the role of input fluctuations inducing transitions between the UP and DOWN 528 network attractors because noise-induced alternations generate periods with large variability as found 529 in the data ( Fig. 2A-B). Adaptation was also necessary to introduce positive serial correlations and to 530 reproduce the observed gamma-like UP-DOWN distributions (compare Fig. 2A using a spiking EI network to produce UP-DOWN alternations shows that to cause noise-driven 543 transitions from a quiescent state synaptic inputs into each neuron need to be correlated and non-544 Gaussian (not shown). Gaussian uncorrelated inputs must generate an unrealistically high DOWN 545 firing rate in order to yield a measurable escape probability to the UP state (because escape requires 546 the synchronous occurrence of multiple independent neuronal discharges). That is why synaptic noise, (Supp. Fig. 2B). This reasoning favors instead synchronous external input kicks as the inducers of Adult, male Sprague-Dawley rats (250-400g) were anesthetized with urethane (1.5 g/kg) and 612 supplemental doses of 0.15 g/kg were given when necessary after several hours since the initial dose. 613 We also used an initial dose of Ketamine (15-25 mg/kg) before the surgery to induce the anesthetized 614 state quickly. We then performed a craniotomy over the somatosensory cortex, whose position was 615 determined using stereotaxic coordinates. Next 32 or 64 channels silicon microelectrodes 616 (2) 674 where: 675 and equivalently for < > and ( ). We controlled whether variability in U i was produced by slow 677 drifts by computing CV 2 a measure of variability not contaminated by non-stationarities of the data The serial correlation between and + , with k setting the lag in the U-D series, e.g. k =0 (k = 1) 680 refers to the immediately previous (consecutive) DOWN period, was quantified with the Pearson 681 correlation coefficient defined as: 682 where the covariance was defined as: 684 Values of and differing more than 3 standard deviations from the mean were discarded from the 687 correlation analysis (circles in Fig 2C). To remove correlations between and produced by slow 688 drifts in the durations we used resampling methods developed to remove slow correlations among 689 spike trains (Amarasingham et al., 2012). We generated the -ℎ shuffled series of U periods { ̂} The autocorrelogram of the instantaneous population rate was defined as: 739 , for > 0 (11) 740 with the sum in t running over the L time bins in a window of size W. The average < ( ) > and 741 variance were performed across time in the same window. To avoid averaging out a rhythmic structure 742 in the instantaneous population rate due to slow drift in the oscillation frequency, we computed ( ) 743 in small windows W=20 s thus having a more instantaneous estimate of the temporal structure. With 744 the normalization used, the autocorrelograms give ( = 0) = 1 and the values with > 0 can be 745 interpreted as the Pearson correlation between the population rate at time t and the population rate at 746 time + (Fig. 1D). 747 and U-offset (Fig 3 and 6), we computed for each neuron the mean of ( ) over the window 750 = (50,200) (U-onset) and the mean of ( ) over the window = (−200, −50) (U-offset). We 751 positioned the windows 50 ms away of the DU and UD transitions in order to preclude the possibility of 752 contamination in the mean rate estimations due to possible misalignments from the U and D period 753 detections. In the averaging we used U and D periods longer than 0.5 s, so that onset and offset 754 windows were always non-overlapping. Equivalent D-onset and D-offset windows were defined in 755 order to compare individual rates during D periods. To make the distribution of mean rates across the 756 cell population Gaussian, we normalized each the rates ( ) and ( ) by the overall time-757 averaged rate of the neuron =< ( ) > finally obtaining onset and offset-aligned normalized 758 averaged rates (e.g. ( )/ ). Despite this normalization, the distribution of the normalized rates in 759 the D-onset and D-offset was non-Gaussian (most neurons fired no spikes). Thus we used the non-760 parametric two-sided Wilcoxon signed rank test to compare onset and offset rates (Fig. 3E). To test 761 the rates changes during U periods in E and I neurons we used a four-way mixed-effects ANOVA with 762 fixed factors onset/offset, E/I and random factors neuron index and animal. We compared the 763 distribution of normalized averaged rate difference at the U-onset minus the U -offset (Fig. 3E right,  764 dark gray histogram) with a distribution obtained from the same neurons but randomly shuffling the 765 onset and offset labels of the spike counts but preserving trial and neuron indices (Fig. 3E  increase that compensate to yield no significant difference on the population averaged rate. The same 770 procedure was followed with the normalized rates in the D-onset and D-offset but the limited number 771 of non-zero spike counts limited the analysis yielding inconclusive results (Fig. 3E left). 772

Computational Rate Model. 773
We built a model describing the rate dynamics of an excitatory ( ) and inhibitory population ( ) 774 recurrently connected that received external inputs (Wilson and Cowan, 1972). In addition, the 775 excitatory population had an additive negative feedback term, ( ), representing the firing adaptation 776 experienced by excitatory cells (McCormick et al., 1985). The model dynamics were given by: 777 The time constants of the rates were = 10 ms and = 2 ms, while the adaptation time constant was from Y to X, were = 5, = 1, = 10, = 0.5 s. Because we are modeling low rates, the adaptation 784 grows linearly with with strength = 0.5 s. The fluctuating part of the external inputs ( ) was 785 modeled as two independent Ornstein-Uhlenbeck processes with zero mean, standard deviation = 786 3.5 and time constant 1 ms for both E and I populations. Because population averaged firing rates 787 during spontaneous activity fell in the range 0-10 spikes/s, we modeled the transfer functions φ X as 788 threshold-linear functions: 789 Fixed points and stability. We start by characterizing the dynamics of the system in the absence of 813 noise. Assuming that the rates evolve much faster than the adaptation, i.e.
, ≪ , one can 814 partition the dynamics of the full system into (1) the dynamics of the rates assuming adaptation is 815 constant, (2) the slow evolution of adaptation assuming the rates are constantly at equilibrium. Thus, 816 the equations of the nullclines of the 2D rate dynamics at fixed a, can be obtained from the 2D system 817 given by Eqs. 12-13. The nullclines of this reduced 2D system are obtained by setting its left hand side 818 to zero: 819 The conditions for this UP state solution to exist are derived from imposing that the right hand side of 836 Eqs. 18-19 is positive. The stability of this solution (Eq. 21 below) implies that the determinant | | is 837 positive and that if is positive, then is also positive. Thus, provided the stability (Eqs. 21-22), the 838 only condition for the solution to exist is that the right hand side of Eq. 19 is positive: 839 Given the separation of time scales described above, this fixed point is stable if the eigenvalues of the 841 matrix of coefficients of Eqs. 16 and 17 without the term a (that we assume is constant) have all 842 negative real part. Because the coefficients matrix is 2 x 2, this is equivalent to impose that the 843 determinant of the matrix has a positive determinant and a negative trace. These conditions yield the 844 following inequalities, respectively: 845 ′ ′ < (21) 846 ( + 1) < ( + 1) (22) 847 Equation 21 is equivalent to the condition that the I-nullclines of the 2D reduced system has a larger 848 slope than the E-nullcline. From the U existence condition in Eq. 20 and D stability condition, it can 849 also be derived that ′ > 0, implying that at fixed inhibition, the E-subnetwork would be unstable (i.e. 850 slope of the E-nullcline is positive). In sum, the conditions for the existence of two stable U and D 851 states imply that the U state would be unstable in the absence of feedback inhibition but the strength 852 of feedback inhibition is sufficient to stabilize it. These are precisely the conditions that define an the ( , β )-plane (Fig. 5A). In the absence of noise, given that θ I ≥ 0, a stable D state exists in the 857 semi-plane (Fig. 5A, violet and red regions): 858 > 0 (23) 859 Provided that our choice of synaptic couplings and time constants hold the stability conditions 860 (Eqs. 21-22), the U state is stable in the semi-plane given by Eq. 20 (Fig. 5A, orange and red regions): 861 In the intersection of these two semi-planes both D and U are stable (bistable region, Fig. 5A red). In 863 contrast, in the complementary region to the two semi-planes, neither U nor D are stable (Fig. 5A  864 white region). There, a rhythmic concatenation of relatively long U and D periods is observed where 865 the network stays transiently in each state until adaptation triggers a transition (see e.g. Fig. 4E). 866 Because of the separation of time-scales, we refer to this stability to the rate dynamics but not to the 867 adaptation dynamics as quasi-stable states. 868 The addition of noise makes that some of the stable solutions now become meta-stable, meaning that 869 the network can switch to a different state by the action of the noise (i.e. the external fluctuations in our 870 model). This is the case of the bistable region (Fig. 5A red) where fluctuations generate stochastic 871 transitions between the two metastable U and D states (Fig. 4D). In the region of D stability > 0, we 872 find a new subregion with noise-driven transitions from a metastable D state to a quasi-stable U state, 873 and back to D by the action of adaptation (Fig. 5A light violet). This subregion is delimited by the 874 condition that U is not stable (i.e. Eq. 24 does not hold) but just because of the existence of 875 adaptation. This can be written by saying that Eq. 24 holds if β=0: 876 < ′ (25) 877 878 Equivalently, within the region of U stability, noise creates a new subregion with noise-driven 879 transitions from a metastable U state to a quasi-stable D state, and back to U by the recovery from 880 adaptation (Fig. 5A light orange). This subregion is given by the condition that there is a negative 881 effective threshold < 0 (i.e. caused by a supra-threshold mean external drive) but the adaptation 882 recruited in the U state is sufficient to counterbalance it: + > 0. This makes the D transiently 883 stable until adaptation decays back to zero. Substituting = (Eq. 14) and by the equilibrium 884 rate at the U state given by Eq. 18, the limit of this subregion can be expressed as (Fig. 5A,