State-dependent pontine ensemble dynamics and interactions with cortex across sleep states

The pontine nuclei play a crucial role in sleep-wake regulation. However, pontine ensemble dynamics underlying sleep regulation remain poorly understood. By monitoring population activity in multiple pontine and adjacent brainstem areas, here we show slow, state-predictive pontine ensemble dynamics and state-dependent interactions between the pons and the cortex in mice. On a timescale of seconds to minutes, pontine populations exhibit diverse firing across vigilance states, with some of these dynamics being attributed to cell type-specific activity. Pontine population activity can predict pupil dilation and vigilance states: pontine neurons exhibit longer predictable power compared with hippocampal neurons. On a timescale of sub-seconds, pontine waves (P-waves) are observed as synchronous firing of pontine neurons primarily during rapid eye movement (REM) sleep, but also during non-REM (NREM) sleep. Crucially, P-waves functionally interact with cortical activity in a state-dependent manner: during NREM sleep, hippocampal sharp wave-ripples (SWRs) precede P-waves. On the other hand, P-waves during REM sleep are phase-locked with ongoing hippocampal theta oscillations and are followed by burst firing in a subset of hippocampal neurons. Thus, the directionality of functional interactions between the hippocampus and pons changes depending on sleep states. This state-dependent global coordination between pontine and cortical regions implicates distinct functional roles of sleep.


Introduction 36
The sleep-wake cycle is a fundamental homeostatic process across animal species (Anafi et al.,37 2019; Aulsebrook et al., 2016;Siegel, 2005). Characterizing these physiological properties is crucial to uncover the roles of brainstem populations 71 in sleep regulation and ultimately the functions of sleep states. 72 In the present study, we adopt several in vivo electrophysiological approaches in mice to investigate 73 state-dependent ensemble dynamics in the brainstem, mainly the pons. We show that on a timescale 74 of seconds to minutes, pontine neurons show state-dependent firing with cell type-specificity. They 75 also have a longer predictive power for vigilance states compared to those in the hippocampus. On 76 a timescale of sub-seconds, we find state-dependent functional interactions between the pons and 77 the cortex, with a focus on P-waves: during NREM sleep, the timing of P-waves is phase-locked with 78 various cortical oscillations and hippocampal SWRs precede P-waves. During REM sleep, P-waves 79 co-occur with hippocampal theta and precede burst firing of hippocampal neurons. These results 80 imply that pontine populations not only play a regulatory role in the sleep-wake cycle, but also 81 contribute to global state-dependent dynamics across brain regions. 82

83
Brainstem population recording across sleep-wake cycles 84 To investigate the state-dependency of brainstem population activity, we inserted a silicon probe into 85 the mouse brainstem in a head-fixed condition, together with simultaneous monitoring of cortical 86 electroencephalograms (EEGs), electromyograms (EMGs) and pupil dilation (Fig. 1). Recorded 87 regions spanned across multiple nuclei, including the sublaterodorsal nucleus, pontine reticular 88 nucleus, medial preoptic nucleus, parabrachial nucleus, pontine central gray, laterodorsal tegmental 89 nucleus and other surrounding areas according to post-mortem histological analysis 90 (Supplementary Fig. 1). Although a majority of neurons were recorded from the pons, we refer to 91 recorded populations as "brainstem" neurons because some cells were located in the midbrain and 92 medulla, but not the hypothalamus. 93 The sleep-wake cycle was classified based on cortical EEGs and EMGs in every 4 second. Based 94 on the classified states, we observed clear state-dependency across measurements (Figs. 1C-F): 95 wakefulness was characterized by high muscle tone (Fig. 1D) and pupil dilation (Fig. 1F) whereas 96 NREM sleep was characterized by higher power of slow oscillations (Fig. 1C) and a wider dynamic 97 range of pupil diameter (Fig. 1F). REM sleep was distinct from the other states, with respect to 98 prominent theta oscillations (Fig. 1C), low muscle tone (Fig. 1D), higher brainstem LFPs power ( Fig.  99 1E) and fully constricted pupil (Fig. 1F). The higher power of brainstem LFPs during REM sleep was 100 preserved across animals (7 animals, 9 recordings) (Supplementary Fig. 2). 101

Figure 1. Population activity in the brainstem across the sleep-wake cycle.
A. A diagram of experimental approaches, showing the insertion of a silicon probe for extracellular recording in the brainstem and a screw for cortical EEG recording. Pupil dilation and EMGs were also monitored in a head-fixed condition.
B. An example of multiple electrophysiological readings across three behavioral states, including local field potentials (LFPs) in the brainstem (locally subtracted LFP signals), brainstem single unit activities (SUAs), cortical EEG, EMG and normalized pupil diameter. REM, rapid eye movement sleep; NREM, non-REM sleep; AW, wakefulness. C and E. Power spectrum density of cortical EEGs (C) and brainstem LFPs (E) across three behavioral states. Spectrum was computed during every 4-sec window. Errors indicate SEM. D and F. Distribution of EMG signals (root mean square) (D) and normalized pupil diameter (F) across three behavioral states. Pupil diameter was normalized as z-score.
Neuronal spiking activity in the 102 brainstem also demonstrated rich state-103 dependent properties (Fig. 1B) To assess state-dependent firing of 122 individual neurons in the brainstem on a 123 timescale of seconds to minutes, we 124 performed in vivo silicon probe recording 125 ( Fig. 2A)   C. Classification of functional classes. Firing rates across three behavioral states were normalized as z-score for individual cells, then a hierarchical clustering was applied. D. A diagram of an experimental approach for fiber photometry-based Ca 2+ imaging from pontine cholinergic neural populations in a freely behaving condition. E. An example of fluorescent signals across sleep-wake cycles. Fluorescent signals (470 nm) were normalized by off-peak (405 nm) signals. red, wakefulness, blue, NREM sleep, green, REM sleep. F. Distributions of fluorescent signals across three behavioral states. G. Group statistics of average signals from 12 recordings from 4 mice (F2,32 = 5.12, p = 0.012, one-way ANOVA). *, p < 0.05 with post-hoc Tukey's honest significant difference criterion. NREM sleep. The largest class (52.6 %) was REM-on neurons, which showed the highest firing rate 148 during REM sleep. Thus, we confirmed highly diverse state-dependent firing in the brainstem. 149 Because the recorded neurons were distributed across various nuclei in the brainstem, it was difficult 150 to determine their state-dependency in each nucleus. However, a subset of neurons was likely 151 recorded from the cholinergic system, namely the pedunculopontine tegmental nucleus and the 152 laterodorsal tegmental nucleus, which show AW-on or REM-on activity (Supplementary Fig. 3). To 153 verify this, we performed in vivo fiber photometry of Ca 2+ signals from pontine cholinergic neurons 154 by expressing GCaMP6s in freely behaving mice (4 animals, 12 recording sessions) (Fig. 2D).

155
Consistent with the data from in vivo electrophysiology, cholinergic populations showed larger 156 activity during REM sleep and wakefulness, compared to NREM sleep (F 2,32 = 5.12, p = 0.012, one-157 way ANOVA) (Figs. 2E-G). Therefore, although state-dependency of individual neuronal firing in the 158 brainstem is diverse, we also confirmed state-dependent and cell-type-specific firing. 159 160

Figure 3. Pupil dilation across the sleep-wake cycle and prediction of pupil dilation by brainstem populations.
A. Mean normalized (z-scored) pupil diameter across the sleep-wake cycle (n = 18, F2, 53 = 220.33, p < 0.0001. one-way ANOVA). ***, p < 0.0001, with post-hoc Tukey's honest significant difference criterion. can predict pupil dilation quantitatively. To address these issues, we analyzed datasets from head-168 fixed mice with either silicon probe recordings from the brainstem (6 animals, 6 recording sessions) 169 or the hippocampus (2 animals, 3 recording sessions), or field potential recording from the brainstem 170 with a bipolar electrode (6 animals, 9 recording sessions). 171 As previously reported (Yuzgec et al., 2018), mice in a head-fixed condition kept their eyes open, 172 allowing us to monitor pupil dilation across states along with cortical EEG and EMG. The effects of 173 behavioral states on pupil diameter was statistically significant (Fig. 3A, F 2, 53 = 220.33, p < 0.0001, 174 one-way ANOVA). More specifically, pupil diameter was constricted during REM sleep and dilated 175 during wakefulness. 176 With respect to pupil dilation dynamics across states ( Fig. 3B and Supplementary Fig. 4), pupil 177 diameter dynamically fluctuated during wakefulness and gradually constricted during NREM sleep.

178
Typically, 10-20 sec before REM sleep, the pupil diameter would further decrease and was fully 179 constricted during REM sleep, with rapid eye movement (Supplementary Fig. 4). 180 Taking advantage of the simultaneous neural population recording and pupil monitoring, we 181 examined how brainstem neurons can predict pupil dilation. First, we predicted pupil dilation based 182 on the activity of individual neurons (Fig. 3C) by applying a linear regression analysis. Because it 183 was expected that the preceding neural activity can better predict pupil dilation, we systematically 184 shifted the temporal relationship between spike trains and pupil diameter (see Methods). As 185 expected, most of the neurons showed asymmetric profiles of R 2 values (Fig. 3C). Although 186 individual profiles were diverse, the average profile showed predictive activity of brainstem neurons 187 for pupil diameter around 10 seconds in advance. We also predicted pupil diameter based on 188 simultaneously recorded brainstem neurons ( Fig. 3D) by applying a multiple linear regression 189 analysis. As with individual neurons, we observed an asymmetric profile of predictability. Thus, 190 changes in brainstem neural activity preceded those in pupil diameter. Thereby, brainstem 191 populations have predictive power for pupil diameter. For example, module 1 represented REM-on activity whereas module 2 was activated at the end of 199 NREM sleep and module 3 was most active during wakefulness. Indeed, the weights in each module 200 were consistent with state-dependency of individual neural firing (Supplementary Fig 5). 201 Besides modules capturing firing patterns across neurons, NMF also yielded activation coefficients 202 of these modules. We noticed that dynamics of these activation coefficients show predictive activity: 203 in the case of Figure 4A, the activation coefficients of modules 1 and 2 gradually built up during REM 204 and NREM sleep, respectively. Therefore, we hypothesized that brainstem population activity 205 exhibits not just state-dependency, but also predictive power for behavioral states (i.e., wakefulness, 206 NREM sleep and REM sleep). To test this, we took the activation coefficients profiles from three 207 modules and classified behavioral states by training a linear classifier, with systematic time shifting 208 (Fig. 4C). Brainstem populations showed predictive activity 10s of seconds before transitions to all 209 three behavioral states. To compare, we also performed the same analysis for hippocampal neural 210 activity (Fig. 4C). Although hippocampal neurons also had predictive power for several tens of 211 seconds, the profile was relatively short-lasting compared to that of brainstem neurons. Thus, 212 brainstem neurons have long-lasting predictive power for behavioral states compared to 213 hippocampal neurons. 214 Brainstem population activity underlying P-waves during NREM and REM sleep 215 On a timescale of seconds to minutes, brainstem neurons show diverse but specific state-dependent 216 firing and have predictive power for pupil dilation and behavioral states. To investigate brainstem 217 neural firing on a sub-second timescale, we focused on P-waves (Callaway et al., 1987;Datta, 1997).

218
Although these sub-second neural events in the brainstem have long been recognized, the 219 underlying neural ensembles and relations to other sleep-related oscillations are fully understood. 220 Taking advantage of our dataset, we first examined whether the mouse pons exhibits P-waves like 221 other mammalian species. We implanted a bipolar electrode in the pons (n = 16 recordings) (Fig.  222 5A) and monitored LFPs by subtracting signals. During REM sleep, we observed large amplitude 223 irregular neural events, which often appeared as a burst (Fig. 5B right). We also observed similar, 224 but isolated neural events during NREM sleep (Fig. 5B left). These neural events appeared more 225 often during REM sleep (p < 0.0001, two-tailed t-test) (Fig. 5C). Intriguingly, the frequency of these 226 events gradually increased during NREM to REM sleep transitions and decreased during REM sleep 227 to wakefulness transitions (Fig. 5D). Because these characteristics generally resemble to those in 228 other species (Callaway et al., 1987;229 Datta, 1997), we concluded that 230 these neural events are P-waves in 231 mice. 232 P-waves can also be seen in silicon 233 probe recordings (Fig. 5E). Similar 234 large-amplitude, irregular neural 235 events were observed in subtracted 236 and filtered LFPs (Fig. 5F) Addressing these questions would 266 provide insight into functions of P-267 waves. To this end, first, we 268 investigated the relationship between 269 P-waves and cortical EEGs (Fig. 6).  F. An example of P-wave, showing LFPs from a shank, filtered, subtracted LFPs, and multiple single unit activity.
G and H. Pooled peri-event time histograms of brainstem single units relative to P-wave timing during NREM (G) and REM sleep (H). Time zero is the timing of P-waves (trough time). Each peri-event time histogram is color-coded by normalizing the maximum firing rate for each cell. The order of single units was sorted by the peak timing in each state. 6B), which consisted of delta (1-4 Hz), theta (~7 Hz) and beta (15~30 Hz) frequencies. To examine 276 this trend further, we assessed the phase relationship between P-wave timing and cortical 277 oscillations (Fig. 6C). We found significant phase preferences of P-wave timing (p < 0.01, Rayleigh's 278 test). We further assessed this phase-locking activity by computing phase modulation index, which 279 is defined as the difference in proportions of P-waves between the preferred phase and opposite 280 phase by dividing phases into four bins (e.g., the higher phase modulation index reflects the larger 281 difference in the proportion of P-waves between two opposing phase bins) (Fig. 6D). We found larger 282 phase modulation at delta and beta ranges. 283 On the other hand, P-waves REM exhibited distinct associations with cortical oscillations (Figs. 6B-D). 284 We observed significant phase modulation at theta range (p < 0.01) (Fig. 6C), indicating that two 285 prominent neural markers in REM sleep, that is, theta oscillations and P-waves, are temporally 286 organized. 287 288 Next, we investigated the relationship between P-waves and hippocampal activity (Fig. 7). We 289 started by assessing the phase relationship between hippocampal LFPs and P-waves across 290 frequency bands (Fig. 7B). While the timing of P-waves was phase-locked strongly at theta range 291 (~7 Hz) in both sleep states, we also observed stronger phase modulations with high frequency 292 components during REM sleep. We also examined underlying spiking activity in the hippocampus 293 (Fig. 7C). Intriguingly, while a subset of hippocampal neurons fired most strongly around the timing 294 of P-waves during both NREM and REM sleep, the temporal order between hippocampal neural 295 firing and P-waves was state-dependent ( Fig. 7C): during NREM sleep, co-firing of hippocampal 296 neurons was followed by P-waves whereas P-waves were followed by burst firing in subset of 297 C. Examples of phase-histograms. Cortical EEGs were filtered at certain frequency bands and the proportion of P-waves elicited in each phase bin was calculated. *, p < 0. 01, Rayleigh's test.

D.
Phase modulations across frequency bands of cortical EEGs. The phase modulation index was defined as the proportion in the preferred bin (the bin with maximal percentage) minus that in the opposite bin (the bin 180° apart). Error, SEM.
E. Effect size of states across frequency bands. hippocampal neurons during REM sleep. To test the hypothesis that co-firing of hippocampal 298 neurons during NREM sleep may reflect sharp-wave ripples (SWRs), we detected high-frequency 299 ripple events based on hippocampal LFPs to assess the temporal relationship between ripples and 300 P-waves NREM (Fig. 7D). We found that ripple events preceded P-waves during NREM sleep. Thus, 301 P-waves are strongly associated with hippocampal activity in both sleep states. However, their 302 associations are state-dependent.

305
Although state-dependent neural ensembles have been intensively characterized in the cortex, little 306 is known about the brainstem. Here, we investigated state-dependent neural population activity in 307 the brainstem, primarily the pons, on two distinct timescales. On a timescale of seconds to minutes, 308 brainstem neurons show diverse state-dependent firing, with cell-type-specificity in pontine 309 cholinergic neurons. Brainstem activity can collectively predict pupil dilation as well as behavioral 310 states. The ability to predict behavioral states is longer lasting compared to hippocampal neurons. On a timescale of sub-seconds, we characterized P-waves in the mouse, with respect to underlying 316 neural firing as well as associated cortical activity. P-waves typically appear during REM sleep and 317 less during NREM sleep. P-waves in both sleep states are accompanied by synchronous firing of 318 brainstem neurons, suggesting that underlying local activity during P-waves is similar between sleep 319 states. However, their relationship to cortical neural events is state-dependent: the timing of P-320 waves NREM are phase-locked to various cortical oscillations and hippocampal SWRs precede P-321 waves NREM , suggesting that P-waves are part of the brain-wide neural events triggered by SWRs. 322 During REM sleep, P-waves are phase-locked most strongly at theta frequency in both the neocortex 323 and hippocampus. Crucially, P-waves precede firing in a subset of hippocampal neurons, suggesting 324 that P-waves may trigger brain-wide neural events. Thus, P-waves are part of the state-dependent 325 coordinated activity across the brain. 326

Technical considerations 327
State-dependent activity in the brainstem has been described over the past several decades by 328 using various methods. The present study utilized a silicon probe to monitor neural activity from 329 multiple neurons simultaneously at a high temporal resolution. This approach allowed us to (1) 330 quantify state dependency of brainstem neural ensemble dynamics on a timescale of seconds to 331 minutes and (2) characterize neural population activity underlying P-waves for the first time. However, 332 because silicon probe recording alone has a limitation to identify cell types, additional approaches, 333 such as Ca 2+ imaging (Fig. 2)  and Dan, 2016). Crucially, the asymmetric profile of the predictability for pupil diameter suggests that 341 the modulation of brainstem activity precedes pupil dilation, rather than simple correlations. 342 The long lasting predictability of brainstem populations for behavioral states is not trivial. Intriguingly still unknown, but we hypothesize that the modulation of neural activity in the brainstem occurs tens 348 of seconds before global brain state transitions from one state to another. In other words, each state 349 emerges from complex interactions across various regions of the brain. 350

P-waves in mice 351
Although PGO or P-waves have been studies since the 1970s in several mammalian species, to the 352 best of our knowledge, we are the first to characterize P-waves in mice. Given the growing 353 importance of the mouse as an animal model for sleep research (Herice et al., 2019), the 354 confirmation of P-waves in mice is important for further interrogation. 355 We have noticed several similarities and differences in P-waves between mice and other species. 356 First, the waveform of P-waves in mice is generally consistent with those in other species, such as 357 cats (Callaway et al., 1987;Jeannerod et al., 1965) and rats (Datta, 1997;Farber et al., 1980), 358 suggesting that underlying neural ensembles may be similar across species. Second, the frequency 359 of P-waves during REM sleep is generally consistent across species (Datta, 1997). However, we 360 have also noticed that the frequency of detected P-waves varied across our experiments. This may 361 be explained by either the variation of REM sleep quality or the variation of electrode positions. 362 Further analysis of P-waves across brainstem nuclei will provide insights into their relationship with 363 sleep homeostasis and the mechanism of P-wave genesis. Third, the temporal evolution of P-wave 364 frequency generally agrees between mice and cats: the frequency of PGO-waves gradually 365 increases before the transition of NREM to REM sleep in cats (Steriade et al., 1989). Although it was 366 weak, a similar tendency was observed in our recordings (Fig. 5D). Rather, the frequency of P-367 waves increases during REM sleep. This subtle difference may be explained by anatomical 368 differences between species (Datta, 2012). Finally, although P-waves appear more frequently during 369 REM sleep, it is important to note that similar neural events also occur during NREM sleep. Given 370 state-dependent interactions between P-waves and cortical oscillations (Figs. 6 and 7), the 371 mechanisms of P-wave genesis are likely distinct. 372

State-dependent global coordination of neural ensembles 373
The temporal correlation between P-waves and hippocampal theta rhythms during REM sleep is 374 consistent with previous studies in cats and rats (Karashima et al., 2004;Sakai et al., 1973). The 375 phase-locked activity with fast gamma (80-110 Hz) in the hippocampus may relate to the recent 376 observation that demonstrated the close association between local hippocampal theta and fast 377 gamma events and brain-wide hyperemic events (Bergel et al., 2018). Because a number of 378 hippocampal neurons fire immediately after P-waves (Fig. 7C), P-waves may play a role in the 379 regulation of hippocampal ensemble dynamics as well as these brain-wide events during REM sleep. 380 On the other hand, the picture during NREM sleep seems to be distinct. Because SWRs precedes 381 P-waves (Fig. 7D) and SWRs are known to be generated within hippocampal circuits (Buzsaki, 2015), 382 SWRs play a leading role in brain-wide sub-second neural events including P-waves. These state-383 dependent brain-wide neural ensembles imply distinct functional roles of sleep states. 384 consisted of wild-type, ChAT-IRES-Cre (JAX006410), or ChAT-IRES-Cre::Ai32 (JAX012569) on 391 C57BL/6 background. ChAT-IRES-Cre::Ai32 mice were used to identify pontine cholinergic neurons 392 post-hoc histological analysis. For brainstem silicon probe recordings, 6 animals were used (9 393 recordings). For hippocampal silicon probe recordings, 2 animals were used (4 recordings). For P-394 wave recordings, 10 animals were used, but 4 were excluded due to electrode mispositioning or lack 395 of histological data. 16 datasets were used for further analysis. For pupil monitoring, 14 animals were 396 used, but one animal was excluded due to their eye closure during recording. 18 datasets were used 397 for further analysis. For fiber photometry, 4 animals were used (12 recordings). The detailed 398 information of genotypes, age, sex, body weight and the number of recordings was provided in 399 Supplementary Table 1. 400

403
Body temperature was maintained at 37°C with a feedback temperature controller (40-90-8C, FHC). 404 Lidocaine (2%, 0.1-0.3 mg) was administered subcutaneously at the site of incision. Two bone 405 screws were implanted on the skull as electrodes for cortical EEGs and twisted wires were inserted 406 into the neck muscle as electrodes for EMG. Another bone screw was implanted on the cerebellum 407 as a ground/reference. 408 For pontine EEG recording, bipolar electrodes were bilaterally implanted in the medial parabrachial 409 nucleus of the pons (5.1 mm posterior, 1.2 mm lateral from bregma, 3.2 mm depth from brain surface).

410
The bipolar electrodes consisted of 75 or 100 µm diameter stainless wires (FE631309, Advent 411 Research Materials and FE205850, Goodfellow, respectively). The tip of two glued wires were 412 separated by 0.5-1.0 mm vertically to differentiate EEG signals. All electrodes were connected to 413 connectors (SS-132-T-2-N, Semtec) and securely attached on the skull with dental cement. A pair 414 of nuts was also attached on the skull with dental cement as a head-post. After the surgery, 415 Carprofen (Rimadyl, 5 mg/kg) was administered intraperitoneously. 416 For brainstem or hippocampal silicon probe recording, in addition to bone screws for cortical EEGs 417 and a ground/reference, a pair of nuts was attached on the skull with dental cement as a head-post.

418
After the head-post surgery, the animals were left to recover for at least 5 days. During the 419 habituation period, the animals were placed in a head-fixed apparatus, by securing them by the 420 head-post and placing the animal into an acrylic tube. This procedure was continued for at least 5 421 days, during which the duration of head-fixed was gradually extended from 10 to 120 min. A day 422 after the habituation period, the animals were anesthetized with isoflurane and a craniotomy to insert 423 silicon probe the brainstem and hippocampus was performed. A craniotomy on the left hemisphere 424 (4.0 mm to 5.5 mm posterior, 1.0 to 1.3 mm lateral from bregma) for the brainstem recording and on 425 the left hemisphere (2.0 mm posterior, 1.5 mm lateral from bregma) for the hippocampus recording 426 were performed, respectively. To protect and prevent the brain from drying, the surface was covered 427 with biocompatible sealants (Kwik-Sil and Kwik-Cast, WPI). In the following day, the animals were 428 placed in the head-fixed apparatus for electrophysiological recording as described below. 429 For fiber photometry experiments, cortical EEG and EMG electrodes were implanted as described 430 above and connected to a 2-by-3 piece connector (SLD-112-T-12, Semtec). Two additional anchor 431 screws were implanted bilaterally over the parietal bone to provide stability and a small portion of a 432 drinking straw was placed horizontally between the anchor screws and the connector. The viral 433 vector (AAV5-CAG-flex-GCaMP6s-WPRE-SV40, Penn Vector Core; titer 8.3x10 12 GC/ml) was 434 microinjected (500 nl at 30 ml/min) (Nanoliter2010, WPI) to target the PPT/LDT area (-4.5 mm 435 posterior, 1 mm lateral from bregma and 3.25 mm depth from brain surface). The micropipette was 436 left in the brain for an additional 10 minutes and then slowly raised up. An optic fiber cannula 437 (CFM14L05, Thorlabs) was then implanted 3 mm deep from the surface of the brain and all 438 components were secured to each other and the skull with dental cement. For pontine EEG recording, cortical EEG, EMG and pontine EEG signals were amplified (HST/32V-445 G20 and PBX3, Plexon), filtered (0.1 Hz low cut), digitized at a sampling rate of 1 kHz and recorded 446 using LabVIEW software (National Instruments). Recording was performed for 5 hrs from 9:00 to 447 14:00. For brainstem or hippocampal silicon probe recording, a 32 channels 4 shank silicon probe 448 (A4 x 8-5 mm-100-400-177 for brainstem recording or Buzsaki32 for hippocampus recording) was 449 inserted slowly with a manual micromanipulator (SM-25A, Narishige) into the brainstem (3.75 mm -450 4.3 mm from the brain surface) or the hippocampus (1.55 mm -1.85 mm from brain surface). Probes 451 were inserted perpendicularly with respect to the brain surface. Broadband signals were amplified 452 (HST/32V-G20 and PBX3, Plexon) relative to screw on the cerebellum, filtered (0.1 Hz low cut), 453 digitized at 20 kHz and recorded using LabVIEW software (National Instruments). The recording 454 session was initiated > 1 hr after the probe was inserted to its target depth, to stabilize neural signals. 455 Recording preparation started from 8:00 and terminated by 15:00. For verification of silicon probe 456 tracks, the rear of the probes were painted with DiI (∼10% in ethanol, D282, Invitrogen) before 457 insertion. 458

459
Pupil monitoring 460 In a subset of in vivo electrophysiological experiments under a head-fixed condition, pupil was also 461 monitored with an off-axis infrared (IR) light source (860 nm IR LED, RS Components). A camera 462 (acA1920-25µm, Basler Ace) with a zoom lens (M0814-MP2, computar) and an IR filter (FGL780, 463 Thorlabs) was placed at ~10 cm from the animal's left eye. Images were collected at 25 Hz using a 464 custom-written LabVIEW program and a National Instruments image grabber (PCIe-8242). The fiber photometry system consisted of two excitation channels. A 470 nm LED (M470L3, 468 Thorlabs) was used to extract a Ca 2+ -dependent signal and a 405 nm LED (M405L3, Thorlabs) was 469 used to obtain a Ca 2+ -independent isosbestic signal. Light from the LEDs was directed through 470 excitation filters (FB470-10, FB405-10, Thorlabs) and a dichroic mirror to the fiber launch 471 (DMLP425R and KT110/M, respectively), The fiber launch was connected to a multimode patch 472 cable (M82L01, Thorlabs) which attached to an implantable optic fiber on the mouse via a ceramic 473 mating sleeve (CFM14L05 and ADAF1, respectively). Light emissions from GCaMP6s expressing 474 neurons were then collected back through the optic fiber, and directed through a detection path, 475 passing a dichroic mirror (MD498) to reach a photodetector (NewFocus 2151, Newport). A National 476 Instruments DAQ (NI USB-6211) and custom-written LabVIEW software was used to control the 477 LEDs and acquire fluorescence data at 1 KHz. LEDs were alternately turned on and off at 40Hz in a 478 square pulse pattern. Electrophysiology signals were recorded at 1 KHz using an interface board 479 (RHD2000, Intan Technologies) and connected to the mouse via an amplifier (RHD2132, Intan 480 Technologies based on the sleep score. This initial step was critical to reflect the dynamic range of pupil 518 dilation/constriction in each experiment. After choosing the file, 20 frames were randomly selected 519 to manually mark the left and right edges of pupil. ImageJ was used for this manual marking. Using 520 these labeled frames, the deep convolutional neural network was then trained, and all video files 521 were processed to detect the left and right edges of pupil. After processing, visual inspection was 522 performed by generating a down-sampled (50-100 times) video clip. The same procedure was 523 applied across all recordings. To compute pupil diameter, the distance between the left and right 524 edges of the pupil was calculated across frames. The profile was filtered (low-pass filter at 0.5 Hz) 525 and z-scored. To compute eye movement, the middle point of pupil was determined and the distance 526 of the middle points between two continuous frames was calculated. The profile was then normalized 527 by the maximal value of the profile. From 20 pupil recordings (14 animals), 2 recordings from a single 528 animal were excluded due to eye closure during most of recording period (Supplementary providing the largest spike amplitude. All subsequent analysis were performed by using custom-546 written codes (MATLAB, Mathworks). 547 To categorize functional classes of single units, average firing rate for each behavioral state was 548 calculated and a hierarchical clustering approach with the Ward's method was applied. 549 To predict the pupil diameter from single unit activity, spike trains were filtered (band-pass filter 550 between 0.5 and 25 Hz) and then a liner regression analysis was performed. To evaluate the 551 goodness-of-fit of the linear model, R-squared value was calculated. The same process was 552 repeated by shifting the time relationship between spike trains and pupil diameter to determine an 553 optimal time window to predict pupil diameter from spike train. Then the sequence of R-squared 554 values were normalized by computing Z-scores. 555 To predict the pupil diameter from simultaneously monitored multiple single unit activities, spike 556 trains were filtered (band-pass filter between 0.5 and 25 Hz) and a linear regression model was 557 trained by using a regularized support-vector machine algorithm with 10-hold cross-validation. Then 558 cross-validated mean squared error (MSE) was computed. The same process was repeated by 559 shifting the time relationship between spike trains and pupil diameter. The sequence of MSEs were 560 normalized. 561 To decompose neural population activity into "space" (neurons) modules and activation coefficients 562 of these modules, we adopted non-negative matrix factorization (NMF) ( all elements of the matrix R and then decomposed the shuffled matrix like we previously 581 decomposed the original R. We quantified the quality of a decomposition using the variance 582 accounted for (VAF) (Clark et al., 2010). Starting with one module, we incremented the number of 583 modules until the VAF of the unshuffled data decomposition did not exceed 3/4 of the average VAF 584 of the decompositions of 100 shuffles. 585 To classify three behavioral states based on the activation coefficients, a linear classifier was trained 586 by fitting a multivariate normal density to each state with 10-fold cross validation. Then classification 587 performance was calculated. This procedure was repeated by shifting the time relationship between 588 the activation coefficients and sleep scores. 589 To detect P-waves, two EEG or LFP signals in the pons were subtracted and filtered (5-30 Hz band-590 pass filter). If the signals cross a threshold, the event was recognized as P-waves. To determine the 591 detection threshold, a 1-min segment of the subtracted signals was extracted from the longest NREM 592 sleep episode for a reliable estimation of stable noise level. The noise level was estimated by 593 computing root-mean-square (RMS) values in every 10 ms time window. The threshold was defined 594 as mean + 5 x standard deviation of the RMS values. The timing of P-waves was defined as the 595 timing of the negative peak. 596 The phase analysis was essentially the same as that described previously (Yague et al., 2017 instantaneous phase of each band was estimated from the Hilbert transform and the phase of P-604 wave occurrence was computed. To quantify the relationship between P-waves and EEG/LFP phase, 605 the percentage of P-waves elicited in each phase bin was calculated. The phase modulation was 606 defined as the percentage in the preferred bin (the bin with maximal percentage) minus that in the 607 opposite bin (the bin 180° apart). Rayleigh's test for non-uniformity of circular data was performed 608 to assess the statistical significance (p < 0.01) of the non-uniformity of the P-wave vs EEG/LFP 609 phase distribution using CircStats Toolbox (Berens, 2009). 610 To detect ripples in the hippocampus, LFP signals from the channel which detected spiking activity 611 were used. Band-limited signals at 140-250 Hz were computed by using a Kaiser finite impulse 612 response filter (see above). Two sequences of RMS values were calculated with two different time 613 window: 2 sec (long RMS) and 8 ms (short RMS). If the short RMS exceeds 4 times larger long RMS 614 for 8 ms, then signals were recognized as a ripple event. 615

Statistical analysis 616
Data was presented as mean ± SEM unless otherwise stated. Statistical analyses were performed 617 with MATLAB. In Figs. 2G and 3A, one-way ANOVA with post-hoc Tukey's Honest Significant 618 Difference (HSD) criterion was performed. In Fig. 4C, repeated measures ANOVA was performed. 619 In Fig. 5C, two-tailed t-test was performed. In Fig. 6C