The study was approved by the Institutional Review Board of the Oregon Health & Science University (OHSU), and the full ongoing RCT was registered at the (ID # NCT03858309). We received verbal and written informed consent from all study subjects prior to all study procedures. We recruited healthy participants from the Portland metropolitan area using OHSU’s study participation opportunities website, Oregon Center for Clinical and Translational Research Institute (OCTRI) research match for recruitment, flyers throughout the OHSU campus and communities in Portland, and social media (Facebook). We aimed to enroll participants 18 to 65 years of age who were able and available for study activities including undergoing non-invasive MRI scans, had no current or previous regular practice of mind–body therapies focusing on breath awareness and/or training (e.g., yoga, meditation, Tai-Chi, Qi-Gong), and were in good health without any history of neurological disorders, sleep disorders, respiratory disorders, problems with heart, circulatory system, and lungs. See Table S1 for a full list of RCT inclusion/exclusion criteria. Of the 65 participants contacted for the study, 57 were phone screened, 26 were enrolled, 21 completed the baseline procedures (September–October, 2019 at OHSU), and 18 were included in final baseline data analysis (N = 18, eight females, ten males; mean age: 34.9 ± 14 (SD) years; age range: 18–61 years). See Fig. 1 for the study flow chart, and Table 1 for the study group characteristics (N = 18).

Figure 1
figure 1

Study flow chart. *Our study utilized physiological data devices to objectively track participants’ home practice during the 8-week interventions. We excluded participants who did not have a compatible electronic device such as smartphone or tablet (see Table S1).

Table 1 Study group characteristics.

Experimental methodology

Each subject’s imaging visit lasted approximately 3-h including study instructions, 1-h MRI scans, and a set of questionnaires (as part of the RCT activities; not reported herein). Upon arrival for their imaging visit, we measured each subject’s temperature, blood pressure, and height and weight for body mass index (BMI). We then transitioned subjects to a mock scanner (0T room) for a ~ 30-min instruction for the breathing practices to be performed during the RT-PCMRI scans. The mock scanner room is designed to prepare research study subjects prior to entering the MRI Instrument suite, and is useful to help acclimate subjects to the enclosed space inside an MR Instrument. We first explained and demonstrated each breathing practice, then guided subjects to perform at their own pace first seated on a chair, and then in supine in the mock scanner to mimic the MRI environment.

MRI breathing protocol

We instructed subjects to perform the following breathing protocol first in the mock scanner for training purposes, and then in the MRI instrument during the ~ 1-min RT-PCMRI measurements, each to be collected twice with 30–60 s between consecutive measurements in the following order: (1) spontaneous breathing (SponB), (2) slow breathing (SlowB), (3) deep abdominal breathing (DAB), (4) deep diaphragmatic breathing (DDB) (5) deep chest breathing (DCB). See Table 2 for the MRI breathing protocol details.

Table 2 MRI breathing protocol.

Rationale for the MRI breathing protocol design

We chose breathing practices that were easily performed in supine in an MRI instrument without any constraints, and were less likely to cause head motion artifacts. We began with spontaneous breathing to observe each subject’s unique resting-state (natural) breathing patterns, and corresponding instantaneous CSF velocity waveforms. Since a key principle of a regular yogic breathing practice — that have been associated with health-benefits59,60,61 — is to make the breath slower and deeper, we included slow and deep breathing practices that were likely to have immediate impact on pulsatile CSF motion, and create larger changes compared to spontaneous breathing in magnitude and frequency of pulsatile CSF motion based on our pilot studies and literature review12,13,14. Based on our investigations, we hypothesized slow breathing would create an impact in between spontaneous and deep breathing. Therefore, we utilized slow breathing in between spontaneous and deep breathing measurements. Three-part breath components were utilized in the order they are performed in a traditional yoga practice: abdominal, diaphragm and chest. We utilized this same order for each subject—instead of randomizing—to prevent any “carrying over” effect from slow and deep breathing patterns to spontaneous breathing. This would allow us (i) to compare each subject’s unique spontaneous versus yogic breathing patterns, and corresponding CSF velocity waveforms (ii) to then quantify changes in magnitude and frequency components of CSF for identifying the primary driving force of CSF (respiratory versus cardiac components) during spontaneous versus yogic breathing practices.

While we herein are interested in the impact of different deep breathing practices on CSF velocities — performed as in traditional yogic practices —, it is important to note the physiology of different breathing practices: a recent study15 provided the influence of deep abdominal versus deep thoracic breathing while noting the two-breathing patterns exerting different muscle groups. During abdominal breathing, the diaphragm is utilized as the inspiratory muscle, and changes in intrathoracic volume and intrathoracic pressure are greater during abdominal versus thoracic breathing. Particularly, more pronounced contracting of the diaphragm during abdominal breathing results in greater opening of costodiaphragmatic recess.

Breathing rate and depth

At the core of yogic breathing practices lies awareness and training of the breath. We designed our RCT yogic breathing interventions from Raja Yoga69 practices in the Himalayan Tradition, in which yogic breathing is suggested to be performed within each person’s own capacity for safety reasons, with inhale/exhale to be extended and expanded with caution through regular long-term practice. With that goal in mind, for the MRI breathing protocol, we specifically avoided enforcing any specific rate or depth for inhale/exhale other than giving a choice of rate (e.g., 3 to 5 counts with a count rate of 1/sec).

Subject preparation in the MR instrument

After being introduced to the breathing techniques in the mock scanner, we transitioned subjects to a 3T MRI instrument (MAGNETOM Prisma, Siemens Healthineers, Erlangen, Germany) for baseline data acquisitions using a 64-channel head and neck coil. We positioned subjects in supine, and provided them with (i) a bolster placed under the knees, (ii) foam pads under the elbows, (iii) pads around head and neck for comfort and minimizing motion artifacts, (iv) blankets for warmth, (v) a wireless finger pulse sensor (Siemens Health) for pulse data collection, and (vi) a respiration bellow (Siemens Health) for respiration data collection during the entire RT-PCMRI data acquisitions. We instructed subjects to lie still in supine during the entire data acquisition.

Data acquisition

We utilized a 1-h data acquisition protocol, similar to our previous work14, consisting of anatomical MRI acquisitions, followed by simultaneous recordings of our previously established RT-PCMRI14 acquisition, respiration and finger pulse acquisitions. Briefly, for consistency across all subjects, we aimed to measure CSF at an angle perpendicular to the spinal cord at the level of the foramen magnum (FM) (Fig. 2A1–2 green lines). To determine the location of FM, we first collected anatomical MR images using a T2-weighted fast spin echo (HASTE, repetition time (TR) 1200 ms, echo time (TE) 80 ms; Fig. 2A1); a 3D T1-weighted gradient echo sequence (MPRAGE; TR 2300 ms; TE 2.32 ms; Fig. 2A2). We then acquired a cardiac gated PCMRI (TR 26.4 ms, TE 9.04 ms; Fig. 2A3) prior to RT-PCMRI to ensure proper slice location and angle for visibility of CSF pulsations, and that CSF was not obstructed. Upon confirming the pulsatile CSF motion (Fig. 2A3), we then acquired ~ 1-min RT-PCMRI (Fig. 2A4–5) at the same slice location and angle when subjects performed each of the breathing practices. RT-PCMRI sequence parameters included: velocity encoding value (VENC) 5 cm/s, temporal resolution ~ 55 ms, flip angle 30 degrees, matrix size 78 × 128, field of view (FOV) 196 × 323 mm (in-plane resolution ~ 2.5 × 2.5 mm), EPI factor 7, slice thickness 10 mm, TR 108.88 ms, TE 8.74 ms). RT-PCMRI has previously been described in detail14. During the RT-PCMRI acquisitions, we simultaneously collected respiration and pulse data with a sampling frequency of fs = 400 Hz.

Figure 2
figure 2

(A1–2) Sagittal anatomical MRI images showing the CSF measurement location at the level of foramen magnum (green lines) of a 37-year-old female. (A3) Axial images for cardiac gated PCMRI and RT-PCMRI velocity distribution. Cardiac-gated PCMRI is first collected for confirming CSF pulsation visibility prior to RT-PCMRI. ~ 1-min RT-PCMRI is then collected at the same location during five breathing patterns. (A5) A detailed image of region delineates the voxels of CSF and spinal cord (orange lines) and surrounding tissue. (B1) RT-PCMRI DICOM phase images (N = 1021) are collected for each breathing pattern repeated twice, resulting a total of 214,410 images processed for the n = 21 subjects, and 183,780 images utilized for the results of n = 18 subjects. (B2) Sample time series of CSF from single voxel RT-PCMRI (2.5 mm × 2.5 mm). Respiration and pulse data were simultaneously recorded, and temporally registered with the RT-PCMRI time series. (CD) Time and frequency domain analysis of five breathing conditions: spontaneous breathing (SponB), slow breathing (SlowB), deep abdominal breathing (DAB), deep diaphragmatic breathing (DDB), and deep chest breathing (DCB) (with the last three forming a specific yogic breathing technique called three-part breath; see Table 2). When compared to SponB, both time domain maximum (positive; cranially directed) instantaneous CSF velocity values (in C), and peak respiration frequency amplitudes (in D) increase during SlowB, DAB, DDB, and DCB. Results are produced with MATLAB ( and presented using Adobe Illustrator (

Data processing

Each RT-PCMRI acquisition series produced 2042 images (1021 magnitude and 1021 phase, Fig. 2B1). In total (for N = 21 subjects, five breathing conditions (SponB, SlowB, DAB, DDB, DCB) each repeated twice) we have acquired 428,820 RT-PCMRI (magnitude and phase) images, and processed the needed 214,410 RT-PCMRI phase images for obtaining CSF velocity time series. We have developed a semi-automated protocol for post-processing all MRI DICOM images, and respiration and pulse data time series using MATLAB software packages [MATLAB, R2019-2020. Natick, Massachusetts: The MathWorks Inc.].

CSF ROI and velocity waveforms

A common method30,32 in conventional PCMRI studies to obtain CSF velocity time series is to average CSF across all voxels within the outlined region of interest (ROI), which may potentially cause spatial noise due to border zone partial volume effects71. Achieving a high temporal resolution for RT-PCMRI further may reduce the spatial resolution compared to conventional PCMRI, for which we previously developed a correlation mapping technique that allowed us to extract and average only highly correlated CSF voxels for an averaged CSF velocity time series. In this study, we are interested in obtaining and comparing true spatial and temporal CSF velocity values [cm/s] for each breathing practice. Therefore, to capture true spatial peak velocities to our best ability, we utilized a 2-step process to evaluate CSF velocity waveforms at a single voxel9 (Fig. 3). We first extracted highly correlated CSF voxels (greater than 0.7 correlation coefficient14) with our previously developed correlation mapping technique14, and then visually compared the CSF ROI voxels on RT-PCMRI images (spatial resolution of 2.5 × 2.5 mm) with the cardiac-gated PCMRI images (higher spatial resolution of 0.625 × 0.625 mm) to confirm the location of a single voxel of interest within CSF ROI. We are interested in maximum capacity of (participant) breathing impact on CSF. Anterior CSF velocities were usually larger than posterior velocities across our study population. We selected a single anterior voxel, within the CSF space, with greater velocity, which was usually among the highest correlated voxels (greater than 0.9 correlation coefficient) obtained from the correlation mapping technique. In addition, to confirm deep breathing practices did not cause artifacts in velocity values, potentially by B0 field changes in the head caused by the motion of torso during deep breathing, we computed CSF velocities in a set of voxels within static tissue, and confirmed there were no respiratory or cardiac frequency components (Fig. S1).

Figure 3
figure 3

RT-PCMRI CSF ROI single voxel selection. (A) Sample conventional cardiac-gated PCMRI image (spatial resolution of 0.625 mm × 0.625 mm) showing the CSF region of interest (ROI) within orange circles. (B) Increasing temporal resolution for RT-PCMRI reduces the spatial resolution (2.5 mm × 2.5 mm). To capture true spatial velocity peak values, we implemented a 2-step process. Utilizing our previously developed correlation mapping technique, we first computed highly correlated CSF ROI voxels (greater than 0.7 correlation coefficient) — e.g., 33, 44, 48, 57 as labeled within each ROI. We then compared PCMRI image (A) with RT-PCMRI correlation map (B) to visually confirm the location of voxels. We are interested in maximum capacity of breathing impact on CSF. Anterior CSF velocities were usually greater than posterior velocities across our study population. Therefore, we have selected a single anterior voxel, within the CSF space and not contaminated by partial volume effects, and with greater velocity, which was usually among the highest correlated voxels obtained from the correlation mapping technique (greater than 0.9 correlation coefficient) such as voxel #44 for this sample. (C) For instance, while voxels #33, 48, and 57 are highly correlated one another, partial volume effects due to spinal cord or outside the subarachnoid space, the velocity waveforms and/or peak values are not truly preserved compared to voxel #44 which is within the CSF space. Results are produced with MATLAB ( and presented using Adobe Illustrator (

Time and frequency domain signals of interest

Previous studies9,13,14,37,48 reported vasomotion, respiration, and cardiac (1st harmonic) components of CSF signals. We observed (Fig. S2) higher order harmonics of cardiac pulsations in our preliminary analysis of frequency domain CSF velocity signals, which provides important information for determining the mechanisms, and their relative contribution to pulsatile CSF velocities. Having observed 1st and 2nd cardiac harmonics but not 3rd or 4th harmonics in all subjects, we have included 2nd cardiac harmonics in our analysis. In short, we are interested in four distinctive CSF velocity time series: instantaneous (iCSF), respiratory (rCSF), and cardiac 1st (c1CSF) and 2nd harmonics (c2CSF), which were below 4 Hz. To remove higher harmonics and high frequency signals, and investigate only the rCSF, respiration, c1CSF and c2CSF, we then low-pass filtered raw CSF velocity time series (see Fig. S3) using a 4th order Butterworth filter with a cut-off frequency of 4 Hz, which provided time domain iCSF velocity signals in Fig. 2C1–5. We then computed frequency domain signals using a fast Fourier transform (Fig. 2D1–5). See Fig. S4 for sufficient resolution of signals in Fig. 2C1–5, D1–5.

To separate and investigate rCSF, c1CSF and c2CSF (Fig. 4), we first computed frequency bands of each of the three components for each subject and breathing condition. This allowed us to take into consideration the individual and unique variations of frequency bands for accurate time domain velocity waveforms as well as frequency domain power calculations for each of the three components. We then filtered instantaneous CSF velocity waveforms (Fig. 4A1), using the individual frequency bands (Fig. 4A2), in (i) respiration frequency band of (estimated as typically f < ~ 0.6 Hz band), and (ii) cardiac 1st harmonic frequency band (estimated as typically ~ 0.6 < f < ~ 1.6 Hz band), and (iii) in cardiac 2nd harmonic frequency band (estimated as typically ~ 1.6 < f < ~ 2.7 Hz band) (Fig. 4A3–5). We repeated the above procedure for each breathing condition of each subject, and obtained time and frequency domain signals for all four distinctive CSF velocity waveforms.

Figure 4
figure 4

Separation of instantaneous CSF (iCSF) velocity waveforms (measured during deep abdominal breathing (DAB) for a 37 y–o female) into three components : respiratory (rCSF), cardiac 1st harmonic (c1CSF) and cardiac 2nd harmonic (c2CSF). (A1) Time domain iCSF velocity waveforms. (A2) Frequency domain iCSF obtained from Fast Fourier Transform (FFT) of iCSF velocity time series presenting peak frequency amplitudes for respiration (rpfa), cardiac 1st harmonic (c1pfa) and 2nd harmonic (c2pfa). We filtered iCSF within the bandwidth of each component — respiratory (magenta), cardiac 1st harmonic (light green) and 2nd harmonic (dark green) — in order to obtain and compare the characteristics of individual pulsatile velocity time series (A3) rCSF, (A4) c1CSF, and (A5) c2CSF. We computed maximum (positive; cranially directed) and minimum (negative; caudally directed) CSF velocity value in two ways: (i) peak value obtained as the highest and lowest point (blue solid line) during the entire time series, and (ii) averaged peak value (orange dash line) obtained by averaging the local maximums (and minimum, respectively) above a threshold set at 97.5th percentile (and 2.5th percentile, respectively). We used averaged peak values in statistical analysis to reduce temporal noise due to any transient events that may cause abrupt peaks (e.g., unexpected deep sigh) or lower peaks (indicating less than maximum capacity) caused by “participant fatigue” among novices during deep breathing conditions. Results are produced with MATLAB ( and presented using Adobe Illustrator (

In parallel, to confirm that estimated frequencies of CSF signals match with the physiological data, we filtered respiratory sensor and pulse sensor data in the same frequency band of respiratory and cardiac components of CSF velocity waveforms. For visualization purposes, we arbitrarily scaled the respiration and pulse data to compare with CSF velocity waveforms (Fig. 2C–D1–5). We confirmed the respiratory component in respiration data, and cardiac (1st and 2nd) harmonic components in pulse data (blue and purple lines in Fig. 2D1–5).

CSF metrics

Time domain CSF metrics

We computed the following metrics from time domain CSF velocity waveforms iCSF, rCSF, c1CSF, c2CSF for each subject and breathing condition: (1) peak13,14,37 maximum (cranially directed value at the highest point) and peak minimum (caudally directed at the lowest point) during ~ 1-min time series; (2) averaged peak maximum and averaged peak minimum obtained from the average of local peak maximums (and peak minimums, resp.) above (and below, resp.) a threshold set at 97.5th percentile72,73 (and 2.5th percentile, resp.); (3) range32 of peak maximum to peak minimum; (4) range of averaged peak maximum to averaged peak minimum; (5) displacement computed from the integration of the CSF velocity time series [cm/s] and converted to [mm]. Lastly, we computed % change in these metrics, from SponB to SlowB, DAB, DDB, and DCB.

Note that the traditional method for computing cranially- and caudally-directed velocities are to compute peak maximum and peak minimum values. In addition to peak maximum and minimum, for this study, we also computed averaged peak maximum and averaged minimum values. See Fig. 4 blue and orange dash lines for a comparison of peak versus averaged peak values. Since our goal during each breathing condition is to capture true maximum capacity of CSF velocity, the use of averaged peak approach allowed us to reduce temporal noise caused by (i) random transient events that are not part of the regular breathing pattern (e.g., unexpected deep sigh) resulting in greater peak values, or (ii) “participant fatigue” experienced — among novice practitioners — while performing slow and/or deep breathing conditions resulting in lower peak values. It is possible for novice practitioners to experience “fatigue” resulting in brief pause and/or reduced capacity of respiratory movements due to novelty of the practices. The longer and more frequently novices perform these breathing practices, the less likely they would experience “fatigue” as in meditation and other mindfulness practices74.

Frequency domain CSF metrics

From the instantaneous CSF velocity waveforms, we have computed (1) peak frequencies, and (2) peak frequency amplitudes (Fig. 4A2). Additionally, to observe individual peak frequency amplitude ratio changes, we computed (3) a peak-to-peak frequency amplitude ratio [r/c1peak] and [r/c2peak] calculated as the ratio of rCSF peak frequency amplitude to the c1CSF, and c2CSF peak frequency amplitudes in the frequency domain (Fig. 5A1–5).

Figure 5
figure 5

Sample datasets from a 37 y–o female presenting all four distinctive CSF signals (black: iCSF, magenta: rCSF, light green: c1CSF, dark green: c2CSF) during SponB, SlowB, DAB, DDB, and DCB. (A1–5) Frequency domain iCSF signals presenting the changes in peak frequencies [x-axis: rpf, c1pf, and c2pf; Hz] and peak frequency amplitudes [y-axis; rpfa, c1pfa, and c2pfa; a.u]. There was a decrease in respiration peak frequency during SlowB compared to SponB, and increase in respiration peak frequency amplitude rpfa during all four yogic breathing techniques (A2–5) due to increased respiratory movement. We computed peak frequency amplitude ratios [r/c1peak] and [r/c2peak] (e.g., [r/c1peak] indicated by sample purple arrows) for observing the changes, and power ratios [r/c1power] and [r/c2power] for testing the relative contribution of each component to instantaneous CSF. (B1–5) Time domain CSF signals presenting an increase in both cranially directed iCSF and rCSF velocities during four yogic breathing techniques compared to SponB; with detailed waveforms presented in [0–10] s time window in (C1–5). Results are produced with MATLAB ( and presented using Adobe Illustrator (

Relative contribution of rCSF, c1CSF, c2CSF signals

To compare the contribution of rCSF, c1CSF, and c2CSF and determine the primary regulatory force(s) for pulsatile CSF, we computed (1) estimated frequency band of rCSF, c1CSF, and c2CSF (Figs. 4A2, 5A1–5), (2) power of rCSF, c1CSF, and c2CSF, and (3) relative contribution13,14,37 of the respiration versus cardiac components by defining power ratio14 [r/c1power] and [r/c2power] calculated as the ratio of the power of the rCSF to the power of c1CSF and c2CSF, respectively.

To compute the power of rCSF, c1CSF, and c2CSF (defined14 as the integral of the square of the amplitude spectrum over the corresponding frequency band), we used trapezoidal numerical integration in frequency given by Eq. (1).

$${\int }_{a}^{b}S\left(f\right)df=\frac{b-a}{2N}\sum_{i=1}^{N}(S\left({f}_{i}\right)+S({f}_{i+1}))$$


where \(a\) and \(b\) are the first and last frequencies of the estimated frequency band (of signal of interest, e.g., rCSF, c1CSF, and c2CSF), (\(b-a)\) is the frequency bandwidth, \(S\) is the square of the signal amplitude spectrum, \(N\) is the number of samples within the frequency band, and (\(b-a)/N\) is the spacing between samples, i.e., \(df\) frequency resolution. We then computed [r/c1power] as the ratio of rCSF power to the c1CSF power, and [r/c2power] as the ratio of rCSF power to the c2CSF power.

Statistical analysis

To test the differences between spontaneous and four yogic breathing techniques, we used mean and standard deviation (SD) for the following time and frequency domain CSF metrics (i) averaged peak maximum and minimum values for iCSF, rCSF, c1CSF, c2CSF; (ii) range of iCSF averaged peak maximum and minimum values, (iii) iCSF displacement, (iv) peak frequencies, (v) peak frequency amplitudes, (vi) power values for rCSF, c1CSF, c2CSF, (vii) frequency peak-to-peak amplitude ratios, and (viii) power ratios for a total of 23 metrics (Table S2).

Independent statistical analysis was conducted by A.H. in R version 4.0.3 (R Core Team, Vienna, Austria). Though data were visually inspected, and extreme values were double checked based on percentiles (with a lower threshold of 5th percentile and a higher threshold of 95th percentile), all values remained in the data to preserve the—diversity of—differences in CSF velocities. Mixed effects linear regression models were built to analyze the associations between each outcome measure and each of the four experimental breathing conditions (SlowB, DAB, DDB, and DCB). Each model included a random subject effect to characterize within-person correlations over repeated measures. Normality assumptions of each model were checked by visual inspection of Q–Q Plots. For multiple comparisons, type I error rate was controlled by using the Benjamini–Yekutieli false detection rate (FDR) procedure75, with an overall FDR of 0.05, using the R function “p.adjust”. Unadjusted p values and FDR corrected p values are provided in Table S2, and p values mentioned within the text are FDR corrected p values.

In addition, associations between demographic covariates and outcomes were inspected visually, and tested for significance by Spearman’s rank-order correlation (continuous covariates) or t-test (dichotomous covariates) if a possible association was seen. We reported associations found to be significant in “Results” section. However, the regression models of the main analyses were not adjusted for these covariates in order to maintain consistency between models and so as not to overfit the models. Note this study was designed and powered to detect changes in CSF metrics after a yogic breathing intervention, rather than the pre-intervention baseline analysis that we present here. For full clinical trial sample size calculations, see Supplemental Materials.

Ethics approval and consent to participate

The study was approved by the Institutional Review Board of the Oregon Health & Science University (OHSU). The clinical trial was registered at the (ID # NCT03858309, February 28 2019, We received verbal and written informed consent from all study subjects prior to all study procedures. All the methods were carried out in accordance with relevant guidelines and regulations.

Source link