Increased very low frequency pulsations and decreased cardiorespiratory pulsations suggest altered brain clearance in narcolepsy.

Background: Narcolepsy is a chronic neurological disease characterized by daytime sleep attacks, cataplexy, and fragmented sleep. The disease is hypothesized to arise from destruction or dysfunction of hypothalamic hypocretin-producing cells that innervate wake-promoting systems including the ascending arousal network (AAN), which regulates arousal via release of neurotransmitters like noradrenalin. Brain pulsations are thought to drive intracranial cerebrospinal fluid flow linked to brain metabolite transfer that sustains homeostasis. This flow increases in sleep and is suppressed by noradrenalin in the awake state. Here we tested the hypothesis that narcolepsy is associated with altered brain pulsations, and if these pulsations can differentiate narcolepsy type 1 from healthy controls.
Methods: In this case-control study, 23 patients with narcolepsy type 1 (NT1) were imaged with ultrafast fMRI (MREG) along with 23 age- and sex-matched healthy controls (HC). The physiological brain pulsations were quantified as the frequency-wise signal variance. Clinical relevance of the pulsations was investigated with correlation and receiving operating characteristic analysis.
Results: We find that variance and fractional variance in the very low frequency (MREGvlf) band are greater in NT1 compared to HC, while cardiac (MREGcard) and respiratory band variances are lower. Interestingly, these pulsations differences are prominent in the AAN region. We further find that fractional variance in MREGvlf shows promise as an effective bi-classification metric (AUC = 81.4%/78.5%), and that disease severity measured with narcolepsy severity score correlates with MREGcard variance (R = -0.48, p = 0.0249). Conclusions: We suggest that our novel results reflect impaired CSF dynamics that may be linked to altered glymphatic circulation in narcolepsy type 1.

Introduction

Narcolepsy is a debilitating, chronic neurological disease characterized by severe daytime sleepiness, sleep attacks, cataplexy, and fragmented nocturnal sleep[1-3]. The diagnosis of narcolepsy relies on symptoms common to other diseases and tests requiring hospital/sleep center resources including mean sleep latency test, polysomnography, and invasive cerebrospinal fluid (CSF) sampling for analysis of the hypocretin level. Up to 60% of patients with narcolepsy are initially misdiagnosed with, for example, other hypersomnias or depression, and, conversely, up to 46% of initial narcolepsy diagnoses are later found to be misdiagnoses[4,5]. Thus, a correct narcolepsy diagnosis is often delayed for years, which can lead to inappropriate treatment, thus compounding the disease burden in patients with narcolepsy or conditions mistaken for narcolepsy[6,7].

Narcolepsy type 1 is hypothesized to arise from specific degeneration or dysfunction of the hypocretin-producing cells in the hypothalamus due to an autoimmune reaction that leads to decreased or absent hypocretin signaling in the ascending arousal network (AAN), hypothalamic tuberomammillary nucleus, and throughout the neocortex[2,8]. Hypocretin 1 and 2 (also known as orexin A and B) are 33 and 28 residue peptide neurotransmitters, which have a direct excitatory effect on cortical and brainstem neurons. The AAN of the brainstem comprises neurons of the midbrain reticular formation (MRF), locus coeruleus (LC: noradrenergic), ventral tegmental area and periaqueductal gray matter (VTA, PAG: dopaminergic), and the raphé nuclei (DR: serotonergic), all of which give rise to ascending projections. These AAN innervations participate in the control of autonomic functions and arousal/sleep via their respective neurotransmitters in dense innervations to the hypothalamus, thalamus, basal ganglia, and neocortex[9]. Earlier imaging studies in narcolepsy have shown structural alterations in cortical and subcortical brain regions including the limbic system and brainstem along its nuclei e.g., LC and neocortex[10,11]. Research conducted with task-fMRI have found altered activation of the amygdala, a part of the limbic system, in tasks with humor/reward paradigms along with increased deactivation of the default mode network (DMN) under cognitive burden in narcolepsy[12-14]. Furthermore, a study by Drissi et al. suggests instability of the DMN even in the resting condition[15]. In line with these, in our earlier study utilizing dynamic lag analysis of fast fMRI data at rest, we found the signal propagation as delayed and monotonic, especially between the DMN and other major networks in narcolepsy type 1[16]. Taken together, it is postulated that the degeneration/dysfunction of hypocretin neurons in narcolepsy type 1 manifests in impaired functioning of the AAN, resulting in arousal/sleep disturbances and imbalance in autonomic control of physiology[3,17,18].

The convection of CSF is driven by physiological brain pulsations, especially cardiorespiratory-related pulsations and vasomotor waves[19]. In a current model of brain fluid homeostasis, CSF water enters the brain parenchyma from arterial perivascular spaces driven by cardiovascular pulsations and exits the brain via paravenous spaces leading to peripheral lymphatic outlets[20,21]. Furthermore, respiration-evoked pressure changes in the thoracic cavity promote venous blood outflow from the brain and reciprocal movement of the CSF from the spinal canal towards the brain, thus facilitating CSF dynamics[22,23]. The CSF flow through the brain parenchyma is an essential aspect of physiology, as the passage of metabolites and immune cells, intracranial fluid dynamics, and brain homeostasis all depend on interstitial CSF flow, which has recently been termed the glymphatic system[24].

The glymphatic system is most active during sleep and is suppressed by noradrenaline transmission during wakefulness[25]. Moreover, hypertension stiffens arterial walls, thus reducing perivascular pumping and impeding the flow of CSF in perivascular spaces and glymphatic clearance[26]. Furthermore, the reduced noradrenergic tonus prevailing during normal sleep facilitates CSF flow increases in brain interstitium and paravascular spaces[25,27,28]. These findings emphasize the importance of hypocretin, AAN brainstem nuclei, and nocturnal sleep for maintaining brain CSF homeostasis. In the present context, narcolepsy type 1 is characterized by daytime sleep attacks, fragmented nocturnal sleep architecture, a high risk of hypertension, and hypothesized inconsistent noradrenaline transmission in the AAN[1,3,29]. Despite these associations, no previous research has considered whether narcolepsy type 1 impacts brain pulsations and related CSF dynamics.

Only recently it has become possible to non-invasively measure physiological brain pulsations driving CSF and blood flow in the human brain using analysis of fast functional magnetic resonance (fMRI) signal dynamics. Such techniques have revealed marked impairments of the pulsations in patients with Alzheimer’s disease and in epilepsy[30-33]. From a clinical perspective, studies have shown that fMRI signal variability is associated with cognitive performance, increases in dementia, and declines with normal ageing[32,34-36].

In the present study, we undertake the first investigation of fast fMRI signal variability in narcolepsy type 1 and discuss its relation to CSF dynamics. We use an ultrafast fMRI sequence (MREG: magnetic resonance encephalography) to investigate whether vasomotion-related very low frequency (MREGvlf) or individual cardiac (MREGcard) and respiratory (MREGresp) frequency-related brain pulsations are disrupted in narcolepsy type 1. This concept is enabled by the 10 Hz temporal resolution of the modern MREG procedure, which provides temporal signals that are free of aliasing and slice time inaccuracies, thus enabling accurate individual separation of physiological signal instances[37-39]. We further extend our analysis to examine the bi-classification accuracy of these pulsations according to receiving operating characteristic (ROC) and test the clinical association by correlation analysis between brain pulsations and disease severity measured with narcolepsy severity score (NSS)[40,41]. Our a priori hypotheses are that (1) brain pulsations driving CSF flow are disrupted in narcolepsy type 1, (2) brain pulsation measurements can be used to differentiate by non-invasive examination of patients with narcolepsy type 1 and healthy controls, and (3) the extent of pulsation changes in narcolepsy type 1 are associated with disease severity.

Here we find that all three hypotheses are met with confirmatory findings. The very low-frequency pulsations are increased while cardiorespiratory-related pulsations are decreased in narcolepsy type 1. Cardiac-related pulsations in the AAN region negatively correlate with increasing disease severity while very low-frequency pulsations show promise as a bi-classificator. These results imply the clinical relevance of our findings. Taken together, we suggest that our results reflect impaired CSF flow in the brain that may be linked to altered glymphatic function in narcolepsy type 1.

Methods

Participants

A registry run from the Oulu University Hospital’s electronic patient records for patients diagnosed with all-type narcolepsy was conducted, resulting in 66 matching diagnoses. All the diagnoses were reassessed with the International Classification of Sleep Disorders (ICSD) 3rd edition[42] to exclude outdated diagnoses and to confirm narcolepsy type 1 diagnoses according to the latest diagnostic criteria. Twenty-three narcolepsy type 1 patients were recruited for the study. Inclusion criteria were (1) confirmed narcolepsy type 1 diagnosis and (2) clear-cut cataplexy (short, usually bilateral loss of muscle tone without loss of consciousness commonly triggered by emotional stimuli). Other brain-related confounding conditions were excluded by a screening of clinical history and examination of structural MRI brain images. All data were collected between 3/2018 and 9/2020.

Data from one patient was excluded due to corrupted MREG data resulting from the failure of the off-resonance correction during imaging. The blood pressure data from one patient and two controls were missing. The final population consisted of 22 patients with narcolepsy type 1 (NT1 group, 12 females, of mean age of 28.1 ± 8.9 standard deviation (SD)). Four of the NT1 patients were unmedicated and 18 were medically treated for daytime sleepiness and cataplexy (Table 1). All subjects in the NT1 group filled out a disease severity form (NSS: narcolepsy severity scale[40,41]), which ranges from 0 to 57 points. Our patients’ NSS scores ranged from 8 to 45 (mean 26.3 ± 10.5 SD). Twenty-two healthy individuals age- (±3 years) and sex-matched controls without continuous medication were recruited by advertisement from the general population as a control group (HC group, 12 females, of mean age 28.2 ± 8.9 SD). Written informed consent was obtained from all participants. This study was approved by the Ethical Committee of Medical Research in the Northern Ostrobothnia District of Finland and was conducted in accordance with the declaration of Helsinki and GDPR regulations.

Table 1

Disease severity by narcolepsy severity score and medication information.

Patients with narcolepsy type 1
ID	NSS	Medication
1	36	Mo, SSRI
2	25	Mo, Me
3	18	-
4	18	Mo
5	11	Mo
6	13	S
7	45	Me
8	31	S, SNRI
9	28	Me, S
10	31	-
11	33	-
12	22	Mo
13	37	Mo, Me
14	8	Me
15	33	Me, S
16	26	Mo, SNRI
17	20	Mo
18	45	Mo, S, SNRI
19	14	Me
20	19¤	Me
21	32	-
22	39	Mo
23	22	Me, S

ID subject, NSS narcolepsy severity scale score, Mo modafinil, Me methylphenidate, S sodium oxybate, SSRI selective serotonin reuptake inhibitor, SNRI serotonin–norepinephrine reuptake inhibitor, ¤ excluded for corrupted magnetic resonance encephalography data.

Data acquisition

NT1 and HC participants were scanned with the fast fMRI sequence MREG using a Siemens Magnetom Skyra 3 T MRI scanner (Siemens Healthineers, Germany) with a 32-channel head coil. MREG is a single-shot three-dimensional (3D) sequence that uses a spherical stack of spirals and undersamples the 3D k-space trajectory[37,43,44]. The following parameters were used for the 3D whole brain MREG sequence: repetition time (TR) = 100 msec, echo time (TE) = 36 msec, flip angle (FA) = 25°, field of view (FOV) = 192 mm, voxel size = 3 × 3 × 3 mm3. MREG data were reconstructed by L2‐Tikhonov regularization with lambda = 0.1, with the latter regularization parameter determined by the L‐curve method[45], giving an effective isotropic spatial resolution of 4.5 mm. MREG includes a dynamic off‐resonance in the k‐space method, which corrects for respiration-induced dynamic field‐map changes in fMRI using the 3D single-shot technique[46]. For registration purposes, the participants also had a T1-weighted Magnetization Prepared Rapid Acquisition with Gradient Echo (MPRAGE) scan with parameters as follows: TR = 1900 msec, TE = 2.49 msec, inversion time (TI) = 900 msec, FA 9°, FOV = 240, and slice thickness 0.9 mm). Simultaneously with MREG scanning, end-tidal CO2 monitoring and photoplethysmogram were recorded to measure each subject’s respiration and heart rate, respectively.

The subjects were instructed to lie still and awake with eyes open and fixating on a cross on the screen for the entire 10-min resting-state scanning. To minimize head motion and to reduce the effects of scanner noise, soft pads were fitted over the study subjects’ ears, along with earplugs.

Preprocessing

A preprocessing pipeline from Järvelä et al.[16] with small changes was implemented here: MREG data were preprocessed with Oxford Centre for Functional MRI of the Brain (FMRIB) software library (FSL, version 5.09) pipeline[47]. The data were high-pass filtered with a cut-off frequency of 0.008 Hz (125 s). To minimize T1-relaxation effects, 180 time points were ignored from the beginning of the data, resulting in 5822 whole brain volumes. Motion correction was carried out using FSL MCFLIRT[47], and all data and motion parameters were visually inspected for spurious signal fluctuations. The motion was further controlled by the exclusion of any motion exceeding the voxel size (no subject had mean relative motion over 0.07 mm or mean absolute motion over 0.6 mm) and by calculating each subject’s mean frame-wise displacement. Brain extraction for 3D MPRAGE volumes was performed with the FSL Brain Extraction Tool using neck and bias-field correction and the following parameters: fractional intensity = 0.20–0.22 and threshold gradient = 0.05–0.25. To obtain optimal quality, the extracted brain images were visually inspected. Images were spatially smoothed with a 5 mm full width and half maximum (FWHM) Gaussian kernel using fslmaths. MREG images were aligned to the 3D anatomical images (full‐search, 12 degrees of freedom (DOF)) and to the Montreal Neurological Institute (MNI152) 4 mm standard space (full‐search, 12 DOF) as a preprocessing step in the FSL multivariate exploratory linear optimized decomposition into independent components tool.

Brain pulsation range estimation

Individual respiratory and cardiac frequencies were extracted from end-tidal CO2 and photoplethysmogram data with the MATLAB version R2019b fast Fourier transformation (FFT) fft function. From the signal frequency spectrum, respiratory and cardiac minimum, maximum, and peak values were obtained (Fig. 1). Minimum values were subtracted from the maximum values to calculate individual cardiorespiratory frequency ranges.

Fig. 1

Analysis pipeline.

Temporal magnetic resonance encephalography (MREG), respiratory (end-tidal CO2), and cardiac (photoplethysmogram) signals are transformed into frequency spectra with fast Fourier transformation (FFT). Individual minimum, maximum, and peak values for respiration (blue) and cardiac (red) frequencies are obtained and the cardiorespiratory ranges are calculated. MREG full band signal is filtered to these physiological ranges and to very low frequency (green). Voxel-wise variance maps are calculated for each MREG frequency band. MREGvlf very low-frequency filtered MREG, MREGresp respiratory frequency filtered MREG, MREGcard cardiac frequency filtered MREG.

In cases with corrupted physiological data (5 end-tidal CO2 and 11 photoplethysmogram data recordings), the MREG data were used to estimate physiological frequency ranges. First, Analysis of Functional NeuroImages’ (AFNI, version 18.0.05)[48] 3dPeriodogram was used to transform the preprocessed four-dimensional MREG data into a voxel-wise FFT spectrum. For respiratory and cardiac frequency estimation, the fourth ventricle and anterior/middle cerebral arteries, respectively, were chosen as reference points, as the investigated physiological events in the MREG FFT spectrum were visually most pronounced in these regions. The corresponding minimum, maximum, and peak values for cardiorespiratory frequencies were then estimated. To test the precision of the estimation of the MREG FFT spectrum, minimum, maximum, and peak frequencies were extracted from 41 available end-tidal CO2 and 35 photoplethysmogram recordings and correlated with their MREG FFT spectrum counterparts using the R software ggscatter function (Supplementary Fig. 1). The MREG FFT spectrum estimations of cardiorespiratory frequencies were found to be accurate when compared with the gold standard estimators, and were thus used in further analyzes in those cases where the peripheral physiological data were corrupted or missing. While estimating respiratory ranges across groups, the minimum value of respiration was found to be lower than 0.1 Hz in two NT1 subjects. To avoid any confounding effects of respiration on very low-frequency variance calculation, the upper limit of the very low-frequency band was set to 0.08 Hz, paralleling previous research[49-51].

Variance analysis

The general workflow is described in Fig. 1. To investigate individual physiological pulsation spectrum, AFNI 3dTproject was used to bandpass-filter all the preprocessed data sets to very low frequency (0.01–0.08 Hz), and individual respiratory and cardiac frequencies. Signal variance is a measure of variability and probability distribution, defined as the expectation of the squared deviation of a random variable from its mean:where χi is the variable (signal amplitude value at a given time), µ is the signal mean, and N is the number of time points.

Fractional variance (Varfrac) is a measure used, for example, in principal component analysis (indicating the percentage of variance explained by a component)[52], which is defined as the variance of a random variable divided by the total variance of the measurement:where Varfrac is fractional variance, Vari is partial signal variance, and Vart is total signal variance. The variance was chosen as the measure of signal variability rather than the standard deviation and coefficient of variation, as used in the previous research[31,32,53] as these latter measures cannot be used to calculate proportional variability of different frequency bands.

Brain voxel-wise variance was calculated with FSL fslmaths for full band (0.008–5 Hz) time series and all physiological pulsation frequencies. Then, all voxel-wise variance data were registered to the MNI152 3 mm standard space with FSL flirt[54] and masked with an MNI152 3 mm standard brain binary mask to remove any residual MREG voxels outside the co-registered standard space. The Varfrac values were obtained by dividing the very low frequency, respiratory, and cardiac frequency variances with the full band variance and then applying the same processing steps described above to acquire brain vowel-wise maps in standard space. Finally, variance and Varfrac maps were compared between HC and NT1 groups with 10,000 iterations of FSL randomize to extract the corrected p value t-statistic maps (p < 0.05)[55]. The results were then displayed on top of the MRIcroGL MNI152 standard brain.

Correlation and receiving operating characteristic analyses

AAN in 1 mm MNI152 standard space mapped by Edlow et al.[8] was used as a hypothesis-driven ROI. Variance results were first registered to 1 mm MNI152 standard space. A whole AAN ROI spanning 2042 voxels/mm3 was generated by creating a binary mask of individual AAN nuclei (for example, LC, MRF, VTA as defined by Edlow et al.) with fslmaths, and then applied to the variance maps. Individual nuclei remain small but at the group level, MREG has been shown to be capable of delineating minimal respiratory centers even in the relatively pulsatile brainstem[56]. The correlation between disease severity (NSS score) and mean variances from the AAN ROI were calculated with the RStudio version 1.3.1093 software ggplot2 library commands. ROC curves were plotted with the R software pROC library commands, and also applied to individual brainstem areas comprising the AAN. The final editing of the figures was done with GNU Image Manipulation Program version 2.10.30.

Statistics and reproducibility

All analyses were conducted between the 22 subjects in the NT1 group and the 22 subjects in the HC group except between the blood pressure measurements where the available 21 NT1 and 20 HC were used for the comparisons, and when correlating end-tidal CO2 and photoplethysmogram derived cardiorespiratory frequencies to the corresponding MREG data (41 end-tidal CO2 vs. MREG pairs and 35 photoplethysmogram vs. MREG pairs). The normality of the data distributions was estimated visually and with Shapiro–Wilk test. Two-tailed Student’s t-test and Wilcoxon rank-sum test were used for hypothesis testing between the study groups in cardiorespiratory frequency estimation (minimum, maximum, peak values, and range), motion, blood pressure (significant p value < 0.05), and individual nuclei variance results that were further corrected for multiple comparisons with the Benjamini–Hochberg method (significant p value < 0.033 or 0.016, please see Supplementary Tables 1, 2). Wilcoxon and Cohen’s d effect sizes were calculated with Rstudio’s wilcoxon_effsize and cohens_d functions. Pearson correlation coefficient was used to determine association between NSS and variance results (significant p value < 0.025). The data describing the study population is available in Supplementary Tables 1, 2. Randomize uses conditional Monte Carlo random permutations implementing family-wise error-corrected threshold-free cluster enhancement correction in both directions (HC > NT1, HC < NT1) separately[55], thus taking into account the multiple comparisons, and was used to test for differences between the study groups in variance and fractional variance maps (significant p value <0.05).