Increased interictal synchronicity of respiratory related brain pulsations in epilepsy.

Respiratory brain pulsations have recently been shown to drive electrophysiological brain activity in patients with epilepsy. Furthermore, functional neuroimaging indicates that respiratory brain pulsations have increased variability and amplitude in patients with epilepsy compared to healthy individuals. To determine whether the respiratory drive is altered in epilepsy, we compared respiratory brain pulsation synchronicity between healthy controls and patients. Whole brain fast functional magnetic resonance imaging was performed on 40 medicated patients with focal epilepsy, 20 drug-naïve patients and 102 healthy controls. Cerebrospinal fluid associated respiratory pulsations were used to generate individual whole brain respiratory synchronization maps, which were compared between groups. Finally, we analyzed the seizure frequency effect and diagnostic accuracy of the respiratory synchronization defect in epilepsy. Respiratory brain pulsations related to the verified fourth ventricle pulsations were significantly more synchronous in patients in frontal, periventricular and mid-temporal regions, while the seizure frequency correlated positively with synchronicity. The respiratory brain synchronicity had a good diagnostic accuracy (ROCAUC = 0.75) in discriminating controls from medicated patients. The elevated respiratory brain synchronicity in focal epilepsy suggests altered physiological effect of cerebrospinal fluid pulsations possibly linked to regional brain water dynamics involved with interictal brain physiology.

Introduction

Epilepsy is the second most onerous neurological disorder, with a lifetime prevalence of about 7.6 per 1000. One third of patients suffer from intractable epilepsy, a debilitating form of the disease, that is resistant to any anti-epileptic medication. Spontaneous unprovoked seizures are the hallmark of epilepsy, but the underlying etiology is highly variable, including genetic conditions or structural abnormalities, such as cortical developmental malformation or gliotic changes occurring after a brain injury. While it is understood that seizures most likely result from an imbalance of excitatory and inhibitory neurotransmission or changes in potassium metabolism leading to net and hyperexcitability, physiological mechanisms that contribute to seizure generation and the factors leading to chronic epilepsy are incompletely understood.

Functional magnetic resonance imaging (fMRI) is a widely used non-invasive research and diagnostic tool, which has also been used to reveal hemodynamic changes in patients with epilepsy. fMRI studies are usually recorded during the inter-ictal period and therefore do not represent acute changes during seizures but rather reflect the more subtle changes prevailing during interictal epileptic discharges. Nevertheless, many studies suggest that hemodynamic changes in epilepsy involve widespread cortical and subcortical brain structures extending beyond the area generating epileptic seizures.

Recent advances have demonstrated the usefulness of fast fMRI scans in epilepsy diagnostics.[9,10] While these studies focus on the blood oxygen level dependent (BOLD) response to epileptic spikes, fast fMRI imaging can also serve to measure physiological brain pulsations (vasomotor, respiratory, and cardiovascular pulsations), which may strongly affect the BOLD signal.[11-13] The latter approach has furnished a better understanding of BOLD variability, but measurements of physiological brain pulsations, might also shed light on the pathomechanisms contributing to epilepsy and other neurological diseases. Indeed, physiological pulsations contain information about different cognitive, vigilance and pathophysiological conditions,[14-17] as they are ubiquitous in brain tissue, blood and cerebrospinal fluid (CSF).[18-22]

Epileptic seizure onset, i.e. ictogenesis, is characterized as an imbalance between the inhibitory and the excitatory neurotransmission. However, according to a novel theory, the glymphatic brain water dynamics can contribute to seizure generation, especially by creating an interstitial environment that facilitates the processes sustaining chronic epilepsy. Alterations in glymphatic convection influence factors like interstitial electrolyte concentrations, inflammation response, and medication delivery,[14,23,24] all of which are important mediators influencing the efficacy of pharmacotherapy for epilepsy. Physiological brain pulsations, especially respiration, are the main drivers of CSF flow,[18-20] and thus are important facilitators of glymphatic brain convection and thus influence epileptic activity.

A recent study has shown correlations between CSF oscillation, and electroencephalography (EEG) activity. Indeed, fast fMRI studies have revealed elevated respiratory brain pulsations in patients with epilepsy.[26,27] Intracranial recordings of intractable epilepsy patients indicate that respiratory brain pulsations drive the local field potentials of neurons,[28,29] which may influence the threshold for spontaneous seizures. Moreover, a recent epilepsy study showed changes in function of neuron-glia structures, the glia limitans, which are aligned spatially with astrocyte endfeet rich in Aquaporin-4 (AQP4) water channels affecting CSF flow in the brain. Taken together these studies suggest altered respiratory drive of glymphatic CSF convection in patients with epilepsy.

In this study, we aimed to characterize the role of respiratory drive in spontaneous brain activity in epilepsy by quantifying differences in respiratory brain synchronization (Respsync) between patients and controls. We used fast fMRI recordings to test the hypothesis that Respsync is increased in epilepsy, and further investigated the relationship of Respsync to disease severity and duration, and tested its diagnostic accuracy.

Methods

Participants

The study included 40 patients with epilepsy (PWE, age 34.0 ± 10.5 years, 25 females), 20 drug-naïve patients with probable focal epileptic seizures (DN, age 36.0 ± 15.3 years, 5 females) and 102 healthy controls (HC, age 37.3 ± 15.4 years, 51 females). We analyzed the PWE (Table 1) and DN groups (Table 2) separately. The PWE group subjects were treated with 1–4 antiepileptic drugs (AEDs, Table 1). The majority of the PWE group (31/40) had intractable epilepsy with continuing seizure occurrences despite their medication. All study participants were imaged between 2012 and 2020. Nine of 20 DN had a confirmed epilepsy diagnosis (Table 2). Seventeen PWE (10 females) and fourteen HC (6 females) were recruited and imaged in the University Medical Center Freiburg. The remainder of the PWE (n = 23), DN (n = 20) and HC (n = 88) groups were recruited from the outpatient clinic at Oulu University Hospital or from different ongoing studies at the University of Oulu.

Table 1.

Clinical characteristics of recruited patients with epilepsy (PWE).

PWE	Age	Sex	Duration from diagnosis	Seizure frequency	Diagnosis	MRI finding	AED
1	21	M	3 years	4	TLE R	MTS R	LEV, LCM
2	32	F	13 years	5	TLE R	Normal	LEV, LTG
3	23	F	18 years	3	TLE L	Normal	LCM, LEV
4	38	M	36 years	4	TLE L	FCD L T	OXC, LEV, ZNS
5	23	F	18 years	5	TLE R	FCD R insula	LTG
6	37	M	10 years	4	TLE bilateral	MTS bilateral	LTG, CBZ, PGB
7	45	F	14 years	4	TLE bilateral	EC T bilateral	OXC, LTG, PB
8	34	M	4 years	4	TLE R	Normal	TPR
9	22	F	6 years	4	TLE bilateral	MTS bilateral	BRV, LCM
10	59	M	38 years	4	TLE L	Cavernoma,FCD L T	LCM, CBZ
11	20	F	5 years	5	TLE R	Glioma R T	GBP, LCM
12	35	F	19 years	4	TLE R	Normal	LEV, PER
13	33	F	28 years	4	TLE bilateral	NF1 suspected	OXC, BRV
14	53	F	26 years	4	TLE R	Normal	OXC, LEV
15	29	M	9 years	4	TLE L	MTS L	LCM, BRV
16	43	F	9 years	5	TLE bilateral	EC T bilateral	ZNS, BRV
17	27	M	26 years	5	FLE R	FCD F R	LTG, LEV, ZNS, PER
18	41	M	3 months	2	TLE L	Normal	OXC
19	40	M	4 months	2	TLE	Normal	OXC
20	32	M	8 years	4	TLE R	Hippocampal edema R	CBZ, PGB
21	38	F	10 years	3	TLE R	HS R	LTG, TPR
22	26	F	19 years	4	FLE	Normal	CLB, LCM, LTG, VPA
23	52	F	4 years	4	TLE L	FCD L T	LCM, ZNS
24	36	M	6 years	3	TLE L	Arachnoid cyst	LEV, OXC, LCM
25	53	M	13 years	3	TLE R	DVA / Telangiectasia R	LTG, LCM, TPR
26	22	F	12 years	3	TLE R	FCD R T	VPA, LTG, CLB
27	25	F	11 years	3	TLE L	FCD L	OXC, ZNS
28	35	F	12 years	5	TLE R	Normal	TPR, LCM, CLB
29	26	F	7 years	3	TLE R	Arachnoid cyst RPituitary microadenoma RFCD R	LEV, LCM, LTG
30	39	F	16 years	3	TLE	Normal	LTG
31	35	F	17 years	3	TLE	Normal	PGB, ZNS
32	26	F	4 years	5	TLE R	Normal	LCM, TPR, BRV
33	43	F	4 years	2	TLE R	Normal	LTG
34	43	F	33 years	4	TLE L	FCD L, Small aneurysm L	CBZ, CLB, LCM
35	17	F	4 years	3	TLE	Normal	OXC, LEV
36	26	M	6 years	4	TLE L	Heterotopia L T	OXC, LEV
37	20	M	2 years	2	TLE R	Normal	OXC
38	49	F	38 years	1	Focal epilepsy	Normal	VPA
39	26	M	11 years	1	TLE R	Normal	OXC, LEV
40	36	F	1 month	1	TLE	Normal	OXC

Seizure frequency grades; Seizure within 1: a year, 2: two months, 3: a month, 4: a week and 5: a day.

AED: antiepileptic drug; M: male; F: female; FLE: frontal lobe epilepsy; TLE: temporal lobe epilepsy; T: temporal; L: left; R: right; MTS: mesio-temporal sclerosis; FCD: focal cortical dysplasia; EC: encephalocele; NF1: neurofibromatosis type 1; HS: hippocampal sclerosis; DVA: developmental venous anomaly; BRV: brivaracetam; CBZ: carbamazepine; CLB: clonazepam; GBP: gabapentin; LCM: lacosamide; LEV: levetiracetam; LTG: lamotrigine; OXC: oxcarbazepine; TPR: topiramate; PB: phenobarbital; PER: perampanel; PGB: pregabalin; VPA: valproic acid; ZNS: zonisamide.

Table 2.

Drug-naïve patients (DN) characteristics.

DN	Age	Sex	Duration of symptoms	Seizure frequency	Diagnosis	MRI finding
1	16	M	4 years	4	TLE L	Normal
2	41	M	10 years	2	TLE L	Normal
3	39	M	1 year	2	TLE	Normal
4	17	M	One seizure	1	No diagnosis	Normal
5	30	M	3 months	1	Focal epilepsy	Normal
6	60	F	One seizure	1	No diagnosis	DNET suspected
7	49	M	2 years	1	TLE	Temporal vascular degeneration L
8	49	F	One seizure	1	No diagnosis	Normal
9	62	M	2 months	1	Focal epilepsy	Normal
10	64	M	2 months	1	TLE	FCD R
11	23	F	No information	1	Focal epilepsy	FCD L
12	22	M	One seizure	1	No diagnosis	Normal
13	57	M	One seizure	1	No diagnosis	Normal
14	25	M	One seizure	1	No diagnosis	Normal
15	25	M	One seizure	1	No diagnosis	Normal
16	53	F	One seizure	1	No diagnosis	Normal
17	37	M	One seizure	1	No diagnosis	MTS L suspected
18	24	M	One seizure	1	No diagnosis	Normal
19	28	M	One seizure	1	No diagnosis	Normal
20	21	F	2 months	2	Focal epilepsy	Normal

Seizure density grades; Seizure within 1: a year, 2: two months, 3: a month, 4: a week and 5: a day.

M: male; F: female; TLE: temporal lobe epilepsy; L: left; DNET: dysembryoplastic neuroepithelial tumor.

The PWE group received a standardized multidisciplinary evaluation including high-resolution structural T1-weighted MRI and T2-FLAIR MRI, electroencephalography and neurological assessment for clinical diagnostics. Additionally, patients diagnosed with intractable epilepsy underwent a neuropsychiatric and neuropsychological evaluation. All patients from this group had diagnosis of focal epilepsy of various etiologies, and had started the treatment with AEDs. DN patients had no earlier diagnosis of epilepsy or other neurological disorders, and were imaged after having had a suspected focal onset epileptic seizure, but had not started AED. The HC group had no history of diagnosed neurological disorders nor other chronic medical conditions or medications. We did not have any IQ data available for our study subjects.

Written informed consent was obtained from all participants according to the Declaration of Helsinki. The research protocol was approved by Ethics Committee of the Northern Ostrobothnia Hospital District, Finland and Research Ethics Committee of the University Medical Center Freiburg, Germany.

Data acquisition

Whole brain fMRI acquisition were performed using a 3D single-shot sequence (Magnetic Resonance Encephalography, MREG) that under-samples k-space to reach a sampling rate of 10 Hz (images per second). The MREG gathers k-space in spherical stack of repeated clock-/anticlockwise spirals, which minimizes air-sinus off-resonance effects. The raw MREG data were reconstructed using MATLAB recon tool with L2-Tikhonov regularization parameter (lambda = 0.1) and the latter regularization parameter determined by the L-curve method resulting in an effective spatial resolution of 4.5 mm anisotropic.

Imaging was undertaken in Oulu University Hospital (Finland) and University Medical Center Freiburg (Germany) with parameters as listed in Table 3. We note that in Freiburg the neck channels from the receive coil were turned off and a smaller slab thickness was used, resulting in exclusion of the inferior parts of the brainstem in nine subjects of whom four were HCs and five in the PWE group.

Table 3.

Scanning parameters at the different sites.

	Oulu	Freiburg
Scanner	Siemens 3 T SKYRA	Siemens 3 T PRISMA
Coil	32-channel head coil	64-channel head-neck coil
Specifications for fMRI imaging (MREG)
TR (repetition time)	100 ms	100 ms
TE (echo time)	36 ms	36 ms
Flip angle	25°	25°
3D matrix	64 × 64 × 64	64 × 64 × 50
Voxel size (isotropic)	3 × 3 × 3 mm3	3 × 3 × 3 mm3
Slab thickness	192 mm	150 mm
Specifications for anatomical reference T1
TR (repetition time)	1900 ms	2000 ms
TE (echo time)	2.49 ms	4.11 ms
TI (inversion time)	900 ms	1000 ms
Flip angle	9°	12°
Field of view	240 mm	256 mm
Slice thickness	0.9 mm	1 mm

The data acquisition lasted up to 20 min. Due to scheduling reasons and/or subject compliance, some acquisitions had to be stopped after 5 min. Therefore, to enable a consistent analysis, only the first 5 min from every dataset was used, which had the benefit of minimizing effects from any potential decline in vigilance of the subjects.[26,31] Study subjects wore earplugs to reduce scanner noise and soft pads were fitted over their ears to further reduce noise and to minimize head motion. Prior to imaging, all participants received instructions to stay awake with their eyes open and fixating at a cross on the screen. No breathing schemes were used in this study and the study subjects were allowed to breath normally ad lib.

Data preprocessing

The preprocessing of the MREG data was done with a standard FSL (Functional Magnetic Resonance Imaging of the Brain’s software library) pipeline, following the procedures and steps described in our previous papers,[26,31] using high-pass filtering with a cut-off frequency of 0.008 Hz (125 s). Motion correction was performed using FSL MCFLIRT, and FSL BET was used for brain extraction. Anatomical T1 images were used for anatomical reference and co-registration with specifications listed in Table 3.

Temporal signal-to-noise ratio (tSNR) is an established way to evaluate fMRI signal, and the global tSNR values between groups were compared. Additionally, it has been shown that MREG only results in a 18.4% reduction in tSNR, which is countervailed by high temporal resolution that increases statistical power. For secondary artifact removal, the FIX ICA (Functional Magnetic Resonance Imaging of the Brain’s independent component analysis -based Xnoiseifier) -method for preprocessed MREG data was used.[32,33] The used FIX classifier was trained on previously collected MREG data from healthy subjects scanned with identical parameters and was used identically for all subject groups.[26,27]

Respiratory synchronicity map analysis

During the respiratory cycle, the CSF has a net in and out of the brain depending on the prevailing respiratory phase as dictated by the Monro-Kellie doctrine (Figure 1(a)).[18,20] This CSF flow effect oscillates reciprocally with venous blood volume, thus pumping the perivenous CSF convection that facilitates fluid solute clearance (Figure 1(b)).

Figure 1.

(a) Cerebrospinal fluid (CSF) flow modulation by respiratory intrathoracic pressure changes (blue arrow). During inhalation the venous return from the brain increases (long white arrow) and CSF flows (yellow) towards the brain to counterbalance the blood volume loss. During exhalation, CSF returns from the brain and venous return declines, which gradually inflates cortical veins (short white arrow). (b) CSF convection (yellow) in the brain tissue is driven by arterial (red) pulsations (cardiac) into periarterial influx channels (light blue surrounding red) and into the interstitium via AQP4 water channels (pink). In the veins (blue) the counterbalancing CSF/venous pulsations (inhalation/exhalation) may move water into perivenous efflux channels (light blue surrounding blue). (c) The mean correlation map of respiratory BELT vs. MREG data (representative signals and raw image data shown) from 30 healthy subjects determined in a region of interest (white area) placed on the fourth ventricle, yielding MREG data that tracks the respiratory phase with high precision. (d) 30 healthy subject's respiration rate in breaths per minute (bpm) determined from respiratory belt and MREG indicates a high correlation coefficient of 0.97 (Pearson, p < 0.0001). Additionally, the differences between the groups (PWE, DN, HC) respiratory frequencies were compared and (e) The individual examples of extracted MREG respiratory signals (mean signal of white area from c) were used for each subject to calculate respiratory pulse synchronization at each voxel across the whole brain as individual Respsync maps (n = 162). Representative maps from one PWE case and one HC, with single brain voxel signals and their Respsync values against the fourth ventricle.

Previous work has shown that physiological signals arising from MREG data strongly correlate with cardiorespiratory pulsations.[16,35] Therefore, a sub-sample of 30 controls (age 43.0 ± 15.7 years, 19 females) with good quality respiratory belt signal was selected to identify brain areas that correlate best with individual respiration signal. The respiratory belt and MREG data signals were bandpass filtered according to individual mean respiratory frequency ±0.025 Hz (Figure 1(c)) and normalized. MATLAB xcorr cross-correlation analysis indicated that the MREG signals in the white area (Figure 1(c)) within the fourth ventricle match most accurately (correlation coefficient >0.9) with individual respiratory pulsations. Additionally, respiratory frequencies calculated from respiratory belt and equivalent MREG data region were compared with Pearson’s correlation to verify similarity between the two respiratory signal sources (Figure 1(d)). Moreover, the respiratory rates of these signals were compared statistically between all groups (Figure 1(d)).

Thus, the respiratory signal for each individual in following analyses was obtained in a data driven way from the MREG respiratory signal region in the fourth ventricle (Figure 1(c)). The measured respiratory frequencies ranged between 0.11–0.51 Hz (∼7–30 breaths per minute, bpm), and this range was used for filtering purposes over all study populations. The normalized mean MREG respiratory signals from the fourth ventricle were cross-correlated voxel by voxel with every other brain voxel, yielding individual whole-brain respiratory synchronization maps (Respsync, Figure 1(e)), which were then statistically compared between groups. We set six seconds (0.1 Hz) as maxima for xcorr temporal sliding lag between the analyzed signals in Respsync and used the sliding lag value of maximal correlation to obtain also respiratory lag maps of the brain. The formula for maximum (i.e. positive) Respsync values was:

Additionally, we calculated minimum (i.e. negative) correlation values similarly, as we wanted to know the direction of measured wave front. Group-wise positive and negative mean Respsync maps and their voxel-wise histograms were calculated to reveal global mean differences and deviations of the synchronization values of brain voxels between the three study groups. The individual Respsync maps were Fisher Z-transformed before calculating the group mean and comparing the histograms statistically. For visualization purposes, the values were back-transformed.

To verify the data, we checked random controls (51 vs. 51 subjects) and Freiburg vs. Oulu controls (14 vs. 14, age and sex matched), to ensure homogeneity without covariates by statistical analysis between these groups. Furthermore, the extracted respiratory signals from MREG data were analyzed for dissimilarities in the standard deviation between all three groups. All calculations were performed in MNI space at 4 mm resolution for computational efficiency. For display purposes, we interpolated the resultant figures in FSL and MRIcroGL.

To analyze the effect of disease severity, we combined the PWE and DN groups and assessed positive Respsync values as a function of the seizure frequency, irrespective of medication status. Patients were stratified into five levels according to their seizure frequency at the time of scanning verified from patient records. The grades were: 1: Seizure within a year, 2: two months, 3: one month, 4: one week and 5: daily. Two patients were included in both datasets, as they were scanned at the baseline and after beginning of medication. Thus, the duplicates were removed from the combined group consisting of 58 seizure patients. Additionally, we investigated the effect of disease duration, i.e., the voxel-wise relationship between whole-brain positive Respsync values and time since diagnosis.

Moreover, we chose the hippocampus as a representative landmark structure for its physiological and diagnostic relevance in epilepsy, to investigate how Respsync values might differ in this a priori region between groups by choosing the FSL’s hippocampus mask (Harvard-Oxford Subcortical Structural Atlas). Furthermore, we calculated the coefficient of variation (CV) and Fast-Fourier Transformation power (FFT) values in this hippocampi mask to afford a comparison of present results with a previous method we published earlier. For this analysis, we calculated the mean value of voxels in the hippocampi showing statistically significant group differences in each respective method (Respsync, CV and FFT, p < 0.05), because the significant voxel clusters were not identical for each method.

To exclude the possibility that observed differences might be due to head motion induced artefacts, we used FSL’s motion correction procedure to inspect for between-group differences in the mean relative head motion. Additionally, a complementary movement analysis of framewise displacement values was used to further control for possible movement differences.

Statistical analysis

The voxel-wise contrasts between the calculated synchronization maps were analyzed with FSL with 50,000 conditional Monte Carlo random permutations implementing family-wise error-corrected threshold-free cluster enhancement corrections separately in both directions (HC>PWE/DN, HC

To assess the potential utility of Respsync as an epilepsy biomarker, the receiver operating characteristic (ROC) curve and the area under the ROC curve (AUC) were calculated as a measure of classification accuracy from previously calculated hippocampi voxels. The ROC curve graphical plot is created by plotting the true positive rate against the false positive rate resulting in a graph showing binary classifier performance of certain disease, e.g. epilepsy. The bootstrap approach was used to estimate the 95% confidence intervals of the AUC. Moreover, we calculated the Spearman correlation coefficient between synchronization values from these same voxels and the frequency of epilepsy seizures. Statistical testing for differences in relation to age, relative head motion, and physiological respiratory rate differences were performed using the non-parametric Kruskal-Wallis one-way analysis of variance method (ANOVA) and post-hoc tests between the three groups. The non-parametric ANOVA was chosen after checking normality and lognormality of the variables between groups with Shapiro-Wilk test, and because not all variables were Gaussian distributed. The ordinary ANOVA was used to evaluate differences between the three histograms of positive and negative correlations (p < 0.05). With histograms, Tukey’s method was used for multiple comparisons in post-hoc tests.

Data availability and used software

The data supporting the findings of this study are available from the corresponding author, upon reasonable request. All the procedures described herein were analyzed and calculated with MATLAB, GraphPad Prism 8, MRIcroGL, Origin Pro 2019 b, Analysis of Functional NeuroImages and FSL, and procedural details are available upon request. For statistical analyses FSL and GraphPad Prism 8 were used. Additionally, Blender version 2.83.4 was used for graphical 3D modeling. Final figures were prepared with GraphPad Prism 8.

Results

Subjects

Between the groups, the relative head motion during the scanning was statistically significant (ANOVA p = 0.048). The PWE had less intra-scan relative head motion compared to HC (0.040 ± 0.023 mm and 0.050 ± 0.025 mm respectively), and it was nearly significant in the post-hoc test (p = 0.0503). The DN group movement (mean = 0.041 ± 0.014 mm) was of the same mean magnitude as in the PWE group, and no motion differences were detected between groups in post-hoc tests (HC vs. DN p = 0.77, PWE vs. DN p > 0.98). With framewise displacement value the results were parallel, with notion that with HC vs. PWE, the post-hoc test was significant (ANOVA p = 0.028, post-hocs HC vs. PWE: p = 0.030, HC vs. DN: p = 0.61, PWE vs. DN: p > 0.99). Furthermore, there were no differences between the calculated global tSNR values (ANOVA p = 0.27).

Additionally, there were no significant differences between groups in age (ANOVA p = 0.79) or respiratory frequencies determined from the power spectrum of MREG data (HC = 15.03 ± 3.84 bpm, PWE = 16.32 ± 4.51 bpm, and DN = 14.47 ± 4.39 bpm. ANOVA p = 0.36). The respiratory rates obtained from the respiratory belt and MREG data had, a correlation coefficient of 0.970 (p-value <0.0001, Figure 1(d)), indicating near perfect agreement between two modalities.

Elevated respiratory brain pulsation synchronicity

The brain area in which the MREG time series exhibited the highest similarity with the respiratory belt signals, was localized in the fourth ventricle, with the highest correlation coefficient being 0.923 (white area, Figure 1(c)). Furthermore, there were no differences in extracted respiratory signals between groups when normalized signal standard deviations were compared (ANOVA p = 0.75).

The group level synchronicity histograms were calculated to compare entire population values with each other. Both of these histograms (Figure 2) from positive and negative Respsync group mean maps showed same distinct statistical differences (ANOVA p < 0.0001, post-hocs HC vs. PWE: p < 0.0001, HC vs. DN: p < 0.0001 and PWE vs. DN: p < 0.0001).

Figure 2.

Distributions of mean positive (a) and negative (b) Respsync maps in each group calculated from maps on the right of the figure. Despite an overlap between the three similar distributions (HC: yellow, DN: orange, PWE: red), the differences between group means were highly significant (p < 0.0001) in both positive and negative Respsync values. Additionally, histograms show an increase of Respsync values as a function of diagnosis occurrence and (c) Mean positive and negative Respsync correlation lag maps. Despite visual dissimilarities, there were no statistical differences between these maps.

The statistical analysis of the PWE and HC groups showed markedly elevated positive Respsync values (p < 0.05, Figure 3(a)). The significant areas are widespread over the whole brain, covering essentially all functional structures, e.g. the upper brain stem respiratory pneumotaxic center, and the midbrain and temporal lobes, including hippocampi and pallida. There were no differences in positive Respsync values between HC and DN or between PWE and DN groups, although a trend was seen in the histograms (HC vs. DN, smallest p-value = 0.055). With negative Respsync there were statistical differences in temporal and occipital structures (Figure 3(b)). There were no differences in lag times between any of the groups (Figure 2(c)). Additionally, there were no differences in the control population homogeneity (51 vs. 51 controls, 14 vs. 14 Freiburg-Oulu age-sex matched), i.e. the control population synchronization values were consistent irrespective of scan time or site.

Figure 3.

(a) Synchronization of positive respiratory related brain pulsations (Respsync) was significantly (p < 0.05) higher in PWE compared to HC group. Significant increases were brain-wide, covering key functional structures, e.g., the upper brain stem respiratory pneumotaxic center and midbrain, along with the hippocampi and pallida. Additionally, the cingulate gyrus and temporal lobes clearly showed increased Respsync in the PWE group. (b) With negative synchronization Respsync values there were significant increases especially in the temporal and occipital structures and (c) The seizure frequency single group analysis showed that the Respsync values increase in the caudate, putamina, superior frontal gyrus and cingulate gyrus areas as a function of seizure density (p < 0.05).

Moreover, we conducted a single group analysis of the effect of seizure occurrence with patient data. The seizure frequency results shows explicitly that the correlation values increased in the caudate, putamen, superior frontal gyrus and cingulate gyrus regions as a function of seizure incidence (p < 0.05, Figure 3(c)). There were no differences in Respsync values as a function of disease duration.

Disease prediction and comparison between different methods

There was a clear and significant monotonic dependency between selected hippocampi area pulsations, as quantified by synchronization values and frequency of epileptic seizures (R = 0.432, p < 0.001). This same area was used to calculate ROC curves between groups, i.e., showing how well the values from different methods predicted disease. Between HC and PWE groups Respsync the AUC was 0.75 (p < 0.001) and between HC and DN groups the AUC was 0.66 (p = 0.025, Figure 4(a) and (b)). Between the patient group data, i.e., for PWE and DN the AUC was 0.61 (non-significant, p = 0.16). Additionally, the ROC and AUC values were calculated based on previously published coefficient of variation and power features from the hippocampi area (Figure 4(c) and (d)). The only significant pairs were the HC vs. PWE, in both coefficient of variation and power (HC vs. PWE: AUC = 0.69, p < 0.001, and AUC = 0.67, p = 0.002 respectively, Figure 4(c) and (d)).

Figure 4.

For ROC analysis, group-level differences in Respsync values of significantly different areas/voxels from the hippocampi were used. (a) The AUC between HC and PWE was 0.745 (p < 0.0001) and (b) with HC against DN the AUC was 0.677 (p = 0.0125). The same analysis was performed with previously published coefficient of variation and spectral power data with the addition of gathered subjects included to the previous population. Here, (c) coefficient of variation provided an AUC of 0.678 between HC and PWE (p = 0.001) and with (d) power data AUC was 0.664 (p = 0.002).

Discussion

In this study we developed a novel way to analyze the synchronicity of respiratory brain pulsations interictally in patients with epilepsy. We show that the high sampling rate of MREG data allows the precise estimation of respiratory signal from the fourth ventricle. The brain-wide synchronicity of this signal was significantly increased in epilepsy patients (treated and untreated), and this increase was a function of epileptic seizure frequency. The respiratory synchronicity quantified in the hippocampi area showed promising preliminary results as a diagnostic classifier.

Brain synchronous respiratory pulsations

As for how brain physiology affects MRI physics, the T2* weighted BOLD signal originates from transverse magnetization signal from coherent para/intravascular water proton spins. Thus, during inhalation, the venous blood drains from brain and CSF flows inwards, countering the perivascular volume loss, whereas in exhalation this net flow reverses resulting in pumping of perivenous convection, c.f. Figure 1(a) and (b), resulting in BOLD signal changes respectively.

Our current study showed significant elevation in synchronization values between HC (n = 102) and PWE (n = 40) groups, which were widespread over the whole brain, encompassing key functional structures known to be related to epileptic pathology in multiple epilepsy types.[37,41-43] Interestingly, these same areas are located around the cavernous sinus which is the main area of venous return from the temporomesial brain structures, and which is the site where CSF flow re-enters the brain. Therefore, in these brain areas, the respiratory tidal effects in perivenous space may be particularly strong, causing a measurable fMRI signal change, and thus is of different origin than the neurogenic signals in fMRI.

Previous intra-cranial needle electrode (iEEG) studies from temporomesial structures have shown that respiration has a strong effect on brain electrophysiological activity.[28,29] In fact, respiration modulates local field potentials and also cognitive functions. Those earlier studies were necessarily conducted only on epilepsy patients, given the invasiveness of the presurgical iEEG technique. Furthermore, the iEEG studies employed DC-EEG amplifiers covering EEG frequencies from 0.01 Hz upwards. Thus, the observed respiration synchronization to fMRI could also be present to some extent in clinical EEGs, since such recordings are usually high-pass filtered just above the respiration frequency range, i.e., over 0.5 Hz. Although iEEG researchers postulate that these electrophysiological effects may be a general brain phenomenon,[28,29] our recent and present results indicate that the increased variation, power and synchronicity of the respiratory brain pulsations are strongly associated with epilepsy in at least fast fMRI data.[26,28]

Importantly, there were no differences in respiratory pulsation signals from the fourth ventricle between any of the groups, denoting that the actual difference is not due to respiratory pulsation in the ventricle per se, but rather is a manifestation of the manner whereby the pulse spreads in the brain. Along the same line, the breathing frequency did not differ between groups, suggesting that the respiratory rates are not the source of the increase interictal synchronicity either.

It has been shown, that some AEDs, affect the cognitive performance and thus fMRI results. Interestingly, there were no significant group-level differences in Respsync maps between the HC and DN (n = 20) groups nor between the PWE and DN groups. This, together with histogram plots (Figure 2), show that the changes were still present also in the DN group. However, with additional drug-naïve patients, the contrast between Respsync maps in HC and DN could be significant also in the spatial analysis. Additionally, these results might have arised due to the more heterogeneous composition of the DN group, which included non-medicated patients with epilepsy (diagnosed during follow-up) and subjects without a final epilepsy diagnosis. This suggests that the seen effect is not AED related.

Earlier studies have shown progressive gray matter atrophy in focal epilepsy and temporal lobe epilepsy,[45-47] and that atrophy progresses in the seizure-free phase of the TLE. Interestingly, our results show increased interictal synchronization in PWE when compared with HC, reflected by the higher synchronization values in the same (e.g., frontocentral and cingulate areas, Figure 3(a)) as in previous studies.[45-48] Furthermore, taking seizure frequency into account, our results correlate to seizure frequency, and they localize to these same brain areas (Figure 3(c)) as in previous studies.[45-48] This indicates that higher synchronization values associated with more severe forms of epilepsy (seizure balance), which may be a factor in the tendency of severely afflicted patients to develop atrophy.

The more detailed inspection of the hippocampi provided further evidence for the aforementioned findings. There was significant monotonic correlation between respiration synchronicity and seizure occurrence in the hippocampi area, proposing that respiratory synchronicity grows stronger as a function of epilepsy severity. However, this correlation should be verified with larger multi-site cohorts where more homogenous sub-groups could be achieved. Additionally, the ROC curve between PWE and HC synchronization values shows potential as a diagnostic marker, especially keeping in mind that there is not always clear hippocampal dysfunction in epilepsy. As such, local maxima in synchronization from specific regions of the cortex, could potentially yield even better discriminatory values.

Pathophysiology behind altered pulsations in epilepsy

The blood–brain barrier (BBB) normally supports interstitial brain electrolyte homeostasis required for unhindered brain functioning. There is evidence of impaired BBB integrity in epilepsy, resulting in altered brain function associated with epileptic seizures examined with contrast-based measures. Additionally, glymphatic system sustains brain homeostasis through metabolite clearance. Therefore, the BBB transport and glymphatic system serve complementary roles. Thus, the impaired function of the BBB also affects the perivascular glymphatic brain clearance system neuron-glia functional unit at the outer edge of the BBB.[24,49,50]

The glymphatic system is driven by cardiovascular brain pulsations that convect CSF from perivascular space into the brain interstitium over the BBB glia limitans and then out of the brain via perivenous CSF.[19,51,52] The cardiorespiratory pulsations are the two main drivers of the intracranial CSF convection,[18-21] and these pressure pulsations extend into perivascular CSF spaces throughout the incompressible cranial vault, as dictated by the Monro-Kellie doctrine. Additionally, during seizures the blood pressure elevates, which have been shown to influence BBB functionality. In addition to blood pressure, the dysfunction seems to be mediated by an increased pinocytotic activity in the endothelial cells. Our group has discovered evidence that, on a macroscopic scale, respiration waves may function as a driver of the glymphatic pulsations in addition to the effects of slow vasomotor waves and faster cardiovascular pulsations.[22,27,31,35,54]

In resected epileptic foci (mesial temporal sclerosis), the perivascular spaces are enlarged, with reduction in the size or even absence of astrocytes. The remaining astrocytes were devoid of AQP4 water channels at their perivascular endfeet but had an ectopic excess in interstitial tissue.[23,24,49,55] The disposition of AQP4 channels has also been linked to impaired water movement which alters osmotic gradients between the neuropil and perivascular space, contracts the extracellular space and leads to increased neuronal excitability, possibly due to volume transmission. AQP4 water channels normally accompany inward-rectifier potassium channels, but in the absence of polarized AQP4 expression (as occurs in epilepsy), the inward-rectifier potassium channels cannot properly balance ion-concentration fluctuations during epileptic discharges. The increased respiratory brain pulsations that have been detected both in intracranial EEG[28,29] and in MREG signals[26,27] could thus reflect increased perivenous water convection over the glia limitans in the BBB.

In our study, we found no group-wise differences between the fourth ventricle signals. The causal reasons behind wide respiratory-related signal synchrony outside of the fourth ventricle in PWE remain unclear. We suggest that this may be due to a compensatory mechanism mediated by increase CSF convection at microscopical level or hindered impulse propagation in the system turning the CSF pulse nearly simultaneous over the brain in PWE. Indeed, further investigations are needed to resolve the issue. Whatever the cause is, it is most likely mediated via (peri)venous pulsations into the whole brain.

Technical considerations, limitations, and future prospects

Previous studies have demonstrated effects of the respiratory cycle on the BOLD signal. In a study of CSF flow oscillations in the fourth ventricle, showed that slow wave activity in the EEG signal correlated with CSF influx into the ventricles. In the present study, the high temporal resolution MREG scans over the whole brain with 3 mm spatial resolution suffice for the separation of respiratory pulsatility without aliasing over very low frequency hemodynamic signals. Combined with previous results, we are now able to show how these aspects are interrelated and how respiratory physiology affects brain pulsatility.

Correlations of the fMRI signal bear a relation to head motion and the spatial distance between compared components, i.e., larger movement is associated with higher correlation values at the local scale. That same study also showed that head motion did not account for the major proportion of the correlation values. In our study, patients with epilepsy had higher correlation values, i.e., greater synchronization than the control group even though their relative head movement was lower. Thus, the obtained brain-wide results are highly unlikely to be explicable by differences in head motion in ANOVA and nearly significant post-hoc test. Additionally, the DN group showed closely similar results to the PWE group.

In this study, the majority of the PWE subjects were diagnosed with TLE, and all had a focal epilepsy diagnosis. Additionally, the large group sizes enhanced our statistical power, thus ensuring reliable results. However, due to the great variety of epilepsy etiologies and affected brain areas, it is a challenging task to acquire a truly homogeneous patient population with matching diagnosis. Despite some heterogeneity in our study population, we consider it sufficiently homogenous for the present analysis. Moreover, the Siemens 3 T scanners used in Oulu and Freiburg were of different models with different scanner bore size and head coils. Nevertheless, the entry of imaging site along with head motion as covariates in the analysis minimized potential confounding effects to ensure the reliability of the results.

Even though we show clear evidence of respiratory related brain pulsations in fast fMRI, there is an obvious need for a more comprehensive set of studies, not only in fMRI, but also using different modalities. Fast fMRI and EEG in conjunction with hyperventilation or even sleep deprivation could provide further insights into the physiological mechanisms of our findings. Additionally, direct measures of CSF convection need to be quantified with respect to breathing in CSF spaces all over the brain.

Conclusions

Our analysis revealed compelling evidence of how respiratory brain pulsations are modulated by actual respiration in patients with epilepsy. The underlying pathophysiological reason for the differences is unclear, thus calling for further studies. The observed group differences were not explicable by head motion or physiological respiratory rates and are therefore most likely due to intrinsic brain pulsations in epilepsy, as was likewise seen in the drug-naïve patient population. Moreover, the fourth ventricle emerges as a significance landmark for different future fMRI, especially when effects of respiration are of interest.