{'fraction_non_alphanumeric': 0.04862630098790952, 'fraction_numerical': 0.017114150421844578, 'mean_word_length': 4.249716088328076, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Entrainment in up and down states of neural populations: non-smooth and stochastic models

We study the impact of noise on a neural population rate model of up and down states. Up and down states are typically observed in neuronal networks as a slow oscillation, where the population switches between high and low firing rates (Sanchez-Vivez and McCormick, 2000). A neural population model with spike rate adaptation is used to model such slow oscillations, and the timescale of adaptation determines the oscillation period. Furthermore, the period depends non-monotonically on the background tonic input driving the population, having long periods for very weak and very strong stimuli. Using both linearization and fast-slow timescale separation methods, we can compute the phase sensitivity function of the slow oscillation. We find that the phase response is most strongly impacted by perturbations to the adaptation variable. Phase sensitivity functions can then be utilized to quantify the impact of noise on oscillating populations. Noise alters the period of oscillations by speeding up the rate of transition between the up and down states. When common noise is presented to two distinct populations, their transitions will eventually become entrained to one another through stochastic synchrony.

a medium with typical extracellular ion concentrations (Sanchez-Vives and Mc-Cormick, 2000). A key feature of these slow oscillations is that they tend to be an alternating sequence of two bistable states, referred to as the up and down states. Up states in networks are characterized by high levels of firing activity, due to depolarization in single cells. Down states in networks typically appear quiescent, due to hyperpolarization in single cells. There is strong evidence that up states are neural circuit attractors, that emerge due to synaptic feedback (Cossart et al., 2003). This suggests up states may be spontaneous remnants of stimulus-induced persistent states utilized for working memory (Wang, 2001) and other network computations (Major and Tank, 2004).
Several different cellular and synaptic mechanisms have been suggested to underlie the transitions between up and down states. One possibility is that the network is recurrently coupled with excitation, stabilizing both a quiescent and active state (Amit and Brunel, 1997;Renart et al., 2007). Fluctuations due to probabilistic synapses, channel noise, and randomness in network connectivity can then lead to spontaneous transitions between the quiescent and active state (Bressloff, 2010;Litwin-Kumar and Doiron, 2012). Alternatively, switches between low and high activity states may arise by some underlying systematic slow process. For instance, it has been shown that competition between recurrent excitation and the negative feedback produced by activity-dependent synaptic depression can lead to slow oscillations in firing rate whose timescale is set by the depression timescale (Bart et al., 2005;Kilpatrick and Bressloff, 2010). Excitatory-inhibitory networks with facilitation can produce slow oscillations, due to the slow facilitation of feedback inhibition that terminates the up state, the down state is then rekindled due to positive feedback from recurrent excitation (Melamed et al., 2008). These neural mechanisms utilize dynamic changes in the strength of neural architecture. However, Compte et al. (2003) proposed that single cell mechanisms can also shift network states between up and down states. The up state is maintained by strong recurrent excitation balanced by inhibition, and transitions to the down state occur due to a slow adaptation current. Once in the down state, the adaptation current is inactivated, and excitation reinitiates the up state. A similar mechanism has been utilized in models of perceptual rivalry, where dominance switches between two mutually inhibiting populations are due to the build up of a rate-based adaptation current (Laing and Chow, 2002;Moreno-Bote et al., 2007).
In this paper, we utilize a rate-based model of an excitatory network with spike rate adaptation to explore the impact that noise perturbations have upon the relative phase and duration of slow oscillations. We find that, as in the spiking model studied by Compte et al. (2003), the interplay between recurrent excitation and adaptation produces a slow oscillation in the firing rate of the network. In fact, for slow timescale adaptation currents, the oscillations evolve as fast switches between a low and high activity state, stable fixed points of the adaptation-free system. Since the timescale and slow dynamics of the oscillation are set by the adaptation current, we mainly focus on the impact of perturbation to the adaptation variable in our model. As we will show, perturbations of the activity variable have much lower impact on the oscillation phase. Introducing noise into the adaptation variable of the population model leads to a speeding up of the slow oscillation, due to early switching between the low and high activity state.
Another remarkable feature of slow oscillations, observed during slow-wave sleep and anesthesia, is that the up and down states tend to be synchronized across different regions of cortex and thalamus (Steriade et al., 1993;Massimini et al., 2004). Specifically, both the up and down states start near synchronously in cells located up to 12mm apart (Volgushev et al., 2006). Such remarkable coherence between distant network activity cannot be accomplished by single cell mechanisms, but require either long range network connectivity or some external signal forcing entrainment (Traub et al., 1996;Smeal et al., 2010). Activity tends to originate from several different foci in the network, quickly spreading across the rest of the network on a timescale orders of magnitude faster than the oscillation itself (Compte et al., 2003;Massimini et al., 2004). The fact that the onset of quiescence is fast and well synchronized means there must be either a rapid relay signal between all foci or there is some global signal cueing the down state. Rather than suggest a disynaptic relay, using long range excitation acting on local inhibition, we suggest that background noise can serve as a synchronizing global signal (Ermentrout et al., 2008). Noisy but correlated inputs have been shown to be capable of synchronizing uncoupled populations of phase oscillators (Teramae and Tanaka, 2004) as well as experimentally recorded cells in vitro (Galán et al., 2006). Here we will show correlated noise is a viable mechanisms for coordinating slow oscillations in distinct uncoupled neural populations.
The paper is organized as follows. We introduce the neural population model in section 2, indicating the way external noise is incorporated into the model. In section 3, we demonstrate the periodic solutions that emerge in the noise-free model, demonstrating it is possible to derive analytical expressions for the oscillation period in the case of steep firing rate functions. Then, in section 4 we show how to derive phase sensitivity functions that describe how external perturbations to the periodic solution impact the asymptotic phase of the oscillation. As demonstrated, the impact of perturbations to the adaptation variable is much stronger than activity variable perturbations, especially for longer adaptation timescales. Thus, our studies of the impact of noise mainly focus on the effects of fluctuations in the adaptation variable. We find, in section 5, that adding noise to the adaptation variable leads to up and down state durations that are shorter and more balanced, so that the up and down state last for similar lengths of time. In section 6, we demonstrate that slow oscillations in distinct populations can become entrained to one another when both populations are forced by the same common noise signal. This phenomenon is robust to the introduction of independent noise in each population, as we show in section 7.

Adaptive neural populations: deterministic and stochastic models
We begin by describing the models we will use to explore the impact of external perturbations on slow oscillations. Motivated by Compte et al. (2003), we will focus on a neural population model with spike rate adaptation, akin to mutual inhibitory models used to study perceptual rivalry (Laing and Chow, 2002;Moreno-Bote et al., 2007).
Single population model. In a single population, neural activity u(t) receives negative feedback due to a subtractive spike rate adaptation term (Benda and Herz, 2003) Here, u represents the mean firing rate of the neural population with excitatory connection strength α. The negative feedback variable a is spike frequency adaptation with strength φ and time constant τ . For some of our analysis we will utilize the assumption τ 1, based on the fact that many form for spike rate adaptation tend to be much slower than neural membrane time constants (Benda and Herz, 2003). The constant tonic drive I initiates the high firing rate (up) state, and slow adaptation eventually attenuates activity to a low firing rate (down) state. Weak but positive drive I > 0 is meant to model the presence of low spiking threshold cells that spontaneously fire, utilized as a mechanism for initiating the up state in Compte et al. (2003). The firing rate function f is monotone and saturating function such as the sigmoid Commonly, in studies of neural field models, the high gain limit of (2) is taken to yield the Heaviside firing rate function (Amari, 1977;Laing and Chow, 2002) which often allows for a more straightforward analytical study of model dynamics.
We exploit this fact extensively in our study. Nonetheless, we have also carried out many numerical simulations of the model for a smooth firing rate function (2), and they correspond to the results we presentfor sufficiently high gain. Note, this form of adaptation is often referred to as subtractive negative feedback, as current is subtracted from the population input. Alternative models of slow neural population oscillations have employed short term synaptic depression (Tabak et al., 2000;Bart et al., 2005;Kilpatrick and Bressloff, 2010), a form of divisive negative feedback.
A primary concern of this work is the response of (1) to external perturbations, acting on the activity u and adaptation a variables. To do so, we will use both an exact method and a linearization to identify the phase response curve of the limit cycle solutions to (1). Understanding the susceptibility of limit cycles (1) to inputs will help us understand ways in which noise will influence the frequency and regularity of oscillations. Stochastic single population model. Following our analysis of the noisefree system, we will consider how fluctuations influence oscillatory solutions to (1). To do so, we will employ the following Langevin equation for (1) forced by white noise where we have introduced the independent Gaussian white noise processes ξ u (t) and ξ a (t) with zero mean ξ u (t) = ξ a (t) = 0 and variances ξ u (t) 2 = σ 2 u t and ξ a (t) 2 = σ 2 a t. Extending our results concerning the phase response curve, we will explore how noise forcing impacts the statistics of the resulting stochastic oscillations in (4). In particular, since we find noise tends to impact the phase of the oscillation more strongly when applied to the adaptation variable, we will tend to focus on the case ξ u ≡ 0.
Stochastic dual population model. Finally, we will focus on how correlations in noise-forcing impact the coherence of two distinct uncoupled populations Thus, the system (5) describes the dynamics of two distinct neural populations u 1 and u 2 , with inputs I. Our main interest lies in the impact the noise terms have upon the phase relationship between the two systems' states. In this version of the model, noise to the activity variables ξ u is totally correlated, as is noise to the adaptation variables ξ a . Thus, all means are zero and ξ 2 u (t) = σ 2 u = D 11 t. Furthermore, ξ 2 a (t) = σ 2 a t = D 22 t. For this study, we assume there are no correlations between activity and adaptation noise, so ξ u (t)ξ a (t) = 0. A more general version of the model (5) would consider the possibility of independent noise in each population Noise terms all have zero mean and variances defined ξ 2 uj (t) = σ 2 uj t = D uj t and ξ 2 aj (t) = σ 2 aj t = D aj t (j = 1, 2, c). To ease calculations, we take D u1 = D u2 ≡ D ul = σ 2 u and D a1 = D a2 ≡ D al = σ 2 a . The degree of noise correlation between populations is controlled by the parameters χ u and χ a , so in the limit χ u,a → 1, the model (6) becomes (5).

Periodic solutions of a single population
We begin by studying periodic solutions of the single population system (1), as demonstrated in Fig. 1A. First, we note that for firing rate functions f with finite gain, we can identify the emergence of oscillations by analyzing the stability of the equilibria of (1). That is, we assume (u,ȧ) = (0, 0), so the system becomes which can be reduced to the single equation Roots of (7), defining fixed points of (1) are plotted as a function of the input I in Fig. 1C,D. Utilizing the sigmoidal firing rate function f given by (2), we can show that there will be a single fixed point as long as φ > α. In this case, we can compute crosses zero once, so (7) has a single root when φ > α. Stability of this equilibrium is given by the eigenvalues of the associated Jacobian We note that the sigmoid (2) satisfies the Ricatti equation f = γf (1 − f ), so we can use (7) to write Oscillations arise when stable spiral equilibria destabilize through a Hopf bifurcation. Hopf bifurcations will occur when complex eigenvalues associated with fixed points (ū,ā) cross from the left to the right half plane. We can require this with the pair of expression: tr(J) = 0 and tr(J) 2 < 4det(J). Thus, a necessary condition for the Hopf bifurcation point is that the equilibrium valueū satisfy Solving this forū yields Thus, Hopf bifurcations will only occur when the timescale of adaptation is sufficiently large τ > [αγ/4 − 1] −1 . Plugging the formula (8) back into the fixed point equation (7) and solving for the input I, we can parameterize Hopf bifurcation curves based upon the equation along with the additional condition tr(J) 2 < 4det(J) which becomes which will always hold as long as φ > α. We partition the parameter space (τ, I) using our formula for the Hopf curve (9) in Fig. 1B. As demonstrated, there tend to be either two or zero Hopf bifurcation points for a given timescale τ , and the coalescence of the two Hopf points is given by the point where τ = [αγ/4 − 1] −1 .
In the limit of slow adaptation τ 1, we can separate the timescales of the activity u and adaptation a variables, finding u will equilibrate according to the equationû and subsequently a will slowly evolve according to the equatioṅ We always have an implicit formula forû(t) in terms of a(t), so the dynamics will tend to slowly evolve along the direction of the a variable. This demonstrates why periodic solutions to (1) are comprised of a slow rise and decay phase of a, punctuated by fast excursions in the activity variable u. In general, it is not straightforward to analytically treat the pair of equations (11) and (12), but we will show how computing solutions of the singular system becomes straightforward when we take the high gain limit γ → ∞. The period of the oscillation T computed from numerical simulations (dots) is accurately approximated by the analytical formula (solid lines) given by (14). Other parameters are α = 0.5 and φ = 1.
Having established the existence of periodic solutions to (1) in the case of sigmoid firing rates (2), we now explore the system in the high gain limit γ → ∞ whereby the firing rate function becomes a Heaviside (3). In this case, fixed points (ū,ā) satisfy the equations Thus, assuming φ > α, thenū = 0 when I < 0 andū = 1 when I > (φ − α). In both cases, the fixed points are linearly stable. When 0 < I < (φ − α), there are no fixed points and we expect to find oscillatory solutions. Assuming τ 1, we can exploit a separation of timescales to identify the shape and period of these limit cycles. To begin, we note that on fast timescaleṡ where a 0 is a quasi steady state. On slow timescales on the order of τ , then u(t) quickly equilibrates and Periodic solutions to (1) must obey the condition (u(t), a(t)) = (u(t + nT ), a(t + nT )) for t ∈ [0, T ] and n ∈ Z, so we focus on the domain t ∈ [0, T ]. Examining (13), we can see oscillations in (1) involve switches between u(t) ≈ 1 and u(t) ≈ 0. We translate time so that u(t) ≈ 1 on t ∈ [0, T 1 ) and u(t) ≈ 0 on t ∈ [T 1 , T ). Subsequently, this means for t ∈ [0, T 1 ] the system (13) becomes u ≡ 1 and (13), and the argument of H(x) must have crossed zero at t = 0. In a similar way, we find on t ∈ [T 1 , T ) that u ≡ 0 and a(t) = (I + α)e −(t−T 1 )/τ . Using the conditions a(T 1 ) = I +α and a(T ) = I, we find that the rise time of the adaptation variable (or the duration of the up state) is and the decay time (or the duration of the down state) is and the total period of the oscillation is Thus, approximate periodic solutions to (1) in the case of a Heaviside firing rate (3) take the form We demonstrate the accuracy of the approximation (15) in Fig. 2A. Furthermore, we show that relationship between the period T and model parameters is well captured by the formula (14). Notice there is a non-monotonic relationship between the period T and the input I. We can understand this further by noting that the rise time T 1 of the adaptation variable a increases monotonically with input . Furthermore, the decay time T 2 of the adaptation variable a decreases monotonically with input . Thus, as I → 0 + , the slow oscillation's period T is dominated by very long decay times T 2 1 and as I → (φ − α) − , it is dominated by very long rise times T 1 1. We can identify the minimal period as a function of the input I by finding the critical point of T (I). To do so, we differentiate and simplify so the critical point of T (I) is I crit = (φ − α)/2, which corresponds to the minimal value of the period T min (I) = 2τ ln [(φ + α)/(φ − α)] as pictured in Fig. 2B.

Phase response curves
We can further understand the dynamics of the slow oscillations in (1) by computing phase response curves for both the case of a sigmoidal firing rate (2) and the Heaviside firing rate (3). As we will show, perturbations of the activity variable u have decreasing impact as the timescale of adaptation τ and the gain γ of the firing rate are increased. Perturbations of the adaptation variable a tend to dominate the resulting dynamics, as it is the evolution of this slow variable that primarily determines the phase of the oscillation. To begin, we derive a general formula that linearly approximates the influence of small perturbations on limit cycle solutions (u 0 (t), a 0 (t)) to (1). Essentially, we utilize the fact that solutions Z(t) to the adjoint equation associated with linearization about the limit cycle solution (u 0 (t), a 0 (t)) provide a complete description of how infinitesimal perturbations of the limit cycle impact its phase (Ermentrout, 1996;Brown et al., 2004). To start, we note that is the linearization of (1) about the limit cycle (u 0 (t), a 0 (t)). Defining the inner product on T -periodic functions in R 2 as F (t), G(t) = T 0 F (t) · G(t)dt, we can find the adjoint operator L * by noting it satisfies F, LG = L * F, G for all L 2 integrable vector functions F, G. We can then compute It can be shown that the null space of L * describes the response of the phase of the limit cycle (u 0 (t), a 0 (t)) to infinitesimal perturbations (Brown et al., 2004). Note that if (u 0 (t), a 0 (t)) is a stable limit cycle then the nullspace of L is spanned by scalar multiples of (u 0 (t), a 0 (t)). Furthermore, appropriate normalization requires that Z(t) · (u 0 (t), a 0 (t)) = 1 along with L * Z = 0 (Ermentrout, 1996). To numerically compute Z(t) = (Z u (t), Z a (t)), we thus integrate the systeṁ backward in time, taking the long time limit to find (Z u (t), Z a (t)) on t ∈ [0, T ], and normalizing (Z u (t), Z a (t)), (u 0 (t), a 0 (t)) = 1 by rescaling appropriately. We demonstrate this result in Fig. 3, showing the relationship between the shape and relative amplitude of the phase sensitivity functions (Z u , Z a ) and the parameters.
Notably, perturbing the activity variable u become less influential as the timescale of adaptation τ is increased (Z u ). Furthermore, there is a sharper transition between phase advance and phase delay region of the adaptation phase response (Z a ) for larger timescales τ . In addition to a general formula for the phase sensitivity functions (Z u (t), Z a (t)), we can derive an amplitude dependent formula for the response of limit cycle solutions (u 0 (t), a 0 (t)) of (1) with a Heaviside firing rate (3), assuming τ 1. In this case, we utilize the formula for the period (14) and limit cycle (15), derived using a separation of timescales assumption. Then, we can compute the change to the variables (u, a) as a result of a perturbation (δ u , δ a ), which we denote (u 0 (t), a 0 (t)) (δ u ,δ a ) −→ (ũ 0 (t),ã 0 (t)). We are primarily interested in how the relative time in the limit cycle is altered by a perturbation δ u -how much closer or further the limit cycle is to the end of the period T after being perturbed. We can readily determine this by first inverting the formula we have for (u 0 (t), a 0 (t)), given by (15), to see how this value determines the time t 0 along the limit cycle Using this formula, we can now map the value (ũ 0 ,ã 0 ) to an associated updated relative time t 0 along the oscillation.
Here, we decompose the impact of perturbations to the u and a variables. We begin by studying the impact of perturbations δ u to the activity variable u. We can directly computeũ Thus, the singular system (13) will be unaffected by such perturbations if sgn(I + α[u + δ u ] − a) = sgn(I + αu − a). This is related to the flatness of the susceptibility function Z u over much of the time domain in Fig. 3D,E,F. However, if sgn(I + Fig. 4 Phase response curves of the fast-slow timescale separated system τ 1. (A,B) Amplitude δu-and δa-dependent phase response curves Gu(θ, δu) and Ga(θ, δu) characterizing phase advances/delays resulting from perturbation of neural activity u and adaptation a. We compare analytical formulae (solid lines) to numerically computed PRCs (dashed lines). (C) Phase response curve associated with perturbations of the adaptation variable a in the small amplitude 0 < |δa| 1 limit. We compare the large amplitude formula (solid line) determined by (24) to the linear approximation (dotted line) given by (25)  α[u + δ u ] − a) = sgn(I + αu − a), thenũ 0 (t) = 1 − u 0 (t), as detailed in the following piecewise smooth map: where (u 0 (t), a 0 (t)) are defined by (15). The formula (18) can then be utilized to compute the updated relative timet 0 := t 0 (ũ 0 ,ã 0 ), finding where a 0 = φ − (φ − I)e −t 0 /τ if u 0 = 1 and a 0 = (I + α)e −(t 0 −T 1 )/τ if u 0 = 0. We can refer to the functiont 0 /T , wheret 0 is defined by (19), as the phase transition curve for u perturbations. Thus, the function G u (θ, δ u ) = (t 0 − t 0 )/T will be the phase response curve, where θ = t 0 /T , and phase advances occur for positive values and phase delays occur for negative values. We plot the function G u (θ, δ u ) in Fig. 4A for different values of δ u , demonstrating the nontrivial dependence on the perturbation amplitude is not simply a rescaling but an expansion of the non-zero phase shift region. Due to the singular nature of the fast-slow limit cycle (15), the size of the phase perturbation has a piecewise constant dependence on the amplitude of the u perturbation. Note, this formulation allows us to quantify phase shifts that would not be captured by a perturbative theory for phase sensitivity functions, as computed for the general system in (17). For perturbations δ a of the adaptation variable a, there is a more graded dependence of the phase advance/delay amplitude on the perturbation amplitude δ a . We expect this, as it was a property we observed in Z a as we varied parameters in Fig. 3. We can partition the limit cycle (u 0 (t), a 0 (t)) into four different regions: two advance/delay regions of exponential saturation and two early threshold crossings. First, note if u 0 (t) = 1 and a 0 (t) + δ a < I + α, theñ Determining the relative time of the perturbed variables (ũ 0 (t),ã 0 (t)) in (20) is straightforward using the mapping (18). However, to determine the relative time described by (21), we compute the time, after the perturbation, untilã 0 (t) = I +α, which will be t w = τ ln [(a 0 + δ a )/(I + α)], sot 0 = T 1 − t w . Second, note if u 0 (t) = 0 and a 0 (t) + δ a > I, theñ In the case of (23), we note that t Combining our results, we find we can map the relative time to the perturbed relative time as where a 0 = φ − (φ − I)e −t 0 /τ if u 0 = 1 and a 0 = (I + α)e −(t 0 −T 1 )/τ if u 0 = 0. Again, we have a phase transition curve given by the functiont 0 /T and phase response curve given by G a (θ, δ a ) = (t 0 − t 0 )/T , where θ = t 0 /T . As opposed to the case of u perturbations, the phase perturbation here depends smoothly on the amplitude of the a perturbation δ a . Furthermore, we can obtain a perturbative description of the phase response curve for the singular system (13) in two ways: (a) Taylor expand the amplitudedependent phase response curve expressions defined by (19) and (24) and truncate to linear order or (b) solving the adjoint equation (17) in the case of a Heaviside firing rate (3) and long adaptation timescale τ 1. We begin with the first derivation, which simply requires differentiating (19) to demonstrate that the infinitesimal phase response curve (iPRC) associated with perturbations of the u variable is zero almost everywhere. However, differentiating (24) reveals that the iPRC associated with perturbations of the adaptation variable a is given by the piecewise smooth function Furthermore, note we could derive the same result by solving the adjoint equations (17) in the case of Heaviside firing rate (3), so thaṫ Note, we have reversed time t → −t, so we can simply solve the system forward. Furthermore, we can use the identity Utilizing the separation of timescales, τ 1, we find that almost everywhere (except where t = 0, T 1 , T ), we have that (26) becomes the systeṁ τŻ a (t) = Z a (t).
As before Z u (t) will be zero almost everywhere, whereas Z a (t) = A(t)e t/τ , where A(t) is a piecewise constant function taking two different values on t ∈ (0, T 1 ) and t ∈ (T 1 , T ), determined by considering the δ distribution terms. This indicates how one would derive the formula (25) using the adjoint equations (26). Note, in previous work (Jayasuriya and Kilpatrick, 2012), we explored the entrainment of slowly adapting populations to external forcing, comprised of smooth and non-smooth inputs to the system (1). In the next section, we explore the impact of external noise forcing on the slow oscillations of (1), subsequently demonstrating that noise can be utilized to entrain the up and down states of two distinct networks.

Impact of noise on the timing of up/down states
We now study the effects of noise on the duration of up and down states of the single population model (1). Switches between high and low firing rate states can occur at irregular intervals (Sanchez-Vives and McCormick, 2000), suggesting internal or external sources of noise determine state changes. This section focuses on how noise can reshape the mean duration of up and down residence times. Due to our findings in the previous sections, we focus on noise applied to the adaptation variable in this section. As we have shown, very weak perturbations to the neural activity variable have a negligible effect on the phase of oscillations. Analytic calculations are presented for the piecewise smooth system with Heaviside firing rate (3), as accurate approximations of the mean up and down state durations can be computed.
Our approach is to derive expressions for the mean first passage times of both the up and down state (T 1 andT 2 ) of the stochastic population model (4). Focusing on adaptation noise allows us to utilize the separation of fast-slow timescales, and recast the pair of equations as a stochastic-hybrid system where ξ a is white noise with mean ξ a = 0 and variance ξ 2 a = σ 2 a t. To begin, assume the system has just switched to the up state, so the initial conditions are u(0) = 1 and a(0) = I. Determining the amount of time until a switch to the down state requires we calculate the time T 1 until the threshold crossing a(T 1 ) = I + α where a(t) is determined by the stochastic differential equation (SDE) which is the well-known threshold crossing problem for an Ornstein-Uhlenbeck process Gardiner (2004). The meanT 1 of the passage time distribution is thus given by defining the potential V (a) = a 2 2τ − φa τ and computing the integral Next, note that the duration of the down state T 2 will be the amount of time until the threshold crossing a(T 2 ) = I given u(0) = 0 and a(0) = I + α, where a(t) obeys the SDE da(t) = [−a(t)] dt/τ + dξ a (t).
Again, defining the potential V (a) = a 2 2τ , we can computē We compare the theory we have developed utilizing passage time problems to residence times computed numerically in Fig. 5C. Notice that increasing the noise amplitude tends to shorten both up and down state durations on average, due to early threshold crossings of the variable a(t). Furthermore, we can examine how noise reshapes the relative balance of up versus down state durations. Specifically, we will explore how the relative fraction of time the up state persistsT 1 /(T 1 +T 2 ) changes with noise intensity σ a and input I. First, notice that, in the absence of noise the ratio .
The up and down state have equal duration when T 1 /(T 1 + T 2 ) = 1/2, or when the input I = (φ − α)/2, as shown in Fig. 5D. Interestingly, this is the precise input value at which the period obtains a minimum, as we demonstrated in section 3. Along with our plot of (29) in the noise-free case (σ a = 0), we also study the impact of noise on this measure of up-down state balance. Noise leads to up and down state durations becoming more similar, so the ratio (29) of the meansT 1 and T 2 flattens as a function of the input I. This is due to the fact that long durations, wherein the variable a(t) occupies the tail of exponentially saturating functions A 0 + A 1 e −t/τ , are shortened by early threshold crossings due to the external noise forcing.

Synchronizing two uncoupled populations
Now we demonstrate that common noise can synchronize the up and down states of two distinct and uncoupled populations. We begin with the case of identical noise and then, in section 7, relax these assumptions to show that some level of coherence is still possible when each population has an intrinsic and independent source of noise. This is motivated by the observation that the neural Langevin equation derived in the large system-size limit of a neural master equation tends to possess intrinsic noise in each population, in addition to an extrinsic common noise term (Bressloff and Lai, 2011). As we will show, intrinsic noise tends to disrupt the phase synchronization due to extrinsic noise. To begin, we recast the stochastic system (5), describing a pair of adapting noise-driven neural populations, as a pair of phase equations: where θ 1 and θ 2 are the phase of the first and second neural populations. As we demonstrate in Fig. 6A, this introduction of common noise tends to drive the oscillation phases θ 1 (t) and θ 2 (t) toward one another. Note that since the governing equations of both populations are the same, then the phase sensitivity function Z(θ) will be the same for both. Furthermore, the synchronized solution θ 1 (t) = θ 2 (t) is absorbing -once the phases synchronize, they remain so. We can analytically calculate the Lyapunov exponent λ of the synchronized state to determine its stability. In particular, we will be interested in how this stability depends on the parameters that shape the dynamics of adaptation.
We next convert the pair of Stratonovich differential equations into a equivalent pair of Ito differential equations: dθ 1 (t) = ω + Z (θ 1 (t)) T DZ(θ 1 (t)) dt + Z(θ 1 (t)) · dξ(t), introducing a drift term due to our change in definition of the noise term. Now, to determine stability of the synchronized state θ 1 (t) = θ 2 (t), we assume we are infinitesimally close to it. We define ψ(t) = θ 1 (t) − θ 2 (t) and require |ψ(t)| 1. Linearizing the system of Ito differential equations (31) with respect to ψ(t) then yields dψ(t) = ψ(t) Z (θ(t)) T DZ(θ(t)) dt + Z (θ(t)) · dξ , where θ(t) obeys either one of the equations in (31). Applying the change of variables y(t) = ln ψ(t), we can rewrite (32) as Notice, on average, the log of the phase difference y(t) tends to decrease over time (Fig. 6B), indicating the phases θ 1 and θ 2 move toward one another. Subsequently, we can integrate equation (33) to determine the mean drift of y(t) λ := lim t→∞ t 0 Z (θ(s))DZ(θ(s)) − Z T (θ(s))DZ (θ(s)) ds, so the phase difference ψ(t) will tend to decay grow if the Lyapunov exponent λ < 0 (λ > 0), and the synchronous state will be stable (unstable). Now, utilizing ergodicity of (33), we can compute λ using the ensemble average across realizations of Z (θ(t)) · ξ(t), so where P s (θ) is the steady state distribution of θ. Since noise is weak (D jk 1, j, k = 1, 2), we can approximate the distribution as constant P s (θ) = 1. Upon applying this to the integrand of (34) and noting the periodicity of Z(θ), we find we can approximate the Lyapunov exponent Assuming noise to the activity variable u and adaptation variable a is not correlated, D will be diagonal. In this case, we can further decompose the phase sensitivity function into its Fourier expansion where a k = (a k1 , a k2 ) T and b k = (b k1 , b k2 ) T are vectors in R 2 so that and we can expand the terms in (35) to yield Thus, as long as Z(θ) is continuous and non-constant, the Lyapunov exponent λ will be negative, so the synchronous state θ 1 = θ 2 will be stable. Note, continuity is Fig. 7 Stationary density M 0 (ψ) of the phase difference ψ = θ 1 − θ 2 for two slowly oscillating neural population driven by both common and independent noise (6). As the degree of noise correlation is decreased from (A) χa = 0.95 to (B) χa = 0.90, the density spreads, but there is still a peak at ψ = 0, the phase-locked state. We focus on noise in the adaptation variable, so σu = 0 and σa = 0.01. Other parameters are α = 0.5, γ = 15, φ = 1, and τ = 20.
not satisfied in the case of our singular approximation to Z(θ). We demonstrate the accuracy of our theory (35) in Fig. 6C,D, showing that λ decreases as a function of τ and is non-monotonic in I. Thus, slow oscillations with longer periods are synchronized more quickly, relative to the number of oscillation cycles. Since the Lyapunov exponent has highest amplitude |λ| for both low and high values of the tonic input I, we also suspect this is related to the period of the oscillation T .

Impact of intrinsic noise on stochastic synchronization
We now extend our results from the previous section by studying the impact of independent noise in each population. Independent noise is incorporated into the modified model (6). Noting, again there is a periodic solution to the noise-free version of this system, phase-reduction methods can be used to obtain approximate Langevin equations for the phase variables (Nakao et al., 2007) where the noise vectors ξ c = (χ u ξ uc , χ a ξ ac ) T and ξ j = ( 1 − χ 2 u ξ uj , 1 − χ 2 a ξ aj ) T (j = 1, 2). We can reformulate the pair of Stratonovich differential equations as Ito stochastic differential equations given by the system where the statistics of the noise terms dζ j (θ, t) = Z(θ j (t)) · dξ c (t) + dξ j (t) (j = 1, 2) are specified by dζ j (θ, t) = 0 and dζ j (θ, t)dζ k (θ, t) = C jk (θ)dt where separating the impact of correlated and local sources of noise. The drift terms can thus be calculated A j (θ) = ω + 1 4 ∂ ∂θ j C jj (θ). The Fokker-Planck equation describing the evolution of the probability density function P (θ, t) of the phases is given Now, we apply a change of variables to the Fokker-Planck equation (39) defined θ j = ωt + ϑ j . Assuming noise is weak, the function Q(ϑ, t) varies slowly compared with the period of the phase oscillators θ j . Thus, we average the drifts A j (θ) and correlation function C jk (θ) over a single period. The resulting Fokker-Planck equation is then where the averaged correlation function is given by the formula where the correlation functions are defined We study the relationship between the phases ϑ 1 and ϑ 2 by substituting the formula for the averaged correlation matrix We can write (40) as a separable equation by employing a change of variables that tracks the average ρ = (ϑ 1 + ϑ 2 )/2 and phase difference ψ = ϑ 1 − ϑ 2 of the original position variables, so Thus, we can solve for the stationary solution of the system (41a) by serving U t = M t ≡ 0 and requiring periodic boundary conditions. We find that the stationary distribution of the position average is U 0 (ρ) = 1. In addition, we can integrate the stationary equation for M (ψ, t) to find M 0 (ψ) = m 0 σ 2 u [(2 − χ 2 u )g u (0) − χ 2 u g u (ψ)] + σ 2 a [(2 − χ 2 a )g a (0) − χ 2 a g a (ψ)] , where m 0 = 1/ 1 0 [g(0) − g(x) + h(0)] −1 dx is a normalization factor and we have simplified the expression using D u1 = D u2 ≡ D ul = σ 2 u and D a1 = D a2 ≡ D al = σ 2 a and defined g j (ψ) = 1 0 Z j (θ)Z j (θ + φ)dθ.
When noise to each layer is independent (χ u , χ a → 0), then M 0 (ψ) = 1 is constant in space. When noise is totally correlated (χ u , χ a → 1), then M 0 (ψ) = δ(φ). The stationary distribution M 0 (ψ) will broaden as the correlations χ u and χ a are decreased from unity, with a peak remaining at φ = 0. We demonstrate the accuracy of the formula (42) for the stationary density of the phase difference ψ in Fig. 7, showing that it widens as the level noise correlation is decreased. Again, we focus on the impact of adaptation noise. Thus, even when independent noise is introduced, there is some semblance of synchronization in the phases of two noise-driven neural populations (6).

Discussion
We have studied the impact of deterministic and stochastic perturbations to a neural population model of slow oscillations. The model was comprised of a single recurrently coupled excitatory population with negative feedback from a slow adaptive current (Laing and Chow, 2002;Jayasuriya and Kilpatrick, 2012). By examining the phase sensitivity function (Z u , Z a ), we found that perturbations of the adaptation variable lead to much larger changes in oscillation phase than perturbations of neural activity. Furthermore, this effect becomes more pronounced as the timescale τ of adaptation is increased. Introducing noise in the model decreases the oscillation period and helps to balance the mean duration of the oscillation's up and down states. When two uncoupled populations receive common noise, their oscillation phases θ 1 and θ 2 eventually become synchronized, which can be shown by deriving a formula for the Lyapunov exponent of the absorbing state θ 1 ≡ θ 2 (Teramae and Tanaka, 2004). When independent noise is introduced to each population, in addition to common noise, the long-term state of the system is described by a probability density for ψ = θ 1 − θ 2 , which peaks at ψ ≡ 0. Our study was motivated by the observation that recurrent cortical networks can spontaneously generate stochastic oscillations between up and down states. Guided by previous work in spiking models (Compte et al., 2003), we explored a rate model of a recurrent excitatory network with slow spike frequency adaptation.
One of the open questions about up and down state transitions concerns the degree to which they are generated by noise or by more deterministic mechanisms, such as slow currents or short term plasticity (Cossart et al., 2003). Here, we have provided some characteristic features that emerge as the level of noise responsible for transitions is increased. Similar questions have been probed in the context of models of perceptual rivalry . In addition, we have provided a plausible mechanism whereby the onset of up and down states could be synchronized in distinct networks (Volgushev et al., 2006).
There are several potential extensions to this work. For instance, we could examine the impact of long-range connections between networks to see how these interact with common and independent noise to shape the phase coherence of oscillations. Similar studies have been performed in spiking models by Ly and Ermentrout (2009). Interestingly, shared noise can actually stabilize the anti-phase locked state in this case, even though it is unstable in the absence of noise. Furthermore, it is known that coupling spanning long distances can be subject to axonal delays. In spite of this, networks of distantly coupled clusters of cells can still sustain zero-lag synchronized states (Vicente et al., 2008). Thus, we could also explore the impact of delayed coupling, determining how features of phase sensitivity function interact with delay to promote in-phase or anti-phase synchronized states.