key: cord-0075170-p4ejdth1 authors: Ioannides, Andreas A.; Orphanides, Gregoris A.; Liu, Lichan title: Rhythmicity in heart rate and its surges usher a special period of sleep, a likely home for PGO waves date: 2022-02-15 journal: Curr Res Physiol DOI: 10.1016/j.crphys.2022.02.003 sha: f9f0f2ae6f571be3d42a79f673674d1137694fc3 doc_id: 75170 cord_uid: p4ejdth1 High amplitude electroencephalogram (EEG) events, like unitary K-complex (KC), are used to partition sleep into stages and hence define the hypnogram, a key instrument of sleep medicine. Throughout sleep the heart rate (HR) changes, often as a steady HR increase leading to a peak, what is known as a heart rate surge (HRS). The hypnogram is often unavailable when most needed, when sleep is disturbed and the graphoelements lose their identity. The hypnogram is also difficult to define during normal sleep, particularly at the start of sleep and the periods that precede and follow rapid eye movement (REM) sleep. Here, we use objective quantitative criteria that group together periods that cannot be assigned to a conventional sleep stage into what we call REM0 periods, with the presence of a HRS one of their defining properties. Extended REM0 periods are characterized by highly regular sequences of HRS that generate an infra-low oscillation around 0.05 Hz. During these regular sequence of HRS, and just before each HRS event, we find avalanches of high amplitude events for each one of the mass electrophysiological signals, i.e. related to eye movement, the motor system and the general neural activity. The most prominent features of long REM0 periods are sequences of three to five KCs which we label multiple K-complexes (KCm). Regarding HRS, a clear dissociation is demonstrated between the presence or absence of high gamma band spectral power (55–95 Hz) of the two types of KCm events: KCm events with strong high frequencies (KCmWSHF) cluster just before the peak of HRS, while KCm between HRS show no increase in high gamma band (KCmNOHF). Tomographic estimates of activity from magnetoencephalography (MEG) in pre-KC periods (single and multiple) showed common increases in the cholinergic Nucleus Basalis of Meynert in the alpha band. The direct contrast of KCmWSHF with KCmNOHF showed increases in all subjects in the high sigma band in the base of the pons and in three subjects in both the delta and high gamma bands in the medial Pontine Reticular Formation (mPRF), the putative Long Lead Initial pulse (LLIP) for Ponto-Geniculo-Occipital (PGO) waves. The hypnogram is a most valuable summary of a night's human sleep. In its original form it is a succession of 30-s periods, with each period assigned to one of seven labels (Rechtschaffen and Kales, 1968) : undefined/movement, awake, rapid eye movement (REM), and four non-REM (NREM) stages that separate into light sleep (NREM1 and NREM2) and deep or slow wave (SW) sleep (NREM3 and NREM4). The definitions of sleep stages and the criteria for their definitions have been continuously updated. However, changes have been minimal over the half century of sleep staging. The main differences between the latest rules (Silber et al., 2007) compared to the original ones (Rechtschaffen and Kales, 1968) are the removal of the movement sleep stage and the amalgamation of the two sleep stages SW sleep (old NREM3 and NREM4) into the new NREM3 sleep stage. The new classification achieves the primary objective of enforcing more uniformity in sleep scoring through visual inspection by different experts. While fully acknowledging the tremendous impact sleep staging had in advancing sleep medicine (Shepard et al., 2005) , we will argue that the improvement made recently are not sufficient. These changes make the hypnogram less informative because the new definitions impose boundaries between sleep stages which cannot be related to fundamental changes in the brain. Also, the established sleep scoring scheme does not take into account abrupt changes in physiology, specifically the activity of the heart that can be objectively and easily captured with today's technology. The assignment of segments into sleep stages is based on the identification of graphoelements on the electroencephalogram (EEG) that are the hallmarks of each sleep stage and of the eyes closed waking (ECW) before sleep. These distinctive features are either highly rhythmic or large amplitude events. The highly rhythmic events include the prominent alpha activity in posterior areas during ECW that drops drastically as we enter NREM1 and the spindles of NREM2. The high amplitude events include the relatively modest vertex sharp waves (VSW) during NREM1 and the saw-tooth waves (STW) of REM and the prominent K-complexes (KC) during NREM2. Here, we grouped high amplitude events across sleep stages into the family of Global-Local Excitations Across the Brain during Sleep (GLEABS) and distinguished single unitary GLEABS (labelled as VSW1 and KC1) from multiple GLEABS appearing as sequences of such events (labelled as VSWm and KCm) . In addition to the EEG graphoelements, features of three other electrophysiological measurements contribute to sleep scoring: the electromyogram (EMG), the electrooculogram (EOG) and the electrocardiogram (ECG). The EMG monitors the muscle tone, and it is usually recorded using bipolar submental electrodes. A drop in muscle tone (EMG amplitude) characterises the transition from ECW to NREM1 and a further drop the transition from NREM sleep to REM. The EOG uses the differential signal of two electrodes for monitoring horizontal (electrodes on the side of the lateral Canthus of each eye) and vertical (electrodes on one eye one immediately above the eye and the other in the infraorbital ridge) eye movements. The onset of NREM1 is usually heralded by slow rolling eye movements on the EOG while fast eye movements, or REM saccades (REMS) are often mentioned as one of the hallmarks of REM. A horizontal REMS appears as a pair of prominent deflections of opposite polarity in the two lateral EOGs, while a vertical REMS shows as a prominent deflection on the bipolar EOG of one eye. A diagonal REMS has contributions from both horizontal and vertical EMs. The ECG uses a combination of electrode pairs (leads) to record the spread of the electrical activity through the body (that acts as a passive volume conductor) generated by the cardiac muscle depolarization. The ECG produces two key outcomes. The first is the sinus rhythm describing each cycle and dominated by the QRS complex. The second is the changes in the repetition of each cycle, usually expressed in terms of the time duration between successive R peaks (the R-R interval) or its inverse, the heart rate (HR). Sleep staging in sleep laboratories uses visual recognition of the rhythmic events and GLEABS to define sleep stages. EMG and EOG provide additional information particularly for the transitions from ECW to NREM1 and from NREM sleep stages to REM. Both measures of cardiac activity show changes in awake state and sleep but in only few cases these changes are well matched to individual sleep stages (Herzig et al., 2018; Kontos et al., 2020) and show some specific changes in only a few cases (Burgess et al., 1999) . The ECG has nevertheless been known to relate to periodic changes of arousal with shorter duration than sleep stages: the cyclic alternating patterns (CAP) for NREM sleep (Terzano et al., 1985) . In general, the ECG has consistently been the least used component of the standard polysomnography measurements, at least for sleep staging. In between periods of high activity, with the distinctive hallmarks of each sleep stages, there are "unremarkable" quiet periods, which inherit the sleep stage of their neighbours. In particular, if the pair of neighbours on either side belong to the same sleep stage, then the quiet segment sandwiched in between them is assigned the same sleep stage label. Also, a quiet segment showing no evidence (hallmark) of another sleep stage or excessive noise inherits the preceding label, if such a label is assigned to the previous 30-s period. This "by default" assignment of sleep stage label has survived throughout as a key element of sleep staging. We further identify within the quiet periods, the periods well away of the hallmarks of each sleep stage and any other high amplitude events as "core" periods and refer to them as NREMnc with the n is a number for the nth NREM sleep stage and REMc for the core period of REM. For example, a quiet period sandwiched between two NREM2 periods is assigned to NREM2, and hence its quiet periods makes up what we define an NREM2c. Tomographic analysis of core states showed distinct changes in spectral powers at widespread brain areas with distinct variations at specific cortical sites (Ioannides et al., 2009) . The most principled change was the increase in gamma band activity from awake state through NREMnc periods, culminating to its highest value during REMc; to the best of our knowledge this remains the only consistent change encountered as sleep progresses in each cycle from light to deep sleep and finally REM. This increase in gamma band activity is specific: it is most prominent in two left paramedial foci of dorsal brain. The two sites are MSRC1 in the left paramedial dorsal frontal cortex and MSRC2 in the midline posterior parietal cortex. For reasons explained below, the two foci are interpreted as the midline self-representation core (MSRC) . A detailed study of the whole brain activity for NREM2 has revealed a systematic change from NREM2c to the periods just before and during spindles and KC1 events (Ioannides et al., 2017 (Ioannides et al., , 2019 , supporting specific interpretations of the roles of KC1 and spindles. A KC1 event has a Janus face (Halász, 2016) because it serves simultaneously arousal and sleep promoting roles. A spindle emerges as a very special event, arising only after extraordinary measures are taken to further dampen influences from the environment, specifically in the areas that have just increased their low frequency activity (increase of NREM2c relative to NREM1c) (Ioannides et al., 2017) . These results together with the continuing relevance of core periods for sleep staging imply that far from being "unremarkable", core periods of sleep represent the "foundational elements" of each sleep stage. They are analogous to the ground state of a physical system (e.g., an atom) from which transitions to different excited states take place. This analogy becomes highly instructive as the changes from the core state of NREM2 are mapped in the periods before spindles and KC1; the changes identified are exclusively increases for theta and higher frequencies. These increases are modest and specific to a small number of areas in the periods just before the hallmarks, with further and much stronger increases for both spindles and KC1 events (Ioannides et al., 2017 (Ioannides et al., , 2019 . The increases in the pre-spindle changes are consistent with a careful preparation for sensitive memory consolidation process, which in the light of many other studies points to changes in the neural representation of self and specifically in the two areas identified with increased gamma band activity during REMc. These two areas and their individual activity are likely to be the closest we can come to a neural representation of self and hence labelled as MSRC1 and MSRC2 (Ioannides, 2018) . In recent years, a great effort was put on the development of automatic classification of sleep stages using biological signals that can be used at home in response to the move for self-monitoring of health and the obvious advantages of regular sleep monitoring at affordable cost for sleep quality monitoring (Scherz et al., 2017) . Most new methods use ECG and particularly heart rate variability (HRV) as a key component (Radha et al., 2019; Sun et al., 2020) . Claims have also been made for successful sleep staging using only heart rate (Mitsukura et al., 2020) . The emerging field of Network Physiology (NWP) focuses on the interaction between distinct networks representing complex physiological systems and particularly the cardiac, respiratory, muscular and central nervous systems. Each of these systems is complex with its own organization and processes and sophisticated regulatory mechanisms. In addition these systems continuously interact to maintain normal physiological function . Our findings should be viewed as part of NWP with good prospects of generalizing the analysis using multi-layer graph theory techniques. In this work, we use the conventional sleep stages as a broad reference to study the HRV changes across sleep and how the systematic changes in HRV and especially heart rate surges (HRS) relate to core states, GLEABS and the other electrophysiological measures used in sleep scoring. Introducing the HRV into the description of conventional sleep stages separates periods of sleep with well-defined features that combine hallmarks of different sleep stages and high arousal. Such periods are often marked as undefined. They may also be assigned to a sleep stage if the hallmarks representing that sleep stage prevail in the standard 30-s periods used for sleep staging. We will refer to these periods of sleep as REM-zero (REM0). A first and rather simplified definition for REM0 is "periods characterized by increase in the variability of both HR and the other main electrophysiological signals EEG/MEG/EMG". We examine the role of REM0 within the existing sleep staging framework using a unique set of polysomnography data with simultaneous whole night MEG recordings that allows a description of REM0 in terms of properties of the corresponding electrophysiological signals that define it. The work of this paper is not an attempt to propose an alternative sleep staging scheme, this will be the subject of future research. We simply propose an improvement: we add to the conventional sleep staging an additional step, which allows well-defined quantitative measures to either reinforce the conventional scoring or to add the labels REM0 or REM, to some of the periods defined as ECW, NREM, REM and to periods that are undefined in the conventional sleep scoring. There are three immediate uses of REM0 as it is proposed here. The first is to define a distinct target period, in addition to REM for identifying Ponto Geniculate Occipital (PGO) waves in humans. This is not surprising since the most consistent correlate of PGO waves are the HRS (Rowe et al., 1999) . The second use is an additional characterization of sleep periods, of possible relevance to normal and pathological sleep along the CAP classification (Terzano et al., 1985) and the phasic and tonic division of REM (Simor et al., 2020) . Finally, the further steps beyond sleep staging can be used in the process of reconciling the sleep staging of two or more sleep experts with minimal human involvement. The original sleep classification scheme (Rechtschaffen and Kales, 1968) has revolutionise sleep research. It is remarkable that this system is still largely intact and serving as the main reference point for sleep medicine. Despite the obvious success of the original sleep classification there were always some problems that remain till today despite the changes introduced so far. One of the problems was that the hallmarks of individual sleep stages sometimes appear together in the same 30 s period to be scored, thus making a unique and unambiguous definition impossible. We now know, that many of the hallmarks used for sleep staging are characterized by "widespread intracranial activity with unexpected synchrony" when studied from either intracranial recordings (Frauscher et al., 2015 Latreille et al., 2020) or tomographic estimates of activity extracted with MFT from MEG data (Ioannides et al., 2017 (Ioannides et al., , 2019 . One may wonder why so much emphasis is placed on obtaining consistency in sleep staging classification which is determined by highly variable and wild excursions in the signal. One good reason is that the high amplitude events stand out in the EEG. The fact that this collective phenomenon arises each time from activity at different brain areas has only recently been realized (Frauscher et al., 2015 Ioannides et al., 2017 Ioannides et al., , 2019 Latreille et al., 2020) . As long as identifying these "hallmarks" leads to some classification of sleep (when there was none before) then agreeing on definitions was a great advance for prescribing treatment. However, in pushing for unanimity on the basis of the high amplitude and highly variable excursions of the signal may lead to adopting arbitrary criteria that may not be optimal in terms of a meaningful description of the underlying physiology. One reason for introducing changes in sleep staging was the concern about how REM can be distinguished from NREM2 periods and periods of large movements. We provide a characteristic extract from the summary of (Silber et al., 2007) , with R standing for REM; it reads: "R sleep commences when chin EMG tone falls, unless K complexes or spindles persist, in which case stage N2 persists until rapid eye movements develop. If chin EMG tone is low in stage N2 as well as REM sleep, the transition to Stage R occurs after the last K complex or spindle. If K complexes or sleep spindles are interspersed among what are otherwise epochs of unequivocal REM sleep (as may especially occur in the first REM period of the night), then stage R should be scored if rapid eye movements are present and stage N2 if eye movements are absent. If epochs of REM sleep are followed by epochs with low amplitude mixed frequency EEG and persistently low chin EMG tone, but no rapid eye movements, K complexes, or spindles, then they should continue to be scored as stage R until there is a transition to stage W, N2, or N3, where an arousal or major body movement is followed by slow eye movements, or an increase in chin EMG tone signifying a change to stage N1." One may be justified in concluding that sleep staging under the new guidelines for visual sleep classification may indeed be a good step forward but at the same time one may wonder whether sleep characterization should consider also a view from a new angle. After all, in the more than half a century since the introduction of sleep staging, there were amazing advances in mass electrophysiology (EEG and MEG) and neuroimaging and these have provided new insights about sleep that have not yet influenced sleep staging. In the work we will describe we offer such a new angle. In the new approach, we view the conventional hallmarks of each sleep stage as wild excursions of the system from equilibrium. This immediately poses a new question: what could this new equilibrium be and how can it be studied? We have already suggested that the equilibrium period for each sleep stage are the core states, as we have originally defined them (Ioannides et al., 2009 ). If core states indeed represent something like a ground state for each sleep stage, then one would expect that a detailed study of the core states of each sleep stage will provide a more intelligible description of sleep than the seemingly chaotic events that are now used to define each sleep stage. Our identification of consistent increase in the gamma band activity in medial dorsal brain areas across the core periods of the sleep stages provided the initial evidence supporting a foundational role for core states (Ioannides et al., 2009) . The follow up detailed study of changes from NREM2 core periods to the periods before and during spindles and KC1s (Ioannides et al., 2017 (Ioannides et al., , 2019 brought together recent suggestions for a role of sleep in memory consolidation, with emphasis on memory consolidation related to the neural representation of self (Ioannides, 2018) . We do not suggest that sleep staging based on the hallmarks of sleep as it is currently defined should be abandoned, at least not yet. In this work, we extend conventional sleep staging a little by paying attention to the critical role played by arousal influences as these are captured by changes in the EMG, EOG and particularly HRV. As already mentioned, sleep staging relies on identification of prominent features in the EEG that can be separated into highly rhythmic and high amplitude events. Examples of highly rhythmic events are bursts in the alpha (8-12 Hz) or sigma (11-17 Hz) bands. The alpha band characterises the awake state with eyes closed before sleep and its disappearance is taken as a sign of the start of the transition to sleep. The rhythmic events in the slightly higher sigma band are known as spindles and they are one of the two main hallmarks of NREM2. The prominent high amplitude NREM2 graphoelement is the KC. As already mentioned, we separate single KC (KC1) from multiple KC (KCm) which corresponds to 3-5 KCs appearing in a sequence one after the other. For reasons that will become apparent later we separate KCms into two types. The first type are KCm events with no strong high frequency content and therefore denoted as KCmNOHF. The other type of KCm has strong high frequency content and it is therefore denoted as KCmWSHF. We will show a strong association between KCmWSHF with HRS. It has been argued that the "KCs are isolated "down-states," a fundamental corticothalamic processing mode already characterized in animals" (Cash et al., 2009) . Our data support this claim but only for KC1 and KCmNOHF events. The MEG measurements were the first major experiment of the newly established laboratory for the human brain dynamics (LHBD-J) performed around the turn of the millennium at the Brain Science Institute (BSI), RIKEN near Tokyo, Japan. After the RIKEN ethical committee approved the protocol subjects were selected according to stringent selection criteria that included passing strict requirements for good sleep habits. As the final test of the entire procedure, a member of the team (AAI), was the subject resulting in the first whole night sleep MEG dataset; since this subject did not pass the selection criteria (for regular sleep habits) he was helped to sleep through sleep-deprivation. After the successful completion of this test and careful inspection of the resulting data, the protocol was internally approved for use with external subjects. For the main experiment, the protocol and the experiment were explained, emphasizing that subjects were free at any point to ask to be released from the recordings and stop the experiment. Each subject gave his informed consent and no subject used the opt-out option. Whole night, high quality data were obtained for three subjects. The test-data from the sleep night of the last test of the procedures (with the sleep deprived subject) were of good quality and they were also saved with the data of the three main subjects. All four subjects were right-handed, healthy males with ages 25, 30 and 31 for the three main subjects and 48 for the sleep deprived subject. Participants were all free of neuropsychiatric illness and medication. Despite the resurgence of sleep research and numerous new studies, whole night sleep MEG recordings are few. For a review of what is available today, in terms of comprehensive recording of electrophysiological data describing neuronal, muscular, eye movement and cardiac activity with good monitoring of head movement, see Supplementary Table 1 of (Brancaccio et al., 2020) . The first night began with training and placement of auxiliary channels. Horizontal and vertical eye movements were detected independently using four electrodes. Two electrodes were placed 1 cm from the lateral canthus of the left and right eye to measure horizontal EOG (EOG-H) and two electrodes were placed 1 cm above and below the left eye to measure the vertical EOG (EOG-V). After the eye movement data collection, the subject returned to the preparation room to sleep in a supine position on a replica of the MEG bed with their head in a replica of the MEG helmet. The eye movement data of the first (acclimatization) night have not been used in the current paper. The instructions to the subject, namely, to sleep keeping his head inside the helmet were repeated. A handset next to the bed allowed the subject immediate contact with a medically qualified member of LHBD-J in an adjacent room (demanded by the RIKEN ethics committee but never needed during the sleep experiments). The subject was awakened and debriefed after an estimated 8 h of sleep if he had not awoken spontaneously. Before the main recording (second night) and while the subject was prepared for the experiment, a noise measurement was taken with the system prepared as in the main experiment (for supine recordings). The data without any subject in the room were used as baseline measurements in some analysis. After the preparation, the subject was brought to the shielded room and told again that at any time he could communicate with the operators, ask to go to the toilet or detach the electrodes and stop the experiment if he had become too uncomfortable. When the subject was relaxed with the head in the helmet-like end of the Dewar, he closed his eyes and prepared to sleep. After final checks for safety and working order of devices were completed, the lights were switched off and the recording began (eyes closed waking, ECW condition). The final checks prevented the definition of sleep onset time with precision, so the sleep expert present at the start of the experiment made a subjective evaluation. The sleep onset time was judged a little longer than usual for one subject (S3) but normal for the other two. Throughout the recording night, at least two experimenters (usually more) were in the laboratory manning (quietly) the control console outside the shielded room or in the adjacent and preparation room; at least one experimenter was a native speaker of the language of the subject and one was a medical doctor. The subject was watched with an infrared camera and the MEG, EEG, and auxiliary channels were continuously monitored. For practical reasons, the entire night of sleep is divided into runs and each run is composed of a variable number of 3-min time periods of continuous recordings of MEG, EEG and auxiliary channels. For each whole night experiment the precise time of each sample of data is known by reference to the system clock and its time from the beginning of the segment it belongs. Each 3-min continuous recording is separated from the next by about twenty (20) seconds which correspond to the time of activation of the head localization coils and the integrity checks for their localization. If the integrity checks fail a new head localization coil is made and hence the separation between periods could be approximately twice as long or more. The start of the recording for each 3-min period is recorded with the data in absolute time (by the system clock). The duration of the run was limited by the capacity of the storage medium used to store the recordings. Unexpected, rare events like the need to deal with some technical problem with the equipment (e.g. to recalibrate electrodes), or a request by the subject for use of the toilet (made by two subjects) could shorten a run and they were used to upload the data from the temporary storage medium. After an interruption that required removing the subject from the shielded room, electrode and head coil connections were detached but replaced and retested on return; the long blank period in Fig. 1 , starting a little after 1.5 h is one such interruption. Both subjects that requested a break to use the toilet readily returned to sleep and slept well for the remainder of the night. Independently of the operational segmentation of sleep as defined in the previous paragraph, sleep is organized into 3 to 5 cycles for an entire night's sleep. A sleep cycle (SC) is a sequence of traversal through light and deep sleep culminating in REM; the nth SC is denoted by SCn, e.g., SC3 denotes the third sleep cycle. In the first part of the preparation, three head coils were attached to each subject's nasion and left and right pre-auricular points in order to record the subject's head position. EEG electrodes were also placed at C3 and C4, referenced A1, EOG vertical and horizontal, and chin EMG (Ioannides et al., 2004) . The EOG's were calibrated using controlled eye movements in the MEG room before the sleep recordings. The recordings were made with the whole-head 151-channel Omega magnetometer (CTF Systems Inc., Vancouver, Canada). The dewar was completely tilted into the horizontal position so the subjects could lie supine on the MEG bed, as is normal for sleep, with the head in the MEG helmet. The MEG data was recorded throughout the night, together with auxiliary EEG and marker channels on a 100 GB disk. The final preparation, before the MEG, five (5) head localization coils (HLC) were used. The location of all 5 HLCs was included in the head outline scan using the Polhemus magnetic stylus device. The head outline was fitted to the MRI of the subject, thus establishing a coordinate system based on each subject's true brain anatomy. As can be seen in Fig. 1 large head movements were present throughout the entire night of the sleep MEG experiment. Such movements are perfectly normal during sleep and allowing them makes sleep in the MEG less unnatural, but at the expense of losing the information necessary for co-registering the head relative to the sensors. To address this problem, the continuous recording was interrupted every 3 min for head-coil localization. For each subject, large head movements occurred intermittently and cannot be predicted, so it is critical to define the period between head localization measurements carefully. The choice of the 3-min period with no further restriction for the head inside the MEG helmet (except for the provision of a soft small pillow or cloth to support the head if the subject wanted one) was made after numerous tests. These tests helped us decide the best compromise between the conflicting requirements of continuous uninterrupted measurement (the choice of 3 min period) and making sure that even during extended periods with head movement a continuous segment with no large head movement could be captured. As can be seen in the results section, even during long periods dominated by large head movement artifacts, there were 3-min intervals of very low head movement and some of these were often REM0 periods (e.g., the 3 min of the first zoom in Fig. 3 ). Large head movements generate high EEG/MEG signal that can swamp the usual neuronal activity in the brain; it is almost impossible to use such data for identifying neural generators, partly because of the high noise and partly because the co-registration between sensors and the brain is lost. For this reason, all large signals above a pre-defined threshold are often eliminated as part of pre-processing. However, neural processes in sleep can also generate relatively large signals. Therefore, a blanket elimination of large signals is likely to remove information, possibly unique and valuable information. A strong artifact that does not involve large head movement will have minimal influence on the MFT solutions, because MFT focuses on activity from brain areas. In the opposite case, i.e. strong activity from a neural generator will be properly and easily localized with great accuracy. The Magnetoencephalographic (MEG) data were recorded using the whole-head Omega bio-magnetometer (CTF Systems Inc., Vancouver, Canada) inside a 3 × 4 × 2.4m shielded room (NKK, Kawasaki, Japan). The raw signals from the 151 primary MEG channels, the 28 reference sensors and from auxiliary channels needed for sleep scoring, eye movement and heart monitoring were recorded continuously on a 100 GB disk after low-pass filtering at 208 Hz and digitization at 625 Hz. The continuous recording mode was interrupted every 3 min for head localization. The same recording parameters were used for the eye movements during the first night but with head localization only at the beginning and end of each block. Noise runs were also taken before and after each full session (evening). In synchrony with the MEG, EEG electrodes C3 and C4, referenced A1, vertical and horizontal EOG and two (2) electrode ECG along with chin EMG were recorded. As already stated, every 3 min during the night experiment, the MEG recording was stopped and the HLCs were activated in turn. From these measurements the location of each coil relative to the sensor array was determined. By comparing the location of the coils at the start and end of each 3-min segment period the size of the segment's head movement was determined. The procedure allowed us to identify large head movements and this information was used for two selections. First, entire 3-min period were excluded from the tomographic analysis if it contained a movement above a set threshold (usually above 4 mm). Secondly, when summarizing signal properties (from the EEG and auxiliary channels) or gross properties across the MEG channels, e.g., Global Field Power (GFP) 3-min segments were excluded if the head movement was very large, typically 5 times or more the threshold for exclusion from tomographic analysis (~2 cm). The head movement thresholds are represented as a dashed black line on some of the presented Figures that follow. A preliminary "on the fly" scoring of the sleep stages was made during the experiment, partly for safety to monitor the subject and partly for on-line check on signal quality of MEG and EEG channels by at least one sleep expert. A more thorough sleep scoring was performed offline. First, each of our two sleep experts made an independent hypnogram using the same data sets (evaluating independently each one of the six 30 s segments of each continuous 3 min recording). The two independent scores were compared, a common one agreed and the final night hypnogram constructed. Sleep stages were noted according to the standardized manual (Rechtschaffen and Kales, 1968) , (W = wake state, MT = subject moving in bed, NREM1-NREM4 = first to fourth sleep stages, REM = REM sleep stage). Standard indicators identified REM: desynchronized EEG, minimization of muscle tone, and the emergence of REMS. We followed our original definition of core states (Ioannides et al., 2009) as, relatively quiescent periods that are similar to the shorter-lasting phase "B" periods of the cyclic alternating pattern in NREM (Ferri et al., 2006) and the tonic periods in REM sleep well clear of phasic eye movements (Simor et al., 2020; Wehrle et al., 2007) . A first-order estimate for what might be "core periods" was derived from a single EEG channel (C3 to A1). The core periods (ECWc, NREMnc and REMc) were then identified for each sleep stage from 3-min segments of data with head movement below the threshold for selection (4 mm). For each sleep stage, quiescent periods of 4 s duration were identified as candidate core periods. The core periods to be used for tomographic analysis were then selected from these candidate core periods with further criteria that they were clearly separated from large graphoelements of each sleep stage or any other prominent EEG features, obvious on any other MEG or EEG record. The GLEABS (and rhythmic) events were identified through visual inspection and were selected to be well away from each other and away from other events and noise. They were also selected to have a clear onset time so that a 4 s GLEABS could be extracted and separated into the first half representing a 2-s pre-GLEABS period and the second half representing the 2-s during-GLEABS period. The technical term for attempting to put similar objects into distinct groups is called clustering. Clustering is a rapidly advancing field receiving input from machine learning and its subdivision deep learning (Karim et al., 2021) . By definition, sleep staging is an attempt to put sleep periods into groups in such a way that each group (sleep stage) contains periods that are more similar to each other and less similar to those of periods belonging to a different group (sleep stage). It is therefore fair to conclude that sleep staging is a form of clustering. The classical sleep staging developed as a rule-based clustering technique based on the electrophysiological data available in the 1960s. The available data then could best discriminate large graphoelements and after further processing (e.g. filtering) periods with special properties (e. g. highly rhythmic activity). Looked at it from this point of view, the development of classical sleep staging in the 1960s and its reliance on the wide excursions of the signal during a rather chaotic period of EEG seems almost inevitable, given the technical capabilities of EEG at the time. Sleep scoring demands human expertise and it is time-consuming; it is therefore costly. Consequently, a lot of effort has been expended recently in developing (semi-)automatic methods of sleep staging based on modern clustering techniques. While methods of automated sleep scoring (using unsupervised learning) are beginning to make an impact (Chriskos et al., 2021) most methods use some form of supervised learning with human pre-labelled classifications as the teaching set. As we already stated, we are not concerned with establishing a new sleep staging procedure. Nevertheless, we will use clustering of core periods to test our claim that the core states of sleep are not just uninteresting periods. We will demonstrate that clustering of spectral density of regional brain activity are as capable of characterising individual sleep stages as the active periods relying on the characteristic hallmarks of each sleep stage. We will further demonstrate that the core states can go beyond classical sleep stages and define sub-clusters of REM. The resulting framework describes not only core periods of all sleep stages but also periods before and during hallmarks of and transitions between sleep stages. The key physiological signals were studied across widely different timescales, ranging from a few milliseconds, to tens of minutes and finally across the entire night of sleep. Some quantities can only be studied at low temporal resolutions. For example, in our analysis, the variability of cardiac activity uses as base unit the R-R interval which is of the order of 1 s, therefore it can only be studied with a resolution not finer than a few seconds, so enough base units are included to give a smooth continuous signal. Almost all other quantities can in principle be studied across the entire range of temporal scales we are interested here, but the measurements at some ranges are difficult, e.g., MEG measured at frequencies below 1 Hz cannot be measured directly because the shielded room does not eliminate ambient noise (which is much higher) at such low frequencies. As we will see, low frequency oscillations, below one Hz can be captured via the changes of GFP (or its variants) derived from the higher frequencies, which show marked variations at much lower frequencies (well below 1 Hz) reflecting cross-frequency coupling and nesting. We use this observation to demonstrate at a signal level analysis cross-system coupling and specifically the coupling of the cardiac rhythms with CNS oscillatory activity. For this signalbased analysis (Figs. 1-8 and 10) we processed each 3-min segments using default functions from the Fieldtrip toolbox (Oostenveld et al., 2011) developed at the Donders Institute for Brain, Cognition and behavior (http://fieldtriptoolbox.org). The MEG, EEG, ECG, EMG and EOG channels were band-stop filtered at 47.5-52.5Hz and 95-105Hz to remove power line interference. For the EEG and MEG signals additional filtering was made to obtain a wide band and signal using a high pass at 0.7Hz and a low pass at 95Hz filters. Signals were also studied in the other standard bands. For much of the work here we will present results for signals filtered in the wideband and in the high gamma (from 55 to 95 Hz) bands. For the EMG we cut the frequencies below 20 Hz; we used the signal from 20 to 95 Hz as the wideband signal and kept the high gamma band in the range (55-95 Hz). Depending on what aspects of eye movements it is to be emphasized different filters can be applied. The Global Field Power (GFP) was calculated separately for the (151) MEG and the (2) EEG channels and the single bipolar EMG channel. GFP is calculated using a moving window of N time-slices and M channels. N can range from 1 to 33 and M can be again just one channel or all channels (151 for MEG, 2 for EEG and 1 for EMG). The GFP is the ratio of the sum of the squares of the channel amplitude across all channels (M) and window time slices (N) to the total number of elements (N × M). The Instantaneous R-R Heart Rate (HR) was computed by first forming the difference between the two ECG leads. This difference shows clearly QRS signal with its prominent R-wave. The instantaneous HR can then be computed as the inverse of the periodicity of R, i.e., the inverse of the R-R interval. We will use the mean periodicity over three cycles in our computations, i.e., the inverse of the average of three successive R-R intervals, R − R . Thus, the heart rate becomes HR = 60 R− R and it is presented in the standard medical units of beats per minute. We will use this smoothed version of HR for all results reported hereafter, simply referred to as HR. Measures of variability can provide distinct measures by removing the mean and showing only how variable each signal is around that mean. The simplest measure of variability is the variance or its square root, the standard deviation. The computation of the variance is made within a window which defines the timescale over which the variability is computed. The variance and SD can be computed for either the signal of each channel or the GFP. The results reported in this paper were obtained using the standard MATLAB functions movvar() and movstd() for the variance and standard deviation respectively, using the formulae for a set of values S i shown below. Boxplots are used to display, in a quantitative way, specific aspects of HRV during phasic and Tonic periods of REM0 ( Fig. 2B ) and during GLEABS events and core periods in Fig. 4 (both parts A and B). Boxplots were generated using the standard MATLAB function 'boxplot'. In these boxplots the median is shown as a red line (inside the box) with the 25th and 75th percentile levels represented by the horizontal blue lines of the lower and upper boundaries of the box. Outliers are represented by crosses and a dotted line shows the range (excluding outliers). For the quantification of the HRV, many methods for HRV representations use time changes in consecutive R-R intervals also known as interbeat intervals (IBIs). In contrast, our analysis uses the measure defined earlier which we dubbed variance of the Heart Rate (V HR ). For each subject, the lower and upper limits are set for the acceptable range of HR, usually from about 40 to 120 bpm. We then remove periods outside this range and periods with unclear QRS complexes. This preprocessing eliminates noisy parts of the ECG and ectopic beats. The standard deviation of the HR, SD HR and its square, i.e. V HR are then computed for the remaining signal. Our method most closely resembles the SDNN method considered as the 'gold standard' measure for medical stratification of cardiac risk when recorded over 24 h period. SDNN measures the standard deviation of accepted R-R intervals hence renamed to Normal-Normal intervals (N-N intervals) in each period (Malik and Electrophysiology, 1996) . Our method provides two (2) key advantages: Firstly, the variance rather than the standard deviation amplifies differences rapidly and in a non-linear fashion within the critical range towards surges where the threshold must be set, making regions with high and low HRV easier to distinguish. Secondly, using bpm (or its square) instead of ms as a unit of measurement stays closer to the familiar clinical units of HR. HRV can also be computed using frequency-domain methods. This include calculating the power in each of the 3 main frequency bands Very Low Frequencies (VLF ∈(0 0.04] Hz), Low Frequencies (LF ∈[0.04 0.4] Hz) and High Frequencies (HF ∈[0.4 ∞] Hz) (Malik and Electrophysiology, 1996) . For frequency domain methods two (2) minute-periods should be used for LF (Shaffer et al., 2014) and at least one (1) minute-period for HF. In this paper we use the HRV for sleep classification using thirty, or even 10 s periods, therefore, frequency domain methods are not an option. Delving into the different mechanisms of HRV, as these may be probed using frequency domain methods for HRV computations, although interesting and worth doing, are beyond the scope of the work presented here. For the time-frequency analysis of Fig. 4B we used the generalized Morse wavelets of Olhede and Walden (Lilly and Olhede, 2012) , with parameters set to β = 27 and γ = 3. We used default functions from the jLab toolbox (http://www.jmlilly.net/jmlsoft.html) made for MATLAB to convolve Morse wavelet power and phase for the computations of the spectral powers within different bands. A Normalized percentage difference was employed to quantify the differences between the two types of events, e.g., the two types of KCm (Fig. 4B ). Time-frequency analysis was run on EEG C3 (C3-A1) channel for the whole trials (180 s) with low head movement. The KCm at HRS (KCmWSHF) and between HRS (KCmNOHF) were identified. For all KCmNOHF power over the time of the second oscillation was averaged over time for each frequency bin to obtain mean KCmNOHF power over frequency for all KCm events between HRS. The power of each KCmWSHF at the second oscillation was subtracted by the mean KCmNOHF power over each frequency bin and then normalized by dividing the result at each frequency bin by the maximum power recorded from both KCmWSHF and KCmNOHF. This analysis gives a normalized difference over frequency for each individual KCmWSHF compared to the mean power of all KCmNOHF. In order to avoid bias from electrode placement and different skin conductivity for each subject, the normalized difference was calculated for each subject separately before the boxplot was computed for each frequency and normalized before the percentage change across subjects was computed. More formally, let us use the labels x, y and z to stand for subject, epoch and frequency, so we label the power of a KCm events by KCmWSHF(x, y, z) and KCmNOHF(x, y, z) . We also denote the corresponding mean power across epochs by KCmWSHF(x, z) and KCmNOHF(x, z) . For normalization purposes, we define the maximum power of all KCm epochs regardless of type and frequency for the subject x, which we denote by K CM(x). With these definitions, we can write the Normalized Power (ND) of the yth epoch of type KCmWSHF of subject x at the frequency z relative to the mean of the KCmNOHF as and finally, the corresponding Percentage Difference (PD) The results are plotted in Fig. 4B using boxplot to show the spectral differences between KCm events at, and between HRS. The REM0 periods of sleep can be approximately defined by introducing a second iteration of sleep scoring/labelling to the sleep stages defined in the classical way (Rechtschaffen and Kales, 1968) . We use V HR and V EEG to refer to the measures of variability (variance) of HR or GFP of the electrophysiological measures, respectively. We compute in exactly the same way corresponding thresholds, T HR and T EEG , using the signal at the beginning of the recordings (with eyes closed and well before the first transition to sleep). After the conventional sleep staging is completed, each 30-s period is re-examined computing two variability measures, one related to the variability of the cardiac rhythm and the other related to measures of the electrical activity of the brain that can be extracted from the EEG, MEG or EMG channels. For each comparison, the expert classification is modified according to the following criteria: The critical element for the cardiac criteria is the presence or absence of surge, it is therefore necessary to use sufficiently long period for the computation of the variance (V) of the HR. We set the lowest duration for the window for V HR computation to be ten (10) seconds. The measure of other variances (i.e. for EEG, MEG and EMG) can be reduced further, and we have used for V EEG windows as small as 200 ms. The HRV defines a candidate REM0 period and its duration and this is used to confirm if within this period the second criterion (V EEG > T EEG ) is also satisfied. We will present results using windows for the V HR ranging from ten (10) seconds to three (3) minutes. The actual tomographic analysis employs the non-linear magnetic field tomography (MFT) algorithm (Ioannides et al., 1990; Ribary et al., 1991) which is designed to extract the maximum information out of real time MEG data. This is achieved through independent MFT computation applied to each snapshot of single trial data MEG data (Taylor et al., 1999) . MFT can also be applied independently to each snapshot of average MEG data, the average performed over a time window (or after filtering) or across selected clusters of STs. The output of the MFT analysis is an estimate of the continuous primary current density vector J(r,t). For each snapshot of data, a linear integral expression relates the signal d m from the mth sensor (set of coils) to the primary current density vector J(r) , where for simplicity the time variable has been dropped since each analysis is performed independently for each time point. The integration is over the source space, Q, which encompasses all regions which can contain active generation of electrical activity (primary currents) which are the generators of the MEG signal of interest; in our case Q is the entire brain: The lead field Φ m (r) is a vector function that is completely determined by the geometric properties of the coils making up each sensor and the conductivity details of the biological medium. The above equation is the solution of the forward problem (FP) and it shows that the FP is linear. However, the inverse problem (IP) is not necessarily linear. To see this let us write the general form of the current density vector as an expansion in terms of the product of the lead fields (thus ensuring that no components can arise in the inaccessible parts of the source space) and a weight function w(r, J(r)), where we have explicitly allowed our modulating weight function to depend on the unknown current density vector. We further simplify the weight function by assuming that its dependence on J(r) is through some power of its modulus, J(r), J(r) = |J(r)|. This leads to the ansatz, Substituting (2) in (1) yields the general family of MFT solutions, with the standard MFT, as it has been consistently used through the last three decades corresponding to the choice p = 0. The optimal choice of p was first investigated through simulation campaign (Ioannides et al., 1990) and then again through a detailed analysis of lead field properties (Taylor et al., 1999) . In both cases the standard MFT (p = 0) was superior: the numerical simulations showed a better performance (Ioannides et al., 1990; Taylor et al., 1999) and the theoretical analysis of the lead field had optimal properties for tomographic analysis (Taylor et al., 1999) . The general family of MFT solutions leads to a highly complicated system of equations except for the case with p = − 1, which is the only one that leads to a linear system of equations that is J-independent. The family of linear solutions (p = − 1) includes the minimum norm (MN), weighted MN (wMN) and LORETTA (where weights are introduced depending on the inverse square Laplacian operator). These methods lead to a linear system of equations' and therefore they are popular. The computational simplicity and efficiency of the family of linear methods is bought at the expense of physics. The choice p = − 1 removes the dependence on J and in doing so it makes it possible to use the same single scalar function to recover any current density without further recourse to the data. The very nature of the non-uniqueness of the inverse problem suggest, and the underlying physics shows, that it is not possible to recover both the direction and strength of the current density anywhere in Q with a kernel that is independent of the current density vector (Taylor et al., 1999) . The standard MFT corresponds to the assumption that only the direction of the unknown J(r) can be expressed as a (weighted) linear sum of lead field functions. The full current density must be obtained from a highly non-linear system of equations for each snapshot of data. Specifically, the strength must be determined more explicitly from the MEG signal itself. The sound theoretical base brings two more consequences that endow standard MFT with properties needed for tomographic analysis. Firstly, the final expression to be solved has no dependence on the modulus of the unknown current density vector, because this appears on both sides of the defining equation (2) (with p = 0), a property that allows for sharp discontinuities in the current density vector with small values of the expansion coefficients. Second, standard MFT satisfies the principle of least sensitivity to both variations of the data and iterations of the non-linear norm constraints; a more detailed mathematical analysis for the two properties can be found in (Taylor et al., 1999) . Therefore, standard MFT is designed for real time (single trial) analysis because, it draws by necessity on all available data: because linearity is lost, the direct appeal to the data must be made on every time slice of the data, so a new non-linear system of equations must be solved each time. It is in this sense that MFT uses all information in the MEG signal. The heavy computational load that used to be a problem in the early years when MFT was first implemented, is less important today and this is partly the reason that the re-analysis of the old data that can be done almost routinely today provides so much new information. Throughout the manuscript we will refer to standard MFT (i.e., the choice of p = 0, amongst the wider MFT family) as MFT. The methods described above were developed specifically for further analysis of the MFT estimates of the instantaneous continuous primary current density vector J(r, t). They can be applied in exactly the same way for the results of any analysis capable of producing an independent estimate of continuous primary current density vector J(r, t) for each time slice of data. We first distinguish analysis using either the time or frequency domains. The continuous vector field extracted from MFT is discretized for storage and further processing as vector time series for each voxel (a point on the grid used to store the MFT solutions) and for each ST of each condition of an experiment. Post-MFT analysis is performed on the resulting time series. Here we will be concerned with statistical comparisons, which are performed voxel by voxel and the results are reported after the conservative Bonferroni correction is applied for multiple voxel comparisons. The current study employs the same methodology and in particular uses the same conceptual framework of using MFT solutions computed in the time domain to obtain statistical comparisons either in the time or frequency domain. The time domain statistical comparisons are suitable for evoked responses, where the onset of stimuli or events provide a time-locking signal that can be used to align individual single trials (ST). The frequency domain statistical comparisons are suitable when state comparisons are to be made, and this is the case for studying sleep. Even in the case of events like spindles and KCs the onset of the event is not precisely defined; even if the onset could be defined accurately, each event has a very different pattern of underlying activity so any pooling of similar events must first attempt to assign very different individual events into clusters showing more within cluster homogeneity. In all our recent sleep studies we have used statistical comparisons confined to specific frequency bands. In our first such analysis, we used sleep MEG data filtered within different spectral bands (Ioannides et al., 2009) to compare different core periods. Specifically, we used the time domain estimates of activity extracted from band-pass signals filtered in two bands: wide (1-200 Hz) and gamma band (25-90 Hz). Each snapshot of the band-passed data was analyzed tomographically. The strength of the resulting estimates, i.e., the modulus J(r) = |J(r)| at each grid point, r, in the brain was used as the basic element of comparison. We used 8 trials, each of 4 s duration, extracted from the core periods of each sleep stage (NREM1c, NREM2c, NREM3c, NREM4c, REM) and for the ECW condition to define the six-core distribution. Since each 4 s trial had 2500 samples, the distribution for each core period had (8 × 2500 = 20000 samples). The distributions for each sleep stage were compared with that of ECW and with that of each other sleep stage, with separate statistical comparisons for the wide and gamma bands. Because of the large number of samples, the comparisons gave highly significant results, typically with p < 0.00001 after Bonferroni correction. In our more recent studies' we have standardized the frequency domain analysis. The details of the analysis are fully described in Fig. 1 of (Ioannides et al., 2017) . This methodology was applied for studying the way the NREM2 hallmarks KC1s and spindles emerge from the NREM2 core periods. The results were reported in (Ioannides et al., 2017 (Ioannides et al., , 2019 and some of its implications in (Ioannides, 2018) . In summary, the new method utilizes the MFT time domain reconstructions as a basic analysis unit: MFT is applied independently to each snapshot of time domain wide-band MEG data. Typically, we cut from continuous 3-min-long data (with low head movement) segments with duration of 2, 4 or 8 s. For the work in this paper, all time domain segments were 4 s-long and were divided into two, 2-s parts for the Fourier analysis. For core periods, the two 2-s data from each 4 s segment were equivalent. For GLEABS the center of each 4 s (zero of time) was at the onset of the first event, e.g., at the onset of the first KC in a KCm sequence. Therefore, the first 2 s segment was from a pre-event period, while the second 2-s cut was from the post-event period, i.e., it contained the activity correlated with the event. Since the MFT reconstructions were performed independently for each sample of data, the collection of samples for each voxel (point on the grid used to store the MFT solutions) described the time evolution of J(r) continuously from the pre-event through the onset and the post-event period for each GLEABS trial. Each 2 s segment of J(r) was Fourier transformed, component by component and the 3D Pythagoras' rule was applied to compute the spectral amplitude at a given frequency. For each 2 s time domain segment we computed the spectral density from 0.2 to 98.4 Hz with step of 0.2 Hz. The usual statistical parametric mapping (SPM) analysis is applied to the spectral density, what we call simple spectral SPM (sSPM) comparing the activity between conditions within a specific frequency band. In all statistical comparisons reported in this work, we contrasted the spectral power in windows of bandwidth 3.2 Hz and used a step of 1.6 Hz to cover the range from 3.2 to 94.6 Hz. The result of a statistical comparison is reported as (Bonferroni corrected) p-value for the frequency at the center of each bandwidth. In our first analysis with the new method (Ioannides et al., 2017) we demonstrated that the new analysis reproduced the results of the 2009 study. The foundation of our results are simple sSPMs applied separately to each subject. For much of the work we will report here, we use the same general analysis steps as in our recent sleep studies (Ioannides et al., 2017 (Ioannides et al., , 2019 , but in one respect we extend the grand sSPM analysis when we widen the comparisons to combine together sSPMs for events from different sleep stages, comparing each one with the same baseline; for consistency we will report here comparisons using as common baseline the NREM2 core (NREM2c) period. We accommodate this generalization by defining two types of grand statistics: the first is obtained from the contrast of the spectra of one condition (usually before the event) and that of NREM2c; the second is obtained from the contrast of two conditions, both from either the pre-event or both from the period during the event. In this work, we will report two composite statistical comparisons; these are obtained by demanding that the statistical significance quoted is reached at the voxel level by all or a subset of subjects and conditions. We will refer to the results of statistical comparison across subjects in one condition as grand sSPM (g-sSPM). A g-sSPM map provides for each voxel the percentage of subjects that have satisfied statistical significance above a given threshold and in a specific direction. A grand-sSPM comparison can then be generated by finding common changes across different sSPM comparisons, for example, different subjects and/or different conditions. Integrating results across subjects requires transforming to a common Talairach space the solutions for each subject and then identifying the voxels showing common activations across subjects and/or conditions either in the pre-event period or in the during-event period. We will refer to an even more generalized result, where statistical comparisons are satisfied across subjects and conditions by Grand sSPM and denote it by Gg-sSPM. In the tomography results reported in the paper the Gg-sSPM was used for Fig. 10A while in the last case (Fig. 11B ) the g-sSPM was employed. Consider the comparison between conditions A and B. Consider the following result for a g-sSPM: +100 for voxel v i and − 75 for a voxel v j . This statement means that the statistical comparison between the distributions for conditions A and B, produced a result satisfying the specified threshold for A > B for all subjects for voxel v i and passed the same threshold for 75% (e.g., 3 of four subjects) for voxel v j but this time for A < B. Next, Consider the Gg-sSPM for N subjects and M pairs of conditions, i.e., NxM statistical comparisons, with the first, A, condition corresponding to the pre-GLEABS period and the second condition NREM2c. Consider the Gg-sSPM results +96 assigned to a voxel v i and − 90 to a voxel v j . This means that for voxel v i , the comparison between the distributions produced a result with A > B at the predefined p-value, in 95% of the cases (e.g., failed for just one condition out of 5 in only one of the 4 subjects). For voxel v j , the result means the same threshold was satisfied (but this time with the second condition higher than the first) for 90% of the cases (e.g., failed for just two conditions in one subject or for one condition in one subject and again for one condition (same or different than before) for another subject). The MEG environment is not ideal for sleep and for this reason strict selection criteria were used to pick the subjects; strict selection of subjects with good sleep habits was instrumental in ensuring that all subjects selected were able to complete the acclimatization night and actual sleep recording night. The details are reported already for the three main subjects who were selected for the experiment and passed the selection criteria (Ioannides et al., 2004) . We copy here the key statistics for the sleep quality for the acclimatization and main sleep experiment nights for the three main subjects. Despite a slight reduction in total sleep time (mean: 369.33 min, SD 41.52 min; efficiency index 0.79) the percentage of REM sleep (mean 26.24, SD 11.07) and the latency to REM sleep (85.33, 25.77) were within normal ranges for young adults (Williams et al., 1974) . Although stage 3 and 4 were slightly lower and stage 1 slightly higher than normal, for each subject all sleep stages were in broadly the expected sequence: stage 1 (mean = 15.06, SD = 3.75), stage 2 (46.12, 6.6), stages 3 and 4 (12.62, 4.34). The hypnograms were judged qualitatively physiological, a fact corroborated by the subjects' report of restful sleep similar to their average night. All subjects completed a debriefing questionnaire and scored their sleep in different aspects compared to their average sleep. The average scores for the three subjects for the sleep during the recording night were overall quality (87%), refreshing (70%) and of normal depth (80%). Based on the subjects' own evaluation the sleep during the recording night was judged better than during the adaptation night (68, 53, and 63%, respectively). The fourth subject (AAI) did not pass the sleep selection criteria and did not had a full acclimatization night before the whole night sleep experiment. The subject was sleep deprived and his participation in the experiment constituted the final test of the procedures before other subjects who satisfied the sleep selection criteria were invited to participate in the experiments. Results from this subject were also analyzed and were found to be consistent with the results from the three main subjects so they have been included in recent publications (Ioannides et al., 2009 (Ioannides et al., , 2017 (Ioannides et al., , 2019 and in some of the results of this paper, but excluded from statistics about duration of sleep and individual sleep stages, including REM0, because these may be sensitive to sleep deprivation. The bulk of detailed individual results in the figures of this paper are mainly from one subject to ensure continuity of presentation, but we emphasize that the results from each one of our subjects show similar patterns and this is demonstrated by the statistics and figures summarizing across subjects' results. Fig. 1 shows for one subject the hypnogram and head movement (first and last rows respectively) with the instantaneous HR and the GFP for the MEG (151 channels) and EEG (2 channels) in between. The Figure covers the entire night sleep showing in one display the variability of these measures across the night; it is a typical variation we see for each one of the four subjects of this study. The increases in HRV are clearly seen in the second row as repeated peaks in HR which are somewhat intermittent in the first run which contains mostly SC1. The HRV increase is stronger in the second run with persistent and longlasting bouts of HR peaks. There is some regularity in the prominent patterns seen for the HR and the global field power (GFP) for the MEG (row 3) and EEG (row 4) signal, but these patterns do not match well the sleep stage periods and boundaries. During NREM sleep, the GFP for EEG and MEG tend to increase at different times relative to the HR: the short lasting HRV rise just before 3 h (system clock) occurs before the rise in GFP (in the middle of a short NREM2 segment). In contrast, the next HRV increase, around 3 and a quarter hour, is when the GFP are already high (in the middle of long NREM4 segment). On entering REM, we see a sharp drop in the GFP for both EEG and MEG while the HRV remains high. During the transition from REM to NREM, the HRV drops sharply while the GFP remains low for a little while, before increasing slowly again. The above changes in HRV and GFP are repeated in similar fashion for each sleep cycle. This pattern is often masked by large head and body movements, which hide some of the variations of the different measures of central and autonomic nervous system. As shown in Fig. 1 , during the first cycle (~first 1.5 h), the sleep is irregular with rapid interchange of sleep stages of short duration and numerous assignments of noisy and awake periods. This irregular pattern cannot be attributed to large head movements, which, as can be seen in the lowest row, is not as big during the first run as in the later cycles of sleep. Although some consistent patterns are seen in the HR and in the GFP of EEG and MEG, there are also abrupt changes in these patterns that do not match with any of the conventional sleep stage boundaries. During REM, HRV is high and the GFP for both EEG and MEG is low. To make these patterns easier to see we mark the boundaries of REM periods by red vertical lines, solid for the start and dash for the end of REM. Note that these REM patterns can be seen also for the very short REM period around 3.5 h into the system clock, when there is no large movement as can be confirmed by visual inspection of the last row of Fig. 1. An important feature of conventional sleep scoring is the identification of cycles, with each cycle containing a sequences of sleep stages, culminating in REM. This is a key feature of the original sleep staging (Rechtschaffen and Kales, 1968 ) and preserved in the revised classification (Silber et al., 2007) . The cycles culminating in REM can also be identified in Fig. 1A , but only after the first sleep cycle. A short REM period after about 1 h and a quarter, is a candidate end of SC1 but little else of SC1 has the canonical form. For the other SCs, the end of each cycle consistently coincides with the end of a REM period, and it is marked by a sudden drop of HRV, followed by a slow increase in GFP of EEG and MEG. There is however a lack of consistent correlation between the features in the different rows and the boundaries of sleep stages within each cycle, especially in the periods before REM. The pre-REM periods are often the ones most difficult to score. Furthermore, there are already suggestions to split REM into two sleep stages: tonic REM (REM-T) and phasic REM (REM-P), based primarily on the absence (REM-T) or presence (REM-P) of REMS (Simor et al., 2019 (Simor et al., , 2020 . The period before REM is often treated separately, sometimes referred to as pre-REM period. This period is often lumped together with REM (Fernández-Mendoza et al., 2009; Lim et al., 2007; McCarley and Ito, 1983; Rowe et al., 1999) or defined as transition REM (T-REM) (Brown et al., 2012) . We also notice that periods like pre-REM are seen in other parts of sleep that do not always lead to REM, which often contain large movement and often marked as noisy, thus eliminated from further consideration. It appears that the description of sleep stages is an incomplete characterization. This concern prompted us to examine closely the variations in the different measures of electrophysiology in general and especially around the REM periods. A cohesive description emerged for periods of sleep across the entire sleep sharing some common features; these periods we group together under the label REM0. All sleep stages and REM0 periods have a short duration during early sleep. Away from these early sleep periods, the start of a sustained REM0 is marked by a rather abrupt rise in HRV during NREM sleep that follows a steady rise of the GFP of EEG/MEG. It is noted that this feature describes the reverse process to what is seen at the end of REM. The REM0 sleep stage is maintained for as long as both HRV is high and peaks are seen in the EEG/MEG GFP. During REM0 periods, HRS are strongly coupled to EEG/MEG/EMG GFP and REMS. In this work we are introducing REM0 as a collective term for periods of sleep with features that are distinctive enough to separate it out and study it in isolation. In all displays we will present, we maintain the conventional sleep classification and superimpose on it the periods marked as REM0 as an additional, and at least for now, independent characterization of sleep. Nevertheless, it is a useful exercise to distinguish REM0, as if it was a new sleep stage so we can estimate in addition to the overall relative proportion of sleep it occupies, what part it would take away from each of the other sleep stages (and the periods that cannot be defined as conventional sleep stages). The way HRV correlates with the other neurophysiological measures would provide a way of distinguishing REM0 and REM: only in REM0 HRS have a regular rhythm that is highly correlated with all other electrophysiological measures. We designed and implemented a semi-automatic method for separating REM, REM0 and NREM, using the most basic elements of the properties of REM0 (see Methods). Applying the criteria to each 3-min period, we made a coarse quantification of the prevalence of REM0 and estimated the head movement under each revised sleep stage entry. The classification produced fractional prevalence 20.69 ± 8.95% for REM and 25.14± 1.56% for REM0. For these identified REM and REM0 periods, head movement above 4 mm was found in 20.96± 10.94% for REM and 35.90± 8.16% for REM0. Excluding the sleep-deprived subject (who had long periods with large movements), we computed the percentages of periods classified automatically by our tool as REM0 for each one of the conventional sleep stages (as classified by human sleep experts) and the results were: NREM2 33.64%, awake or not specified 22.73%, NREM4 15.45%, NREM3 11.83%. NREM1 10.00% and 6.36% REM. We explore the similarities and differences between REM0 and REM, by computing the mean and standard deviation of HR during periods when saccades in REM (i.e., REMS events) are present or absent. For the purposes of this comparison, we follow (Simor et al., 2019 (Simor et al., , 2020 and define phasic (P) or tonic (T) periods by presence or absence of REMS respectively; we make this distinction separately for REM sleep stage and REM0 periods. For this analysis we used periods with head movement below 2 cm because for the separation into phasic and tonic we employed only electrophysiological signals and not tomographic solutions (where we usually employ as threshold 4 mm). The results of this comparison are displayed in Fig. 2 , for each of the four subjects on the left ( Fig. 2A ) and the collective results across subjects in the boxplot on the right (Fig. 2B) . The results show that the relationship between the presence and absence of REMS (phasic and tonic distinction) with HRV during REM0 is opposite to that during REM: the phasic part of REM0 is like the tonic part of REM, while the tonic part of REM0 is similar to the phasic part of REM. Fig. 3 displays the second (and longest) continuous run of the entire night's recording (Fig. 1) . The top display includes REM0 periods defined automatically as described in the Methods section using a single sleep stage assignment for each 3-min window. The results of the new iteration after classical sleep staging (using the criteria above) are superimposed in between parts A and B, showing the additional labels REM or REM0 that the second iteration produces: a solid and dash vertical line marks the onset and offset of REM (red) and REM0 (green) with the duration of the re-classified segment printed above a solid line of the corresponding color drawn between the start/end vertical lines. Note that the classification used very broad and simple tools: we used the variance of the HR and the variance of the EEG channels. Similar results were obtained using the variance of the MEG channels and the variance of the EMG. We deliberately chose the most elementary classification (just two channels of EEG) because it is easy for others to apply the same procedure to their sleep EEG data. We maintained the 3-min intervals in the second sleep staging iteration because we only wanted to show first an approximate but robust estimate for the prevalence of REM and REM0. Two successive zooms highlight details easier to see in different timescales. The entire 3-min period of Zoom 1 has low head movement and is within a wider 36-min REM0 period. It shows an underlying infra low frequency (ILF) oscillation with periodicity around 20 s (frequency 0.05 Hz), seen particularly well in the HRSs of Fig. 3E . This ILF oscillation is also seen in the high amplitude peaks in the GFP of the MEG (Fig. 3C) and EMG (Fig. 3D ) traces, which begin just before each HRS (Fig. 3E ). The set of Zoom 1 suggests a coupling between HRS and other brain electrophysiology (i.e., increases in the variance of the EEG, MEG and EMG) during REM0. Zoom 2 goes down to fifty (50) seconds (almost twice the 30-s conventional sleep staging), showing well the graphoelements used for visual sleep scoring, particularly KCs and REMS. For easier reference to the previous results, each of the four traces of Zoom 2 contains the vertical dash lines marking each HRS and down oriented triangles marking the peaks of the EMG of Zoom 1. Fig. 3F , shows the two horizontal EOG channels; clear eye movements can be seen in each KCm with amplitude that does not seem to correlate with the KCm strength, or the HRS. The EEG trace (wide band) shows a sequence of high amplitude features that we interpreted as KCms. Some of these KCm events are clearly associated with HRS. Other KCm events are identified in between HRS, with very similar EEG amplitudes as the ones at HRS. The separate displays of the moduli of the EEG signal show that each KCm linked to HRS has a clear increase in both wide (Fig. 3H ) and high (Fig. 3I ) gamma band (55-95 Hz) ranges. Since the wide band signal in Fig. 3H is dominated by low frequencies, it is clear that high frequencies are only present for the KCm events during HRS. Consequently, as described in methods section, we label the KCm events close to HRS as KCmWSHF. Correspondingly, KCm events in between HRS have REMS but no increases in either high gamma EEG or EMG signals and are therefore labelled as KCmNOHF. Fig. 4 and Table 1 shows the clear dissociation: KCm close to HRS have significant HF increase and KCm far from HRS do not. Therefore, the same clusters will be obtained if we group KCm events according to how close they are to HRS or presence/absence of high frequencies. Fig. 4B unfolds the spectral power difference between KCmWSHF and KCmNOHF in the frequency domain. The results confirm the increase in gamma band identified in Fig. 3H and show that the increase is consistent across our subjects. The spectral unfolding in Fig. 4B shows distinctively different variance and average patterns for low (up to the beta band) and gamma band frequencies. In the low frequencies the difference at, and between, the HRS attains high values but with no preponderance of increases or decreases thus averaging around zero. For high frequencies there is a consistent positive excess (higher spectral power at the HRS than between) which rises a little in the low gamma band and becomes more prominent in the high gamma band; the dip around 50 Hz is due to the notch filter around the mains frequency. The results presented so far were obtained using smooth assignments, i.e., using a single sleep stage label for each 3-min period (corresponding to the duration between head coil activations). Such a long window is appropriate for stable periods of sleep like the sleep period of run2 which contains the second and third cycles of the night's sleep. Such long windows are not suitable for unsettled sleep or even for stable periods if finer details of transitions between stages are to be studied. We will present next computations of sleep stages and the variations of the EEG signal and heart rate using shorter windows for the computation of the variances, first for the stable periods of run2 and then for run1. Fig. 5 displays the first two sleep cycles of run2 (SC2 and SC3), each lasting about ninety (90) minutes. Each SC is represented in 3 successive rows, showing in turn the output of sleep staging, i.e., the hypnogram (rows 1 and 4, for SC2 and SC3 respectively), and the two inputs: the variances of the GFP of the EEG (rows 2 and 5, for SC2 and SC3 respectively) and the variances of the HR (rows 3 and 6, for SC2 and SC3 respectively). In each row representing the input a dash horizontal line sets the threshold for the application of the criteria for (re-)defining REM and REM0. The threshold chosen is simply the value of the mean variance of either the EEG GFP or of the HR at the beginning of sleep (i.e., the beginning of run1 well before the start of SC1). As can be seen the changes in the variance are huge, so there is a clear demarcation of values well above and well below the threshold. Clearly, the results will not change much if we had chosen for each thresholds the values at the beginning of run2, i.e., well before the start of SC2. The hypnogram rows (rows 1 and 4) show both the human experts' definition (in the low resolution of 3-min as used so far). The additional REM and REM0 labels produced by the second iteration are marked as red and green blocks respectively, placed at the top of the hypnogram (the noise level of the old hypnogram). Different windows were used for the variance of the EEG GFP (ten (10) seconds) and for the variance of the HR (thirty (30) seconds). The thirty (30) second assignments set the overall limit for invoking the use of the criteria of iteration 2 and it is determined by the HRS. Within this thirty (30) second window a finer resolution can be Fig. 3 . Key features of electrophysiological signals at different time scales within the 4 h long second continuous segment of data (run 2); same subject as in Fig. 1 . The revised automatic marking for REM and REM0 are included as red and green bars respectively. Vertical lines represent the start (solid) and end (dotted) of REM (red) and REM0 (green) periods. The first zoom (Zoom 1) focuses on the marked 3min period that fulfils the requirements for REM0 and also has a very small head movement, less than 1 mm. The second zoom (Zoom 2) focuses on the marked 50 s period containing three clear HRSs and is sufficiently detailed to relate to the graphoelements that are the hallmarks of each sleep stage. Table 1 shows the value of the p-value when carrying a two-sample student's ttest of the percentage difference in HR as described in the Methods section using the standard MATLAB function ttest2. Using null hypothesis that the percentage difference in HR for the two events for the specified subject(s) (Rows) populations are independent random samples following a normal distribution with equal means and equal but unknown variances. The alternative hypothesis is that the percentage difference for the two event populations have unequal means. achieved through recourse to the variance of the EEG GFP computed over ten (10) seconds. Note that the gap between 3-min sessions is largely artificial; it is forced by the interruption of recordings due to the activation of the head localization coils (lasting for twenty (20) or more seconds) to which fifteen (15) seconds must be added (half of the window used for the HR variance computation). Direct comparison of Figs. 3 and 5 shows the small changes made using the shorter window of thirty (30) seconds (for Fig. 5 ) rather than three (3) minutes (for Fig. 3 ). The display in Fig. 5 demonstrates how the criteria for the second iteration are applied, and as already mentioned it demonstrates the robustness of the procedure. In fact, a case can be made for increasing the threshold for the variance of the HR to 1.5 or 2 times the value chosen (i.e., 1.5 or 2 times the value average of the HR variance at the first few minutes of run1). The results with this new threshold will be almost identical to the ones presented in the figure, with only exception the elimination of the first few REM periods of SC2 which fall within ECW and NREM1 periods. The way the results are presented in Fig. 5 does not show well the yoking of the HRV and especially the surges to the rapid changes seen in the EEG, MEG and EMG. To allow for this yoking of different electrophysiological measures to the cardiac rhythm and especially the HRS we return to the presentations showing the HRV and the GFP in the next two displays. Fig. 6 , shows again the two main sleep cycles of run2 this time in pairs of rows, starting with the hypnograms (rows 1 and 2) including REM and REM0 definitions from the second iteration and ending with the head movement (rows 9 and 10). We show also the GFP for the EEG (rows 3 and 4) and EMG (rows 7 and 8); we note that there were some large movements during run2 which are likely the cause of a partial detachment of the EMG electrode early in run2 that continued throughout this run; this showed as an increase in low frequency noise that nevertheless is not affecting the location of the EMG peaks that are the element we wish to highlight in Fig. 6 . Finally, the HR is displayed in two middle rows (5 and 6). The peak of each HRS is identified, and its position marked by a pair of upward and down facing arrows at the top Fig. 4 . Physiological and spectral differences between KCmWSHF and KCmNOH-FAcross subjects changes relating HRV for GLEABS and core periods and changes in spectral power for KCm periods in REM0. In the boxplots, the median is shown by a red line, 75th/25th percentile shown by upper/lower blue lines, a cross represents an outlier and a dotted line represents the range (excluding outliers). We compute separately the changes for the two types of KCm, the ones close to HRS, which have strong high frequency spectral power (KCmWSHF) and the ones in between or well away from HRS, which have NO strong high frequency spectral power (KCmNOHF). (A) Percentage change in maximum computed heart rate (computed for each event separately) during the first 4 s starting at the onset of GLEABS or start of core period relative to the mean heart rate of the 2 s before. This percentage difference is computed separately for each member of the GLEAB family and the quiet periods of each sleep stage. The statistical significance of the result for each subject and across subjects is listed in Table 1 (B) Boxplot of normalized difference for each subject. The power of the KCm periods at the HRS relative to the mean power of KCm in between HRS for each frequency bin. In order to avoid bias from electrode placement and different skin conductivity across each subject the percentage normalized difference was calculated for each subject separately before the boxplot was computed for each frequency. More details and formulae for this case (B) can be found in methods. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) Sleep classification by the experts and after our additional iteration step (run2). The figure shows the results for the second run of the same subject as in Fig. 1 , classifying with each method succesive 30 s windows (six such windows in each 3 min continuous recording). For the second iteration we used the variances of the EEG and Heart rate. New sleep classification is shown on the top row of the hypnogram (corresponding to the NS level in the conventional hypnogram) with red and green blocklines indicating REM and REM0 respectively. In the rows displaying the EEG and HR variances, dotted black lines mark the thresholds used by the automatic scoring routine. The plot is composed of 2 triplets (A,B,C) and (D,E,F) each presenting the same measurements for the 2nd and 3rd sleep cycle of HBD005 respectively. Top and bottom triangles represent Heart rate surges as identified from the program. A green set of triangles is shown if the identified surge corresponds to an increase to the power of both EEG and EMG, a blue set is shown if the surge corresponds to an increase in either the power of EEG or EMG while a red set of triangles shows no corresponding increase in the power of EEG or EMG within a minute window of the identified surge. (A, D) Hypnogram; (B, E) Variance of Global Field Power of the two EEG channels (C3-A1, C4-A1) between 0.75Hz and 95Hz and displayed using a log scale. (C, F) Variance Instantaneous Heart rate. The representation of Fig. 5 shows the quantities used in the criteria of the second iteration. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) and bottom of the display. The period from 1 min before to 1 min after each HRS is searched for peaks in the EEG and EMG GFP and the result is quoted in the color of the corresponding downward and upward facing arrows of each peak: (green → if a peak is found in both EMG and EEG), (blue → if a peak is found in either the EMG or EEG, but not in the other) and (red → if no peak is found in either the EMG or the EEG). The distribution of arrows shows clearly that the periods defined as REM0 correspond to the green triangles showing convincingly the yoking of the EEG and EMG to the HRS. The display of Fig. 6 shows well the highly regular HRS sequence with corresponding regularity in the peaks of HR (row E) for the duration of the long REM0 period (row A) between 200 and 230 min of the system clock, followed by an abrupt dissolution of this regularity as REM proper begins. While the regular HRS lasts, there is a strong association of peaks in the EEG and EMG with HRS, but this becomes weak or non-existent abruptly with the end of REM0 and remains so for the duration of REM. Run1 contains the first cycle of sleep which, as is often the case, is characterized by many transitions between awake state and sleep stages of light sleep with brief duration. The first sleep cycle of the night usually ends with a brief REM period and lasts between sixty (60) to one hundred (100) minutes a somewhat shorter duration than the other sleep cycles that usually last between ninety (90) and one hundred and twenty (120) minutes. The three subjects that were not sleep deprived in our study all had a first sleep cycle lasting between sixty-five (65) and eighty (80) minutes, ending with a short REM period lasting between one and 3 min, under the conventional sleep staging procedures. Fig. 7 shows the SC1 for HBD005 in the same format as that used in Fig. 6 for the first two sleep cycles of run2 (SC2 and SC3). Using a single sleep stage classification, i.e., forcing the single sleep stage assignment for each 3-min continuous recording, shows only glimpses of the canonical form, starting with light sleep and apparently jumping into a brief REM session at the conclusion of the short (~75 min) sleep cycle. The introduction of the new iteration after the sleep staging process (using thirty-and 10-s windows for the computations of the variance of HR and EEG GFP respectively) brings forth numerous occurrences of intermittent short periods with the additional second iteration labels of either REM0 or REM. The additional new labels are displayed in Fig. 7 , in exactly the same way as in Fig. 6 : as red (REM) and green (REM0) blocks placed at the top part of the first row. The third row in Fig. 7 shows the HRV and uses the arrows to mark the presence of peaks within 1 s of the HRS; as before, a green arrow is used to mark the presence of a peak in both the GFP of EEG and EMG while a blue arrow marks the presence of only one GFP peak, either for the EEG or the EMG and finally a red arrow when there is no GFP peak. These arrows show that during the first cycle of sleep of HBD005, the peaks of the GFP for EMG and EEG are tightly yoked with HRS during REM0 and more weakly, usually to only one of the two, during REM. This yoking is established early during sleep and it is already evident in the newly labelled REM0 and REM periods defined during what the human experts identified as awake (ECW). For early sleep in particular, we resist the temptation to proceed with proposing these newly labelled periods as putative REM0 and REM periods, because more work is needed to clearly distinguish these periods from each other and from the awake state. Nevertheless, it is worth noting that, if our REM0 classification for SC1 is counted as sleep then the duration of SC1 increases from 65 to 80 min approaching the canonical SC length of the later parts of sleep. The REM0 periods are often associated with head movements and strong EMG activity and for this reason they often correspond to periods marked as noisy by experts. Careful inspection of the data however Fig. 6 . The regular Heart Rate Surge rhythmicity of REM0 in run2. The figure shows again the results for the second run as in Fig. 5 , classifying with each method succesive 30 s windows. For this figure we keep the same conventions for the different displays as in Fig. 5 , but we show in successive rows the same measure for the SC2 (rows A, C, E, G and I) and SC3 (rows B, D, F, H and J), with successive pairs as follows. Rows A and B show the hypnograms. Rows C and D show the actual EEG GFP (C3-A1,C4-A1) between 0.75Hz and 95Hz using a window of 3 time-slices (4.8 ms). Rows E and F show the actual instantaneous HR. Rows G and H show the GFP of chin EMG electrode between 20Hz and 95Hz using a window of 3 time-slices (4.8 ms). Rows I and N show the head movement within each 3 min. The representation in this Fig. 6 shows well the regular HRS sequence during REM0 periods, which stand out prominently during the extended REM0 periods from 200 to nearly 240 min. shows that this is not always the case as there are clear REM0 periods corresponding to very small head movements well below the threshold of 4 mm we have used for acceptance for tomographic analysis. Specifically, SC1 has relatively low head movement except for a short period just after sixty (60) minutes and at the very end of the REM period (as defined by experts) marking the end of SC1. The large head movement at the end of run1 seen prominently in Fig. 1 is beyond SC1, just outside the period displayed in Fig. 7 . During the sleep cycles of run2 (SC2 and SC3) the REM0 and eventually REM periods coalesce into longer continuous periods, justifying the use of single sleep label for most of the intervals of 3 min continuous recordings. In contrast, during run1 sleep stages, and the additionally Fig. 7 . Using REM0 to assign sleep stages for hard to classify early sleep in run1 The figure shows the results for the first run of the same subject as in Fig. 1 , classifying with each method succesive 30 s windows (six such windows in each 3 min continuous recording). For the second iteration we used the variances of the EEG and Heart rate. The conventions and displays are the same as in the last Figure ( labelled periods as REM0 and REM by the second iteration, maintain separate and distinct presence for only short time intervals. Fig. 8 displays a zoom on ninety (90) seconds of data (roughly from 33.8 to 35.4 min in Fig. 7) which contains 3 successive thirty (30) second periods which experts classified as (NREM2, NREM2 and NREM1); these periods were re-assigned by our second iteration to the labels (REM, REM and REM0). The remaining rows demonstrate the difficulty one faces when one tries to do sleep staging during periods that our second iteration assigns to either REM or REM0, like the periods used in Fig. 8 . The presence of REMS in the horizontal EOG channels in all three thirty (30) second periods would justify the REM label; on the other hand, the EEG (C3) channel has graphoelements that would also pass as graphoelements of NREM2 (spindles and KCs) or NREM1 (VSWs). There are also high frequency activities consistent with muscular artifacts, especially around the end of the period. While in the conventional classification the conundrum is difficult to resolve, the second iteration uses the presence of surges in the HR and the simultaneous presence or absence of high variance in the EEG and EMG GFP to determine whether these periods should have the additional new labels of either REM0 or REM, respectively. The distinction of REM or REM0 is further confirmed by the presence of a peak in the EEG and MEG GFP for REM0 and only in one of them in REM (not shown in Fig. 8 ). We present here the simplest form of clustering analysis using feature extraction with PCA analysis of the sleep staging as these are classically defined. The starting data are regional spectra of 8 or more exemplars of core periods for each sleep or awake category of data to be analyzed. In addition to actual measurements from the brain, the measurements with the empty room are also analyzed with MFT, as if a subject was in place inside the helmet, and 8 exemplars of baseline, or ghost, regional spectra are obtained and used in the PCA analysis. Two distinct computations are presented. For both computations the same set of ROIs is used; this set has 29 ROIs, which are spread throughout the cortex, the length of the cingulate cortex and brainstem; these ROIs have been identified during earlier sleep studies. In Fig. 9 , we show the results of the two computations using displays in the feature space of the reduced 3D space, defined by the three eigenvectors in the ranked list of PCA eigenvalues. For each computation, the exemplars of each category (noise, each sleep stage and/or pre-and post-periods of KC1 and spindles and EC, EO and EF conditions) were collapsed to a representative center point. The center point describes the centroid of the 8 exemplars (defined by human experts according to the classical sleep staging process), in the reduced space of the PCA analysis. The full set of exemplars and the centers were then clustered using the k-Means algorithm (Kanungo et al., 1992) , providing a data driven partition of the data. Finally, the centers representing core sleep stages were expanded one by one and if more than one clusters were evident, separate centers were made for each sub-cluster. Evidence for sub-clusters was strongest in the case of REM and we will represent them, together with the nodes of the 8 exemplars of REM in the displays of Fig. 9 . We use the feature space for the displays of the results so that both centers and individual exemplar nodes can be represented in a way that preserves the natural metric of the data. Note that the k-Means analysis was made without any influence exerted on sleep stage membership, so the placement of exemplars, clusters and sub-clusters produced by k-Means was unbiased. The node naming defines a priori definition of category, e.g. sleep stage defined by classical sleep scoring by human experts. The color of each node identifies the k-means cluster. The centers of a set of exemplars is shown as big filled circles with the prefix C in the label. For the case of REM three centers are shown: the one using all exemplars (C REM) and the centers for the two sub-clusters (C REM-a and C REM-b). In addition, the display includes the nodes for individual REM exemplars (with the label rem). The first computation uses only the noise, ECW and sleep core periods, so there are 56 exemplars. The maximum frequency resolution is used for the regional spectra (i.e. from 0.2 to 100 Hz with step of 0.2 Hz making a total of 50 frequency points). The PCA analysis is therefore performed with approximately 80,000 points. Fig. 9A shows the results of the PCA in the 3D space defined by the first three eigenvectors. The kmeans clustering with the number of clusters set to four (4) shows the following as we traverse the display from C NOISE to C ECW. The first k-Means cluster (blue) is composed of C NOISE, C REM-b and two of its three exemplars. The second k-Means cluster (red) contains C S2 and the two centers of deep sleep (C S3 and C S4), with C3 and C4 very close to each other. The third k-Means cluster (cyan) has the center of all REM exemplars plus the remaining exemplar of C REM-b and one of the five exemplars of C REM-a. Finally, the fourth k-Means cluster (green) is composed of The C ECW, C REM-a and the remaining four of its exemplars. For the second computation two more distinct groups of data are In addition the centroid of the noise exemplars (that can be represented as proxy origin of the coordinate system) and the centroids of the two subclusters of REM are also shown as large filled circles. The display also shows the actual eight (8) exemplars of REM as smaller filled circles. The labels and centroids represent the classical sleep staging output of the human experts. The color of the nodes is the independent clustering result of the k-means algorithm using 4 clusters. B) The same as (A) but with the addition of centroids for the EC, EO and EF conditions (36, 34 and 34 exemplars respectively) and the pre-and during KC1 and spindle periods (8 exemplars each). The three principal axes in (B) are not related to the three principal axes of (A). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) added. The first group consists of 4 sets, each with 8 exemplars: a set from the period before and a set for the period after, for each of the two main hallmarks of NREM2, KC1s and spindles, selected to be well away from NREM2c periods. The second group consists of exemplars from a different experiment with the same subject, performed some 8 years after the sleep study. In this experiment the MEG recording was used Fig. 10 . The Relationship of REM and REM0 through the night sleep Cycles of REM0 across subjects. The display shows the duration of REM0 as this was defined by the automatic scoring routine. The criteria of the second iteration were applied to 30 s windows for EEG and Heart rate variances for the whole sleep for the 3 subjects excluding the sleep-deprived subject. The results are displayed in the following order. with eyes closed (EC; 36 exemplars), and two conditions with eyes open, one just looking straight ahead but not focusing on anything in particular (EO) and the other fixating on a dot placed at the center of gaze (EF), each with 34 exemplars. The second computation uses as input the same 29 ROIs but only the averages within each of the eight (8) standard frequency bands from delta (1.5-4 Hz) to high gamma (55-95 Hz) for the augmented set with more than three times the exemplars (192) . The PCA analysis is therefore performed with approximately 45,000 points. Fig. 9B shows the results of the PCA analysis using the augmented data, in the 3D space defined by the first three eigenvectors. The k-Means clustering with the number of clusters set to five (5) shows the following as we traverse the display from C NOISE to C S4. The first k-Means cluster (cyan) is composed of C NOISE, C REM-b all its four exemplars and all the centers of the EC, EO and EF. The second k-Means cluster (red) is composed of C ECW, CS1, C REM, C REM-a and all its four exemplars. The third k-Means cluster (yellow) groups together C S2 and the period just before KC1 events (C S2-KC1-bef). The fourth k-Means cluster (green) puts together the pre-spindle centers (S2-SP-bef) and the periods during the NREM2 hallmarks (S2-KC1-dur and S2-SP-dur) with the first of the two deep sleep centers (C S3). Unlike part A of the Figure the center for the second deep sleep center NREM4 (C S4) is far away and it is the only member of the final fifth k-Means cluster. Fig. 10 shows the variation of REM0 properties across sleep cycles for each subject. The first row shows the mean duration of REM0, on the left in absolute terms (in minutes) and on the right as a fraction of each sleep cycle. The general trend is for REM0 as a percentage of the total duration of each SC to be high in the first two cycles and to reduce in the last three cycles. The next three rows show the ratio of REM0 to REM across SCs for each subject. We show the result for each subject on a different subplot because the range is very different across the subjects. The overall pattern though is similar across subjects, with high values (more REM0 rather than REM) in the first two cycles and much reduced ratio in the last SCs. In summary the results in Fig. 10 show greater periods spent in REM0 in the first two cycles which reduces to near parity (HBD002 and HBD012) or even well below unity (HBD005). We finally turn to the tomographic analysis to explore which areas are activated in the 2 s before the start of GLEABS. We first ask whether one or more areas are activated for each one of the four subjects, before each one of the three types of KC events (KC1, KCmWSHF and KCmNOHF). We then narrow the question, for the same pre-event period, to a direct contrast between KCmWSHF and KCmNOHF. Since KCmWSHF occur only at the HRS, i.e., at the peaks of HR and KCmNOHF in between HRS, i.e., when HR is at its lowest, the [KCmWSHF vs KCmNOHF] contrast probes also the contrast at the HRS relative to the low HR periods between HRS. All results reported in section 3.7 are from 3-min segments with very low head movement (<4 mm) and they are for common increases in activity identified in each one of the four subjects, with some additional cases seen in three of the four subjects. For both cases, i.e., the ALL (4 subjects) and 75% (3 of 4 subjects), statistical significance was demanded at p < 0.05 after the conservative Bonferroni correction for multiple voxel comparisons is applied. Fig. 11A shows the only Gg-sSPM single cluster, surviving all 12 simple sSPMs contrasts between the pre-KC and NREM2c. This cluster is identified for each of the four subjects in each of the three contrasts [KC1 vs NREM2c] , [KCmNOHF vs NREM2c] and [KCmWSHF vs NREM2c] with higher spectral power in the alpha band for the pre-KC1 period compared to NREM2c with p < 0.05 after the conservative Bonferroni correction. The cluster center has Talairach coordinates (TC) at: x = − 7, y = 6 and z = − 6. The cluster encompasses the Nucleus Basalis of Meynert (NBM) in its caudal part, while its rostral part is in the overlap of NBM and the ventral striatum (VS) (Li et al., 2014) . The pre-event investigation was concluded with the single condition contrast [pre-KCmWSHF vs pre-KCmNOHF] . This g-sSPM had only one common increase across all four subjects in two successive bands (of width 3.2 Hz) centred at 16 and 17.6 Hz with two foci both on the left paramedian sagittal plane (Fig. 11B) , one in the mid-brainstem and the other in the adjacent anterior part of the cerebellum. The dorsal area was on base of the pons (TC: 8 -33 -31) consistent with the Raphe Magnus nucleus; the second focus on the left posterior uvula and its adjacent left tonsil of the cerebellum (ULTC) (TC: 7 -51 -44). Relaxing the criterion to three of the four subjects the g-sSPM identified more foci across the frequency spectrum. In the left paramedian locations the increases seen in all four subjects in the high sigma/low beta range (16 and 17.6 Hz) were expanded for three subjects in the two lowest bands, centred at 3.2 and 4.8 Hz. Increases were also identified in companion locations in the mid-sagittal slice at medial Pontine Reticular Formation (mPRF) (TC: 1-31 -33) and uvula (TC: 3-54 -48), but only in the delta band (3.2 Hz). In more lateral areas, common increases for three subjects, were identified in low frequencies (3.2-6.4 Hz) in the left (TC: − 31 0-28) and right (27 0-30) amygdala. Increases common to three subjects were also identified in the high gamma band in the same areas as above, most prominently in the mPRF (60.8 and 62.4 Hz and again from 83.2 to 94.4 Hz), uvula (60.8-67.2 Hz), uvula and ULTC (73.6-94.4 Hz) and the areas immediately in front of the left NBM (from 80 to 94.4 Hz) that included NBM in the high end of this band (89.9-94.4 Hz). During a similar period as the mPRF increase, a group of areas along the right lateral aspect of Anterior Cingulate Cortex (ACC) and frontal cortex were also activated: right paramedial anterior lobe from 60.8 to 64 Hz (TC: 7 52-10); right rostral ACC (rACC) from 60.8 to 62.4 Hz; (TC: 6 46 8); right pre-genual ACC from (TC: 6 27 22). In this work we analyzed whole night MEG sleep data together with the auxiliary electrophysiological signals that are recorded during polysomnography, across a wide range of temporal scales. Critically for our analysis, the data included a recording of the head location every 3 min, allowing the clear identification of continuous 3-min segments of data with very low head movement (below 4 mm). Focusing our investigation within segments with small head movement, we identified time periods for which sleep staging was problematic when the conventional sleep staging criteria are employed (Rechtschaffen and Kales, 1968) . Even if the new sleep staging criteria (Silber et al., 2007) were applied to force agreement, this agreement will simply hide the inconsistency of finding distinct graphoelements, each used already as hallmarks of a different, single sleep stage; this seems arbitrary with no obvious explanation in terms of underlying mechanisms or processes. In line with the Research Domain Criteria (RDoC) framework (Insel et al., 2010) we introduced objective measures as an addition to the conventional sleep staging: a second iteration is added, which uses quantitative properties that can be objectively computed from the raw electrophysiological signals, particularly the HRV and specifically the presence of HRS. In an effort to maintain continuity with the well-established and vast sleep-scoring literature, we introduced the second step as a refinement of sleep staging. The new iteration is particularly effective in labelling periods of sleep that conventional sleep staging either has difficulty or even completely fails to unambiguously assign to a single sleep stage; the second iteration either leaves the assignment the same (hence no new label), or adds an additional label (either REM or REM0) to that already assigned by human experts. With the head movement under control, we could retain periods of moderate even fairly high signals that nevertheless are likely to have neural origin. Focusing on these low head movement periods allows us to describe in some detail the properties of a range of neurophysiological signals during REM0 periods and define its prevalence during each of the classically defined sleep stages. We finally proceeded to compare the tomographic estimates of activity before GLEABS, focusing here on comparisons between each KC type with NREM2c. We followed in this analysis the methodology we used in recent studies separately for spindles (Ioannides et al., 2017) , and for spindles and KC1 (Ioannides et al., 2019) . In the analysis presented here, we focused on KCm events that are the prominent high amplitude feature of REM0 that is often eliminated by threshold setting ostensibly to eliminate noisy segments. We were able to include KCm events with confidence that the high amplitude signal was not due to head movement by analyzing only events from 3-min continuous segments with head movement below 4 mm. We first point out four general results generated by the analysis at the level of signals. Firstly, the comparison of the results of our signal-based analysis, particularly the way HRS relate to other electrophysiological signals, with conventional sleep stages highlights periods where conventional sleep staging is challenging or even problematic. These are periods early in the night's sleep during SC1 or periods around REM sleep. Our interest in HRS was motivated by one of the original aims of this work: the identification of PGO waves in humans. In animal experiments, it has been demonstrated that the generation of hippocampal theta waves and PGO activity is highest during REM and pre-REM periods peaking around HRS (Rowe et al., 1999) . As already discussed in subsection 3.3.1, we found prominent HRS in periods before REM, which have been referred to as pre-REM and transition REM in previous studies, but also in other parts of sleep that do not always lead to REM. Introducing the label REM0 gathers together these periods into one concept, and simply adding this label to conventionally defined sleep stages, makes sleep staging less ambiguous and subjective, particularly in early sleep and around REM. The periods carrying the label REM0 are characterized by the presence of HRS and they have in general increases in activity in all electrophysiological measures exceeding the levels recorded during ECW periods. The second iteration after conventional sleep staging uses two simple quantitative and objective criteria to easily and cleanly separate out REM0 periods. The criteria require high variances in the HR and one other (usually all) electrophysiological measures, suggesting that REM0 is a period of sustained high arousal. REM0 periods are not easily accommodated by conventional sleep staging, partly because they contain a combination of features with each of these features individually used as hallmarks of distinct sleep stages. REM0 is also distinct from REM in terms of the coupling between HRS with EMG and the EEG/MEG GFP. Finally, REM and REM0 phasic and tonic periods, defined purely in terms of REMS (Simor et al., 2016 (Simor et al., , 2019 , have opposite attributes as shown in Fig. 2 . Secondly, the duration of REM0 is not negligible, we find it to occupy a significant part of sleep, about 25%, with about one third of the original classification as NREM2 and slow wave sleep (NREM3 and NREM4) switching to REM0. The corresponding switch is less for NREM1 10.00% and lowest (but not zero) for REM with 6.36%. Using head movement up to 2 cm as threshold for accepting segments for the second iteration, we find about 23% of periods originally marked as awake or not specified to be assigned the extra label REM0 by the second iteration step; this may well be an underestimation, since much of the remaining classification, as awake or not specified, is because of large movements which may also contain a good fraction of REM0. REM0 absorbs the "difficult to score periods" under one umbrella with clearly defined and objective characteristics, leaving for classification periods that are easier to classify to one of the classical sleep stages. Thirdly, during REM0 the HRS are tightly coupled to periodic increases in other electrophysiological signals and the head movement data. This coupling reveals an ILF oscillation around 0.05 Hz. We note that an ILF oscillations in the range 0.01-0.2 Hz is also observed in the fMRI BOLD signal, with the most stable amplitude seen in the slow (0.027-0.073 Hz) band as defined by Buzsáki (Buzsáki and Draguhn, 2004; Zuo et al., 2010) . These quasiperiodic patterns in fMRI have been attributed to changes in neuromodulation initiated in the reticular formation and they were shown to modulate sensory processing (Belloy et al., 2020) . Our results point to an additional contribution to the fMRI slow oscillations due to increased blood flow during rhythmic surges generated during REM0 periods. This hypothesis can be verified or falsified with fMRI. In a recent EEG/fMRI sleep study, the spectra of the BOLD response was separately studied for the awake state and for light and deep sleep, see Fig. 1 of (Song et al.) . It was found that the mixed frequency rather low-amplitude BOLD activity of the awake state showed two peaks a modest peak at low frequency during NREM1 (0.05 Hz) and stronger peaks at low and high (0.17 Hz) frequencies during NREM2, with the low frequency reducing a little and the high frequency increasing further during deep sleep. In the (Song et al.) no significant correlation was found between cardiac (or respiratory) activity with either the low or high frequency Bold oscillation power, but the outcome may be different if HRS or the variance of HR is tested against the variation of BOLD signal. Fourthly, REM0 makes unsettle and early sleep easier to handle. There is an ongoing debate about the nature of falling asleep (Ogilvie, 2001) with suggestions of separating the transition from awake state to sleep into sub-stages. In Hori's nine-stage EEG system the conventional separation of awake and light sleep has about 3 sub-stages for each one of ECW, NREM1 and NREM2. For early sleep both the standard (30 s) and the longer 3-min scoring fail to describe large parts of the signal, usually labelling much of the signal as noise or awake. In contrast, our description of early sleep, as seen in Fig. 5 , shows for the same period many transitions between short sleep stages and periods labelled as REM0. For later sleep cycles succession of sleep stages is smoother, so much so, that many of the 3-min of continuous sleep with low head movement coalesce to the same sleep stage. Further characterization by the additional second iteration introduced in section 2.5.4 identifies long REM0 periods during the middle part of sleep. Inspection of identified REM0 periods on the signal level at finer temporal resolution (Figs. 3 and 8) identifies periods of HRS containing large graphoelements that we identified as KCm events, which cluster either just before or between successive HRS. These two types of KCm events have similar durations and amplitudes and they are similar in the low frequencies but distinctly different at high frequencies: KCm events just before HRS have strong high frequency oscillations and hence the label KCmWSHF, while those in between HRS do not and hence the label KCmNOHF. Further analysis revealed a strong coupling between KCms and eye movements identified through the EOG channels with only KCmWSHF being found near HRS and EMG excursions. One may well ask why not use clustering techniques (based on either supervised or unsupervised machine learning techniques) to objectively set the borders of sleep stages. This is indeed what we are doing in separate research projects and a glimpse of the results already obtained presented in section 3.6.1 and Fig. 9 . Two key results are obtained. Firstly, a vindication of the hypothesis that the core periods are foundational states of each sleep stage: using the spectral properties of the core states representing sleep stages obtained from the hallmarks of each place preserves the identity of each sleep stage and provide further support for separation of REM into more than one sleep stage as already suggested (Simor et al., 2020) . Secondly, the results also show that a new framework is established as soon as we allow into a clustering scheme other data (e.g. awake resting state measurements or periods just before and during large graphoelements). A clustering approach can be useful for exploring further aspects of REM0, but without proper inclusion of the HRV into the clustering scheme it is unlikely that a complete picture could emerge. The ongoing graph theoretical work suggests that REM0 may be an important conduit period facilitating the transitions between the sleep stages as these are defined by the classical sleep staging conventions. We used post MFT sSPM analysis to ask what is common across all three types of KCs, and what makes KCmWSHF different from, other GLEABS and specifically KCmNOHF. The pre-KC Gg-sSPM identified a single area encompassing the NBM and its overlap with the VS (Fig. 11A) ; in this area activity in the alpha band was increased for all subjects for all pre-KC events (12 simple sSPM comparisons). The results are consistent with a cholinergic model for the generation of GLEABS (Everitt and Robbins, 1997; Irmak and de Lecea, 2014; Obermayer et al., 2017) : alpha band activity in the NBM sends an arousal wave across the cortex that promotes all KC types and possibly most if not all VSWs (results not shown). The different types of GLEABS may then emerge as different manifestations of this arousal wave. The appearance of a specific type of GLEABS is likely to depend on the background neuromodulation levels in each sleep stage (Hobson, 2009) , especially the availability of other neurotransmitters like noradrenalin (Lelkes et al., 2013) . All the areas identified before GLEABS, the ones we reported here and the ones with lower threshold including the VSW events and the ones in our earlier studies focusing on NREM2 sleep stage with KC1 and spindles (Ioannides et al., 2017 (Ioannides et al., , 2019 are known to be connected with the NBM, as recently documented using resting state fMRI measurements (Li et al., 2014) : one of the strongest connections reported was with the Dorsal caudal part of ACC (dcACC)/Supplementary Motor Area (SMA)/pre-SMA complex, almost exactly matching our identification of the key generator area of KC1 (Ioannides et al., 2019; Voysey et al., 2015) . In Fig. 3 of the recent EEG/fMRI study of (Song et al.) the strongest correlation between both the low (0.05 Hz) and high (0.17 Hz) frequency BOLD oscillation is seen in the subgenual area, which is amongst the areas showing the strongest positive correlation in the resting state functional connectivity BOLD activity with the NBM (Li et al., 2014) . As expected, a wider set of activations is seen when the VSW were included and the threshold for acceptance was reduced to 3 of the four subjects. These activation reflect the struggle between balance excitation bouts with a sleep promoting action, which is seen as suppression of brain regions that respond to reward in a circuit involving the ACC, the OrbitoFrontal Cortex (OFC) and the ventral striatum (Haber and Knutson, 2010 ). An alternative view of this tussle is activity at the interface of cognitive control and decision making (Botvinick and Braver, 2015) , along the lines of our recent observation that the dcACC1, that generates KC1 events, coincides with the generation of the error related negativity (ERN), thus suggesting that the same area (dcACC1) performs the same task (monitoring the environment) during sleep (through KC1) and in awake state (through ERN) (Ioannides et al., 2019) . The critical comparison between the two types of KCm produced common increases across all subjects in the pre-event contrast [KCmWSHF vs KCmNOHF] in two places: in the midbrain, possibly the Raphe Magnus nucleus and the other at ULTC (Fig. 11B ). Lowering the threshold to three subjects revealed more areas notably in the midsagittal slice at mPRF (TC: 1 -31 -33) and uvula (TC: 3 -54 -48), but only in the delta band (3.2 Hz) and more laterally in low frequencies (3.2-6.4 Hz) in the left (TC: − 31 0 -28) and right (27 0 -30) amygdala. Increases common to three subjects were also identified in the high gamma band (56-94.4 Hz) in the same areas as above, most prominently in the mPRF (60.8 and 62.4 Hz and again from 83.2 to 94.4 Hz). The increases in the pre-KCm period are consistent with a conflict between inhibitory tendencies (increases in delta band) presumably acting in a sleep promoting role to prevent wakening up by the excitatory influences in the lower brainstem (in the high sigma and low beta bands); this is a recurring theme; it is also seen before spindles in the rACC (Ioannides et al., 2017) and in the dcACC1 area in the pre-KC1contrast with NREM2c (Ioannides et al., 2019) . Our results show that KCmWSHF events occur when there is high amplitude EMG, EOG and ECG activity; it is therefore legitimate to ask whether the KCmWSHF events are really of neural origin or are they mostly collections of muscular, eye movement and cardiac artifacts? In particular, the strong EMG activity raises the possibility that the high frequency content of the KCmWSHF is entirely due to muscle artifacts (Muthukumaraswamy, 2013) . Of course, correlation does not imply causation. As we already stated REM0 has all the hallmarks of high arousal; therefore, a high muscular activity is indeed to be expected during REM0. However, a closer look at the features of EMG, ECG, EOG and KCm events, e.g., as these are exhibited in Fig. 3 , shows higher activity in these channels but no evidence that any one of them is necessarily tightly yoked to any other or to the heartbeat. For KCm events in particular, both the signal-based and the tomographic analyses show neural origin. At the signal level, Figs. 3 and 8 the overall morphology of the wide band signal (dominated by low frequencies) is very similar for KCmWSHF and KCmNOHF, with KCmNOHF clearly away from HRS and with no strong association to high EMG, EEG, MEG or EOG. Also, after entering REM, a decoupling is observed between eye movements, muscular activity, HRS and EEG/MEG graphoelements especially in the low frequency spectrum. HRS in REM differ greatly from HRS during REM0 as they are not regular, they are instead disorganized, leaving no significant artifact to the EMG or EEG electrodes and they do not occur in coincidence with eye movements. The same is true for eye movements as well, there are no KCm events near REM saccades. It is therefore difficult to interpret the KCmWSHF during REM0 as a reflection of other, mainly muscular activity, since no such reflection is seen in REM around REMS which are known to generate one of the strongest EEG and MEG artifacts. The more parsimonious explanation is that the precursor activity is a general arousal wave (marked by an alpha band oscillation in the NBM) that alerts areas of the brain involved in cognition and motivation, possibly reaching the muscles but not enough to generate a huge activation as seen near HRS. The same stimulation that gives rise to increases in brain and muscle activity during REM0 fails to do so during REM (proper) possibly because muscular activity is suppressed by a range of mechanisms leading to atonia which are still poorly understood (Brooks and Peever, 2008) . The increase in EOG and EMG just before HRS suggest that increasing HR has a strong and fast influence on muscular excitability. In this way symmetry in heart-muscle interaction is restored, because fast increase in HR (faster than corresponding blood pressure increase) has been shown to be caused by increase in muscle (Voluntary isometric contraction of triceps Surae) activity (Gladwell and Coote, 2002) . The results presented in Fig. 11 relate to the 2 s before the start of the GLEABS events. As already discussed, the common focus of alpha band increases identified from each and every KC event, KC1, KCmWSHF and KCmNOHF is the NBM in the basal forebrain, the key provider of cholinergic input to the cortex and limbic system. It is therefore necessary to turn our attention to this property of NBM as the main provider of the neurotransmitter Acetylcholine (ACh). ACh is a neurotransmitter involved in both sympathetic and parasympathetic autonomic nervous system and the only one that continues to exert an influence when all others stop during REM. In this effort we are hampered by the scarcity of evidence and the different models used in the detailed studies available. Much of the related work is performed in rats and mice including experiments where optogenetic techniques are used to identify and record from cholinergic neurons in the basal forebrain. Bridging the homology gap is not easy in general, and particularly so in our case because rats and mice do not have alpha activity, at least not the alpha activity as defined in humans. It is therefore difficult to explain the NBM arousing alpha wave we identified in all pre-KC periods for each one of the KC types by drawing information from data obtained from rats, even from data performed separately in humans and rats by the same group (Corsi-Cabrera et al., 2000 , 2001 . The arousal generated by KC1 and KCmNOHF events is not enough to increase muscle excitability (as this is seen in the EMG). To identify the extra push for the increase in muscle tone we used the direct comparison [KCmWSHF vs KCmNOHF] which, for the pre-event period identified increases in activity in the caudal Raphe complex and uvula/cerebellum in the high sigma/low beta band (16 and 17.6 Hz), just a little higher (by 2-3 Hz) than the frequency of posterior mid-parietal increases during spindles during NREM2. From the other areas identified when the criterion was reduced to identifying increases common to three subjects, we note the increases in the delta (sign of inhibition) and high gamma (sign of excitation), at the mPRF the putative Long Lead Initial Pulse (LLIP) for the generation of PGO waves. If the increase in mPRF is indeed the LLIP for PGO generation then this adds another possible cognitive role of KCmWSHF activity: a PGO-initiated memory consolidation process (Gott et al., 2017; Walker and Stickgold, 2004) played out as an elaboration of the corresponding operation by spindles in NREM2. Many of the areas identified in our tomographic analysis are known to be involved in autonomic control, e.g. Ventromedial OrbitoFrontal Cortex (vmOFC), dcACC, Amygdala, medullary nuclei (Green and Paterson, 2020) inferior colliculus (Iigaya et al., 2012) , uvula (Nisimaru, 2004) , ACC and anterior insula (Medford and Critchley, 2010) . This is consistent with the prominent role played by the HRV in our analysis and the demonstrated ability to identify brain correlates of autonomic modulation using HRV in fMRI (Napadow et al., 2008; Park and Thayer, 2014) , including the role of vagal tone in top-down and bottom-up modulations of cognition and attention (Thayer et al., 2012) . The main limitation of our work is the low number of subjects. It is this limitation which makes us cautious to suggest REM0 as a future candidate for a putative new sleep stage. The addition of the second iteration to sleep staging allows REM0 periods to be displayed as an additional label to the conventionally defined sleep stages. We feel this is a fair beginning; it paves the way for the many further studies that need to be completed to characterize REM0 in normal sleep and pathology. At the time the experiment was designed it was a compromise between a study with more subjects but severely limited by the constraints of MEG and the detailed study we did with few subjects and all preparatory activity to ensure acclimatization and making as comfortable sleep environment as possible for MEG, using all polysomnography channels and particularly introducing the first capability for whole night sleep with regular head localization. Getting right the recording of head movement throughout the night proved to be the critical element for our analysis to produce useful results. Even today, the expense and difficulty of performing a whole night sleep MEG study is formidable (Brancaccio et al., 2020) and the additional difficulties introduced by the COVID-19 pandemic make it even harder. At the theoretical level, the strong association we found between HRV and central and motor system activity is very much in line with reports of interaction between distinct networks representing complex physiological systems and particularly the cardiac, respiratory, muscular and central nervous systems in the emerging field of NWP . Sleep offers itself as an excellent natural laboratory for NWP: sleep allows one to study how the dynamics and topology of the overall network changes within and across the different sleep stages which are the exemplars par excellence of primary physiological states. Already pioneering studies have revealed relationships between network topology and physiological function (Bashan et al., 2012; Liu et al., 2015) . It is likely that REM0 will play a particularly strong influence in future NWP studies and that NWP results may eventually force restructuring of sleep staging. The areas identified in the tomographic analysis reported here are targets for varied interventions with deep brain stimulation. The NBM is already a candidate target for reducing the loss of memory and cognitive function in dementia (Kumbhare et al., 2018) while areas related to PGO waves, the Pedunculopontine Nucleus (PPN) and Sub-Thalamic Nucleus (STN), which receive the LLIP from mPRF, are already targeted for Parkinson's Disease (PD) patients (Stein, 2009 ). Other critical applications where our findings may be useful are cases where sleep staging is difficult either because of pathology or because of the conditions, such as situations when brain activity monitoring is needed during highly demanding tasks or missions, for example long space missions which was a prime motivation for carrying out this work. In summary, we demonstrated the non-invasive identification of precursors of GLEABS in the vmOFC and the NBM for each and every KC type, while the direct contrast of KCmWSHF with KCmNOHF identified precursor activity deep in the brainstem. KCm events within REM0, particularly those just before HRS, are candidates for PGO-related activity (Rowe et al., 1999) . The clear differences between REM0 and REM, implies that PGO waves during periods carrying the REM0 label may be generated by different mechanisms and have different roles than the PGOs in REM without the REM0 label. A recent multi-structure recordings in macaque monkeys during conditions of anesthesia as analogues of sleep, provided evidence for two distinct types of PGO waves, corresponding to opposite hippocampal spike-field coupling (Ramirez-Villegas et al., 2020). The current paper was confined to a description of REM0, focusing in particular the tomographic analysis on the pre-KCm period where putative precursors of PGO waves may be found. It is left to future works to search for PGO waves, separately in REM0 and REM periods. Ongoing work in our team focuses on understanding how the precursor activity of KCmWSHF develops into activity in deep brainstem structures during KCmWSHF in two midline bands one including the mPRF and the other including what we provisionally identified as the Raphe Magnus of the lower pons and medulla. During KCmWSHF new activations are also seen in the amygdala and wider amygdala/hippocampus area. Accessing such information non-invasively is significant not only for demonstrating unambiguously the existence of human PGO waves but also for sleep medicine and psychiatry (Arnulf et al., 2010; Gott et al., 2017) . Lead Contact: Further information and reasonable requests for resources should be directed to the Lead Contact, Andreas A. Ioannides (AAI) (a.ioannides@aaiscs.com). Freely available software and algorithms used for analysis are listed in the Key Resource Table. All Custom scripts for signal processing and processed data contained in this manuscript are available upon reasonable request from the Lead Contact. The MFT suite of programs for MEG tomography and post-MFT analysis is based on software developed over the last three decades using originally FORTRAN and then C++ and IDL modules; in its current form it is only available for processing at the AAI Scientific Cultural Services Ltd (AAISCS), but efforts to make the software available to the research community as part of widely used open-source toolboxes are in progress. The original raw data is also only available for analysis at AAISCS because of stipulations under the Material Transfer Agreement between the Brain Science Institute, RIKEN (where the experiments were performed), and AAISCS. Researchers and competent students are welcome to use the software to analyze the anonymized raw data at AAISCS using either the native software of AAISCS or their own software under reasonable requests and arrangements ensuring health and safety procedures prevailing at the time, national and international laws for work practices and agreed procedures for the dissemination of the results. The cost for both the visitor(s) and AAISCS support will be covered by either the institution of the visitors or through joint research grants provided for relevant research and academic collaborations. In rare cases, AAISCS will consider using its own research budget for full or partial funding of proposals submitted by individuals or organizations that demonstrate exceptional promise and naming appropriately qualified individuals committed to doing the work at AAISCS, including work to further develop the MFT analysis and make it available for the wider community as modules of open-source software toolboxes. Sleep induced by stimulation in the human pedunculopontine nucleus area Network physiology: how organ systems dynamically interact Network physiology reveals relations between network topology and physiological function Resting brain fluctuations are intrinsically coupled to visual response dynamics Motivation and cognitive control: from behavior to neural mechanism Cortical source localization of sleep-stage specific oscillatory activity Unraveling the mechanisms of REM sleep atonia Control of sleep and wakefulness Cardiac activity during sleep onset Neuronal oscillations in cortical networks The human K-complex represents an isolated cortical down-state A review on current trends in automatic sleep staging through bio-signal recordings and future challenges EEG bands during wakefulness, slow-wave and paradoxical sleep as a result of principal component analysis in man EEG bands during wakefulness, slow-wave, and paradoxical sleep as a result of principal component analysis in the rat Central cholinergic systems and cognition Evidence of subthalamic PGO-like waves during REM sleep in humans: a deep brain polysomnographic study The time structure of the cyclic alternating pattern during sleep Scalp spindles are associated with widespread intracranial activity with unexpectedly low synchrony REM sleep sawtooth waves are associated with widespread cortical activations Heart rate at the onset of muscle contraction and during passive muscle stretch in humans: a role for mechanoreceptors Towards a functional understanding of PGO waves Using deep brain stimulation to unravel the mysteries of cardiorespiratory control The reward circuit: linking primate anatomy and human imaging The K-complex as a special reactive sleep slow wave -a theoretical update Reproducibility of heart rate variability is parameter and sleep stage dependent REM sleep and dreaming: towards a theory of protoconsciousness Synchronized activation of sympathetic vasomotor, cardiac, and respiratory outputs by neurons in the midbrain colliculi Research domain criteria (RDoC): toward a new classification framework for research on mental disorders Neurofeedback and the neural representation of self: lessons from awake state and sleep MEG identifies dorsal medial brain activations during sleep Continuous probabilistic solutions to the biomagnetic inverse problem MEG tomography of human cortex and brainstem activity in waking and REM sleep saccades Using MEG to understand the progression of light sleep and the emergence and functional roles of spindles and K-complexes The emergence of spindles and Kcomplexes and the role of the dorsal caudal part of the anterior cingulate as the generator of K-complexes Basal forebrain cholinergic modulation of sleep transitions Psyllium husk II: effect on the metabolism of apolipoprotein B in African green monkeys Deep learning-based clustering approaches for bioinformatics The inconsistent nature of heart rate variability during sleep in normal children and adolescents Nucleus basalis of Meynert stimulation for dementia: theoretical and technical considerations The human K-complex: insights from combined scalp-intracranial EEG recordings Cholinergic basal forebrain structures are involved in the mediation of the arousal effect of noradrenaline Resting state functional connectivity of the basal nucleus of Meynert in humans: in comparison to the ventral striatum and the effects of age Generalized Morse wavelets as a superfamily of analytic wavelets Characterization of REM-sleep associated ponto-geniculo-occipital waves in the human pons Major component analysis of dynamic networks of physiologic organ interactions Heart rate variability: standards of measurement, physiological interpretation, and clinical use Intracellular evidence linking medial pontine reticular formation neurons to PGO wave generation Conjoint activity of anterior insular and anterior cingulate cortex: awareness and response Sleep stage detection using only heart rate High-frequency brain activity and muscle artifacts in MEG/EEG: a review and recommendations Brain correlates of autonomic modulation: combining heart rate variability with fMRI Cardiovascular modules in the cerebellum Cholinergic modulation of cortical microcircuits is layer-specific: evidence from rodent, monkey and human brain The process of falling asleep From the heart to the mind: cardiac vagal tone modulates top-down and bottom-up visual perception and attention to emotional stimuli Sleep stage classification from heart-rate variability using long short-term memory neural networks Coupling of hippocampal theta and ripples with pontogeniculooccipital waves A Manual of Standardized Terminology, Techniques and Scoring System for Sleep Stages of Human Subjects Magnetic field tomography of coherent thalamocortical 40-Hz oscillations in humans Heart rate surges during REM sleep are associated with theta rhythm and PGO activity in cats Heart rate spectrum analysis for sleep quality detection A healthy heart is not a metronome: an integrative review of the heart's anatomy and heart rate variability History of the development of sleep medicine in the United States The visual scoring of sleep in adults EEG spectral power in phasic and tonic REM sleep: different patterns in young adults and children The paradox of rapid eye movement sleep in the light of oscillatory activity and cortical synchronization during phasic and tonic microstates The microstructure of REM sleep: why phasic and tonic? BOLD signatures of sleep Akinesia, motor oscillations and the pedunculopontine nucleus in rats and men Sleep staging from electrocardiography and respiration with deep learning Mathematical analysis of lead field expansions The cyclic alternating pattern as a physiologic component of normal NREM sleep A meta-analysis of heart rate variability and neuroimaging studies: implications for heart rate variability as a marker of stress and health Electrical stimulation of the anterior cingulate gyrus induces responses similar to K-complexes in awake humans Sleep-dependent learning and memory consolidation Functional microstates within human REM sleep: first evidence from fMRI of a thalamocortical network specific for phasic REM periods Electroencephalography (EEG) of Human Sleep: Clinical Applications The oscillating brain: complex and reliable FieldTrip: Open Source Software for Advanced Analysis of MEG, EEG, and Invasive Electrophysiological Data Andreas A. Ioannides: conceived initiated and directed the study, for both the experimental phase at the RIKEN BSI and together with Dr. Lichan Liu the analysis in Cyprus. The main MFT analysis was performed, Formal analysis, Dr. Ioannides adapted analysis methods to specific needs of the study, some together with Mr Gregoris A. Orphanides who also implemented them and performed most of the signal based and post-MFT analysis in Cyprus, Writingoriginal draft, All three authors contributed to the writing of the paper. Gregoris A. Orphanides: Formal analysis, Dr. Ioannides adapted analysis methods to specific needs of the study, some together with Mr Gregoris Orphanides who also implemented them and performed most of the signal based and post-MFT analysis in Cyprus, Writingoriginal draft, All three authors contributed to the writing of the paper. Lichan Liu: The main MFT analysis was performed, Writingoriginal draft, All three authors contributed to the writing of the paper. The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. Both authors were employed by the company AAI Scientific Cultural Services Ltd (GO on a part-time basis). The whole night MEG data were recorded at the Laboratory for Human Brain Dynamics (1998-2009), RIKEN Brain Science Institute, Japan. After the closure of the MEG laboratory, the data were anonymized and transferred under a material transfer agreement to the Laboratory for Human Brain Dynamics at AAI Scientific Cultural Services Ltd. (AAISCS) in Nicosia, Cyprus for follow up research and data analysis. Dr. Evangelos Papaefthymiou contributed to theoretical discussions and helped the development of the PCA approach with Mr. Constantinos Kourouyiannis who continuously updated some of the software needed for the analysis. The work reported here was mainly supported by the work for European Space Agency (contract CY2_02 (4000127142) NEA2RS). The clustering and PCA analysis was performed using methods and software developed and partly supported by the European Regional Development Fund and the Republic of Cyprus through the Research & Innovation Foundation (GRATOS EXCELLENCE/1216/0207 and CAVER COMPLEMENTARY/0916/0174). The time frequency analysis benefitted in its early development by work by AAI and Mr. Matthew Fenech during reciprocal visits to Aston University UK and AAISCS as part of exchanges and secondments under the EU Horizon project HOPE (grant number: 823958). The opinions expressed herein belong solely to the authors. The European Space Agency, the European Commission and the Cyprus Research & Innovation Foundation were not involved in the study design, collection, analysis and interpretation of data, the writing of this paper and the decision to submit this paper for publication.