# Experimental and computational analysis of REM sleep distributed cortical activity in mice

Mathias Peuvrier, Laura Fernandez, Sylvain Crochet, Alain Destexhe, Paul Salin

<figure><img src="image_1.png" /></figure>

## To cite this version:

Mathias Peuvrier, Laura Fernandez, Sylvain Crochet, Alain Destexhe, Paul Salin. Experimental and computational analysis of REM sleep distributed cortical activity in mice. Journal of Neurophysiology, 2025, 134(5), pp.1378-1395. 10.1152/jn.00344.2024. hal-05288972

HAL Id: hal-05288972

https://hal.science/hal-05288972v1

Submitted on 6 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Experimental and computational analysis of REM sleep distributed cortical activity in mice

Mathias Peuvrier<sup>a</sup>, Laura Fernandez<sup>c</sup>, Sylvain Crochet<sup>d</sup>, Alain Destexhe<sup>a</sup> and Paul Salin<sup>b</sup>

## ARTICLE INFO

**Keywords:**

REM sleep; Delta oscillations; Cerebral cortex; Computational models; electrophysiology

**New & noteworthy :**

Frequency and shape of REM sleep slow oscillations (SO) are different from SWS and wakefulness and are predominantly observed in the primary and secondary areas of the somatosensory cortex.

REM sleep SO alternate with periods of desynchronized activity of low amplitude.

REM sleep SO are poorly synchronized between cortical areas.

Cortical modeling suggests that spike frequency adaptation differs between cortical areas.

Modeling results are consistent with a cortical heterogeneity of cholinergic signaling during REM sleep.

## ABSTRACT

Although classically Rapid-Eye Movement (REM) sleep was thought to generate desynchronized activity similar to wakefulness, it was found that some brain regions could express Slow Wave activity (SWA), a pattern which was normally typical of slow-wave sleep. To investigate possible underlying mechanisms, we analyzed experimental recordings and introduced a computational model of mice cerebral cortex in REM sleep. We characterized the patterns of slow-wave activity across somatosensory and motor areas, and found that the most prominent REM-related SWA was present in the primary (S1) and secondary (S2) somatosensory areas, more rarely seen in motor cortex, and absent from prefrontal cortex or hippocampus. The SWA also tended to be synchronized in S1 and S2. We next investigated possible mechanisms by using a computational model of the mouse brain consisting of adaptive Exponential (AdEx) mean-fields connected together according to the mouse connectome. To compare with experimental data, the local field potential was calculated in each mouse brain region. To reproduce the experiments, we had to assume a heterogeneous level of adaptation in different cortical regions during REM sleep. In these conditions, the model reproduced some of the experimental observations in the somato-motor areas and the other cortical areas. We then used the model to test how the presence of SWA affected cortical responsiveness. Indeed, we found that the areas expressing SWA had diminished evoked responses, which may have participated in a diminished responsiveness during REM sleep.

## 1. Introduction

Classically, brain states such as wakefulness (W) and Rapid-Eye Movement (REM) sleep are thought to generate essentially desynchronized activity, while slow-wave sleep (SWS) typically generates synchronized slow oscillatory waves [1]. However, studies have shown that slow-wave activity (SWA) can also be seen during quiet wakefulness in some brain areas [2, 3, 4, 5]. SWA can also be observed in REM sleep, in human [6, 7], and in mice[8, 9]. Funk et al demonstrate local SWA in the superficial layers of primary areas of the mouse sensory cortex, while the activity in secondary sensory and associative cortical areas appears desynchronized [8]. These results suggest that SWA is present in all primary sensory areas and filters sensory information during REM sleep.

### Experimental recording of slow-wave activity

The slow waves of SWS are composed of delta oscillations of 2-4 Hz and the slow oscillations (SO <2 Hz). This slow oscillatory activity corresponds to the coordinated spiking activity of a large population of neurons, alternating be-

tween active (upstates) and silent (downstates) periods. Several studies demonstrate that the slow and delta oscillations are generated in the neocortex [10, 11, 12, 13]. Taken together, these oscillations make up the slow wave activity (SWA) which reflects the homeostatic processes of sleep. Thus with an increase of the sleep pressure, the power of SWA increases and conversely, in the course of SWS when the sleep pressure decreases, the SWA power decreases [14?].

In contrast, REM activity is associated with low amplitude higher frequency theta (5-10 Hz) and gamma (35-80Hz) activity [15]. Studies have observed phasic events of faster theta with strong gamma coupling in the hippocampus [16, 17, 18]. Thus the brain states of SWS and REM examined with EEG recordings appear entirely opposite. It should also be noted that cholinergic neurons in the brainstem and the basal forebrain that innervate the cortex are active during REM sleep, as in wakefulness, but not in SWS [19, 20, 21, 22, 23, 24, 13].

Advances in in vivo electrophysiological techniques (intracellular recordings, multisite local field potential (LFP) and spike recordings and high-density EEG recording) have revealed that in the neocortex W is characterized by substates alternating SWA and desynchronized activity during quiet and active wake respectively. These SWA in quiet wake share some common features with delta oscillations in SWS

<footnote><sup>a</sup>Paris-Saclay University, CNRS, Paris-Saclay Institute of Neuroscience (NeuroPSI), 91400 Saclay, France</footnote>

<footnote><sup>b</sup>University Claude-Bernard Lyon1, Forgetting processes and cortical dynamics' team, Centre de Recherche en Neurosciences de Lyon (CRNL), CNRS UMR 5292, INSERM U1028, F-69500, BRON, France.</footnote>

<footnote><sup>c</sup>Université de Lausanne, UNIL, Département des neurosciences fondamentales (DNF), Lausanne, Switzerland</footnote>

<footnote><sup>d</sup>Laboratory of Sensory Processing, Brain Mind Institute, Faculty of Life Sciences, Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland</footnote>

<footnote>peuvrierm@gmail.com (M. Peuvrier); Laura.Fernandez@unil.ch (L. Fernandez); sylvain.crochet@epfl.ch (S. Crochet); alain destexhe@cnrs.fr (A. Destexhe); paulantoine.salin@outlook.fr (P. Salin)
ORCID(s): 0000-0002-7942-3369 (L. Fernandez); 0000-0002-2965-8899 (S. Crochet); 0000-0001-7405-0455 (A. Destexhe)</footnote>
and are generated in sensory, motor and associative cortical areas [2, 3, 4, 5]. Intriguingly, SWA is also observed in REM in the superficial layers of primary areas of the mouse sensory cortex [8]. In contrast, the activity of the secondary sensory, motor and associative cortical areas would be desynchronized as in EEG and hippocampal recordings [8, 16]. This finding seems to be opposed to EEG work showing that REM activity is characterized by rapid oscillations of low amplitude. However, the results of Funk's study may help explain the REM "paradox": an EEG activity close to wakefulness that is associated with high sensory threshold due to synchronized slow activities in the primary sensory cortices. If SWA are present in these primary cortical areas they could filter sensory perception during REM like during SWS. Thus, this finding leads to reconsider the importance of SWA in REM as compared to SWS and W [8, 6, 25].

Advances in in vivo electrophysiological techniques (intracellular recordings, multisite local field potential (LFP) and spike recordings and high-density EEG recording) show that in the neocortex W is characterized by sub-states alternating SWA and desynchronized activity during quiet and active wake respectively. These SWA in quiet wake share common features with delta oscillations in SWS and are generated in sensory, motor and associative cortical areas [2, 3, 4, 5]. Intriguingly, SWA also occurs in REM in the superficial layers of primary areas of the mouse sensory cortex [8]. In contrast, the activity of the secondary sensory, motor and associative cortical areas is desynchronized as in EEG and hippocampal recordings [8, 16]. This finding is opposed to EEG work showing that REM activity is characterized by rapid oscillations of low amplitude. However, the results of Funk's study may help explain the REM "paradox": an EEG activity close to wakefulness that is associated with high sensory threshold due to synchronized slow activities in the primary sensory cortices. When SWA is present in these primary cortical areas it may filter sensory perception during REM like during SWS. Thus, this finding leads to reconsider the importance of SWA in REM as compared to SWS and W [8, 6, 25].

Very few studies have characterised the differences between SWA in REM and those in SWS and W. A detailed analysis of cortical SWA in REM versus SWS and wakefulness remains to be performed. Are REM slow waves different in amplitude and shape compared to SWS and W? Is the temporal dynamic of this slow activity similar in the areas where it is present during REM? What makes this slow activity possible? By recording LFP activity in several cortical areas using multi site recordings it becomes possible to directly address these questions and to determine similarities and differences in the scale of the individual slow wave as well as in the temporal and spatial distribution of SWA in REM compared to SWS and W.

Using multisite LFP recordings, we first analyze REM SWA in comparison to the SWA of SWS and W in the mouse somato-motor cortex. Since this SWA is considerable in this cortex, we examine whether it is correlated between these

cortical areas on the one hand, but we also determine whether it is correlated to other areas where faster oscillations are predominant such as the hippocampus. We demonstrate that REM SWA is characterized by Delta waves in primary and secondary cortical somatosensory (S1 and S2 respectively), as well as in primary motor areas (M1), and that SWA alternates with faster low voltage oscillations during REM episodes. This indicates a coordination in the different cortical regimes in place, and provides us new insights into the brain capacity for information integration in REM.

## Computational models of sleep

Until now, most modeling efforts regarding REM aimed at building models concerning the genesis, maintenance or regulation of this sleep state in relation to the other states of vigilance [26, 27]. Surprisingly, little has been done to model cortical activity during REM, compared to the numerous computational works concerning SWS. Models of cortical dynamics of sleep (and wakefulness) with fine temporal precision rely on neuronal network models with single neurons, of various degrees of detail. Those models, and the thalamocortical models in particular, were very fruitful in describing the slow oscillatory Up and Down regime of the SWS [28, 29, 30, 31]. In such models the neocortex is composed of a few thousand neurons that would represent one to two cortical columns, with a specific connectivity, often inspired by biological observations [32, 33]. However, expanding such models to the whole cortex would require impractical amounts of computational power to both run and analyze. This is why we considered a population models, abstracting the neuron activity to drastically reduce the computational cost, while remaining completely relevant and very detailed in regard to experimental LFP recordings. It is a mesoscopic approach that describes the average collective properties of neuronal populations, instead of simulating each individual neuronal dynamics [34].

Until now, most modeling efforts regarding REM have aimed at building models concerning the genesis, maintenance or regulation of this sleep state in relation to the other states of vigilance [26, 27]. Surprisingly, little have been done to model cortical activity during REM, compared to the numerous computational works concerning SWS. Models of cortical dynamics of sleep (and wakefulness) with fine temporal precision rely on neuronal network models with single neurons, of various degrees of detail. Those models, and the thalamocortical models in particular, have been very fruitful in describing the slow oscillatory Up and Down regime of the SWS [28, 29, 30, 31]. In such models the neocortex is composed of a few thousand neurons that would represent one to two cortical columns, with a specific connectivity, often inspired by biological observations [32, 33]. However, expanding such models to the whole cortex requires impractical amounts of computational power to both run and analyze. We therefore use population models, abstracting the neuron activity to drastically reduce the computational cost, while remaining completely relevant and very detailed in regard to experimental LFP recordings. It is a mesoscopic
approach that describes the average collective properties of neuronal populations, instead of simulating each individual neuronal dynamics [34].

Furthermore, this modeling approach is inspired by the work of [35] that produces important features of global cortical dynamics comparable to W and SWS in human brain-scale simulation. This approach uses The Virtual Brain (TVB) simulator (https://www.theyirtualbrain.org/tvb) to integrate population models (mean-fields) as the nodes of a network in which the edges (connections) are based on brain connectivity as defined from diffusion imaging experiments. In our model we will use the TVB mouse counterpart, The Virtual Mouse Brain (TVMB; [36]). TVMB can construct individual mouse brains from dMRI experiments or from tracer experiments realized by the Allen Institute. This allows us to use a detailed connectome of the mouse cortex that fits well the level of description we have in our biological dataset.

To model neuronal populations, we will use an AdEx mean-field model, as described in [37]. This biologically realistic mesoscale model is built to correctly predict the level of spontaneous activity observed in AdEx spiking neural networks with Asynchronous Irregular (AI) or Up and Down regimes. Two network regimes that are the computational analogs of biological Wakefulness and SWS sleep respectively [38]. In those models, the spike-triggered adaptation parameter is sufficient to switch from one regime to the other. It is interesting to note that this adaptation parameter is strongly related to the type of adaptation produced by high-threshold calcium-gated current [39]. This current is known to mediate the spike frequency adaptation in cortical pyramidal neurons but also to be deactivated by acetylcholine [40]. Therefore it is not surprising that in models with strong spike-triggered adaptation, the adaptation mediated by this current is strong, such as in a cortex without acetylcholine signaling, such as in SWS.

Using this model of the whole mouse cortex, we can simulate the transition from a Wake model (W-model) to a SWS model (SWS-model) by changing only the spike-triggered adaptation parameter. Yet there is no REM-model. Here we present a whole-brain model that we tune to fit LFP recording of mouse REM sleep. In particular, we focus on cortical dynamics in the somatosensory primary and secondary cortices (S1 and S2), the primary motor cortex (M1), the medial prefrontal cortex (mPFC) and the hippocampus (dCA1), with a particular interest regarding the REM SWA as recently observed [8, 7, 41]. With these considerations to build our REM-model, we propose a simple topological organization of the spike-triggered adaptation parameter in the model, a gradient of adaptation similar to the gradient of acetylcholine proposed by Nazari et al [25]. This organization allows the simulation of different patterns of activity across the cortex that reproduce a part of the diversity of activities we observe experimentally. Finally we use this new model to predict the potential cortical response to stimulation.

## 2. Methods

### 2.1. Electrophysiological signal analysis

The data analyzed in the present study come from a database acquired in unanesthetized head restrained mice recorded by L.F. and S C. Part of this dataset has been recently analyzed and published to characterize the cortical slow waves of wakefulness and SWS [5].

#### Animal preparation and surgery

Briefly, nineteen C57BL/6 male mice (Janvier SAS, St. Berthevin, France), aged of 6 to 8 weeks at the time of surgery, were used for the study. All procedures were approved by the University of Lyon 1 Animal Care Committee (project DR2013-4) and were conducted in accordance with the French and European Community guidelines for the use of research animals.

We analyzed the data of seven mice out of the nineteen available, as these were the only mice with recordings at the specific locations of interest: the barrel field of the primary Somatosensory cortex (S1), the secondary Somatosensory cortex (S2), the primary Motor area (M1), the prelimbic area of the medial prefrontal cortex (mPFC), and the area CA1 of the dorsal hippocampus (dCA1).

In addition to the LFP microelectrodes, two surface EEG electrodes were implanted over the parietal (AP 2.0, lat 1.5) and frontal areas (AP 5.3, lat 1.5) onto the dura mater of the contralateral hemisphere, two EMG electrodes were inserted bilaterally in the neck muscles, and two silver wires were inserted on both sides of the cerebellum for reference and grounding. Then, to allow for painless head-fixation during recording sessions, a light-weight metal head-post was cemented to the skull.

#### Habituation and recording

Habituation and recording sessions were performed in the second part of the rest period of the mice, when the sleep pressure was limited, between 12AM and 5PM, during the light period in a semidark light condition, with ambient room temperature (22±2°C), in an electrically shielded recording chamber.

After recovery from surgery, mice were progressively habituated to head-restraint conditions with repetitive daily training. The recording sessions began when the first SWS episodes were observed.

#### Sleep scoring

Epochs of wakefulness, SWS and REM were scored using EEG and EMG. REM sleep was identified by theta oscillations on the EEG and neck muscle atonia on the EMG. Artifacts, drowsiness periods and intermediate sleep were discarded. Therefore, each episode of vigilance states (W, SWS, REM) showed a clear separation.
<table><caption>**Table 1**</caption><thead><tr><th>Area</th><th>S1</th><th>S2</th><th>M1</th><th>A1</th><th>V1</th><th>PtA</th><th>mPFC</th><th>dCA1</th></tr></thead><tbody><tr><td>Antero-Posterior</td><td>2.8</td><td>3.1</td><td>5.1</td><td>1.3</td><td>0.5</td><td>1.8</td><td>5.7</td><td>1.3</td></tr><tr><td>Lateral</td><td>3.2</td><td>4.0</td><td>1.8</td><td>4.0</td><td>2.5</td><td>1.8</td><td>0.3</td><td>2.0</td></tr><tr><td>Depth</td><td>0.5</td><td>0.8</td><td>0.4</td><td>0.6</td><td>0.4</td><td>0.4</td><td>1.85</td><td>1.3</td></tr></tbody></table>

LFP microelectrodes implantation coordinates. In bold is highlighted the areas that will be studied in this article.

## LFP dataset

For each mouse we had sets of brain state episodes, determined with the sleep scoring, corresponding to specific brain states with the following data: EEG, EMG and a set of 5 to 6 LFP (all sampled at 1kHz). Electrodes implanted in each mouse corresponded to the configuration presented in table 1. In this study, we used only the first configuration, with 7 mouse recordings in: S1, S2, M1, mPFC, dCA1. Recording sites in areas S1, S2 and M1 were mainly located in layers 2 and 3. This allowed the study of interactions between S1 and other areas of the somatomotor cortex. Quantitatively this database included 341 episodes of W (total duration of 14866 s, mean episode duration: 43.6 ± 36.05), 318 episodes of SWS (total duration 12274.08s, mean episode duration: 38.6 ± 25.1s) and 120 episodes of REM (total duration 7418.68s, mean episode duration: 61.82 ± 40.1s).

## Spectral analysis

To compute the spectrogram of LFP, from the dataset or from the model, we used the Welch method. The power spectral density was estimated using a Hamming window of 4 seconds with 2 seconds overlapping, with the exception of 1000ms windows with 300ms overlap for the phasic theta analysis. The power spectrum was expressed in $V^2$. One spectrogram was computed for each brain state episode recorded, and they were used to compute band power and to find frequency peaks in the delta band.

The maximum frequency peak was the highest peak in the frequency range. The absolute delta power was the integral of the area, between the frequency bands and below the periodogram line. We approximated this area using the composite Simpson's rule, which decomposes the area into multiple parabolas and then sums their areas.

## Delta activity in the signal

In order to understand the delta activity, we first detected the individual delta cycles in the signals. Then we detected delta bouts, defined as short periods (from 1 to 10 seconds) with high density of delta cycles. The different steps of this analysis are shown in Supplemental Fig S1 (10.5281/zenodo.15063947), where we have the raw signal, the delta cycle identification and the delta bout determination.

## Bycycle algorithm

To detect and analyze Delta cycles in our data, we used a method inspired and based on the Bycycle algorithm (for

cycle-by-cycle analysis) and python library [42]. In this method the LFP signal was segmented into cycle periods using a multi-step process. The signal was lowpass filtered, and the zero-crossings were indexed, then, on the raw data peaks and troughs were identified as the minima and maxima between two zero-crossing time points. With these, a cycle was defined from trough to trough, consisting of an initial trough, a peak, and a final trough. Using these points, we computed the rise duration as the time between the initial trough and the peak, and the decay duration as the time between the peak and the final trough. The rise-decay symmetry was then computed as the rise duration divided by the cycle period (the total duration between the two troughs). In the original use of the bycycle algorithm, it was proposed to validate a putative oscillatory cycle when it was surrounded by similar cycles, taking part into a larger, sustained oscillatory bout. However, delta cycles were frequently interspersed with faster activities, and we observed isolated delta cycles. Thus, we used permissive settings from the original library, but added another threshold as compensation and redefined the detection of Delta bouts.

## Determination of delta cycle

To find delta cycles we used the algorithm to compute cycles, from trough to trough, centered on the peak. For the initial determination of the zero-crossings, a FIR filter with a passband between 2 and 4.5 Hz and filter length of three cycles of the lower cutoff frequency (2 Hz) using a Hamming window was applied. Then, using a simple amplitude threshold, cycles with amplitude high enough (above median amplitude) compared to the rest of the episode signal, were classified as putative delta (Supplemental Fig S1 (10.5281/zenodo.15063947), blue line) without regard to other detected cycles. This threshold was not restrictive at all and captured many false positives. But it was only used to calculate the actual amplitude threshold, as the median amplitude value for all putative delta cycles found in S1 and S2 during SWS for the specified individual This second decision parameter (Supplemental Fig S1 (10.5281/zenodo.15063947), purple line) restricted the delta cycles to those above 80% of the amplitude of the Slow wave sleep Delta cycles recorded in the primary and secondary somatosensory areas for each mouse. The cycle-by-cycle analysis of delta oscillations was progressively optimized and verified by direct observation of the dataset by 2 authors (MP, PS). While this method was effective at detecting a significant number of Delta cycles in the signal, it still missed few (usually isolated) delta cycles.

## Delta bouts

To understand the repartition of delta cycle in the signal, we categorized the LFP signals in delta and non-delta bouts. Those bouts are short periods (usually about 2 to 4 seconds) of signal characterized by prominent Delta activity, or its absence. To determine the temporal distribution of Delta, we describe LFP signals in delta and non-delta bouts. Validated delta cycles that were close enough (less than 800ms) were gathered together, adding non detected cycles if they
were surrounded by delta cycles and that at least 75% of the bout was made of validated delta cycles. Then those periods were validated as delta bouts if they lasted at least 1 second and were constituted of at least 3 delta cycles (Supplemental Fig S1 (10.5281/zenodo.15063947), bold green). Non-delta bouts were defined as periods of at least 1 second without any validated delta cycle. Thus, isolated delta cycles and very short non-delta portion of signal were excluded of this categorization as bouts.

To assess Delta bouts co-occurrence across S1, S2, and M1, we quantified their temporal overlap. For each area pair, we created boolean arrays (1 for Delta bout presence, 0 for absence) spanning the episode duration. We then calculated the percentage of time Delta bouts in one area overlapped with those in another by summing these boolean arrays. The overlap was expressed as the proportion of concurrent Delta activity relative to the total active time of the first area. We analyzed all area combinations across 120 REM episodes from multiple mice.

## Statistical analysis

Statistical analyses were performed using R software and the "psych" and "rstatix" libraries. For figure 3B and C, an ANOVA was performed for S1, S2 and M1, followed by a post hoc analysis with Bonferroni correction for multiple comparison. Complete statistical tables can be found in supplementary. Statistical comparisons were based on individual sleep state episodes as the unit of analysis. For figures 1 (B,E,F), 2 (B,D,E) and 3 (B,C,D,E) we grouped animals together, but for all statistical tests we took into account the inter-animal variability.

Estimation plots: Data were visualized using estimation statistics created with the DABEST python library [43]. These plots quantified effect sizes and their precision, rather than relying solely on null-hypothesis testing. In these plots (Fig. 1C,D and Fig. 2C), individual data points show mean values per mouse, with group mean ± SD indicated by gapped lines. The lower panel displays the mean difference between conditions, with its 95% confidence interval (black dot and vertical line) and sampling distribution (colored distribution).

## 2.2. Whole brain modeling

The whole brain models we discuss were built at the neuronal population resolution and implemented in The Virtual Mouse Brain. It consisted of a network of AdEx mean fields connected together following Connectome-based information from the brain Allen institute (https://alleninstitute.org/). To keep this section concise, we omitted some mathematical details regarding the mean-field model or the TVB implementation. The corresponding information can be found in [37] and [35, 44] respectively.

All simulations presented in the paper were conducted for ten seconds, with the first two seconds (the simulation initiation) being systematically excluded. The results presented for the brain state model in the figure 4 were results from a single simulation, while results concerning the stimulation part (Fig.5) were the average of 10 simulations (See the related part of the method).

## The Virtual Mouse Brain

We used the simulation framework The Virtual Mouse Brain (https://www.thevirtualbrain.org/tvb) to perform detailed connectome-based mouse brain simulations. With the Allen Connectivity Builder pipeline, experimental connectivity data from dMRI or tracer experiments were used to build a precise mouse connectome. The connectome was a square matrix (origin to target structure and vice versa) informed with connection strength and tract length. The Allen Connectivity Builder also provided a 3D matrix which represented the volume of the mouse brain to create a region volume mapping [36]. This framework allowed us to virtualize a mouse brain with a population model at the nodes corresponding to a cortical column or a brain area, and edges corresponding to connection strength between regions. We used a connectome composed of 104 cortical nodes [45]. We selected this connectome over more detailed ones that included thalamus and other subcortical nuclei for two principal reasons. First we did not have simultaneous recordings including subcortical areas in our dataset. Second, the population model we decided to use was built to match the dynamics of cortical spiking neuronal network models. Therefore more complex connectomic dataset were not more relevant regarding our biological dataset while it could have induced some errors in the network dynamics.

## Mean field model

The population model used was an AdEx mean-field model as described by [37]. This mean-field model, based on a previous Master Equation formalism [46] was built to match the spontaneous activity of networks of 10000 Regular and Fast Spiking AdEx (Adaptive Exponential) neurons. It took into account the nonlinear effects of the conductance-based synaptic interactions as well as the spike-frequency adaptation of excitatory regular spiking neurons. The mean-field had two neuron populations with different intrinsic properties. One was excitatory, called the “regular spiking” (RS), and was affected by the spike-frequency adaptation, as cortical pyramidal neurons of the somatomotor cortex. The other was an inhibitory population called “fast spiking” due to its usually higher firing frequency as cortical parvalbumin GABAergic interneurons. The dynamics of the two population firing rate and their covariance were computed as a function of the transfer function (TF; the firing rate of one population as function of the inputs it received), which was itself dependent on the population adaptation variable W. This population adaptation had two important parameters: the sub-threshold adaptation *a* and the spike-triggered adaptation *b*. It should be noted that this parameter *b* affected only the excitatory population. Another difference in the dynamics of the two populations was the inputs considered in their transfer function. Indeed, only the excitatory population was affected by external noisy drive, while both populations were affected by external inputs (such as thalamic inputs).
This mean field described the population dynamics of a single, isolated cortical column. The parameters used in our study were the same as in previous ones. They produced biologically realistic and sustained brain dynamics. With those parameters we were able to switch from an AI to an Up and Down model regime changing only the spike triggered adaptation parameter *b*. Therefore, except this single parameter *b*, all parameters were the same in every node and every simulation presented.

## Networks of mean-field models

In order to model large brain regions, we represented interconnected cortical columns as a network of mean-field models. The inter-cortical mean-fields connections were done by excitatory interactions only, while inhibitory connections remained locally within the column.

Therefore, once in the network, the dynamics of the single mean-field was changed as it was now impacted by the excitatory input exerted from multiple excitatory populations. The integration of all excitatory inputs was computed using the distance between the origin to target node and the axonal propagation speed to account for the delay of axonal propagation. Thus it was at this step that the connectomics information impacts and shaped the whole model dynamics.

## REM model spike-triggered adaptation parameter

The REM model tried to account for the diversity of activity observed experimentally. We aimed to produce a simulation in which some areas were marked by a Delta activity intertwined with periods of faster oscillations (in different proportions according to the area), while other areas would be closer to an 'Awake-like' activity. We hypothesized that regional variation of adaptation could be a key to get such a model as it was the parameter used to switch from Wake like to SWS like in the reference model.

In order to keep the model easy to use and describe, we limited the variability of adaptation to only three values. The different adaptation values used had a conductance of 5 nS (as for our reference W-model), 30nS and 60nS (as for our reference SWS-model). The topological organization of this parameter was a decreasing gradient, close to the observation of acetylcholine suggested in [25], but centered on the S1 barrel cortex to match our observation of the decreasing prevalence of delta activity in the LFP going from S1 to S2 and to M1. We placed the maximum adaptation level (60nS) centered at S1 barrel cortex. At 12mm from this center (encompassing S2) the adaptation was decreased to 30nS, and then, 12mm further away (encompassing M1) it was decreased to 5nS. We added some exceptions to this method of application of the *b* parameter regarding hippocampal (field CA and dentate gyrus) and prefrontal (anterior cingulate and prelimbic areas) areas that were parameterized with *b* = 5nS disregarding their position in the model. We thought that this setup about the hippocampus and the mPFC was more in accordance with the literature [15], although we did not look further into optimizing the activity of those regions.

## REM shuffled as control

With our REM model we introduced a mirror model that was used as a control: the REM shuffled. This model was built to have exactly the same number of nodes with the different values of adaptation. However in this model the topological arrangement was randomized, losing the spatial features of the model for the adaptation. The goal of this model was to look at the importance of this spatial arrangement of adaptation parameter *b*. We used it to question if the phenomenon needed the organized recruitment of multiple nodes, or if any node with high enough adaptation could have produced the delta activity disregarding (or with randomized) the average input of the other cortical areas. Most results concerning this model were not shown to keep the results clear, although it was noted that the results of those REM shuffled models were always more variable than any other model.

523

524

525

526

527

528

529

530

531

532

533

534

535

536

537

538

539

## LFP measurement

The estimate of the LFP signal of each node was computed with a Kernel-based method based on previous work [47, 48]. First demonstrated in a neuron network and then extended to mean-field models, it related the spiking of neurons to the LFP through efferent synaptic connections. The LFP was calculated as the convolution of the unitary LFP waveforms modeled from the spiking activity of the field. This method also took into account the relative contribution of the neuronal population to the signal according to the distance between the origin of the unitary LFP and the recording site.

540

541

542

543

544

545

546

547

548

549

550

551

552

553

554

555

556

557

558

559

## Analysis of the delta activity in the model

To perform a similar analysis on the model LFP data as we did with the experimental data, we ran a one minute long simulation of the REM model. To determine the delta cycle we used a single threshold step with manually adapted parameters to fit the model data (frequency band = 0.5-4Hz, amplitude fraction = 0.35, amplitude consistency = 0.45, period consistency = 0.5, monotonicity = 0.65). Analysis was performed on all 104 nodes independently and the results were aggregated for nodes with the same adaptation value.

560

561

562

563

564

565

566

567

568

569

## Stimulation

With TVMB we used an external input in order to stimulate our model. The stimulus was produced by an external node and consisted of a single pulse of 20ms starting 3 seconds after the beginning of the simulation and repeated every 1.25sec (±250ms) and Supplemental Fig S5 (10.5281/zen-

570

571

572

573

574

575
odo.15063947)).

The connection strength between the external node and its target was modulated to control the intensity or strength of the stimulus. A range of connection strengths was evaluated by stimulating S1 under various state conditions. Analysis of the evoked potential kinetics revealed that beyond the connection strength of $1 \times 10^{-3}$ nS the peak response latency remained consistently minimal. Consequently, this connection strength was used for all results presented in the main article.

The stimulation experiment consisted of a set of ten simulations in which the stimuli was provided to one node only (i.e. simulating the electrical stimulation of neurons in a cortical area), from the right hemisphere, with 5 or 6 stimulus repetition. For each we therefore had a set of above 50 Stimulus Evoked Potentials (SEP). We performed this stimulation experiments for the 5 areas we focused on this article (S1, S2, M1, mPFC (its equivalent was the anterior cingulate in the model), dCA1) in each model state (AI, UpDown, REM-model, REM-shuffled). The results presented were the average of the SEP obtained in each condition.

## 3. Results

We began by analyzing the occurrence and distribution of slow-waves in REM sleep and then introduced a computational model to account for these results.

### 3.1. Analysis of slow-waves in REM sleep

#### Primary somatosensory area signal in different brain states

Figure 1 shows the three states of the wake-sleep cycle in mice, wakefulness (W), slow-wave sleep (SWS) and rapid eye movement (REM) sleep observed in local field potentials (LFP) recorded by microelectrodes in the mouse primary somatosensory area (S1). W and SWS classically displayed two drastically different activity states, asynchronous and synchronized slow-waves respectively as in EEG recordings (data not shown). However, during REM sleep, figure 1 shows that LFP activity is much less homogeneous than during SWS and W. While asynchronous activity usually dominated in EEG and in some LFP recordings of cortical areas, other brain areas exhibited local slow-wave activity, in middle and superficial layers of primary cortices, as pointed out previously [8]. The differences between these states was also apparent in the power spectral density (PSD) of the LFP signals (Fig. 1B). Here, the broad low-frequency peak, culminating in the Slow Oscillation (SO) band (<1.5Hz), during SWS contrasted with the absence of clear peak in W, while in REM, a low-frequency peak appeared around 2-4 Hz (Fig. 1C). The REM slow-waves seemed thus faster than the classic slow wave complexes of SWS. This difference in frequency was characterized further in respect to the Delta component of the slow activity. Traditionally considered characteristic of SWS, the power in the delta frequency band measured very similar, if not higher, in REM (Fig. 1D).

However, the shape of the Delta cycles detected on S1 LFP signals appeared different in REM compared to SWS (Fig. 1E). In regard to the average rising and decay duration of the Delta cycles, the delta rise time was longer for REM oscillations than for SWS oscillations (N delta cycles and densities in S1 during SWS = 15963 cycles, 1.3 cycle/second; REM = 11171 cycles, 1.5 cycle/second. Delta mean rise time: REM = 185.43±13.06ms, SWS = 167.66±12.35ms student t value = 13.17, p ≈ 1e-33), while delta decay was shorter for REM oscillation than for SWS oscillations (Delta mean decay time: REM = 131.09±11.9ms, SWS = 146.19±11.80ms student t value = -11.89, p ≈ 2e-28). Therefore there was a clear difference in the ratio between rise and decay time of Delta cycles in REM compared to SWS. As illustrated in figure 1F, the distribution of rise and decay time of Delta cycles in the two sleep states overlapped poorly, and REM delta was more biased toward high rising time and low decay time. Rise and decay time of delta cycles that were occasionally observed during W (i.e. quiet wake see [5]) were shorter than those of the SWS and REM. It should be noted that our algorithm was not tuned to optimize the detection of W delta, and thus might be less relevant than what we found for REM sleep.

#### Delta waves in different cortical areas in REM

To investigate the spatial distribution of REM slow-waves, we next compared LFP recordings in the different cortical areas during REM sleep.

Figure 2 shows the typical LFP activity observed during REM in simultaneously recorded cortical areas. Major differences in the amplitude as well as in the oscillation rhythms were observed. The area S1 and S2 presented bouts of large amplitude delta oscillations, intermingled with periods of smaller and faster signal fluctuations. In the area M1 we also found some large amplitude Delta waves, but sparser than those of the somatosensory areas, and most of the signal was of small amplitude with no clear prevalent oscillation. Finally, the LFP of the dCA1 area of the hippocampus showed a sustained theta oscillation along the REM episodes, which is a typical signature of REM in rodents [49]. Almost no delta oscillations were found in the dCA1 LFP. The spectral analysis of LFP confirmed these differences between the activities of the cortical areas recorded (Fig 2.B). The spectra of S1 and S2 were very similar with a clear peak in the Delta frequencies. In the area M1 there was also a delta peak above the 1/f trend of the LFP spectrum but its amplitude was considerably smaller than that of the somatosensory areas. LFP spectral analysis of the area mPFC did not have a delta peak and therefore the analysis of delta oscillations in this area was not pursued. On panels D and E of the figure 2 were shown the averaged Delta cycles centered by their peak, and the distribution of the rise and decay duration of those Delta cycles. From those, it appeared that the delta cycles found in S1 and S2 were extremely similar, with a bias toward long rising and short decay time, while Delta cycles from M1 were less stereotyped. They had a tendency to have a longer time decay, which might be due to the fact
<figure><img src="image_3.png" /><figcaption>**Figure 1:** Delta in the somatosensory cortex across brain states (Figure colors: blue for W, red for SWS and green for REM). (A) Twenty seconds of LFP examples recorded in area S1 during sleep and wakefulness. (B) Mean power spectrum is computed from all brain state episodes recorded in S1 for all mice considered in this study. The band shown around the mean is the standard error of the mean. The grey area corresponds to the delta band, between 1.5 and 4.5 Hz. (C) Frequency peak in the SWA band (0.5-4.5 Hz) in S1 across states, shown as estimation plots. Upper: individual values and group statistics. Lower: mean difference from REM with 95% confidence interval. (D) Delta band (1.5-4.5 Hz) power amplitude in S1 across states, shown as estimation plots. Upper: individual values and group statistics. Lower: mean difference from REM with 95% confidence interval. (E) Average Delta cycle, from all individual cycles in S1 across states. Cycles are aligned at the beginning. Lines show mean ± SEM, with line opacity indicating the proportion of cycles contributing (minimum at <30% of total cycles). (F) Scatter plot of Delta cycles rise and decay time in S1 across states: Each dot corresponds to the mean rise and decay time of delta cycles detected during one episode of the three brain states. This figure shows the dispersion of Delta cycles rise-decay symmetry.</figcaption></figure>

that Delta events in M1 were much more isolated and rarely immediately followed by another event.

### Temporal distribution of REM-sleep delta bouts

Next, to determine the temporal distribution of REM-sleep delta waves, we detected successive delta cycles and found that periods of delta waves alternate with periods of asynchronous activity which we characterize in figure 3 as

the delta and non-delta bouts.

The detection of Delta bouts allowed us to discriminate periods (of at least a second) that were dominated by Delta activity from periods without any Delta activity, while excluding periods that would be considered as mixed-activity. With this approach we could describe at least 80% of any of our REM signals. This was observed in figure 3A as the thick green and red lines, that represented Delta bouts and

692

693

694

695

696

697

698

699
A

<figure><img src="image_4.png" /></figure>

B

<figure><img src="image_5.png" /></figure>

C

<figure><img src="image_6.png" /></figure>

D

<figure><img src="image_7.png" /></figure>

E

<figure><img src="image_8.png" /><figcaption>**Figure 2:** REM sleep delta oscillation in the cortex (Figure colors: Purple: S1, dark blue: S2, light blue: M1, brown: dCA1, orange: mPFC)</figcaption></figure>

A) 20 seconds of simultaneous LFP recording in 5 different areas during REM. Signals are shown with a common scale. Delta cycles detected in the signal are highlighted in green. (B) Average power spectrum for the different areas for all REM episodes. Grey square shows the 1.5-4.5Hz delta band. (C) Delta band (1.5-4.5 Hz) power amplitude during REM sleep across brain areas, shown as estimation plots. Upper: individual values and group statistics. Lower: mean difference from REM with 95% confidence interval. (D) Representation of the mean delta cycles, centred around their peak, over all detected cycles in the somatosensory and motor area. Line intensities represent the portion of delta cycle from the total pull recorded at specific area used for each datapoint. (Approximate number of cycles considered in each area: S1: 10000, S2: 8000, M1: 5000). (E) Scatter plot of REM Delta cycles rise and decay times in the somatosensory and motor area. Each dot corresponds to the mean rise and decay time of delta cycles detected during REM episodes recorded in areas exhibiting Delta activity in REM.
Non-Delta bouts respectively, occupied the major part of the signals. While in S1 the Delta bouts occupied almost half of REM episodes and lasted about 2.9 seconds, it described only 7.36% of the LFP signal of M1 in REM and lasted only 2 seconds on average. The result was intermediate for S2, where Delta bouts accounted for approximately 30% of the recording duration during REM sleep. In S2, the average duration of Delta bouts was 2.49 seconds.

Next, we determined whether Delta bouts occurred simultaneously in these different cortical areas (Fig. 3B and Supplementary Tables S1-1 and S1-2 (10.5281/zenodo.15063947)). We thus assessed the level of temporal overlap of the delta bouts in the area S1, S2 and M1. First, we computed at the overlap of S1 on S2, the distribution peak was around 0.4, meaning that for about 36.1%(±17.9%) of S1 Delta bouts duration, we simultaneously observed delta bouts in S2. This was notable considering that S1 Delta bouts were longer and occupied a larger part of the REM episodes. The overlap of S2 on S1 was even higher (with 51.9% overlap (±20.5%), student t-test p value = 1e-09), indicating that more than half of the time S2 spent in Delta bouts coincided with S1 Delta bouts. In contrast, the distribution of overlaps from either S1 or S2 on M1 Delta bouts was much smaller (respectively 9.6% (±16.4%) and 7.9% (±13.9%), and we rarely observed Delta bouts in M1 during Delta bouts of S1 or S2. However, M1 Delta bouts showed higher temporal overlap with S1 and S2 (26.4% (±34.8%) and 21.5% (±32.4%) respectively). This asymmetry reflected that while M1 had fewer Delta bouts overall, they tended to co-occur with Delta bouts in S1 and S2 (Bonferroni test, p < 0.001, see supplemental table S1-2 (10.5281/zenodo.15063947)).

Since S1 and S2 are characterized by numerous long delta bouts, we performed an additional analysis (results not shown) to verify that the overlaps observed between delta bouts are not entirely due to chance. We therefore compared the observed overlaps between delta bouts with a null distribution obtained by mixing the positions of the delta episodes in each cortical area for 1,000 iterations, while keeping their number of delta bouts and their duration constant. We then computed the z-score relative to this null distribution, as well as the associated p-values using a two-tailed normal distribution test. This procedure revealed a significant deviation from chance, with a z-score ranging from 2.69 to 6.40 (p-value: 10e-10 to 7.1e-3) for all pairs of cortical areas, demonstrating a synchronization of delta bouts between these areas rather than random co-occurrence.

Next, using spectral analysis, we examined the changes in the delta band power during Delta and non-Delta bout periods. Figures 3C-E also showed the change in delta power in the different areas (two way ANOVA comparing Delta band power according to the Delta/non-Delta bout "state" and the area, see supplemental table S1-3 and S1-4 (10.5281/zenodo.15063947)).

Figure 3C, as expected from our previous results (Fig 2B and C), shows a significant difference in Delta power for the different recorded areas (p=1.34e-05, ges=0.464). In the

<table><caption>Table 2</caption><thead><tr><th></th><th>S1</th><th>S2</th><th>M1</th></tr></thead><tbody><tr><td>Delta bouts</td><td>47.33%</td><td>30.373%</td><td>7.362%</td></tr><tr><td>Non-Delta bouts</td><td>41.441%</td><td>53.25%</td><td>84.283%</td></tr><tr><td>Mean Delta bout duration</td><td>2.951s</td><td>2.494s</td><td>2.006s</td></tr></tbody></table>

Percentage of somatomotor REM sleep LFP signal characterized by Delta bouts and Non-Delta bouts.

In the spectral analysis revealed a significant change in Delta power in Delta bout periods compared to non-Delta bout (p=9e-03, ges=0.176), which clearly showed the efficiency of the method to distinguish periods dominated by delta activity in S1 in the REM signals. We also observed a significant interaction between these two factors for the area S1 (i.e. cortical areas and Delta versus non-Delta bout discrimination, p=2.8e-2, ges=0.180; after Bonferroni correction, the interaction was significant for S1 (p=2.31e-4), but not for S2 (p=0.056) and for M1 (p=0.1)).

In figure 3 panels D and E, we found significant difference in Delta band power according to the areas (p=2.56e-5, ges = 0.444 and p=7.12e-4, ges=0.331) and the Delta/non-Delta bout discrimination performed on S2 or M1 respectively (p=7e-3, ges=0.184 and p=0.4e-2, ges=0.209). The interaction between the two factors was no longer significant for the areas S2 and M1 (i.e. cortical areas and Delta versus non-Delta bout discrimination, p=1.12e-1 and p=0.46). For the Delta versus non-delta bouts found on S2 (panel D), we observed significant differences in delta band power for the area S1 (p=0.038) and the area S2 (p=0.003). For Delta versus non-Delta bouts in the area M1 (panel E), we observed significant differences for the three areas (S1 p=9.16e-4, S2 p=0.03, M1 p=3.45e-4).

Altogether these results showed that our method was efficient to detect periods dominated by delta activity in the areas S1, S2 and M1. They also indicated differences in Delta power across areas in the distribution of Delta and non-Delta bouts, highlighting nuanced spatial variations of delta activities in the different cortical areas.

## 3.2. Computational models of REM sleep in mice

To investigate possible mechanisms to explain these local delta waves during REM sleep, we employed computational models. We used The Virtual Brain simulation environment (https://www.thevirtualbrain.org/tvb) to simulate the whole brain of the mouse (see Methods). The model was inspired from previous publications [35, 44] who had used mean-field models with conductance-based synaptic interactions to simulate human wakefulness and slow-wave sleep using the TVB platform. Like the human model, the TVB mouse model produced the asynchronous dynamics reminiscent of Wakefulness (Supplemental Fig S2 B, D b=5 (10.5281/zenodo.15063947)), and the synchronized slow waves reminiscent of SWS (Supplemental Fig S2 C, D b=60 (10.5281/zenodo.15063947)). In this model, the transition was driven
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>Figure 3: Temporal presence of Delta in REM.</figcaption></figure>

A) Figure representing the results of the delta bout analysis using a color code. Isolated delta cycles are represented in light green with horizontal black bars above them, while delta bouts are represented in bold green with a green overlay. Non-delta periods are indicated in red. Non-scored periods are represented in blue. Same LFP recordings as those presented in figure 2A. (B) Temporal overlap of Delta bouts between brain areas. Temporal overlap is calculated for each REM episode as the portion of time where Delta bouts in the origin and target areas coincide. Each violin plot shows the distribution of overlap percentages (y axis from 0, no overlap, to 1 (100%) overlap) across the 120 REM episodes. (C) Power spectrum of S1 Delta bouts and Non-Delta bouts. For each area (with the same color code as in figure 2), the power spectrum shown is the average power spectrum of this area, for the Delta bout periods of S1 in plain line, and for the non-Delta periods of S1 in dashed line. The thickest line highlights S1, the area from where the Delta / non-Delta periods are determined for this figure. (D) Same as C centred on S2. (E) Same as C centred on M1.

by modulating the level of adaptation in excitatory populations, which mimicked the action of neuromodulators such as acetylcholine, noradrenaline or serotonin, which strongly reduced adaptation in cortical neurons (reviewed in [40]).

In the wakefulness model (W-model), with low adaptation (b=5), the neuronal populations were constantly firing with higher firing frequencies for inhibitory cells as observed experimentally [50]. We also saw that the different mean-fields were not synchronized together and did not exhibit clear patterns, thus causing asynchronous dynamics. While in the SWS-model with high adaptation (b=60), the

neuronal populations alternated between periods of high frequency with synchronized firing (up to 50Hz) and periods of silence where the mean firing rate fell to 0 (Hz), inducing a shift to inverted synchronized slow waves in the LFP with Up and Down state dynamics. We also presented the average LFP and spectrum of two intermediate adaptation values for this model. At an intermediate low value (b=15nS) most nodes engaged in an oscillatory dynamics but it was less consistent and of smaller amplitude than the oscillations of the SWS-model and there is no real period of silence within the neuronal population. Then, with b at 30nS, a value that will
be used in our REM model, the simulation was already in a regime similar to our SWS-model but with shorter periods of silence it oscillates faster.

## 826 Simulating REM sleep

In order to simulate REM sleep activity, we hypothesized that the level of adaptation (or equivalently the action of neuromodulators) was inhomogeneous in the cortex. Figure 4 shows an example of the heterogeneous distribution of the adaptation parameter b. Motivated by our experimental analysis, we assumed higher values of b in areas such as S1 and S2, and low values in association cortex.

This type of distribution produced the spatial patterns of local delta waves in the simulated REM sleep. In figure 4B, the nodes with high adaptation displayed spectrums dominated by the power between 2 and 6 Hz, which could correspond to the model equivalent of the Delta band. It should be noted that there was a considerable contrast between this relatively broad power in the Delta band and the strong and sharp peak that the SWS-model produced (Supplemental Fig S2 (10.5281/zenodo.15063947)). For nodes with intermediate level of adaptation, although the PSD peaks at a faster frequency, about 6Hz, there was more power in the low frequencies (below 6Hz) compared to nodes with low adaptation.

In the simulated LFP example of the node corresponding to S1, with high adaptation (Fig 4C), we could see some slow fluctuations in the signal. This activity was not sustained along the whole simulation but few consecutive cycles were observed around the 7th second of simulation. Those slower oscillations were of higher amplitude than the portions of the signal that oscillated faster. This slow activity was the model equivalent of the Delta activity we described in the previous analysis. It was harder to find such well defined Delta activity in the simulated LFP for S2. However we could see that the signal of nodes with intermediate level of adaptation, such as the node corresponding to S2, oscillated on average slightly slower and in a slightly more synchronized way than nodes with low adaptation (Fig 4.D,E).

To better describe the presence of Delta activity according to the adaptation parameter in the model, we reproduced a Delta cycle and Delta bout analysis on a one minute long simulation (Supplemental Fig S4 (10.5281/zenodo.15063947)). It showed that the number of detected Delta cycles at each node increased with higher adaptation, and their shape is more consistently skewed toward a long rising time and a short decay time. Then, the decomposition of the signals in Delta and non-Delta bouts showed that with higher node adaptation the average number and duration of Delta bout increased, while we observed the inverse relation for the non-Delta bout (Supplemental table S2 (10.5281/zenodo.15063947)). Finally, the power spectrum of Delta bout versus non-Delta bout revealed a large difference in Delta power, which was more pronounced for nodes with high adaptation than it was for nodes with lower adaptation.

Our REM-model, with heterogeneous parametrization of

the adaptation presented some sporadic slow activity in the nodes with high adaptation that was close to the Delta oscillation observed experimentally. Here, areas with an intermediate level of adaptation oscillated slower than nodes with low adaptation of the reference model of wakefulness, and were prone to produce a few delta events.

In order to grasp the importance of the topological setup for the spike-triggered adaptation we proposed a control model: REM-shuffle, with spatially randomly distributed parameter b. The control model indicated that the observed patterns were not solely determined by the presence of nodes with varying adaptation values or the absolute level of adaptation. Instead, the essential factors were the spatial arrangement and connectivity between nodes for the described dynamics to occur. An example was shown in Supplemental Fig S3 (D,E,F) (10.5281/zenodo.15063947). As expected each REM-shuffle presented different global activity and the nodes with high b did not always present high slow frequency power (<5Hz). In the example we presented, nodes with high b peaked at faster frequencies, but some other REM-shuffle (data not shown) could produce dynamics more biased toward slower oscillations than our REM-model. The main point of this control was to show that the adaptation parameter alone was not sufficient to generate the dynamics we wanted to produce.

We explored alternative control models to test the robustness of the REM model. A "reversed model" that mirrored the REM model by featuring high adaptation in areas exhibiting asynchronous irregular activity and low adaptation in regions marked by Delta activity. With this configuration, nodes with low adaptation presented a signal with smaller amplitude and ratio of high frequency over low frequency power higher than nodes with high adaptation, but we observed very clear up and down in all nodes. Thus, the control model did not have the same dynamic as the REM model, since the difference between local and global dynamics was very slight in the control.

We also considered variations of our REM model, preserving the central high-adaptation and peripheral low-adaptation gradient but shifted the central focus to different nodes beyond the S1 barrel cortex. This resulted in the creation of 51 control models, with various patterns of activity. Compared to our REM model, we categorized the different results into three categories: those with stronger slow activity in nodes with low or average adaptation, those similar to our REM model, and those with less slow activity. These control models demonstrated that our modeling approach, which involved applying a gradient of adaptation, produced a variety of cortical dynamics. The REM model configuration was not the only one producing local slow activity. This suggested the importance of applying the high adaptation to the correct nodes in the model. For future work, it would be interesting to relate those different dynamics to the connectivity pattern of the nodes in the central focus of the adaptation gradient.
Stimulation in REM model 932

Finally, to investigate the possibility that Delta in pri-
mary sensory cortex impair information integration and thus
decrease responsiveness to sensory stimuli [51, 8], we used
our model to compare the responsiveness of different brain
states with respect to external inputs (Fig. 5; see also Sup-
plemental Fig S6 (10.5281/zenodo.15063947)).

Stimulus evoked response could be described as com-
posed of two main parts. The first part was a large amplitude
and rapid negative deflection that occurred in less than 50ms
after the stimulation and corresponded to the recruitment of
a large population of neurons. In figure 5A and Supplemen-
tal Fig S6 (10.5281/zenodo.15063947), we observe that this
initial deflection was slightly slower and of higher amplitude
in the SWS-model compared to the W-model. Then in a sec-
ond time the LFP returned to the baseline level. In the W-
model this return was fast and accompanied with a second
much smaller negative peak in the LFP. In the SWS-model
this second part was slower and was associated with a pos-
itive peak around 280ms, corresponding to a post-stimulus
suppression of the neural activity. Although the amplitude
of the first peak response was exaggerated in our model and
the response was a bit slower, the general kinetics were very
similar to what could be observed experimentally [52, 53].
In our REM model, the response to stimulation in areas with
high adaptation such as S1 was a bit different to what we saw
for the W or SWS models. The amplitude of the first peak
was between the two other model levels, although closer to
W-model. Then, with a slow hyperpolarized return to the
baseline, although without a clear positive peak, the second
phase of the response was closer to what we observed in the
SWS-model.

In figure 5.B we saw in our different model states how
the response to S1 stimulation propagated across the cor-
tex. We observed that the response in the non-stimulated
nodes was much smaller, when there was a clear response.
For S2 we noticed that the response followed at a smaller
scale the dynamics seen in S1 for the W-model and SWS-
model. However for the REM-model the second phase of
the response in S2 was closer to a slower version of that ex-
hibited by the W-model, with no clear post-stimulation hy-
perpolarization. For M1 and the other areas, the response
was not clear except in the SWS-model because of ampli-
tude issues. However in the SWS-model the response of M1
peaked later than for S2 and with a smaller amplitude. Such
observation could have been guessed for the two other states
and it followed the response propagation pattern observed in
the barrel cortex of awake mice [53]. We did not show stim-
ulation responses in our control REM-shuffle for clarity, but
when averaged over multiple simulations (therefore multiple
random parametrization) the results were very close to what
we observed in the W-model.

Finally we investigated whether our REM-model responded
to stimulation in the same way when nodes with different
adaptation values were stimulated. As we observed in the
figure 5.D, in homogeneous models such as the Wm or the

SWS-model, the response to stimulation in the different areas was extremely similar. On the other hand, in our REM-model we saw that each of our adaptation values elicited a slightly different response. The high adaptation values had the largest first peak and the slowest second phase with marked hyperpolarization. The nodes with intermediate adaptation such as S2 had a faster and more pronounced return to the base line in the secondary part of the response, but it was still relatively slow and there was no negative peak there. And then, the nodes with low adaptation responded to stimulus similarly to the W-model, with the second part of the response as a fast return to the baseline accompanied with a second negative peak. Those results exhibited the variability in stimulation responses one could expect when stimulating during REM, with some areas closely matching the W-model while others had a larger initial peak response and an hyperpolarized second phase as observed in SWS.

4. Discussion

In this paper, we have analyzed and modeled cortical
REM sleep activity in mice, and in particular the occurrence
of slow wave activities. The analysis showed that the REM
SWA was mostly constituted of delta oscillations in the 2-
4 Hz frequency range. This analysis also showed that REM
SWA differed from the SWA of SWS. On a spatial scale,
REM sleep expressed Delta waves only in some cortical ar-
eas but not in others, in contrast with SWS. On a tempo-
ral scale, again in contrast to SWS, REM Delta bouts alter-
nated with periods of asynchronous activity. We next in-
vestigated a computational model of the mouse brain build
assuming that the neuromodulatory drive, expressed with
the adaptation parameter, was heterogeneous in REM sleep.
The model further showed that assuming such a heteroge-
neous neuromodulatory drive, there was a coexistence be-
tween Delta activity and asynchronous activity in the same
brain. Below we relate these findings to previous work and
delineate possible perspectives to follow-up on this work.

Experimental analysis of REM sleep slow-waves

The occurrence of local slow-waves during REM has been
reported in humans and mice [6, 7, 8, 41? ], but our work
reveals considerable differences with the work of Funk et
al. First, REM slow-wave activity is significantly faster than
that of SWS and is mostly made of Delta activity (2-4Hz),
and not SO (<2Hz). It reinforces the idea that SWS SWA is
the association of SO and Delta waves [54, 55] and proposes
that, contrary to SWS, REM Delta (2-5Hz) is relatively dis-
sociated from slower activity (<2Hz). Second, it has been
claimed that REM slow-waves occur only in primary sen-
sory areas and therefore not in secondary or association ar-
eas. The clear occurrence of delta oscillations in S2 during
REM (Fig. 2) appears to contradict this statement. More-
over, while we find clear evidence of Delta activity in S1 and
M1, we do not observe it in the primary visual cortex (V1) or
and in the primary auditory cortex (Au1) (see Supplemental
Fig S8 (10.5281/zenodo.15063947)). Thus, it seems that the
<figure><img src="image_12.png" /><figcaption>Figure 4: Cortical model with heterogeneous B setting as model of REM.</figcaption></figure>

(A) Representation of the Mouse TVB network with nodes for each mean-field. The color of each node corresponds to the Spike-Adaptation B parameter. Nodes surrounded by a color circle (purple, dark blue, light blue, green and yellow) correspond to the areas studied from the biological dataset (see Supplemental Fig S7 (10.5281/zenodo.15063947)). (B) Mean power spectrum for each B parameter value in the model. Power spectrum is computed for each node LFP (with a time window of 4 sec and an overlap of 2 sec), and then averaged across nodes with the same b value. Averaged power spectrum are shown with the SEM. (C) Example 8 seconds of LFP signals from the REM model nodes highlighted in A. (D) Mean LFP signal and its deviation obtained at nodes with b=5nS. This corresponds to the adaptation value of the AI reference model. (E) Same as D for nodes with b=30nS. (F) Same as D for nodes with b=60nS. It corresponds to the adaptation value of the Up and Down reference model

results of Funk et al. in visual cortex cannot be generalized to all sensory cortical areas. In these areas theta oscillations are prominent during REM sleep. This topological organization is in line with the multi-site LFP recordings of [41]. These results suggest that LFP recordings in the neocortex result from a competitive combination of two oscillation generators, one close to the hippocampal theta and the other to the delta prominent in some cortical areas. A recent two-photon experiment by [56] even proposes two sub-stages for REM, a quiet and an active REM sleep, which might correspond to what we characterized as periods with prominent cortical delta activity in REM, or its absence. Such differences of Delta power in V1 compared to S1 is also observed in human in the REM episodes at the end of the night [6]. It is also likely that one of the multiple experimental differences (chronic vs acute recording, silicon probes vs tungsten

micro-electrodes, the detection methods...) could explain those apparently opposite results.

A major structure not yet studied is the thalamus which is known to play a key role in SWS oscillations [54]. It is indeed likely that the thalamus, connected to the different sensory modalities, dictates in an essential way the alternation between delta and desynchronized activities. Knowing about thalamic dynamics in REM would helps us to understand the origin of this activity, but would be also essential to produce a better, more complete, computational model. Recent work [57], using intra- and extracellular recordings, shows that the activity of relay neurons in thalamic nuclei connected to S1 and S2 areas changes during REM episodes in a coordinated manner with the LFPs of these cortical areas, as previously shown in wake [58] and NREM sleep [59]. Interestingly, in epileptic humans recorded by intracerebral
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>Figure 5: Model response to stimulation in the different states.</figcaption></figure>

(A) S1 barrel response when stimulated in the different states (Red corresponds to up and down (SWSm), Blue to Asynchronous Irregular (Wm), Green to REM-model (REMm) and olive is REM-shuffled (REMms)). (B) Cortical response to S1 stimulation in the different states. The average response to stimulation in S1 is observed in the different areas of interest across the model states (i: W-model, ii: SWS-model, iii: REM-model). (C) Representation of the model network. Surrounded in color the five nodes subjected to the stimulation experiment. (D) Response of cortical areas when stimulated across the model states. The color code corresponds to the areas highlighted in panel C. Each line is the average response observed in an area when it's stimulated. Each sub-panel corresponds to a different state of the model (i: W-model, ii: SWS-model, iii: REM-model)

EEG, the thalamic medial pulvinar nucleus is characterized by stable, continuous delta waves in REM, sometimes interrupted by short periods of rapid activity [60] and is shown to have a specific profile of thalamo-cortical connectivity between the thalamus and insula during REM sleep [61].

One of the essential highlights of Funk et al study is the large neuronal silence observed in the superficial layers of the sensory cortex during REM SWA. This silence period would participate in the decreased responsiveness of REM by making the cortical areas less prone to integrate inputs [8], but we see that REM delta is not constant and alternates with periods of asynchronous activity which may reflect a cortical state of increased neuronal excitability. The hippocampal activity might be coordinated with the alternation of Delta bouts periods of different activities. Our results suggest that phasic theta often occurs at the beginning and/or end of Delta bouts and might thus participate in the switch from Delta/non-Delta bouts. It has been proposed that phasic theta may be a window of intense synchronization activating populations of neurons in the hippocampus and possibly in the neocortex during brief periods of REM [16]. Thus, phasic theta periods during REM sleep could also correspond to short windows of increase in neuronal activity in cortical sensory areas.

The presence of REM SWA in some cortical areas could also be related to a heterogeneous temporo-spatial distribution of neuromodulatory signaling in the neocortex or thalamus. In classical studies the unique neuromodulation of REM has been associated with elevated Acetylcholine (Ach) [62, 19]. But in [63], they showed that SWA occurs when cholinergic signaling in S1 is blocked during wakefulness. Therefore if we consider neuromodulation as global, slow and tonic, in REM, the high Ach should prevent SWA. Also, recent technical advances [64], suggest in contrast that neuromodulatory signaling is phasic, short duration and very local (synaptic). This conceptual shift helps to relate the overall high Ach with the Delta bouts we observe. Moreover, in a recent study, Nazari and colleagues have observed a varying innervation over the neocortex of the cholinergic neurotransmission in REM, and correlated it with the areas exhibiting slow activity [25]. These findings support our model constructed on the hypothesis that in REM the cholinergic signalling toward the somatosensory and motor cortex is low and of short duration compared to other brain regions.

Our findings suggest several directions for future research. A first study is to determine the neocortical-hippocampal coordination during REM sleep, particularly in relation to SWA. Another critical direction involves understanding that-
lamic contributions, as they represent an influential missing link in determining whether the hippocampus or thalamus drives the temporo-spatial distribution of cortical activity during REM. Since sensory activation threshold is higher during REM, especially for the sense of touch, one may wonder whether SWA in S1 and S2 act as a sensory filter. A related question is how thalamic transmission to sensory cortices differs from W. This would help determine whether cortical delta indeed serves to decrease responsiveness.

Another follow-up would be to understand the relation between the different transient events that occur during REM, which concern the whole body. REM sleep remains a global phenomenon affecting the whole individual. Multiple organs work together and many factors, such as natural whisking, respiratory or cardiac rhythm could potentially affect brain oscillations [2, 65, 66, 67, 68]. Also, it would be interesting to relate the occurrence of Delta activity we observe to the minute-scale fluctuation that was recently described in [69]

More experiments could be planned to better characterize the cortical activity in REM, which has received little attention since the transition from cat to rodent as a biological model and is often discarded as just 'similar to wake'. Such experimental insights would be particularly valuable and could help to guide the development of computational models of REM sleep cortical dynamics. While models of wake and NREM states have become increasingly sophisticated across multiple scales, REM models are still lacking. Such models would generate testable predictions to guide future experiments, driving progress in both domains. One such model is discussed below.

## Computational models of REM sleep slow waves

The occurrence of local SWA has been reported before in humans [6, 7] and mice [8, 41]. Based on those observations we propose a new model of the whole cortical dynamics of REM sleep. This model is designed to relate the Delta waves observed in our analysis as well as the few previous observations of the phenomenon, and to this regard it is the first of such models. We also provide a randomized control showing that the emergence of REM-sleep features is not only due to the presence of high adaptation in the cortex, but is an emergent property of its spatial distribution in the network. In the model, we explain the delta waves by assuming a moderate neuromodulatory drive (using an heterogeneous adaptation) in some areas such as S1, S2, which then produce SWA at around 4 Hz. But those Delta waves we observe in S1 are of small amplitude. We have only little Delta activity in M1, meaning that with such parametrization, the drive from the somatosensory cortex is not sufficient to produce large bouts of Delta waves there. The temporal pattern of the Delta activity in S1, with Delta bouts alternating with faster activity, seems to be relevant compared to the experimental data. This could be further investigated by examining the delta cycles more closely, focusing on the Delta activity of individual nodes not explored in our work and delving

into the specifics of the model (see supplemental material (10.5281/zenodo.15063947)).

1177

1178

We use our model to predict the response to sensory stimulation in the different states. While there are some issues with the responses we observe in the model (addressed later in discussion), the general dynamics of the response across the cortex is consistent with experimental data [52, 53, 70, 71]. In the REM-model, the responsiveness is modulated with the level of adaptation, evolving from similar to W-model for areas with low adaptation to somehow closer to SWS-model responses for the low to high adaptation levels respectively. In those areas, such as in S1, the peak response is smaller than in the SWS-model which might indicate that a smaller portion of the neuronal population is recruited for the response. But the slow hyperpolarized second phase of the response is very close to what is characteristic of the SWS-model second phase. It supports the idea that in REM, when we stimulate an area that can present Delta oscillation such as S1, we should observe large peak response followed by an hyperpolarization that disrupts the propagation of activity. This supports the idea that REM SWA leads to an impairment of information integration and thus decrease responsiveness to sensory stimuli [51, 8].

1179

1180

1181

1182

1183

1184

1185

1186

1187

1188

1189

1190

1191

1192

1193

1194

1195

1196

1197

1198

1199

Such a prediction can be tested experimentally by applying whisker stimulation during REM sleep while recording multiple areas. But it should be noted that the stimulus we use in the model is not as natural as a whisker stimulation or a noise presentation done experimentally, but is closer to TMS experiments such as those performed in humans [72, 73]. In the latest article, authors perfomed TMS on the rostral part of the right premotor cortex during REM sleep. They report a response composed of two peaks close to what they see in wakefulness, but the first peak having a larger amplitude and the second peak a slower amplitude [73]. In our data, M1 behaves very similarly in W-model and REM-model. However when compared to the mPFC or dCA1 we notice a slight tendency toward a slower dynamics in the second phase. We already reported that in our model M1 lack of delta activity compared to the biological data, but it seems likely that with some minor tuning for this node our model could align with such human results.

1200

1201

1202

1203

1204

1205

1206

1207

1208

1209

1210

1211

1212

1213

1214

1215

1216

1217

The model suffers from a number of limitations. First, it is a population model that cannot explain and detail the activity at the single-neuron level. Therefore such an approach cannot be easily related to experimental data at the cellular level but is more appropriate to relate to population-level variables, such as LFP, EEG or imaging. Second, we identify some differences between the model and the experimental analysis. The low Delta activity in S2 and its absence in M1 should be overcome with a larger area that encompasses M1 for the application of adaptation and some small tuning of the adaptation values used (using more different values seems an interesting approach from some testings, however it makes the model complex to examine). Third, we also point to some limits regarding the responses to stimulation. Indeed we have a very large response amplitude and

1218

1219

1220

1221

1222

1223

1224

1225

1226

1227

1228

1229

1230

1231
a slightly slower dynamics. While the first issue might be addressed with a refinement of the stimulus used, the second would require adjusting the parameter of the temporal integration of the model, which would in turn require a complete verification of the different parameters used in the model.

The REM model we proposed can be improved in many ways. The topological parametrization of adaptation is kept simple and some issues can be addressed with finer tuning. But one can also think about more complex parametrization of the model that would take into account other aspects of REM sleep such as the characteristic theta activity [49, 16] in the hippocampus or the very active prefrontal areas [74, 75]. Also, the connectome we use in this model remains relatively limited. We use it for multiple practical reasons (easi-er to test, faster to run), but it lacks multiple brain regions, some of which, like the thalamus, might be critical to under-stand cortical delta oscillations or stimulus-response. Larger connectomes are now available with more regions and finer parcellation, but correctly integrating new regions accord-ing to this newly characterised phenomenon will require a large amount of research and the perspectives of this work are multiple.

Our findings suggest several promising research direc-tions that could extend and validate the present modeling across different species and experimental paradigms. First, one could perform a similar approach based on the human brain. One of the problems to overcome is that local slow-waves may not appear in the EEG, and will ideally require intracranial recordings [60]. However, such recordings are usually rare in primary sensory areas, so one could use the model to simulate local slow-wave activity, use models to calculate the EEG and determine how such local slow-waves affect the EEG. The aim would be to find identifiable EEG signatures of local slow-waves in REM sleep.

Another possible outcome is to use intracranial record-ings in monkey to confirm the presence of slow-waves in some sensory areas, and test the response to stimuli. A TVB model of the macaque monkey brain is currently available based on the Cocomac database [76]. One could use this model constrained by intracranial recordings in monkey to further investigate the role of local slow-waves with respect to sensory or other stimuli in REM sleep.

A final, but not least follow-up is to test the hypothesis behind our REM model: that the occurrence of Delta waves in REM is related to an heterogeneous level of neuromodu-lation. Because techniques now exist to measure cholinergic release locally [77], this question could be addressed exper-imentally. Experiments of this kind have been conducted at the same time of our work [25], and the results of Nazari study mostly agree with our work. It is also possible to use such data to construct a model where the level of adaptation (or cholinergic drive) is included not only spatially but also temporally, constrained by such data, leading to more realis-tic models of state transitions in brain activity during REM sleep.

**Supplementary materials**

1288

Supplementary material can be found at https://doi.org/10.5281/zenodo.15063947

1289

1290

**Code Availability**

1291

The code used for the simulations presented in this paper is available under request and will be available online after publication.

1292

1293

1294

**Acknowledgments**

1295

We want to thank the researchers involved in the ANR PARADOX for their various inputs to the project.

1296

1297

**Grants**

1298

Research supported by the CNRS, the ANR (PARADOX grant) and the European Union (Human Brain Project, H2020-785907, H2020-945539).

1299

1300

1301

**Disclosures**

1302

No conflicts of interest, financial or otherwise, are de-clared by the author.

1303

1304

**References**

1305

[1] E Niedermeyer. The electrocerebellogram. *clinical EEG and neuro-science*, 35:112–115, 2004.

1306

1307

[2] Sylvain Crochet and Carl C.H. Petersen. Correlating whisker behav-ior with membrane potential in barrel cortex of awake mice. *Nature Neuroscience*, 9:608–610, 5 2006.

1308

1309

1310

[3] Vladyslav V. Vyazovskiy, Umberto Olcese, Erin C. Hanlon, Yuval Nir, Chiara Cirelli, and Giulio Tononi. Local sleep in awake rats. *Nature*, 472:443–447, 2011.

1311

1312

1313

1314

1315

1316

[4] Edward Zagha, Amanda E. Casale, Robert N.S. Sachdev, Matthew J. McGinley, and David A. McCormick. Motor cortex feedback in-fluences sensory processing by modulating network state. *Neuron*, 79:567–578, 2013.

1317

1318

[5] Laura M.J. Fernandez, Jean Christophe Comte, Pierre Le Merre, Jian Sheng Lin, Paul A. Salin, and Sylvain Crochet. Highly dynamic spatiotemporal organization of low-frequency activities during behav-ioral states in the mouse cerebral cortex. *Cerebral cortex* (New York, N.Y.: 1991), 27:5444–5462, 2016.

1319

1320

1321

1322

[6] Benjamin Baird, Anna Castelnovo, Brady A. Riedner, Antoine Lutz, Fabio Ferrarelli, Melanie Boly, Richard J. Davidson, and Giulio Tononi. Human rapid eye movement sleep shows local increases in low-frequency oscillations and global decreases in high-frequency os-cillations compared to resting wakefulness. *eNeuro*, 5, 7 2018.

1323

1324

1325

1326

1327

[7] Giulio Bernardi, Monica Betta, Emiliano Ricciardi, Pietro Pietrini, Giulio Tononi, and Francesca Siclari. Regional delta waves in human rapid eye movement sleep. *Journal of Neuroscience*, 39:2686–2697, 2019.

1328

1329

1330

1331

[8] Chadd M. Funk, Sakiko Honjoh, Alexander V. Rodriguez, Chiara Cirelli, and Giulio Tononi. Local slow waves in superficial layers of primary cortical areas during rem sleep. *Current Biology*, 26:396–403, 2016.

1332

1333

[9] Igor Timofeev and Sylvain Chauvette. *Neuronal Activity During the Sleep-Wake Cycle*, volume 30. Elsevier B.V., 1 edition, 2019.

1334

1335
[10] Maria V. Sanchez-Vives and David A. McCormick. Cellular and network mechanisms of rhythmic recurrent activity in neocortex. *nature neuroscience*, 3, 2000.

[11] I. Timofeev, F. Grenier, M. Bazhenov, T.J. Sejnowski, and M. Steriade. Origin of slow cortical oscillations in deafferented cortical slabs. *Cerebral cortex*, 10:1185-1199, 1 2000.

[12] Mircea Steriade and Igor Timofeev. Neuronal plasticity in thalamocortical networks during sleep and waking oscillations. *Neuron*, 37:563-576, 2003.

[13] Maria V. Sanchez-Vives. Origin and dynamics of cortical slow oscillations. *Current Opinion in Physiology*, 15:217-223, 2020.

[14] Alexander A Borbély, Fritz Baumann, Daniel Brandeis, Inge Strauch, and Dietrich Lehmann. Sleep deprivation: Effect on sleep stages and eeg power density in man. *Electroencephalography and Clinical Neurophysiology*, 51(5):483-493, 1981.

[15] Kenji Mizuseki, Sebastien Royer, Kamran Diba, and György Buzsáki. Activity dynamics and behavioral correlates of ca3 and ca1 hippocampal pyramidal neurons. *Hippocampus*, 22:1659-1680, 8 2012.

[16] SM Montgomery, A Sirota, and G Buszaki. Theta and gamma coordination of hippocampal networks during waking and rem sleep. *The Journal of ...*, 28:6731-6741, 2008.

[17] Jurij Brankačk, Claudia Scheffzük, Valeriy I. Kukushka, Alexei L. Vyssotski, Adriano B.L. Tort, and Andreas Draguhn. Distinct features of fast oscillations in phasic and tonic rapid eye movement sleep. *Journal of Sleep Research*, 21:630-633, 2012.

[18] Robson Scheffer-Teixeira, Hindael Belchior, Fábio V. Caixeta, Bryan C. Souza, Sidarta Ribeiro, and Adriano B L Tort. Theta phase modulates multiple layer-specific oscillations in the ca1 region. *Cerebral Cortex*, 22:2404-2414, 2012.

[19] Seung hee Lee and Yang Dan. Review neuromodulation of brain states. *Neuron*, 76:209-222, 2012.

[20] Akihiro Karashima, Mitsuyuki Nakao, Norihiro Katayama, and Kazuki Honda. Instantaneous acceleration and amplification of hippocampal theta wave coincident with phasic pontine activities during rem sleep. *Brain Research*, 1051:50-56, 7 2005.

[21] J Lerma and E Garcia-Austt. Hippocampal theta rhythm during paradoxical sleep. effects of afferent stimuli and phase relationships with phasic events. *Electroencephalogr Clin Neurophysiol*, 60:46-54, 1985.

[22] Terry E Robinson, R C Kramis, and C H Vanderwolf. Two types of cerebral activation during active sleep: relations to behavior. *Brain Research*, 124:544-549, 1977.

[23] Katharine Rowe, Ricardo Moreno, T Rern Lau, Umesha Wallooppil-lai, Bruce D Nearing, Bernat Kocsis, James Quattrochi, J Allan Hobson, Richard L Verrier, and Bernat Koc-sis. Heart rate surges during rem sleep are associated with theta rhythm and pgo activity in cats. *American Physiological Society*, 277:R843-9, 1999.

[24] H Sei and Y Morita. Acceleration of eeg theta wave precedes the phasic surge of arterial pressure during rem sleep in the rat. *Neuroreport*, 7:3059-3062, 1996.

[25] Mojtaba Nazari, Javad Karimi Abadchi, Milad Naghizadeh, Edgar J. Bermudez-Contreras, Bruce L. McNaughton, Masami Tatsuno, and Majid H. Mohajerani. Regional variation in cholinergic terminal activity determines the non-uniform occurrence of cortical slow waves during rem sleep in mice. *Cell Reports*, 42, 5 2023.

[26] R.W. McCarley. Mechanisms and models of rem sleep control. *Arch Ital Biol*, 142(4):429-67, 7 2004.

[27] Charlotte Héricé, Amisha A. Patel, and Shuzo Sakata. Circuit mechanisms and computational models of rem sleep. *Neuroscience Research*, 140:77-92, 2019.

[28] Maxim Bazhenov, Igor Timofeev, Mircea Steriade, and Terrence J Sejnowski. Model of thalamocortical slow-wave sleep oscillations and transitions to activated states. *Journal of Neuroscience*, 22:8691-8704, 2002.

[29] Alain Destexhe, Diego Contreras, and Mircea Steriade. Spatiotemporal analysis of local field potentials and unit discharges in cat cerebral cortex during natural wake and sleep states. *Journal of Neuroscience*, 19:4595-4608, 1999.

[30] Albert Compte, Maria V. Sanchez-Vives, David A. McCormick, and Xiao Jing Wang. Cellular and network mechanisms of slow oscillatory activity (<1 hz) and wave propagations in a cortical network model. *Journal of Neurophysiology*, 89:2707-2725, 5 2003.

[31] Sean Hill and Giulio Tononi. Modeling sleep and wakefulness in the thalamocortical system. *Journal of Neurophysiology*, 93:1671-1698, 3 2005.

[32] Henry Markram, Maria Toledo-Rodriguez, Yun Wang, Anirudh Gupta, Gilad Silberberg, and Caizhi Wu. Interneurons of the cortical inhibitory system. *Nature Reviews Neuroscience*, 5:793-807, 10 2004.

[33] Joachim Lübke and Dirk Feldmeyer. Excitatory signal flow and connectivity in a cortical column: Focus on barrel cortex. *Brain Structure and Function*, 212:3-17, 7 2007.

[34] Gustavo Deco, Viktor K. Jirsa, Peter A. Robinson, Michael Breakpear, and Karl Friston. The dynamic brain: From spiking neurons to neural masses and cortical fields. *PLoS Computational Biology*, 4, 8 2008.

[35] Jennifer S Goldman, Lionel Kusch, Bahar Hazal Yalcinkaya, Damien Depannemaecker, Trang-Anh E Nghiem, Viktor Jirsa, and Alain Destexhe. Brain-scale emergence of slow-wave synchrony and highly responsive asynchronous states based on biologically realistic population models simulated in the virtual brain. *bioRxiv*, page 2020.12.28.424574, 2020.

[36] Francesca Melozzi, Marmaduke M. Woodman, Viktor K. Jirsa, and Christophe Bernard. The virtual mouse brain: A computational neuroinformatics platform to study whole mouse brain dynamics. *eNeuro*, 4:1-14, 2017.

[37] Matteo di Volo, Alberto Romagnoni, Cristiano Capone, and Alain Destexhe. Biologically Realistic Mean-Field Models of Conductance-Based Networks of Spiking Neurons with Adaptation. *Neural Comput*, 31(4):653-680, 4 2019.

[38] Alain Destexhe. Self-sustained asynchronous irregular states and updown states in thalamic, cortical and thalamocortical networks of nonlinear integrate-and-fire neurons. *Journal of Computational Neuroscience*, 27:493-506, 2009.

[39] Moritz Augustin, Josef Ladenbauer, and Klaus Obermayer. How adaptation shapes spike rate oscillations in recurrent neuronal networks. *Frontiers in Computational Neuroscience*, 7, 2 2013.

[40] David A McCormick. Neurotransmitter actions in the thalamus and cerebral cortex and their role in neuromodulation of thalamocortical activity. *Progress in Neurobiology*, 39:337, 1992.

[41] Sara Soltani, Sylvain Chauvette, Olga Bukhtiyarova, Jean Marc Lina, Jonathan Dubé, Josée Seigneur, Julie Carrier, and Igor Timofeev. Sleep-wake cycle in young and older mice. *Frontiers in Systems Neuroscience*, 13:1-14, 2019.

[42] Scott Cole and Bradley Voytek. Cycle-by-cycle analysis of neural oscillations. *Journal of neurophysiology*, 122:849-861, 2019.

[43] Joses Ho, Tayfun Tumkaya, Sameer Aryal, Hyungwon Choi, and Adam Claridge-Chang. Moving beyond p values: data analysis with estimation graphics. *Nature Methods*, 16:565-566, 7 2019.

[44] Jennifer S. Goldman, Lionel Kusch, David Aquilue, Bahar Hazal Yalcinkaya, Damien Depannemaecker, Kevin Ancourt, Trang-Anh E. Nghiem, Viktor Jirsa, and Alain Destexhe. A comprehensive neural simulation of slow-wave sleep and highly responsive wakefulness dynamics. *Frontiers in Computational Neuroscience*, 16, 2023.

[45] Giovanni Rabuffo, Jan Fousek, Christophe Bernard, and Viktor Jirsa. Neuronal cascades shape whole-brain functional dynamics at rest. *eNeuro*, 8, 9 2021.

[46] Sami El Boustani and Alain Destexhe. A master equation formalism for macroscopic modeling of asynchronous irregular activity states. *Neural computation*, 21:46-100, 1 2009.

[47] Bartosz Telenczuk, Maria Telenczuk, and Alain Destexhe. A kernel-based method to calculate local field potentials from networks of spiking neurons. *Journal of Neuroscience Methods*, 344:108871, 2020.

[48] Federico Tesler, Núria Tort-Colet, Damien Depannemaecker, Malory Carlu, and Alain Destexhe. Mean-field based framework for forward modeling of lfp and meg signals. *Frontiers in Computational
1474 Neuroscience, in press, 2022.

1475 [49] G Buzsáki, D L Buhl, K D Harris, J Csicsvari, B Czéh, Czéh, and A Morozov. Hippocampal network patterns of activity in the mouse. Neuroscience, 116:201–211, 1 2003.

1476 [50] Nima Dehghani, Adrien Peyrache, Bartosz Telenczuk, Michel Le Van Quyen, Eric Halgren, Sydney S. Cash, Nicholas G. Hatsopoulos, and Alain Destexhe. Dynamic balance of excitation and inhibition in human and monkey neocortex. Scientific Reports, 6, 3 2016.

1477 [51] Giulio Tononi and Marcello Massimini. Why does consciousness fade in early sleep? Annals of the New York Academy of Sciences, 1129:330–334, 2008.

1478 [52] Yuval Nir, Vladyslav V. Vyazovskiy, Chiara Cirelli, Matthew I. Banks, and Giulio Tononi. Auditory responses and stimulus-specific adaptation in rat auditory cortex are preserved across nrem and rem sleep. Cerebral Cortex, 25:1362–1378, 2015.

1479 [53] Pierre Le Merre, Vahid Esmaeili, Eloïse Charrière, Katia Galan, Paul A. Salin, Carl C.H. Petersen, and Sylvain Crochet. Reward-based learning drives rapid sensory signals in medial prefrontal cortex and dorsal hippocampus necessary for goal-directed behavior. Neuron, 97:83–91.e5, 1 2018.

1480 [54] F. Amzica and M. Steriade. Electrophysiological correlates of sleep delta waves. Electroencephalography and Clinical Neurophysiology, 107:69–83, 1998.

1481 [55] Jaekyung Kim, Tanuj Gulati, and Karunesh Ganguly. Competing roles of slow oscillations and delta waves in memory consolidation versus forgetting. Cell, 179:514–526.e13, 2019.

1482 [56] Yufan Dong, Jiaqi Li, Min Zhou, Yihui Du, and Danqian Liu. Cortical regulation of two-stage rapid eye movement sleep. Nature Neuroscience, 11 2022.

1483 [57] Flore Boscher, Katlyn Jumel, Tereza Dvořáková, Luc J. Gentet, and Nadia Urbain. Thalamocortical dynamics during rapid eye movement sleep in the mouse somatosensory pathway. Journal of Neuroscience, 44(25), 2024.

1484 [58] Nadia Urbain, Paul A. Salin, Paul Antoine Libourel, Jean Christophe Comte, Luc J. Gentet, and Carl C.H. Petersen. Whisking-related changes in neuronal firing and membrane potential dynamics in the somatosensory thalamus of awake mice. Cell Reports, 13:647–656, 10 2015.

1485 [59] Nadia Urbain, Nicolas Fourcaud-Trocmé, Samuel Laheux, Paul A. Salin, and Luc J. Gentet. Brain-state-dependent modulation of neuronal firing and membrane potential dynamics in the somatosensory thalamus during natural sleep. Cell Reports, 26:1443–1457.e5, 2019.

1486 [60] Michel Magnin, Hélène Bastuji, Luis Garcia-Larrea, and François Mauguière. Human thalamic medial pulvinar nucleus is not activated during paradoxical sleep. Cerebral Cortex, 14:858–862, 2004.

1487 [61] Isabelle Lambert, Nicolas Roehri, Julie Fayerstein, Bernard Giusiano, Bruno Colombet, Christian George Bénar, and Fabrice Bartolomei. Cortico-cortical and thalamo-cortical connectivity during non-rem and rem sleep: Insights from intracranial recordings in humans. Clinical Neurophysiology, 143:84–94, 11 2022.

1488 [62] Herbert H Jasper and Jacques Tessier. Acetylcholine liberation from cerebral cortex during paradoxical (rem). American Association for the Advancement of Science, 172:601–602, 1971.

1489 [63] Emmanuel Eggermann, Yves Kremer, Sylvain Crochet, and Carl C.H. Petersen. Cholinergic signals in mouse barrel cortex during active whisker sensing. Cell Reports, 9:1654–1661, 2014.

1490 [64] Martin Sarter and Cindy Lustig. Forebrain cholinergic signaling: Wired and phasic, not tonic, and causing behavior. Journal of Neuroscience, 40:712–719, 2020.

1491 [65] Detlef Heck, Samuel McAfee, Yu Liu, Abbas Babajani-Feremi, Roozbeh Rezaie, Walter Freeman, James Wheless, Andrew Papanicolaou, Miklós Ruszinkó, and Robert Kozma. Cortical rhythms are modulated by respiration. NEURON, pages 1–30, 2016.

1492 [66] Antoine Bergel, Thomas Deffieux, Charlie Demené, Mickaël Tanter, and Ivan Cohen. Local hippocampal fast gamma rhythms precede brain-wide hyperemic patterns during spontaneous rodent rem sleep. Nature Communications, 9, 2018.

1493 [67] Adriano B.L. Tort, Simon Ponsel, Jakob Jessberger, Yevgenij

Yanovsky, Jurij Brankačk, and Andreas Draguhn. Parallel detection of theta and respiration-coupled oscillations throughout the mouse brain. Scientific Reports, 8, 12 2018.

1494 [68] Adriano B.L. Tort, Jurij Brankačk, and Andreas Draguhn. Respiration-entrained brain rhythms are global but often overlooked. Trends in Neurosciences, 41:186–197, 2018.

1495 [69] Lezio S. Bueno-Junior, Maxwell S. Ruckstuhl, Miranda M. Lim, and Brendon O. Watson. The temporal structure of rem sleep shows minute-scale fluctuations across brain and body in mice and humans. Proceedings of the National Academy of Sciences of the United States of America, 120, 2023.

1496 [70] Manon Rampon, Julien Carponcy, Mégane Missaire, Romain Bouet, Regis Parmentier, Jean-Christophe Comte, Gael Malleret, and Paul A. Salin. Synapse-specific modulation of synaptic responses by brain states in hippocampal pathways. The Journal of Neuroscience, 43:1191–1210, 2 2023.

1497 [71] Paul Marchal, Paul A. Salin, Mégane Missaire, Manon Rampon, Julien Carponcy, Régis Parmentier, Gina Poe, Gaël Malleret, and Jean-Christophe Comte. Dynamics of evoked responses in hippocampal pathways are encoded by the duration of vigilance states. Nature Communications, 16:2973, 3 2025.

1498 [72] Marcello Massimini, Fabio Ferrarelli, Steve K Esser, Brady A Riedner, Reto Huber, Michael Murphy, Michael J Peterson, and Giulio Tononi. Triggering sleep slow waves by transcranial magnetic stimulation. Proc Natl Acad Sci U S A, 104(20):8496–501, 2007.

1499 [73] M. Massimini, F. Ferrarelli, M. J. Murphy, R. Huber, B. A. Riedner, S. Casarotto, and G. Tononi. Cortical reactivity and effective connectivity during rem sleep in humans. Cognitive Neuroscience, 1:176–183, 2010.

1500 [74] Leslie Renouard, Francesca Billwiller, Keiko Ogawa, Olivier Clément, Nutabi Camargo, Mouaad Abdelkarim, Nadine Gay, Céline Scoté-Blachon, Rouguy Touré, Paul Antoine Libourel, Pascal Ravassard, Denise Salvert, Christelle Peyron, Bruno Claustrat, Lucienne Léger, Paul Salin, Gael Malleret, Patrice Fort, and Pierre Hervé Luppi. The supramammillary nucleus and the claustrum activate the cortex during rem sleep. Science Advances, 1, 4 2015.

1501 [75] Renato Maciel, Risa Yamazaki, Dianru Wang, Anna De Laet, Sébastien Cabrera, Claudio Agnorelli, Sébastien Arthaud, Paul Antoine Libourel, Patrice Fort, Hyunsook Lee, Claudio Queiroz, and Pierre Hervé Luppi. Is rem sleep a paradoxical state?: Different neurons are activated in the cingulate cortices and the claustrum during wakefulness and paradoxical sleep hypersomnia. Biochemical Pharmacology, 191, 9 2021.

1502 [76] Kelly Shen, Gleb Bezgin, Michael Schirner, Petra Ritter, Stefan Evlerling, and Anthony R. McIntosh. A macaque connectome for large-scale network simulations in the virtual brain. Scientific Data, 6, 12 2019.

1503 [77] Ricardo M. Santos, João Laranjinha, Rui M. Barbosa, and Anton Sirota. Simultaneous measurement of cholinergic tone and neuronal network dynamics in vivo in the rat brain using a novel choline oxidase based electrochemical biosensor. Biosensors and Bioelectronics, 69:83–94, 2015.