A dynamic link between respiration and arousal

Viewing brain function through the lense of other physiological processes has critically added to our understanding of human cognition. Further advances though may need a closer look at the interactions between these physiological processes themselves. Here we characterise the interplay of the highly periodic, and metabolically vital respiratory process and fluctuations in arousal neuromodulation, a process classically seen as non-periodic. In data of three experiments (N = 56 / 27 / 25 women and men) we tested for covariations in respiratory and pupil size (arousal) dynamics. After substantiating a robust coupling in the largest dataset, we further show that coupling strength decreases during task performance compared with rest, and that it mirrors a decreased respiratory rate when participants take deeper breaths. Taken together, these findings suggest a stronger link between respiratory and arousal processes than previously thought. Moreover, these links imply a stronger coupling during periods of rest, and the effect of respiratory rate on the coupling suggests a driving role. As a consequence, studying the role of neuromodulatory arousal on cortical function may also need to consider respiratory influences. Significance statement We characterise the interplay of the respiratory rhythm and pupil diameter dynamics as a well-known proxy for arousal. Although we consistently find respiratory modulation of pupillary changes, they were most pronounced during periods of rest (compared to during task performance) and dependent on respiratory rate (deep vs. normal breathing).


Introduction
Until recently, our understanding of human cognition was largely dominated by an isolated brain-centric perspective of its underlying physiological mechanisms.Challenging this view, an increasing number of current findings embed the brain in a flurry of other non-cerebral physiological processes.Cardiac (Al et al., 2020;Galvez-Pol et al., 2020), gastric (Rebollo et al., 2021(Rebollo et al., , 2018)), respiratory (Kluger et al., 2023(Kluger et al., , 2021;;Park et al., 2020), and dynamics linked to arousal neuromodulation (Schneider et al., 2016;Groot et al., 2021) have all been shown to impact human cognitive function.Understanding these brain-body interactions has afforded overwhelming support for a more integrative view on brain function that acknowledges the meaningful interplay of central and peripheral physiological signals.It remains unclear however how physiological dynamics of the periphery interact with each other (see Criscuolo et al., 2022).The breathing rhythm has attracted particular interest as a pacemaker because it can be voluntarily controlled (Allen et al., 2023;Braendholt et al., 2023).A recent systematic review by Schaefer and colleagues (Schaefer et al., 2023) showed that respiratory coupling to arousal dynamics, as measured by pupillometry, remains critically understudied.Consequently, here we investigate the potential co-variations of endogenous fluctuations in respiratory and pupil-linked arousal dynamics.
From a neurophysiological perspective, the respiratory rhythm itself is initiated by pattern generators in the brain stem, particularly the preBötzinger complex (Del Negro et al., 2018).From here, efferent respiratory signals project to suprapontine nuclei (via locus coeruleus and olfactory nuclei) as well as to the central medial thalamus, which is directly connected to limbic and sensorimotor cortical areas (Yang and Feldman, 2018).Complementing this cascade in a top-down fashion, cortical areas in turn evoke changes in the primary respiratory network, e.g. to initiate specific motor acts (e.g., respiratory coordination during speech) or transitions between brain states (e.g., heightened arousal during a panic attack).
Arousal describes a global physiological preparedness to process and respond to sensory stimulation (Aston-Jones and Cohen, 2005).Physiologically, arousal levels are under the control of brainstem and basal forebrain nuclei.They determine momentaneous wakefulness and change when transitioning between behavioural states (Harris and Thiele, 2011;McGinley et al., 2015b;Schwalm and Rosales Jubal, 2017): The locus coeruleus (LC; brainstem) and the nucleus basalis of Meynert (basal forebrain) respectively control the release of norepinephrine (NE) and acetylcholine (ACh) to widespread cortical regions (Hasselmo, 1995;Lee and Dan, 2012;Steriade, 1996).Coincidentally, release of these neurotransmitters also affect pupil size therefore making pupillometry an effective readout of arousal (McGinley et al., 2015a;Reimer et al., 2014;Vinck et al., 2015).Importantly, variations in pupil-linked arousal have been shown to influence cortical activity (Pfeffer et al., 2022;Podvalny et al., 2021;Radetz and Siegel, 2022) and cognitive function (Waschke et al.,201); Kosciessa et al., 2021).In contrast to breathing, pupil dynamics, linked to arousal neuromodulation, do not usually present with a strong rhythmic component.However, a resting rhythm is known to emerge during prolonged periods of low arousal, the so-called Hippus (Mathôt, 2018).The frequency of this rhythm, typically around 5 seconds per cycle (Bouma and Baghuis, 1971), roughly coincides with the average pace of the breathing rate.
The relative neighbourhood of respiration-and arousal controlling brainstem structures, as well as the existence of concrete synaptic connections from the preBötzinger complex to LC suggests an interaction between both functions (Del Negro et al., 2018) and may explain earlier findings (Ohtsuka et al., 1988).In fact, interrupting the connection between both cores in the rodent brain led to chronic hypoarousal and "lethargic" behaviour (Yackle et al., 2017).Critically, these changes in balancing calm and arousal behaviour were entirely independent from respiratory changes and thus underscore the importance of propagating the respiratory rhythm from brainstem generators to upstream areas.Although a recent review on potential links between respiratory and pupil dynamics in the human body failed to find conclusive evidence (Schaefer et al., 2023), we hypothesise that the anatomical and functional links between the respiratory and arousal system are specifically reflected in respiratory modulation of pupil size.
This study takes a new look at the respiration-arousal interplay, measured as changes in pupil diameter and tidal volume.Considering the projections between their main orchestrators in the brainstem (Fig. 1a), we first tested whether the emerging Hippus would indicate a stronger coupling, or entrainment, of arousal to the strong periodicity of the breathing rhythm (Nakamura et al., 2019;Ohtsuka et al., 1988) in N = 56 healthy volunteers resting for 5 min.In a second dataset (N = 27), we tested how the coupling changed when participants engaged in voluntary deep (vs normal) breathing.Finally, for a subset of participants of the first study (N = 25) we also compared the spectral make-up of respiratory and pupil-linked arousal dynamics between rest and during task performance, and tested whether a potential coupling would depend on their behavioural state.Below, we show evidence for a coupling of respiratory and pupillary dynamics during rest.Coupling characteristics changed with the respiratory rate and coupling strength decreased when participants were engaged in a task.

Results
First, we analysed continuous 5-min long measurements of pupil diameter and respiratory tidal volume from N = 56 resting volunteers who took part in a series of MEG studies.Analyses that investigated the relationships of MEG-recorded cortical activity with either pupil or respiratory dynamics have been published elsewhere (Kluger et al., 2023(Kluger et al., , 2021;;Pfeffer et al., 2022).Here, we focus on the link between respiratory dynamics and changes in pupil size as an indicator of arousal neuromodulation.

Coupling between respiration and pupil-linked arousal at rest
Separate spectral analyses of pupil-and respiration time series showed that respiratory tidal volume fluctuated periodically at the typical breathing rate of about 0.25 cycles per second (Fig. 1B).Further, smaller spectral peaks can be observed at harmonic frequencies due to the non-sinusoidal waveform of the tidal volume fluctuations.Although less pronounced when compared with the respiration spectrum, pupil dynamics also contained a distinct periodicity that peaked at 0.19 Hz.This is consistent with the previously described Hippus, a resting rhythm of the pupil (Bouma and Baghuis, 1971;Pomè et al., 2020).Respiratory pacemaker cells in the preBötzinger complex (pBC) project to locus coeruleus (LC), which controls global release of noradrenaline (NE), and to centromedial thalamus (CM).In a feedback loop, the incoming airstream of each breath triggers receptors reaching from the epithelium into the nasal cavity, propagating breathing-locked activity to the piriform cortex (PC) and mediodorsal thalamus (MD).B Spectral composition (power spectra) of 5-min resting pupil (red) and respiration (blue) time series.Shaded area shows SEM.C Mag-squared coherence spectrum between pupil and respiration time series.Observed coherence (purple) plotted against coherence computed on surrogate pupil time series (grey).Black straight line shows frequency range with significant differences (p < .05,cluster-based permutation test).D Power envelope correlations (Pearson's correlation coefficient) between respiration and pupil time series.Marked areas denote significant differences at p < .05(uncorrected; but relaxed cluster criterion applied -see Methods).E Distribution of individual correlation coefficients (pupil-respiration power correlations) averaged across the largest significant negative cluster in D. Observed correlations (purple) pitted against correlations based on surrogate pupil time series (grey).F Pupil power spectra for three consecutive 2-min time windows.Darker colour = later time window.Black straight line shows frequency range with significant changes between time windows (p < .05,cluster-based permutation test).G Individual power values per time window at the Hippus peak frequency (0.19 Hz) in F. H Respiration power spectra, otherwise same as F. I Coupling (magnitude-squared coherence spectra) between pupil and respiration, otherwise same as F and H.No significant monotonic changes observed between time windows.
Next, we tested for the coupling of both signals by means of spectral cross-coherence.Individual coherence values, quantified as the logarithm of the amplitude of the coherence spectrum were subjected to a cluster-based permutation procedure that tested against surrogate coherence spectra based on reconstructed pupil time series with identical powerbut scrambled phase spectra (IAAFT approach; Schreiber and Schmitz, 1996; also see Kluger et al., 2023).This analysis returned a cluster between 0.17 and 0.37 Hz in which coherence in the original data systematically exceeded the coherence expected by chance (Tsum = 41.57,pcbpt < .001;Fig. 1c).This cluster contained the peak coherence at 0.29 Hz.A smaller adjacent cluster was found between 0.09 and 0.14 Hz (Tsum = 19.03,pcbpt = .002;Fig. 1c) We also tested whether arousal-and respiratory signals covary in magnitude by subjecting wavelet-decomposed time series, i.e., power envelopes, to a frequency-frequency correlation analysis.Again, this was tested against correlations with a surrogate pupil time series generated by means of IAAFT (Fig. 1D).This analysis produced a cluster of negative correlation that coincided with the peak frequencies of respective power spectra, peaking at 0.18 Hz (respiration) and 0.28 Hz for the pupil (cluster peak Pearson's r(54) = -0.19,puncorr < .05for a minimum of five adjacent pixels in the frequency-frequency plane; a cluster-based permutation test did not return any significant clusters; also see Figure 1E for the distribution of individual correlation values averaged across that cluster).A similar smaller cluster at a higher respiration frequency mirrored the effect and was likely due to the strong harmonics in the spectral composition of the respiration time series (see Figure 1B).As the peak power of the respiratory rhythm is a measure of the tidal volume, or depth of breathing, this suggests that shallower breathing coincides with a stronger pupillary Hippus.An interpretation of this negative association may be that, during rest, episodes of lower metabolic requirements reduce respiratory depth, yet increase the coupling of arousal to respiration producing a stronger Hippus.

Hippus increases during rest
To explore changes in the coupling of respiratory and pupil dynamics over time, we split the 5 min resting state recordings into consecutive overlapping segments of 150 sec (3 segments, 50% overlap).Power spectra of each segment were then submitted to a regression analysis using cluster-based permutation testing to identify frequency ranges with monotonous changes across segments.These analyses identified a cluster centred on the Hippus frequency range for pupil dynamics (Tsum = 174.52,p = .001)and ranging from 0.05 to 0.58 Hz (Fig. 1F-G).We also found an increase in power in the frequency ranges between 0.05 and 0.1 Hz in the respiration data (Tsum = 54.95,p = .015).However, this effect markedly excluded the spectral peak indicating that the magnitude of breathing remained constant over time (Fig. 1H).Finally, no significant changes were registered in the coupling between the two measures in each individual time window (Fig. 1I).These results also suggest that the negative power-power correlation reported above is likely emerging from a source of variability that is not explained by time passing (see e.g.Benwell et al., 2018 for time on task as a mediator in a different context), as respiratory tidal volume does not seem to change monotonically (decrease) across the three time windows.

Respiration-pupil coupling follows breathing rate
We explicitly tested the hypothesised influence of breathing depth on respiration-pupil coupling in a separate data set (Kluger et al., 2023).Participants were instructed to rest while breathing normally vs deeply in two separate 5-min measurements.Looking at the power spectra of pupil size and tidal volume time series, deeper breathing led to an overall decrease in respiratory rate, from the typical 0.27 Hz to 0.17 Hz on average (Wilcoxon signed rank test, Z = 3.75, p < .001;see Fig. 2A).Interestingly, pupillary dynamics did not exhibit a strong Hippus in this dataset.The typical peak at around 0.20 Hz seemed more pronounced in the normal breathing condition, however without producing measurable power differences between both conditions.Extracting Hippus peak frequencies therefore required removing the 1/f trend in the spectra (explained in detail in the Methods) but showed a pronounced difference between conditions (Wilcoxon signed rank test, Z = 3.27, p = .001;see Fig. 2B).Finally, we also observed substantial phase coupling between both measures (normal breathing: Tsum = 14.75, pcbpt = .005;deep breathing: Tsum = 19.46,pcbpt = .007;also see Figure 2C).Importantly, the peak frequency of the coupling also showed a substantial decrease from 0.25 to 0.17 Hz, suggesting that the momentary respiratory rate strongly influences the coupling (Wilcoxon signed rank test; Z = 3.77, p < .001).

Respiration-pupil coupling reduced during task performance
For a subset of the resting state data (N = 25) we were able to compare the above findings with similar data recorded while participants were performing a visual spatial detection task (Kluger et al., 2021).As quiescence and task performance have been shown to involve different levels of arousal (Podvalny et al., 2021;Reimer et al., 2014), we tested for a difference in coupling between respiration and pupil dynamics.Specifically, we compared the last of six blocks of task performance in an MEG experiment (length = 446 ± 34s [M ± SD]) with the immediately following resting state recording.We found evidence for decreased power of respiration during rest in frequency ranges above (0.38 -0.52 Hz, Tsum = -6.81,pcbpt = .044;0.58 -0.80 Hz, Tsum = -6.80,pcbpt = .044;1.12 -6 Hz, Tsum = -42.51,pcbpt < .001),but not including the respiratory peak frequency (see Figure 2D).These differences may point at changes in oscillatory waveform of the respiratory rhythm that do not impact power at the fundamental frequency but only at higher harmonics (an index of non-sinusoidal properties of the time series).For the pupil, we observed systematically higher power in the Hippus frequency range (0.19 -0.28 Hz; Figure 2E) during rest (Tsum = 9.96, pcbpt = .033),as well as a difference in the range between 0.90 -2.46 Hz (Tsum = 11.11,pcbpt = .026).Ultimately, we contrasted the coupling in both conditions and found that coupling was stronger during rest within 0.16 to 0.31 frequency range (Tsum = 27.17,pcbpt = .001),i.e., largely commensurate with the frequency range we observed the coupling in the resting state analysis above.

Discussion
In this investigation of respiration-arousal interplay, we demonstrate substantial phase coupling of pupil dilation dynamics to the breathing rhythm during rest.In an independent data set, changing the respiratory pattern towards voluntary deep breathing substantially decreased the peak frequency of respiration-pupil coupling.Finally, coupling decreased during a visual perception task (compared to rest) in a third data set, which opens up instructive avenues for comparing effects of contextual and respiratory interventions on phase coupling between respiration and pupil-linked arousal.

Strong evidence for coupled rhythmic pupillary and respiratory dynamics
In their recent review of the literature on links between respiration and pupil dynamics, Schaefer and colleagues (Schaefer et al., 2023) graded the overall strength of evidence for the effect of breathing phase on pupil dynamics as "low", and for the effects of breathing depth and rate as "very low".Notably, their sample included studies with various measures of pupil dynamics (e.g., tonic pupil size, light reflex) and their correlations with respiration parameters.Focussing on oscillatory coupling in our investigation, we identified the rhythmic Hippus as a measure of interest, due to its prominent peak at around 0.2 Hz in spectra of pupil time series but also due its potential role as an indicator of momentary brain state (McLaren et al., 1992).More specifically, quantifying oscillatory coupling by means of phase coherence (e.g.Gross et al., 2021), an approach that had not been used before (see Schaefer et al., 2023), reliably produced a peak coherence between 0.2-0.3Hz in three datasets (during normal breathing).Given that this is the typical range of the respiratory rate (Fleming et al., 2011), as well as the emerging Hippus during rest (Bouma and Baghuis, 1971;Turnbull et al., 2017), we provide strong evidence in favour of an effect of breathing phase (sensu Schaefer et al., 2023) on the Hippus aspect of pupillary dynamics.
Our approach does not allow a direct inference of whether respiration is the driver of the periodic arousal modulation linked to the Hippus.However, given the relative strength and constancy of the respiratory rhythm, a respiratory drive seems likely, and aligns with earlier ideas about the neural circuitry of this link (Borgdorff, 1975;Ohtsuka et al., 1988) as well as its demonstration in the animal model (Yackle et al., 2017).It is therefore possible that the rhythmic Hippus can be explained by the entrainment of locus coeruleus neurons by the respiratory pacemaker in the pre-Bötzinger complex, a phenomenon known to occur between distinct but connected neuronal populations with the ability to produce concerted rhythmic activity (Thut et al., 2011;Herrmann et al., 2016;Thut et al., 2011).
Phase entrainment can explain the coupling, yet additional mechanisms need to be considered to explain the increase in Hippus magnitude at the peak frequency over the course of a 5 min resting state recording.Although this effect aligns with the idea that an increased Hippus indicates decreased vigilance or alertness as relaxation sets in (McLaren et al., 1992), it cannot simply be attributed to a similar increase in the magnitude of the fundamental respiratory rate, or a change in coupling strength because we found these to remain constant.Moreover, our results suggest that peak-frequency respiratory and Hippus magnitude may share an inverse relationship during rest, whereas an account in terms of entrainment (or resonance) would predict that increased magnitude in respiration leads to an increased Hippus.A more comprehensive understanding of the observed effects requires exploring factors that affected coupling.

Breathing mode and behavioural state alter pupil-respiration coupling
Changes we observed in pupil-respiration coupling based on manipulations of breathing mode (normal vs deep breathing) and behavioural state (rest vs task) shed further light on the interplay between respiratory and Hippus rhythms.Absent a direct measure of respiration, Pomè and colleagues (2020) have recently argued that observed changes in pupillary dynamics during mindfulness meditation may have been a consequence of altered breathing patterns.However, evidence for pupil size variations as a function of breathing depth has remained inconclusive (see Schaefer et al., 2023 and, to date, only few studies have contrasted the effects of different breathing modes directly.Debnath and colleagues (2021), for instance, investigated relationships between several physiological measures during a variety of physical exercises and stress tests that affected breathing and reported no overall effect.In contrast, Schumann and colleagues (2020) do report an increase in pupillary unrest during deep breathing, a measure for which they aggregate spectral power below 1.5 Hz.Similar to Schumann et al., we observed a marked decrease in breathing rate when participants were instructed to breathe deeply.Along with the change in breathing frequency, both the peak frequency of respiratory and Hippus signals decreased, as did the peak frequency of pupil-respiration coupling.This strong dependence of coupling frequency on breathing rate adds further support to the idea that the observed pupil-respiration link is primarily driven by the respiratory pacemaker.It may also provide an explanation for the inverse relationship between peak respiratory and Hippus power observed in the resting state data: Normal fluctuations in respiratory rate (peak frequency) between deeper and shallower (i.e., slower and faster) breathing lead to varying resonance in the neuromodulatory circuitry that influences pupil size.
In contrast to the frequency-changing effect of breathing mode, we observed a frequencyspecific reduction in peak coupling strength (at around 0.2 Hz) when participants engaged in a task vs resting.Notably, this effect occurred in the absence of any peak frequency changes in respiration itself (although power differences in harmonic frequencies may point at waveform changes during task engagement), while the pupillary Hippus was effectively abolished during task performance.Together with its gradual increase during rest, this effect underlines its potential role as an index of current alertness, momentaneous internal vs external focus, and more generally, arousal levels where a higher Hippus indicates a state of lower arousal (Bouma and Baghuis, 1971;Mathôt, 2018;McLaren et al., 1992;Bouma and Baghuis, 1971;Mathôt, 2018).Taken together, the malleability of this indicator to differences in breathing mode suggest a strong respiratory drive, and the decrease in respiratory-coupling during task performance suggests that the respiratory drive may be more permeable in states of low arousal.

Pupil-respiration coupling suggests interplay of breathing and arousal
Physiological arousal is primarily controlled via the noradrenergic neurotransmission as regulated by the locus coeruleus (LC; (Aston-Jones and Cohen, 2005;Mather and Harley, 2016).As the primary source of noradrenergic release in the mammalian brain, the LC projects to subcortical and cortical regions (Aston-Jones and Cohen, 2005).In concert with the release of acetylcholine (ACh) in the basal forebrain, the LC-NE axis profoundly influences cortical activity (Pfeffer et al., 2022;Podvalny et al., 2021;Radetz and Siegel, 2022) with consequences for cognitive function (Waschke et al., 2019;Kosciessa et al., 2021).An intriguing possibility then is that respiratory rhythms constitute another, potentially earlier node in a sequence of physiological processes that modulate cortical activity by controlling noradrenaline release, creating a causal loop between interoceptive sensation, respiration, and arousal.In line with this, a recent study by Yackle and colleagues (2017) identified a subpopulation of neurons in the preBötzinger Complex, the primary respiratory rhythm generator of the brainstem, with direct projections to noradrenaline expressing LC neurons.Ablating these connections further eliminated the breath-by-breath control of noradrenaline release, causing mice to exhibit altered arousal responses to exteroceptive stimuli.
In the light of these connections in the underlying neural and neuromodulatory circuitry, and our finding that pupil-respiratory coupling is stronger during a restful, introspective behavioural state, we can speculate that the brain may adopt a metabolically optimal resting mode that is characterised by a rhythmic intermittent release of neuromodulators, instead of a tonic, more continuous release during active performance that requires a higher sustained level of cortical activation for the processing of task-relevant sensory input (see Schroeder and Lakatos, 2009 for a similar suggestion regarding the functional significance of cortical electrophysiological oscillations).This account allows testable predictions: Cortical markers of sensory processing may underlie stronger periodic fluctuations with frequencies in the range of the coupling observed here when participants are in a more restful than alert state.
In summary, we show evidence for an oscillatory coupling of respiratory and pupil dynamics that critically depends on breathing phase and rate, and also changes depending on behavioural state.These results contrast with the results of a recent survey which showed low confidence overall in pupil-respiration links (Schaefer et al., 2023).However, the survey considered a range of possible links, whereas we focussed on oscillatory coupling, possibly explaining the divergence, and pointing out a promising direction to take to explore these links further.The close link between pupil dynamics and neuromodulatory influences of the arousal system on cortical function on one side, and the potential drive of pupil-linked arousal by respiration, poses new exciting challenges for our understanding of how physiological processes influence cortical activity, and hence, human cognition.

Participants and procedures
Respiratory and pupil data of N = 57 right-handed volunteers (31 women, age 25.3 ± 3.1 years [M ± SD]) were originally recorded for MEG studies published elsewhere (Kluger and Gross, 2020); Kluger et al., 2023;Pfeffer et al., 2022).All participants reported having no respiratory or neurological disease and gave written informed consent prior to all experimental procedures.The studies were approved by the local ethics committee of the University of Münster (approval numbers 2018-068f-S and 2021-785-f-S).During all procedures, participants were seated upright while we simultaneously recorded respiration with a respiration belt transducer around the chest (BIOPAC Systems, Inc, Goleta, USA) and pupil area of the eye with an EyeLink 1000 plus eye tracker (SR Research).Both, respiration and eye tracking data, were sampled at 600 Hz.All N = 57 participants completed a 5 min resting-state recording during which they were to keep their eyes on a fixation cross centred on a projector screen placed in front of them.A subset of n = 27 participants (12 women, age 25.0 ± 2.8 years) took part in a second study (see Kluger and Gross, 2020) in which they were instructed to breathe either naturally or voluntarily deeply through the nose for two separate 5 min runs.A second, independent sample of n = 25 participants (13 women, age 25.5 ± 2.7 years) took part in a third study (see Kluger et al., 2021) in which they performed a simple visual detection task: Gabor patches were presented for 50 ms at near-threshold contrast either to the left or to the right of a central fixation cross displayed on a projector screen in front of them.After a short delay of 500 ms, participants were to report via button press whether they had seen the target on the left, the right, or no target at all.For comparison with the resting condition, we report results from the last of six task runs (120 trials, mean duration 446 ± 34 sec).

Data preprocessing
Eye tracking and respiration time series of all three datasets used in analyses underwent largely similar preprocessing steps.Analyses made use of the Matlab toolbox fieldtrip (Oostenveld et al., 2011) in combination with custom-written code.All analysis scripts can be found here: OSF | Spontaneous covariation of pupil size and respiration.
Pupil area traces were converted to pupil diameter to linearize our measure of pupil size.Using routines available from https://github.com/anne-urai/pupil_preprocessing_tutorial(also see Urai et al., 2017), blinks were identified by an automatic (and visually validated) procedure and linearly interpolated.We slightly modified the original procedure in that blinks were detected in a first pass with a standard criterion of z = 3 SD.Then, blink-interpolated pupil time series were subjected to the procedure again using a relaxed criterion of z 6 SD to capture remaining artifacts.Next, canonical responses to blinks were estimated and removed from pupil time series (Hoeks and Levelt, 1993;Knapen et al., 2016;Wierda et al., 2012).To that end, pupil time series were band-pass filtered (pass band: 0.01-10 Hz, second-order Butterworth, forward-reverse two-pass).Due to its robustness against artifacts, Respiration data underwent the same filtering procedure but no other preprocessing.For analysis, all time series were converted to zscores using a robust procedure (MATLAB function 'normalize' with options set to 'zscore' and 'robust').

Data analysis
Spectral analysis.Power spectra of full-length recordings for all three datasets were obtained by subjecting pupil-and respiration time series to spectral decomposition by means of the multi-taper method as implemented in fieldtrip.Spectral smoothing was set to 0.05 Hz.All time series were zero-padded to a length of 600 sec to match the spectral resolution irrespective of the original length of the recordings.We used a logarithmically spaced frequency axis, a frequency range of interest from 0.065 to 6Hz, avoiding the influence of filter artifacts (see Data preprocessing).As a final step, all power spectra were converted to dB scale by taking the decadic logarithm and multiplying by 10.For pairwise comparisons between power spectra we generally used the cluster-based permutation approach based on dependent-samples T-tests, and clustering results along the frequency dimension (n = 5000 permutations).Generally, cluster-corrected p-values had to satisfy a threshold of ɑ = .05.
Respiration-pupil coupling -coherence.Coupling between respiration and pupil time series was quantified using the phase-based magnitude-squared coherence metric (Carter et al., 1973).To this end, we used the same approach as described above, but retained complex fourier spectra, and subjected these to a coherence analysis as implemented in fieldtrip (ft_connectivityanalysis, option 'method' set to 'coh ', and 'complex' to 'complex').Additionally, this analysis was run on a combination of the original respiratory time series and a surrogate pupil time series.Surrogates were generated based on the spectral characteristics of the original pupil time series, retaining power-but scrambling phase information (IAAFT approach; Schreiber and Schmitz, 1996) also see Kluger et al., 2023).The logarithmized magnitudesquared coherence spectra for respiration coupling with original and surrogate pupil data were then subjected to the cluster-based statistics as described for the power spectra above.For the visualisation of the distribution of individual data points in Figure 1E (as well as 1G, and insets in Figure 2A and B) we used the raincloud plot functionality for MATLAB (Allen et al., 2019; as implemented here: https://github.com/RainCloudPlots/RainCloudPlots/tree/master/tutorial_matlab).Identical approaches were taken when testing respiration-pupil coupling (coherence) in rest-vs-task and normal-vs-deep breathing datasets.
Respiration-pupil coupling -power correlation.In the resting-state only dataset (N = 57), we also tested for potential covariations in power of respiration and pupil dynamics over the of the 300 sec of the recording.To this end, time series were subjected to a spectrotemporal decomposition using the multitaper approach as implemented in fieldtrip (function 'ft_freqanalysis', with 'method' set to 'mtmconvol').Again, data were zero-padded to 600sec, and spectral smoothing was set to 0.2 Hz.For 36 frequencies, spanning a range from 0.065 to 2 Hz logarithmically, spectral representations were computed for frequency-dependent time windows of nine cycles each, and with a temporal step size that had consecutive windows overlap by 50%.The following time-frequency representations (power envelopes) for respiration data were correlated (Pearson's correlation coefficient) with power envelopes of the pupil data, after dB-scaling them.Correlation coefficients were transformed to Fisher's z before subjecting to statistical testing.As a contrast, we correlated original respiration power envelopes with surrogate pupil power envelopes, again derived from IAAFT (Schreiber and Schmitz, 1996;see Kluger et al., 2023; also see section 'Respiration-pupil couplingcoherence').
Initially, z-transformed original correlation coefficients were then compared against the surrogate data by means of cluster-based permutation testing using pairwise comparisons ('ft_statfun_depsamplesT') in the two-dimensional frequency-frequency plane (resulting from correlating respiration power envelopes of every frequency with pupil envelopes of every frequency, and vice versa).Prior to the testing a blurring filter was applied to individual frequency-frequency planes (original and surrogate to reduce the influence of small interindividual differences in frequencies (MATLAB function 'imgaussfilt', default settings).This analysis did not produce any clusters that survived correction for multiple comparisons.We then applied a relaxed cluster thresholding to a map obtained by applying paired t-tests to data at each frequency-frequency pair.In this map a cluster was only considered significant when consisting of at least five adjacent frequency-frequency correlations below the uncorrected ɑ threshold.Correlation values (not z-transformed) were averaged across frequency-frequency correlations of the largest cluster (and displayed in Figure 1E).Note that this relaxed criterion increases the chance of false positives.Time-window analyses.For analysing changes over time in the resting-state only dataset (N = 57), preprocessed time series (all 300 sec in length), were re-epoched into three overlapping time windows of 150 sec each, with 50% overlap between time windows.Windowed data underwent the same spectral decomposition approach as described in section 'Spectral analysis'.To test for monotonous trends, increases or decreases, in pupil and respiratory power, as well as in pupil-respiratory phase coupling (magnitude-squared coherence) across the time windows, we used the cluster-based permutation approach based on dependentsamples regression, as implemented in fieldtrip ('ft_statfun_ depsamplesregrT').
Peak frequency shifts in normal vs deep breathing comparisons.When analysing the normal-vs-deep breathing dataset it became clear that the major effect in respiratory, pupil and coupling dynamics was a change in peak frequencies rather than spectral power (or magnitude-squared coherence when looking at coupling).To capture this effect, we refrained from the general approach to testing spectral differences described above and instead extracted individual peak frequencies as follows: First we increased the resolution of individual spectra by a factor of 100 using spline interpolation (MATLAB function 'interp1', option 'spline').This allowed capturing finer differences in the exact individual peak frequencies.Then, we detected peak frequencies in leave-one-out subaverages using 'findpeaks' (MATLAB) in the frequency range between 0.1 and 1 Hz.For peak detection in the pupil power spectra, we had to remove the 1/f trend first using a procedure described in Keitel et al. (2018); see Supplementary Materials.This approach applied a Savitzky-Golay filter of length 9 and order 3 to the sign-flipped second-order gradient of the pupil power spectra to accentuate peaks.We reduced the spectral range for the peak detection to 0.2 to 0.6 Hz in this case because the second-order cirgradient also accentuates spectral noise.For peak frequencies derived from all spectra, subaverages were corrected for deflated variance following Smulders (2010).This procedure restored estimates for individual peak frequencies but also produced strong outliers in some cases, when the exclusion of a given participant had a strong effect (leverage) on the group average.These outliers have been excluded from depictions in Figure 2A-C, but have been left in the data for statistical comparisons.To reduce their influence on overall outcome we used non-parametric Wilcoxon signed-rank statistics to test differences in peak frequencies for respiration and pupil power spectra, as well as respirationpupil coherence spectra.

Figure 1 .A
Figure 1.A Neuroanatomy of respiratory and pupillary rhythms.Respiratory pacemaker cells in the preBötzinger complex (pBC) project to locus coeruleus (LC), which controls global release of noradrenaline (NE), and to centromedial thalamus (CM).In a feedback loop, the incoming airstream of each breath triggers receptors reaching from the epithelium into the nasal cavity, propagating breathing-locked activity to the piriform cortex (PC) and mediodorsal thalamus (MD).B Spectral composition (power spectra) of 5-min resting pupil (red) and respiration (blue) time series.Shaded area shows SEM.C Mag-squared coherence spectrum between pupil and respiration time series.Observed coherence (purple) plotted against coherence computed on surrogate pupil time series (grey).Black straight line shows frequency range with significant differences (p < .05,cluster-based permutation test).D Power envelope correlations (Pearson's correlation coefficient) between respiration and pupil time series.Marked areas denote significant differences at p < .05(uncorrected; but relaxed cluster criterion applied -see Methods).E Distribution of individual correlation coefficients (pupil-respiration power correlations) averaged across the largest significant negative cluster in D. Observed correlations (purple) pitted against correlations based on surrogate pupil time series (grey).F Pupil power spectra for three consecutive 2-min time windows.Darker colour = later time window.Black straight line shows frequency range with significant changes between time windows (p < .05,cluster-based permutation test).G Individual power values per time window at the Hippus peak frequency (0.19 Hz) in F. H Respiration power spectra, otherwise same as F. I Coupling (magnitude-squared coherence spectra) between pupil and respiration, otherwise same as F and H.No significant monotonic changes observed between time windows.

Figure 2 .
Figure 2. A Contrasted normal vs deep breathing conditions.Marked decrease in respiration frequency.Inset raincloud plot shows individual distributions of peak frequencies in both conditions (jackknife-based individual estimates -see Methods for further details).B Same as in A but for pupillary dynamics.No significant differences in Hippus peak frequencies between normal and deep breathing.C The coupling (magnitude-squared coherence) peak frequency is markedly reduced during deep breathing.Inset raincloud plot shows individual distributions of peak frequencies in both conditions.Straight lines at the bottom signify systematic coherence differences when compared with surrogate data (not shown).D Contrasted rest vs task conditions (different sample than A-C).Respiration power spectrum during rest (dark blue) and task performance (light blue).Straight lines at the bottom indicate frequency ranges with significant changes between both behavioural states (p < .05,pairwise t-statistic, cluster-based permutation test).E Same for power spectra of pupil time series, dark red = rest, light red = task.F Comparison of respiration-pupil coupling, measured as magnitude-squared coherence.Dark purple = rest, light purple = task condition.