- 1Department of Bioengineering, University of Pittsburgh, Pittsburgh, PA, United States
- 2Department of Psychiatry, University of Pittsburgh, Pittsburgh, PA, United States
- 3Department of Radiology, University of Pittsburgh, Pittsburgh, PA, United States
- 4Department of Anesthesiology and Perioperative Medicine, University of Pittsburgh, Pittsburgh, PA, United States
The development of innovative non-invasive neuroimaging methods and biomarkers is critical for studying brain disease. Imaging of cerebrospinal fluid (CSF) pulsatility may inform the brain fluid dynamics involved in clearance of cerebral metabolic waste. In this work, we developed a methodology to characterize the frequency and spatial localization of whole brain CSF pulsations in humans. Using 7 Tesla (T) human magnetic resonance imaging (MRI) and ultrafast echo-planar imaging (EPI), in-vivo images were obtained to capture pulsations of the CSF signal. Physiological data were simultaneously collected and compared with the 7 T MR data. The primary components of signal pulsations were identified using spectral analysis, with the most evident frequency bands identified around 0.3, 1.2, and 2.4 Hz. These pulsations were mapped spatially and temporally onto the MR image domain and temporally onto the physiological measures of electrocardiogram and respiration. We identified peaks in CSF pulsations that were distinct from peaks in grey matter and white matter regions. This methodology may provide novel in vivo biomarkers of disrupted brain fluid dynamics.
1 Introduction
Clearance and exchange of brain fluids promote brain health by removing neurotoxic metabolic byproducts from the brain such as amyloid beta and tau (Nedergaard, 2013; Xie et al., 2013). Clearance of brain fluids is driven by pulsations of the perivascular spaces from autonomic nervous system (ANS) activity, which may vary as a function of brain states such as sleep and wakefulness (Herculano-Houzel, 2013; Xie et al., 2013; Hauglund et al., 2020) as well as brain diseases such as Alzheimer’s disease (Ramanathan et al., 2015; Peng et al., 2016) and major depressive disorder (Hock et al., 1998; Pomara et al., 2012).
Concurrently, the glymphatic system facilitates the convection of cerebrospinal fluid (CSF) between the peri-arterial and peri-venous spaces. This convective flow is thought to at least partially be driven by cardiac-induced blood flow pulsations along the arteries (Adolph et al., 1967; Schroth and Klose, 1992; Martin et al., 2012; Iliff et al., 2013). Although the specific mechanisms of clearance remain under debate, there is evidence that water is propelled by the arterial pulsations through aquaporin channels and supports solute transport from extracellular interstitial spaces, through perivascular spaces, and into CSF spaces. CSF and waste products from the brain are then pushed from parenchyma to subarachnoid spaces and may eventually be cleared via arachnoid granulations and dural and nasal lymphatic vessels (Rennels et al., 1990; Iliff et al., 2013; Kiviniemi et al., 2016; de Leon et al., 2017), or as recently discovered, the meningeal lymphatic vessels (mLVs) (Louveau et al., 2015; Ahn et al., 2019).
Development of quantitative CSF imaging methods is critical to understand factors that influence the CSF dynamics and brain fluid clearance. T1-weighted magnetic resonance imaging (MRI) with intrathecal injection of a gadolinium (Gd)-based contrast agent has been used to characterize CSF flow in human participants with idiopathic normal pressure hydrocephalus (iNPH) and dementia (Ringstad et al., 2017, 2018; Eide and Ringstad, 2019). This technique has afforded fully quantitative, high-resolution imaging of CSF and interstitial fluid (ISF) flow throughout the whole head but is highly invasive, as it requires a lumbar puncture. Gd may also remain deposited in the brain, limiting its longitudinal research utility (Gulani et al., 2017).
Fast acquisition functional magnetic resonance imaging (fMRI) paradigms have also been used to characterize CSF dynamics in iNPH and Alzheimer’s disease (AD) patients and in healthy control participants during sleep (Fultz et al., 2019; Shanks et al., 2019; Yamada et al., 2020; Yang et al., 2022). A resting-state fMRI study has evidenced a potential coupling between CSF pulsations and the global blood oxygen level-dependent (BOLD), which is notably reduced in patients with AD-related diseases, where the CSF region of interest (ROI) was considered at the bottom slice of the acquisitions and delineated using masks obtained from T2*-weighted fMRI and further confirmed with T1-weighted images (Han et al., 2021). While non-invasive, these sequences have relatively poor signal-to-noise ratio (SNR) and spatial resolution and have been limited to narrow fields of view encompassing only the 4th ventricular, bottom edge slice of the acquisitions or cerebral aqueduct. Moreover, similar fMRI techniques were employed comparing CSF pulsations in different regions such as edge slices and 4th ventricle, revealing distinct behaviors depending on the region selected for assessing the CSF dynamics (Kim et al., 2022).
The bulk changes in blood volume at the capillary level could cause widespread fluctuations of measured signal intensity with the cardiac cycle. Furthermore, large vessel pulsatility may cause tissue and CSF movement and production of an influx of unsaturated blood into the slice of interest affecting the measured signal intensity in the areas adjacent to the vessels. This leads to a signal variation when using echo-planar imaging (EPI) acquisitions (Dagli et al., 1999). Hence, fMRI and other techniques have been used to characterize different sources of pulsations in the brain (Poncelet et al., 1992; Biswal et al., 1995; Purdon and Weisskoff, 1998; Dagli et al., 1999; Kiviniemi et al., 2000). Thus, using MRI of CSF dynamics can inform the study of brain diseases and the role of sleep–wake states (Xie et al., 2013; Fultz et al., 2019).
Based on 3 Tesla (T) MRI magnetic resonance encephalography (MREG), it has been demonstrated that AD patients experience abnormalities in cardiovascular brain impulses, which can manifest as slow, fast, or even in a reverse direction of propagation (Rajna et al., 2021). Similarly, patients with epilepsy (Kananen et al., 2018, 2020, 2022) and narcolepsy (Järvelä et al., 2022) also have exhibited altered physiological brain pulsations. MREG imaging have also been used to identify changes in physiological brain pulsations during nonrapid eye movement (NREM) sleep (Helakari et al., 2022), and to identify the spatial location of brain physiological pulsations as well as the multiple sources of BOLD signal (Raitamaa et al., 2021). Despite being a non-invasive technique with high temporal resolution used to characterize alterations in brain pulsations dynamics based on fMRI measurements, MREG’s limited accessibility stands in contrast to the more broadly available fast EPI methods.
Alternative techniques for ultrafast EPI acquisitions have been proposed. These include inverse imaging (InI) (Lin et al., 2012), generalized inverse imaging (GIN) (Boyacioglu et al., 2013), and multi-slab echo-volumar imaging (multi-slab EVI or MEVI) (Posse et al., 2012, 2013). These techniques offer a fast-sampling rate and reduced sensitivity to physiological noise. However, they come with the trade-off of potential loss of spatial resolution or introduction of geometrical distortions (Lin et al., 2012; Posse et al., 2012, 2013; Boyacioglu et al., 2013). Moreover, these studies have primarily utilized 3T MRI scanners.
Ultra-high field MRI (≥ 7 T) provides a major advantage of increased SNR, which can be used either to increase the resolution of the images or to decrease the scanning time (with the use of higher acceleration factors). The 7 T field strength also has higher sensitivity to BOLD signal and better vasculature conspicuity (Moser et al., 2012; Santini et al., 2021b).
Using ultrafast EPI we acquired the CSF MR signal in real-time. We report CSF pulsation patterns through spectral analysis. We applied the same methodology across datasets of seven different participants to validate the observed results. To identify whether spectral peaks in CSF pulsatility aligned with changes in physiological measures thought to drive CSF flow, we simultaneously collected physiological measures of electrocardiogram (ECG) and respiration in one participant.
2 Methods
The overall design of this study is based on two main steps: (1) In-vivo 7 T image acquisition, concurrent with physiology measurements, and (2) image processing with spectral analysis. The processing and analysis of the power-frequency spectrum and its corresponding spatial mappings were fine tuned for detection of bandwidth and peak span, thresholding levels for masks, and smoothing degrees for filtering.
Participants provided informed consent as approved by the University of Pittsburgh’s Institutional Review Board (identification number PRO17030036). Seven healthy volunteers (two male, age range 26–30 years old and five female, age range 21–25 years old) were scanned to obtain EPI data. From one volunteer, we obtained two axial and one sagittal slices. Whole brain EPI data, including the cerebellum, were collected from five participants. For the last participant, we collected simultaneous EPI and physiological data.
2.1 Image acquisition
Images were acquired using a whole-body 7 T MRI system (Siemens 7 T MAGNETOM) and with the human-connectome EPI multiband MR sequence, which can be obtained from https://www.cmrr.umn.edu/multiband/index.shtml (Moeller et al., 2010; Uğurbil et al., 2013). The sequence obtains fast acquisitions, high signal contrast of the CSF pulsation, and high sensitivity to BOLD signal, thus being well-suited for studies of sleep and neurodegenerative and psychological disorders. The imaging was acquired with an in-house developed and fabricated head coil with a 16-channel Tic-Tac-Toe transmit array with a 32-channel receive insert (Santini et al., 2018; Krishnamurthy et al., 2019; Santini et al., 2021b) that is load insensitive (Ibrahim et al., 2008; Kim et al., 2016; Santini et al., 2021b) and capable of whole brain homogenous imaging at 7 T (Ibrahim et al., 2013). By using this coil design, we were able to acquire signal from the entire brain with minimal regions of significant excitation losses and using the single transmit mode of the 7 T scanner.
The acquired EPI images yield a real-time visualization of the brain pulsations. The sequence was optimized to perform fast brain imaging. For Volunteer 1, during the development of the protocol, the EPI acquisition was broken into two sequences: a 2-slice axial view, with echo time (TE) of 17 ms, repetition time (TR) of 102 ms, acceleration factor of 3, field of view (FOV) of 220 mm × 220 mm, resolution 1.5 mm × 1.5 mm × 3 mm, total acquisition time of 1 min and 8 s; and a 1-slice sagittal view, with TE of 18 ms, TR of 100 ms, acceleration factor of 3, field of view (FOV) of 216 mm × 216 mm, resolution 1.5 mm × 1.5 mm × 4.4 mm, total acquisition time of 1 min and 3 s. In both cases, a total of 600 volumes were sequentially acquired.
For five volunteers (Volunteers 2 to 6), whole-brain imaging was acquired. The main data acquisition was done with TE of 20 ms, TR of 155 ms, isotropic resolution of 2 mm, and acceleration factor of 2. The FOV was 192 mm × 192 mm. The acquisition was broken into 19 slabs of 3 axial slices each for a total of 57 slices, where each slice has a thickness of 2 mm, totaling 114 mm, providing a whole-brain coverage. A total of 600 volumes were sequentially acquired per slab in a single sequence run for an acquisition time of 1 min and 36 s per slab. For Volunteer 2, the EPI acquisition was done using TR of 152 ms but only 15 slabs of 3 axial slices for a total of 45 axial slices. Another EPI acquisition on the same volunteer was also performed using TR of 51 ms and a single slice.
In all cases, the TE values were chosen for future potential BOLD analysis. Furthermore, all EPI images were acquired from volunteers in a resting state condition. Moreover, spin-echo EPIs were also acquired for B0 field distortion correction with the same phase encoding (PE) direction of the EPI acquisition (PA) and with the opposite PE direction (AP). The sequence parameters included a TE of 39.4 ms, TR of 6,000 ms, with the other parameters – such as field of view, resolution, number of slices, echo spacing, and position - matched to those of the EPI sequence.
A T1-weighted imaging (MPRAGE) sequence was used for proper localization of the EPI field of view and as a structural scan for the image processing. This acquisition was done using 0.75 mm isotropic resolution, TR of 3,000 ms, TE of 2.17 ms, and 256 slices for a coverage of 240 mm x 173 mm x 192 mm in total time of acquisition of ~5 min.
2.2 Physiological measurements
Electrocardiogram (ECG) and respiration signals were collected for one participant (Volunteer 7) inside the MR scanner using MR compatible ECG leads and an expansion belt attached to the chest to track inflation and deflation of the chest during respiration. Acquisition was digitalized using BIOPAC system (ECG: Cardiology | Research | BIOPAC, n.d.; Respiration Transducer for MRI | TSD221-MRI | Research | BIOPAC, n.d.). The simultaneously collected data allowed temporal signal analysis of both MR and physiologic signals. The imaging data acquired in conjunction with the physiological data used a TR of 75 ms, TE of 28 ms, and 3 axial slices (4 mm-thick each) including the lateral ventricles.
2.3 Image processing
The processing pipeline was developed based on MATLAB (MATLAB - MathWorks, n.d.), ANTs (Avants et al., 2009), and FSL (Jenkinson et al., 2012) packages. It consisted of denoising, distortion correction, bias correction, and skull stripping of each dataset volume. The initial step was loading the slabs and merging them into a single dataset (necessary only for whole-brain coverage images). Next, denoising was performed using a noise estimation tool with variance stabilization transformation (VST) for Rician-distributed noise (Foi, 2011). The Rician heteroscedastic noise was converted to a homoscedastic noise after the forward VST. The block-matching 4D (BM4D) denoising algorithm (Maggioni et al., 2013) could then be applied and the denoised image is obtained after the inverse VST. This tool has been used for other MRI applications (Santini et al., 2021a) and yields a good result when applied to EPI data. Distortion correction was performed using the estimated B0 maps derived from the spin-echo sequence using the topup tool (Andersson et al., 2003) from FSL software package. The generated map was used for correction of the EPI data. Then, the images were bias corrected using the N4 (Tustison et al., 2010) tool from the ANTs software package with spline distance parameter of 200. The final skull stripping was performed using the FSL brain extraction tool (BET). We note the merging of images captured at different times, as each slab is acquired within distinct time frames. However, the transition of the analysis to the frequency domain, with the fast Fourier transform (FFT) executed for each voxel, mitigated these temporal differences. For the subsequent spectral analysis, the whole-brain coverage images were segmented using SynthSeg (Billot et al., 2023a,b), allowing the assessment of the brain pulsations in three different regions: cortical gray matter (cGM), cerebral white matter (cWM) and CSF.
2.4 Spectral analysis
The frequency analysis was performed for each dataset individually and resulted in both a frequency power spectrum and a mask for brain localization of specific frequency bands. A frequency spectrum was calculated for arbitrarily selected points for validation of the findings across different brain regions, for single slices (performed for images not covering the entire brain), and also for the regions described in the previous section (whole-brain, cGM, cWM and CSF).
After processing each individual EPI data, the frequency processing and analysis were performed in MATLAB and Python. The time series of each voxel was used to generate a frequency spectrum using FFT. With the 600 volumes of 152/155 ms TR acquisition, the frequency resolution of the frequency spectrum is 0.011 Hz and the maximum frequency is 3.289/3.226 Hz. The same frequency analysis for the 51 ms TR data produces a much larger frequency spectrum of up to 9.8 Hz. Therefore, frequency components higher than 3 Hz can be observed and analyzed. The analysis was done using both individual points manually selected as well as the average across the 2D or 3D space. For acquisitions encompassing a 2D/3D region of interest (ROI), the FFT was calculated for each voxel, and subsequently the spectral intensity values across all voxels were averaged for each frequency, resulting in a single intensity value per frequency for each 2D/3D ROI.
2.4.1 Spectral peaks identification
To identify the relevant peaks of each spectrum, an algorithm was developed in Python as depicted in Figure 1. Initially, each spectrum plot was smoothed using a Savitzky–Golay filter of order 2. The filter window was defined as 0.133 Hz, corresponding to the normal range of 12 to 20 breaths per minute in adults (Sapra et al., 2020). For a 51 ms TR, the number of samples of this window was doubled.
Figure 1. Procedure to identify peaks for each spectral analysis: (A) original FFT obtained, (B) smoothed FFT after the application of Savitzky–Golay filter, (C) obtaining the inferior envelope, (D) adjusted FFT obtained after dividing each value of the smoothed FFT by the inferior envelope, (E) identification of the most evident peaks, with a height higher than a threshold defined as the sum of the mean and standard deviation of the adjusted FFT, and (F) peaks identified in the adjusted FFT with a minimum height given by the sum of the mean and standard deviation of the regions not excluded in the previous plot.
The lower envelope observed across all spectra, which is a characteristic of the acquisition as evidenced by tests with a phantom (provided in the Supplementary material), was modeled by identifying and accumulating the minimum values along the frequency axis. Whenever the observed values exceed these minimums, indicating an increase in the signal, linear interpolation was employed to bridge between the identified minimum value and the subsequent minimum value observed when the signal decreased again. The smoothed spectrum was then divided by this envelope to reflect the relative increase. By subtracting one unit from this normalized spectrum, we obtained an adjusted spectrum optimized for peak identification.
To establish a baseline level, the most evident peaks in each adjusted spectrum were identified first. We calculated the mean and standard deviation of the adjusted spectrum, and their sum was used as the minimum threshold for detecting primary peaks (findpeaks implementation in Python). An exclusion zone of 0.667 Hz was applied around each identified peak, where this frequency range corresponds to the normal cardiac range of 60 to 100 beats per minute in healthy adults (Sapra et al., 2020). Regions in the adjusted spectrum outside the exclusion zones were assessed to establish a baseline level.
A spectrum peak was identified as those with a minimum height equal to the sum of the mean and standard deviation of the baseline level. To prevent small oscillations above this height from being classified as peaks, a minimum prominence of 1.5 dB was required.
2.4.2 Spectral spatial analysis
Spatial analysis was done by creating image masks based on the localization of voxels with significant signal in each frequency band. It was also performed for single-slice acquisitions, covering the 2D space acquired. Moreover, due to whole-brain coverage, a similar approach was employed for the entire brain images and specific brain regions (cGM, cWM and CSF).
Regarding the voxel analysis, the power map of a given frequency band was determined voxel-wise by averaging respective power values within the frequency band. For better visualization, each power map was then binarized with a chosen threshold (75% of the peak amplitude of the corresponding frequency band) and spatially smoothed using a Gaussian filter (sigma of 1.6), generating the final masks for each frequency band. These masks were overlaid on the original EPI and T1 weighted acquisitions for anatomical reference. The T1 weighted image was registered with the average EPI image of the dataset using SPM12.
For single-slice, whole-brain and brain regions (cGM, cWM and CSF), peaks were identified as shown in the previous section. For each of these peaks, the bandwidth corresponding to a 3 dB drop ( of the peak magnitude) and the area under the spectrum within this bandwidth were determined, aiming to identify potential biomarkers.
3 Results
A video was created based on the image series of the fast EPI data after acquisition and processing in Volunteer 1 (Figure 2). The video visually indicates the presence of periodical signal from the CSF flow. For this first volunteer, the acquired image covers only a few slices. However, it was still possible to apply the developed methodology to identify relevant frequencies, demonstrating the potential of using this approach even with single-slice acquisitions.
Figure 2. Fast EPI acquisition (TR = 102 ms for axial view, TR = 100 ms for sagittal view) for Volunteer 1 showing signal changes due to CSF flow; axial slices with a spatial resolution of 1.5 × 1.5 × 3 mm and a sagittal slice with spatial resolution of 1.5 × 1.5 × 4.4 mm. The blue arrows point to regions of large variation in signal over time. A video showing these pulsations in real time is available at doi: 10.6084/m9.figshare.24022932. For each view, a single-slice FFT spectrum is shown, where the peaks were identified following the procedure proposed in this study.
To confirm the presence of physiological signals such as respiration and cardiac motion, CSF temporal data was aligned with measurements from the electrocardiogram and respiration belt in Volunteer 7 for visual comparison of similarity between the physiological activities and the change in signal intensity from CSF regions (Figure 3A). The frequency spectrum of the datasets was also aligned following the same comparison as the time series data (Figure 3B). The two major signal bands were highlighted between the ECG and CSF data (around 1.1 Hz) and the respiration belt and CSF data (around 0.3 Hz).
Figure 3. MRI data acquisition along the cardiac cycle and respiration for Volunteer 7. Echo-planar imaging acquisition performed with concurrent physiological measurement of electrocardiogram in the 7 T scanner. (A) Time series of ECG, respiration belt, and EPI signals. The EPI signal was temporally aligned with the physiological data using an external trigger signal from the scanner. The red lines represent the R-peaks of the ECG. (B) Frequency spectrum of the ECG, respiration belt, and EPI signals. The purple region highlights the common frequency band between the respiration and EPI signals, and the orange region highlights the common frequency band between the ECG and EPI signals.
To verify that various points of the brain contribute differently on the frequency spectrum, Figure 4 represents the frequency spectrum for 9 arbitrary points throughout the brain. The position of each point is described by the brain anatomy that it belongs to as shown on the top-right corner of each spectrum graph. Most of the points show frequencies around 1.2 Hz. Depending on the position, the signal shows the 0.3 Hz and/or the 2.4 Hz bands.
Figure 4. Frequency spectrum for nine selected points (A–I) throughout the brain for Volunteer 2. Some points show higher intensity on the 1.2 and 2.4 Hz bands (points A,B,E,F,G–I) whereas other show more on the 0.3 Hz band (points C,D,G). The labels describe the brain anatomy where the data was obtained.
A frequency analysis was obtained from the EPI data for the single-slices acquisitions, whole-brain coverage, and in the following brain regions: cGM, cWM and CSF. Based on a developed methodology to identify relevant peaks (Figure 1), the most significant frequency bands were highlighted, and for each of them a peak magnitude, bandwidth for a 3 dB drop and area under this range were obtained (Figure 5). For all volunteers with whole-brain coverage, bands with similar center frequencies of approximately 0.3, 1.2, and 2.4 Hz can be identified. Table 1 shows the center frequency for the bands calculated for each volunteer. The frequency bands with centers at 0.3 and 1.2 Hz closely approximate the respiration and cardiac frequencies of a human adult (around 18 breaths per minute and 72 heart beats per minute, respectively). These bands can be identified as Band 1 and Band 4 on Table 1. In terms of magnitudes discovered within these frequency bands for CSF regions, their averages are 0.594 ± 0.378 (a.u.) and 1.747 ± 0.475 (a.u.), respectively. For the areas, their averages are 0.027 ± 0.017 (Hz) and 0.119 ± 0.037 (Hz), respectively. Finally, for the bandwidths, their averages are 0.053 ± 0.006 (Hz) and 0.078 ± 0.016 (Hz), respectively.
Figure 5. Frequency spectra and identified peaks in different participants and brain regions: (A) frequency spectra for the whole brain, cGM, cWM, and CSF regions for Volunteer 2, (B) frequency spectra of CSF across five different volunteers, (C) frequency bands identified in CSF regions for Volunteer 2, (D) frequency bands identified in the whole-brain coverage for Volunteer 2, (E) frequency bands identified in cGM for Volunteer 2, and (F) frequency bands identified in cWM for Volunteer 2.
Table 1. Values of identified frequency bands and respective magnitudes, areas under the curve and bandwidths for 3 dB drop in whole-brain coverage imaging.
The masks created per frequency band (Figure 6) show a spatial localization for the frequency band centered at lower frequencies, e.g., 0.387 Hz (Figure 6A) and 1.151 Hz (Figure 6C) overlapping with brain regions with larger volume of CSF (the main ventricles and cerebral aqueduct). Similar patterns were observed for all volunteers as shown in Figure 7 as the mask for the heart rate band is demonstrated in each of the volunteer’s data.
Figure 6. Frequency spectrum of Volunteer 4 with spatial localization of the signal from 4 separate frequency bands. For each band, a spatial mask was applied to the T1 weighted image. The acquired data can be obtained from the inferior region of the brain (cerebellum) up to the middle of the brain. The bandwidth for each band is defined as those that result in a 3 dB drop. The center frequencies are (A) 0.387 Hz, (B) 0.763 Hz, (C) 1.151 Hz, and (D) 2.355 Hz. The acquisition was done using an EPI sequence with TR = 155 ms with 19 slabs of 3 slices each for a total of 57 slices.
Figure 7. Visualization of the mask created for each volunteer at approximately the same position in the brain (bottom of the brain and top of the cerebellum and at approximately the same frequency band); (A) for Volunteer 2 at 1.151 Hz; (B) for Volunteer 3 at 1.301 Hz; (C) for Volunteer 4 at 1.151 Hz; (D) for Volunteer 5 at 1.129 Hz; and (E) for Volunteer 6 at 1.226 Hz.
For the larger frequency spectrum (dataset with TR of 51 ms), extra bands can be identified, and the center of one of the most prominent band was calculated at around 3.5 Hz (Figure 8).
Figure 8. Frequency spectrum for Volunteer 2 done using an EPI sequence single slice with TR = 51 ms. Maximum frequency of 9.8 Hz and frequency bands were highlighted based on the procedure developed in this work.
4 Discussion
We demonstrated a method to analyze the physiological brain pulsations in the human brain in-vivo using ultrafast 7 T EPI acquisitions, and an analytic approach to examine the dynamics of regional CSF signal pulsations. The raw visualization of the real-time signal (Figure 2) shows in-vivo CSF pulsations. The flow of CSF within the ventricles and in the subarachnoid space can be clearly visualized with changes in signal intensity. The time series and the frequency spectrum comparison between the collected physiological data and the EPI data (Figure 3) shows a direct alignment between the two types of data where the cardiac and respiration cycles can be observed in the EPI MRI data. The frequency analysis also shows consistent results across multiple volunteers, with similar frequency spectrums are observed.
Compared to previous studies (Fultz et al., 2019; Shanks et al., 2019; Yamada et al., 2020; Yang et al., 2022), we were able to achieve a greater frequency range, improved frequency resolution, and full brain coverage. This was made possible by utilizing ultrafast acquisition times (ranging from 51 to 155 ms), which allowed for whole-brain spectral analysis up to 9.8 Hz. Additionally, we employed high SNR and homogeneous images by using a 7 T MRI with a customized radiofrequency (RF) coil system (Santini et al., 2018; Krishnamurthy et al., 2019; Santini et al., 2021b). To optimize sequence parameters, we tailored the flip-angle to maximize the signal of the CSF flow and adjusted other parameters (e.g., acceleration factor) to minimize susceptibility-related distortions. We also selected a TE that could potentially capture the BOLD signal if functional connectivity data are warranted.
Based on the methodology implemented in this work, we can compare physiological brain pulsations across different tissue types and brain regions. In the frequency domain, there is evidence of more well-defined peaks in CSF regions than in other regions, with cGM following as the second most defined, and cWM as the third. This suggests a potential relationship between these measurements and CSF dynamics. Moreover, when comparing these peaks in terms of their respective magnitudes, areas under the curves, and bandwidths for a 3 dB drop, we observe low variability across different volunteers for similar segments, indicating potential biomarkers for brain pulsations analysis.
The creation of frequency masks allowed for an analysis of the spatial localization of each frequency band. The presence of heart rate frequencies (1.2 Hz) in the ventricles validates the analysis as the arterial pulse wave in the choroid plexus, for instance, is known to influence the CSF motion (Bilston et al., 2010; Martin et al., 2012; Iliff et al., 2013). Additionally, the presence of high frequencies (over 2 Hz) responses can suggest a more turbulent flow that also aligns with regions of the main cerebral aqueduct. This work provides a basis for identifying new biomarkers for brain fluid dynamics. For example, the frequency spectrum can be analyzed for different brain diseases. The lower frequency bands (below 1 Hz) contain physiological signals that corelate with the heart rate and breathing, so brain conditions that affect those variables can be analyzed directly from the MRI data. The magnitude of each band may also provide insights into the coupling between the heart and breathing rates with the CSF pulsations, which may correlate with clearance rate. On the other hand, the higher frequency bands (above 1.8 Hz), can be correlated with sleep cycles and potential sleep studies (Xie et al., 2013).
5 Conclusion
The development of non-invasive neuroimaging methods and biomarkers of brain fluid dynamics is essential for studying brain diseases. This work presents a novel methodology to characterize the frequency spectrum and spatial localization of CSF pulsations in the human brain. The use of ultrafast EPI in conjunction with 7 T human MRI and simultaneous collection of physiological data enabled the identification of primary components of CSF pulsations and their mapping spatially and temporally onto the MR image and physiological domains. The methodology showed low variability and repeatability in-vivo, making it a promising tool for potential studies of brain fluid dynamics and CSF flow. Future studies will explore this methodology in clinical studies to determine its implications for the diagnosis and treatment of brain diseases.
Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.
Ethics statement
The studies involving humans were approved by University of Pittsburgh’s Institutional Review Board (identification number PRO17030036). The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.
Author contributions
TM: Conceptualization, Formal analysis, Investigation, Methodology, Software, Writing – original draft, Writing – review & editing. BD: Investigation, Methodology, Software, Writing – review & editing, Formal analysis, Visualization. MW: Conceptualization, Methodology, Software, Writing – original draft, Writing – review & editing. KW: Conceptualization, Writing – original draft, Writing – review & editing. DM: Writing – original draft, Writing – review & editing, Conceptualization. JI: Conceptualization, Methodology, Investigation, Writing – original draft, Writing – review & editing, Resources. HA: Conceptualization, Funding acquisition, Writing – original draft, Writing – review & editing. TS: Formal analysis, Supervision, Writing – original draft, Writing – review & editing, Conceptualization, Investigation, Methodology, Software. TI: Conceptualization, Funding acquisition, Investigation, Methodology, Resources, Supervision, Writing – original draft, Writing – review & editing.
Funding
The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This work was supported by NIH R01AG063525, R01MH111265, R56AG074467, and P01AG025204. The author TM was partially supported by the CAPES Foundation, Ministry of Education of Brazil, 88881.128222/2016-01. This research was supported in part by the University of Pittsburgh Center for Research Computing, RRID:SCR_022735, through the resources provided. Specifically, this work used the H2P cluster, which is supported by NSF award number OAC-2117681.
Acknowledgments
The authors would like to thank Dr. Marcelo AC Vieira and Vinicius P Campos for their support in the imaging denoising process. Additionally, we express our appreciation to the reviewers for their insightful feedback.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins.2024.1305939/full#supplementary-material
References
Adolph, R., Fukusumi, H., and Fowler, N. (1967). Origin of cerebrospinal fluid pulsations. Am. J. Physiol. Legacy Content 212, 840–846. doi: 10.1152/ajplegacy.1967.212.4.840
Ahn, J. H., Cho, H., Kim, J.-H., Kim, S. H., Ham, J.-S., Park, I., et al. (2019). Meningeal lymphatic vessels at the skull base drain cerebrospinal fluid. Nature 572, 62–66. doi: 10.1038/s41586-019-1419-5
Andersson, J. L. R., Skare, S., and Ashburner, J. (2003). How to correct susceptibility distortions in spin-echo echo-planar images: application to diffusion tensor imaging. NeuroImage 20, 870–888. doi: 10.1016/S1053-8119(03)00336-7
Avants, B. B., Tustison, N., and Song, G. (2009). Advanced normalization tools (ANTS). Insight J 2, 1–35. doi: 10.54294/uvnhin
Billot, B., Greve, D. N., Puonti, O., Thielscher, A., Van Leemput, K., Fischl, B., et al. (2023a). SynthSeg: segmentation of brain MRI scans of any contrast and resolution without retraining. Med. Image Anal. 86:102789. doi: 10.1016/j.media.2023.102789
Billot, B., Magdamo, C., Cheng, Y., Arnold, S. E., Das, S., and Iglesias, J. E. (2023b). Robust machine learning segmentation for large-scale analysis of heterogeneous clinical brain MRI datasets. Proc. Natl. Acad. Sci. 120:e2216399120. doi: 10.1073/pnas.2216399120
Bilston, L. E., Stoodley, M. A., and Fletcher, D. F. (2010). The influence of the relative timing of arterial and subarachnoid space pulse waves on spinal perivascular cerebrospinal fluid flow as a possible factor in syrinx development: laboratory investigation. J. Neurosurg. 112, 808–813. doi: 10.3171/2009.5.JNS08945
Biswal, B., Yetkin, F. Z., Haughton, V. M., and Hyde, J. S. (1995). Functional connectivity in the motor cortex of resting human brain using echo-planar MRI. Magn. Reson. Med. 34, 537–541. doi: 10.1002/mrm.1910340409
Boyacioglu, R., Beckmann, C. F., and Barth, M. (2013). An investigation of RSN frequency spectra using ultra-fast generalized inverse imaging. Front. Hum. Neurosci. 7:156. doi: 10.3389/fnhum.2013.00156
Dagli, M. S., Ingeholm, J. E., and Haxby, J. V. (1999). Localization of cardiac-induced signal change in fMRI. NeuroImage 9, 407–415. doi: 10.1006/nimg.1998.0424
de Leon, M. J., Li, Y., Okamura, N., Tsui, W. H., Saint-Louis, L. A., Glodzik, L., et al. (2017). Cerebrospinal fluid clearance in Alzheimer Disease measured with dynamic PET. J. Nucl. Med. 58, 1471–1476. doi: 10.2967/jnumed.116.187211
ECG: Cardiology | Research | BIOPAC. (n.d.). BIOPAC systems, Inc. Available at: https://www.biopac.com/application/ecg-cardiology/ (Accessed April 20, 2021).
Eide, P. K., and Ringstad, G. (2019). Delayed clearance of cerebrospinal fluid tracer from entorhinal cortex in idiopathic normal pressure hydrocephalus: a glymphatic magnetic resonance imaging study. J. Cereb. Blood Flow Metab. 39, 1355–1368. doi: 10.1177/0271678X18760974
Foi, A. (2011). Noise estimation and removal in MR imaging: the variance-stabilization approach. IEEE Int. Symposium Biomed. Imaging 2011, 1809–1814. doi: 10.1109/ISBI.2011.5872758
Fultz, N. E., Bonmassar, G., Setsompop, K., Stickgold, R. A., Rosen, B. R., Polimeni, J. R., et al. (2019). Coupled electrophysiological, hemodynamic, and cerebrospinal fluid oscillations in human sleep. Science 366, 628–631. doi: 10.1126/science.aax5440
Gulani, V., Calamante, F., Shellock, F. G., Kanal, E., and Reeder, S. B. (2017). Gadolinium deposition in the brain: summary of evidence and recommendations. Lancet Neurol 16, 564–570. doi: 10.1016/S1474-4422(17)30158-8
Han, F., Chen, J., Belkin-Rosen, A., Gu, Y., Luo, L., Buxton, O. M., et al. (2021). Reduced coupling between cerebrospinal fluid flow and global brain activity is linked to Alzheimer disease–related pathology. PLoS Biol. 19:e3001233. doi: 10.1371/journal.pbio.3001233
Hauglund, N. L., Pavan, C., and Nedergaard, M. (2020). Cleaning the sleeping brain – the potential restorative function of the glymphatic system. Curr. Opin. Physio. 15, 1–6. doi: 10.1016/j.cophys.2019.10.020
Helakari, H., Korhonen, V., Holst, S. C., Piispala, J., Kallio, M., Väyrynen, T., et al. (2022). Human NREM sleep promotes brain-wide vasomotor and respiratory pulsations. J. Neurosci. 42, 2503–2515. doi: 10.1523/JNEUROSCI.0934-21.2022
Hock, C., Golombowski, S., Müller-Spahn, F., Naser, W., Beyreuther, K., Mönning, U., et al. (1998). Cerebrospinal fluid levels of amyloid precursor protein and amyloid β-peptide in Alzheimer’s Disease and major depression – inverse correlation with dementia severity. Eur. Neurol. 39, 111–118. doi: 10.1159/000007917
Ibrahim, T. S., Hue, Y.-K., Boada, F. E., and Gilbert, R. (2008). Tic Tac Toe: highly-coupled, load insensitive Tx/Rx array and a quadrature coil without lumped capacitors. 16th Edn, Toronto, 438.
Ibrahim, T. S., Krishnamurthy, N., Wood, S., Raval, S., and Kim, H. (2013). 20-To-8 channel Tx array with 32-channel adjustable receive-only insert for 7T head imaging. In Proceeding of the 21st International Society of Magnetic Resonance in medicine annual meeting.
Iliff, J. J., Wang, M., Zeppenfeld, D. M., Venkataraman, A., Plog, B. A., Liao, Y., et al. (2013). Cerebral arterial pulsation drives paravascular CSF–interstitial fluid exchange in the murine brain. J. Neurosci. 33, 18190–18199. doi: 10.1523/JNEUROSCI.1592-13.2013
Järvelä, M., Kananen, J., Korhonen, V., Huotari, N., Ansakorpi, H., and Kiviniemi, V. (2022). Increased very low frequency pulsations and decreased cardiorespiratory pulsations suggest altered brain clearance in narcolepsy. Commun. Med. 2:122. doi: 10.1038/s43856-022-00187-4
Jenkinson, M., Beckmann, C. F., Behrens, T. E. J., Woolrich, M. W., and Smith, S. M. (2012). FSL. NeuroImage 62, 782–790. doi: 10.1016/j.neuroimage.2011.09.015
Kananen, J., Helakari, H., Korhonen, V., Huotari, N., Järvelä, M., Raitamaa, L., et al. (2020). Respiratory-related brain pulsations are increased in epilepsy—a two-Centre functional MRI study. Brain Commun. 2:fcaa076. doi: 10.1093/braincomms/fcaa076
Kananen, J., Järvelä, M., Korhonen, V., Tuovinen, T., Huotari, N., Raitamaa, L., et al. (2022). Increased interictal synchronicity of respiratory related brain pulsations in epilepsy. J. Cereb. Blood Flow Metab. 42, 1840–1853. doi: 10.1177/0271678X221099703
Kananen, J., Tuovinen, T., Ansakorpi, H., Rytky, S., Helakari, H., Huotari, N., et al. (2018). Altered physiological brain variation in drug-resistant epilepsy. Brain Behav 8:e01090. doi: 10.1002/brb3.1090
Kim, J.-H., Im, J.-G., and Park, S.-H. (2022). Measurement of CSF pulsation from EPI-based human fMRI. NeuroImage 257:119293. doi: 10.1016/j.neuroimage.2022.119293
Kim, J., Krishnamurthy, N., Santini, T., Zhao, Y., Zhao, T., Bae, K. T., et al. (2016). Experimental and numerical analysis of B1+ field and SAR with a new transmit array design for 7T breast MRI. J. Magn. Reson. 269, 55–64. doi: 10.1016/j.jmr.2016.04.012
Kiviniemi, V., Jauhiainen, J., Tervonen, O., Pääkkö, E., Oikarinen, J., Vainionpää, V., et al. (2000). Slow vasomotor fluctuation in fMRI of anesthetized child brain. Magn. Reson. Med. 44, 373–378. doi: 10.1002/1522-2594(200009)44:3<373::AID-MRM5>3.0.CO;2-P
Kiviniemi, V., Wang, X., Korhonen, V., Keinänen, T., Tuovinen, T., Autio, J., et al. (2016). Ultra-fast magnetic resonance encephalography of physiological brain activity – glymphatic pulsation mechanisms? J. Cereb. Blood Flow Metab. 36, 1033–1045. doi: 10.1177/0271678X15622047
Krishnamurthy, N., Santini, T., Wood, S., Kim, J., Zhao, T., Aizenstein, H. J., et al. (2019). Computational and experimental evaluation of the tic-tac-toe RF coil for 7 tesla MRI. PLoS One 14:e0209663. doi: 10.1371/journal.pone.0209663
Lin, F.-H., Tsai, K. W., Chu, Y.-H., Witzel, T., Nummenmaa, A., Raij, T., et al. (2012). Ultrafast inverse imaging techniques for fMRI. NeuroImage 62, 699–705. doi: 10.1016/j.neuroimage.2012.01.072
Louveau, A., Smirnov, I., Keyes, T. J., Eccles, J. D., Rouhani, S. J., Peske, J. D., et al. (2015). Structural and functional features of central nervous system lymphatic vessels. Nature 523, 337–341. doi: 10.1038/nature14432
Maggioni, M., Katkovnik, V., Egiazarian, K., and Foi, A. (2013). Nonlocal transform-domain filter for volumetric data denoising and reconstruction. IEEE Trans. Image Process. 22, 119–133. doi: 10.1109/TIP.2012.2210725
Martin, B. A., Reymond, P., Novy, J., Balédent, O., and Stergiopulos, N. (2012). A coupled hydrodynamic model of the cardiovascular and cerebrospinal fluid system. Am. J. Phys. Heart Circ. Phys. 302, H1492–H1509. doi: 10.1152/ajpheart.00658.2011
MATLAB - MathWorks. (n.d.). Available at: https://www.mathworks.com/products/matlab.html (Accessed April 20, 2021).
Moeller, S., Yacoub, E., Olman, C. A., Auerbach, E., Strupp, J., Harel, N., et al. (2010). Multiband multislice GE-EPI at 7 tesla, with 16-fold acceleration using partial parallel imaging with application to high spatial and temporal whole-brain fMRI. Magn. Reson. Med. 63, 1144–1153. doi: 10.1002/mrm.22361
Moser, E., Stahlberg, F., Ladd, M. E., and Trattnig, S. (2012). 7-T MR—from research to clinical applications? NMR Biomed. 25, 695–716. doi: 10.1002/nbm.1794
Nedergaard, M. (2013). Garbage truck of the brain. Science 340, 1529–1530. doi: 10.1126/science.1240514
Peng, W., Achariyar, T. M., Li, B., Liao, Y., Mestre, H., Hitomi, E., et al. (2016). Suppression of glymphatic fluid transport in a mouse model of Alzheimer’s disease. Neurobiol. Dis. 93, 215–225. doi: 10.1016/j.nbd.2016.05.015
Pomara, N., Bruno, D., Sarreal, A. S., Hernando, R. T., Nierenberg, J., Petkova, E., et al. (2012). Lower CSF amyloid Beta peptides and higher F2-Isoprostanes in cognitively intact elderly individuals with major depressive disorder. Am. J. Psychiatry 169, 523–530. doi: 10.1176/appi.ajp.2011.11081153
Poncelet, B. P., Wedeen, V. J., Weisskoff, R. M., and Cohen, M. S. (1992). Brain parenchyma motion: measurement with cine echo-planar MR imaging. Radiology 185, 645–651. doi: 10.1148/radiology.185.3.1438740
Posse, S., Ackley, E., Mutihac, R., Rick, J., Shane, M., Murray-Krezan, C., et al. (2012). Enhancement of temporal resolution and BOLD sensitivity in real-time fMRI using multi-slab echo-volumar imaging. NeuroImage 61, 115–130. doi: 10.1016/j.neuroimage.2012.02.059
Posse, S., Ackley, E., Mutihac, R., Zhang, T., Hummatov, R., Akhtari, M., et al. (2013). High-speed real-time resting-state FMRI using multi-slab echo-volumar imaging. Front. Hum. Neurosci. 7:479. doi: 10.3389/fnhum.2013.00479
Purdon, P. L., and Weisskoff, R. M. (1998). Effect of temporal autocorrelation due to physiological noise and stimulus paradigm on voxel-level false-positive rates in fMRI. Hum. Brain Mapp. 6, 239–249. doi: 10.1002/(SICI)1097-0193(1998)6:4<239::AID-HBM4>3.0.CO;2-4
Raitamaa, L., Huotari, N., Korhonen, V., Helakari, H., Koivula, A., Kananen, J., et al. (2021). Spectral analysis of physiological brain pulsations affecting the BOLD signal. Hum. Brain Mapp. 42, 4298–4313. doi: 10.1002/hbm.25547
Rajna, Z., Mattila, H., Huotari, N., Tuovinen, T., Krüger, J., Holst, S. C., et al. (2021). Cardiovascular brain impulses in Alzheimer’s disease. Brain 144, 2214–2226. doi: 10.1093/brain/awab144
Ramanathan, A., Nelson, A., Sagare, A., and Zlokovic, B. (2015). Impaired vascular-mediated clearance of brain amyloid beta in Alzheimer’s disease: the role, regulation and restoration of LRP1. Front. Aging Neurosci. 7:136. doi: 10.3389/fnagi.2015.00136
Rennels, M. L., Blaumanis, O. R., and Grady, P. A. (1990). Rapid solute transport throughout the brain via paravascular fluid pathways. Adv. Neurol. 52, 431–439.
Respiration Transducer for MRI | TSD221-MRI | Research | BIOPAC. (n.d.). BIOPAC systems, Inc. Available at: https://www.biopac.com/product/respiration-transducer-for-mri/ (Accessed April 20, 2021).
Ringstad, G., Valnes, L. M., Dale, A. M., Pripp, A. H., Vatnehol, S.-A. S., Emblem, K. E., et al. (2018). Brain-wide glymphatic enhancement and clearance in humans assessed with MRI. JCI Insight 3:e121537. doi: 10.1172/jci.insight.121537
Ringstad, G., Vatnehol, S. A. S., and Eide, P. K. (2017). Glymphatic MRI in idiopathic normal pressure hydrocephalus. Brain 140, 2691–2705. doi: 10.1093/brain/awx191
Santini, T., Koo, M., Farhat, N., Campos, V. P., Alkhateeb, S., Vieira, M. A. C., et al. (2021a). Analysis of hippocampal subfields in sickle cell disease using ultrahigh field MRI. NeuroImage 30:102655. doi: 10.1016/j.nicl.2021.102655
Santini, T., Wood, S., Krishnamurthy, N., Martins, T., Aizenstein, H. J., and Ibrahim, T. S. (2021b). Improved 7 tesla transmit field homogeneity with reduced electromagnetic power deposition using coupled tic tac toe antennas. Sci. Rep. 11, 1–11. doi: 10.1038/s41598-020-79807-9
Santini, T., Zhao, Y., Wood, S., Krishnamurthy, N., Kim, J., Farhat, N., et al. (2018). In-vivo and numerical analysis of the eigenmodes produced by a multi-level tic-tac-toe head transmit array for 7 tesla MRI. PLoS One 13:e0206127. doi: 10.1371/journal.pone.0206127
Schroth, G., and Klose, U. (1992). Cerebrospinal fluid flow. Neuroradiology 35, 10–15. doi: 10.1007/BF00588271
Shanks, J., Bloch, K. M., Laurell, K., Cesarini, K. G., Fahlström, M., Larsson, E.-M., et al. (2019). Aqueductal CSF stroke volume is increased in patients with idiopathic Normal pressure hydrocephalus and decreases after shunt surgery. Am. J. Neuroradiol. 40, 453–459. doi: 10.3174/ajnr.A5972
Tustison, N. J., Avants, B. B., Cook, P. A., Zheng, Y., Egan, A., Yushkevich, P. A., et al. (2010). N4ITK: improved N3 bias correction. IEEE Trans. Med. Imaging 29, 1310–1320. doi: 10.1109/TMI.2010.2046908
Uğurbil, K., Xu, J., Auerbach, E. J., Moeller, S., Vu, A. T., Duarte-Carvajalino, J. M., et al. (2013). Pushing spatial and temporal resolution for functional and diffusion MRI in the human connectome project. NeuroImage 80, 80–104. doi: 10.1016/j.neuroimage.2013.05.012
Xie, L., Kang, H., Xu, Q., Chen, M. J., Liao, Y., Thiyagarajan, M., et al. (2013). Sleep drives metabolite clearance from the adult brain. Science 342, 373–377. doi: 10.1126/science.1241224
Yamada, S., Ishikawa, M., Ito, H., Yamamoto, K., Yamaguchi, M., Oshima, M., et al. (2020). Cerebrospinal fluid dynamics in idiopathic normal pressure hydrocephalus on four-dimensional flow imaging. Eur. Radiol. 30, 4454–4465. doi: 10.1007/s00330-020-06825-6
Keywords: cerebrospinal fluid (CSF), magnetic resonance imaging (MRI), ultra-high field (7 Tesla), brain clearance, echo-planar imaging (EPI), physiological brain pulsations
Citation: Martins T, de Almeida B, Wu M, Wilckens KA, Minhas D, Ibinson JW, Aizenstein HJ, Santini T and Ibrahim TS (2024) Characterization of pulsations in the brain and cerebrospinal fluid using ultra-high field magnetic resonance imaging. Front. Neurosci. 18:1305939. doi: 10.3389/fnins.2024.1305939
Edited by:
Harald E. Möller, Max Planck Institute for Human Cognitive and Brain Sciences, GermanyReviewed by:
Sung-Hong Park, Korea Advanced Institute of Science and Technology, Republic of KoreaVesa J. Kiviniemi, University of Oulu, Finland
Copyright © 2024 Martins, de Almeida, Wu, Wilckens, Minhas, Ibinson, Aizenstein, Santini and Ibrahim. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Tales Santini, santini.tales@pitt.edu; Tamer S. Ibrahim, tibrahim@pitt.edu
†These authors have contributed equally to this work and share first authorship