Coupled Flip-Flop Model for REM Sleep Regulation in the Rat

Recent experimental studies investigating the neuronal regulation of rapid eye movement (REM) sleep have identified mutually inhibitory synaptic projections among REM sleep-promoting (REM-on) and REM sleep-inhibiting (REM-off) neuronal populations that act to maintain the REM sleep state and control its onset and offset. The control mechanism of mutually inhibitory synaptic interactions mirrors the proposed flip-flop switch for sleep-wake regulation consisting of mutually inhibitory synaptic projections between wake- and sleep-promoting neuronal populations. While a number of synaptic projections have been identified between these REM-on/REM-off populations and wake/sleep-promoting populations, the specific interactions that govern behavioral state transitions have not been completely determined. Using a minimal mathematical model, we investigated behavioral state transition dynamics dictated by a system of coupled flip-flops, one to control transitions between wake and sleep states, and another to control transitions into and out of REM sleep. The model describes the neurotransmitter-mediated inhibitory interactions between a wake- and sleep-promoting population, and between a REM-on and REM-off population. We proposed interactions between the wake/sleep and REM-on/REM-off flip-flops to replicate the behavioral state statistics and probabilities of behavioral state transitions measured from experimental recordings of rat sleep under ad libitum conditions and after 24 h of REM sleep deprivation. Reliable transitions from REM sleep to wake, as dictated by the data, indicated the necessity of an excitatory projection from the REM-on population to the wake-promoting population. To replicate the increase in REM-wake-REM transitions observed after 24 h REM sleep deprivation required that this excitatory projection promote transient activation of the wake-promoting population. Obtaining the reliable wake-nonREM sleep transitions observed in the data required that activity of the wake-promoting population modulated the interaction between the REM-on and REM-off populations. This analysis suggests neuronal processes to be targeted in further experimental studies of the regulatory mechanisms of REM sleep.

In the coupled flip-flop model, firing rates (in Hz) of the wake-promoting, sleep-promoting, REM-on and REM-off neuronal populations, f W, f S, f R on , and f R of f , respectively, are governed by the following equations: where indicates differentiation with respect to time (in s). Concentration levels of neurotransmitters expressed by each population, cW, cS, cR on , cR of f (dimensionless with values between 0 and 1), are described by the following equations: The parameters g X,Y for X, Y = W, S, R on , R of f indicate the strength of the effect of the neurotransmitter expressed by population X on the activity of population Y . The time constants τ X for X = W, S, R on , R of f describe the time scale of evolution of firing rate changes in population X, and the time constants τ cX describe the time scale of evolution of neurotransmitter expression by population X. Population firing rate dynamics are determined by the sigmoidal steady state firing rate functions X ∞ (c) given by: where α X dictates the slope of the sigmoidal curve and β ∞,X gives the half-activation threshold. For X = S, R of f , the half-activation threshold (in square brackets) varies with the homeostatic sleep drive h and the REM sleep homeostatic drive stp, respectively, through a linear transformation with constants k 1 X and k 2 X . Neurotransmitter expression dynamics are determined by the saturating steady state neurotransmitter release functions cX ∞ (f ) for X = W, S, R on , R of f given by: where γ X for X = W, S, R on , R of f dictates the slope of the saturating curve. The homeostatic sleep drive, h, increases towards a maximum value h max during wake (f W above a threshold value θ W ) and decreases during sleep states (f W ≤ θ W ) to a minimum value h min with time constants τ h,up and τ h,down , respectively: The REM sleep homeostatic drive, stp, increases towards a maximum stp max during NREM sleep states (f W ≤ θ W and f R on ≤ θ R on ) and decreases towards a minimum stp min during REM sleep states (f W ≤ θ W and f R on > θ R on ) with time constants τ stp,up and τ stp,down , respectively. During wake (f W > θ W ), stp approaches the reset value stp r with time constant τ stp,W : There are five sources of noise in the model. The first four modulate neurotransmitter release by the W, S, R on , and, R of f populations; this is accomplished by randomly updating the multiplicative scaling parameters ξ W , ξ S , ξ R on , and ξ R of f respectively. Each parameter is updated at independent and random times, where each time is drawn from an independent Poisson distribution. For ξ W and ξ S we used Poisson distributions with λ = 0.02, while we used λ = 0.01 for ξ R on , ξ R of f . When the time arose to update ξ W , ξ S , ξ R on or ξ R of f , a new value was chosen from a uniform distribution between 0.5 and 1.5.
The final source of noise is the external random excitatory stimuli to the wake population, δ(t). δ(t) is updated at random times drawn from a Poisson distribution where the λ parameter is variable. When δ is updated, it is reset to δ = 10, and then it decays to 0 with timescale τ stim according to: The time of the next δ update is calculated using λ = ω where ω approaches a maximum value ω max during REM sleep (f R on > θ R on ) and approaches a minimum value ω min during all other states, with a time constant τ ω : Parameter values for the baseline sleep condition are listed in Table S1. The following parameters were changed for the post-REM sleep deprivation condition: h max = 1.1, stp max = 1.3, stp min = −0.6, k 2 R of f = 5.25, ω max = 0.005 and ω min = 0.0007.

Effect of variability on the dynamics of a single flip-flop model: Analysis
Here, we describe in more detail our analysis of how physiologically motivated sources of noise, which we include in our models, perturb hysteresis loop dynamics of a single flipflop. First, we consider the effects of variability in neurotransmitter expression governed by the randomly varying terms ξ X (t) X = W, S) which multiplicatively scale the steady state neurotransmitter expression functions cX ∞ (X = W, S). As ξ W , ξ S take on different values, in the interval around 1, they affect the Z-shaped fixed point curve, and thus the hysteresis loop, by changing the values of RK W and LK W . As shown in Fig. S1, the distance between RK W and LK W decreases for ξ W and ξ S values less than 1 and increases for values greater than 1. This effect on the width of the hysteresis loop is a direct result of decreases and increases in inhibition between populations. Specifically, variations in ξ W have a large effect on the value of RK W by changing the level of inhibition from the wake-promoting population to the sleep-promoting population. When this inhibition is high (ξ W > 1), the transition from the wake state (f W activated/f S deactivated) to the sleep state requires higher values of h to overcome this increased inhibition and induce activation of f S. Thus, RK W increases. Conversely, when ξ W is lower than 1, the sleep-promoting population receives less inhibition and can activate at lower values of h. Thus, RK W decreases. Variations in ξ S have a greater effect on LK W by influencing the level of inhibition from the sleep-promoting population to the wake-promoting population. Higher levels of inhibition to the wake population during the simulated sleep state (f W deactivated/f S activated) require more complete deactivation of the sleep-promoting population to allow the wake population's release. Thus, h must decrease to lower values for the sleep-to-wake transition to occur and LK W decreases. When either ξ W or ξ S are very low, the knees of the Z-shaped curve coalesce and hysteresis loop dynamics disappear. In this case, inhibition between populations is too weak to support a state where one population is fully activated and the other population is fully deactivated. Durations of wake and sleep bouts randomly vary due to the changing shape of the hysteresis loop. Given that variations in ξ W and ξ S induce both lengthening and shortening of the hysteresis loop, one might expect that the distribution of bout durations takes on a generally Gaussian shape centered around the duration of the deterministic model where the hysteresis loop remains fixed. This can occur when the time scale of bout durations in the deterministic model are similar to the time scale of the variability (i.e. the average frequency of ξ W and ξ S random variations). However, features of the bout duration distributions change when the time scales of the variability and the deterministic dynamics are different, as in our model where we assume that the neurotransmitter variability occurs at a higher frequency than wake and sleep transitions. In this situation, model simulations show that bout distributions are more positively skewed than Gaussian, with a tail of longer durations, and that the distribution peaks are at shorter durations than that of the deterministic model (Fig. S2 A,B). The majority of bout durations are shorter than the deterministic durations due to the higher frequency of variability compared to the bout duration. Namely, if we consider the movement of the trajectory around the hysteresis loop during either a wake bout or a sleep bout, the trajectory must evolve to the right or left knee, respectively, before the bout is terminated. To obtain a long bout, ξ W and ξ S must remain at appropriate values that define a long distance between RK W and LK W , or set RK W (LK W ) near h max (h min ) until the trajectory reaches the knee. When variability frequency is high, ξ W or ξ S will most likely change values before that time, thus modulating the hysteresis loop and forcing an earlier bout termination. However, occasionally ξ w and ξ S variability do allow a long bout to occur as evidenced by the distribution tail of long bout durations.
As an example for the introduction of asymmetry in bout duration distributions, we note that for the parameter values we consider here, the wake bout duration distribution has a longer positive tail than the sleep bout distribution, indicating that a small number of very long wake bouts do occur. Such very long wake bouts are possible because, high values of ξ W move RK W to high h values that approach h max (=1.4 in these simulations). Thus, it is possible that ξ W and ξ S could take on values where RK W > h max which would result in the maintenance of the wake state with h saturated at h max until the next random change of ξ W or ξ S , thus promoting a very long wake bout. Such a mechanism for long sleep bouts is not possible for the parameter values considered here since LK W does not approach h min = −1.6 for any values of ξ W or ξ S . In the Discussion, we provide possible physiological mechanisms that could support these model dynamics.
The second source of physiologically motivated variability we include in our model is brief, random excitatory stimuli input to the wake-promoting population, δ(t), representing external inputs from brain areas not included in the model such as sensory or cortical areas. These stimuli have qualitatively different effects on the distributions of wake and sleep bout durations (Fig. S2 C,D). These perturbations generate a bi-modal distribution of wake bout durations, with a sharp peak at very short durations indicating the brief wake bouts induced by the stimuli, and a broader peak at durations shorter than in the deterministic model. Here, wake bouts are generally shorter than the deterministic duration since stimuli-induced transitions to wake result in the trajectory being closer to the right knee at the beginning of the wake bout. The distribution of sleep bout durations takes on a more exponential-shape because sleep bouts are interrupted at random times by either brief wake bouts or early transitions to wake. The maximum possible sleep bout duration is the deterministic duration since the stimuli to the wake population cannot extend sleep bout durations.

Differences between baseline and post-REM sleep deprivation sleep patterning: Analysis
To lengthen the REM-NREM hysteresis loop by increasing the distance between the knees of S, we adjusted the influence of stp on f R of f activity by decreasing its effect on the halfactivation threshold of the REM-off population with a decrease in the parameter k 2 R of f (see Eq. S12). In Fig. S3, the vertical distance between the curves indicates the stp interval between the left and right knees of S, LK S and RK S , respectively, as a function of k 2 R of f . Smaller values of k 2 R of f increase the distance between the knees with a larger effect on LK S which governs the REM to NREM sleep transition.