Niko Huotari1,2*
Lauri Raitamaa1,2
Heta Helakari1,2
Janne Kananen1,2
Ville Raatikainen1,2Aleksi Rasila1,2
Timo Tuovinen1,2
Jussi Kantola1,2
Viola Borchardt1,2
Vesa J. Kiviniemi1,2
Vesa O. Korhonen1,2- 1Oulu Functional NeuroImaging Group, Research Unit of Medical Imaging, Physics and Technology, University of Oulu, Oulu, Finland
- 2Department of Diagnostic Radiology, Medical Research Center, Oulu University Hospital, Oulu, Finland
Low image sampling rates used in resting state functional magnetic resonance imaging (rs-fMRI) may cause aliasing of the cardiorespiratory pulsations over the very low frequency (VLF) BOLD signal fluctuations which reflects to functional connectivity (FC). In this study, we examine the effect of sampling rate on currently used rs-fMRI FC metrics. Ultra-fast fMRI magnetic resonance encephalography (MREG) data, sampled with TR 0.1 s, was downsampled to different subsampled repetition times (sTR, range 0.3–3 s) for comparisons. Echo planar k-space sampling (TR 2.15 s) and interleaved slice collection schemes were also compared against the 3D single shot trajectory at 2.2 s sTR. The quantified connectivity metrics included stationary spatial, time, and frequency domains, as well as dynamic analyses. Time domain methods included analyses of seed-based functional connectivity, regional homogeneity (ReHo), coefficient of variation, and spatial domain group level probabilistic independent component analysis (ICA). In frequency domain analyses, we examined fractional and amplitude of low frequency fluctuations. Aliasing effects were spatially and spectrally analyzed by comparing VLF (0.01–0.1 Hz), respiratory (0.12–0.35 Hz) and cardiac power (0.9–1.3 Hz) FFT maps at different sTRs. Quasi-periodic pattern (QPP) of VLF events were analyzed for effects on dynamic FC methods. The results in conventional time and spatial domain analyses remained virtually unchanged by the different sampling rates. In frequency domain, the aliasing occurred mainly in higher sTR (1–2 s) where cardiac power aliases over respiratory power. The VLF power maps suffered minimally from increasing sTRs. Interleaved data reconstruction induced lower ReHo compared to 3D sampling (p < 0.001). Gradient recalled echo-planar imaging (EPI BOLD) data produced both better and worse metrics. In QPP analyses, the repeatability of the VLF pulse detection becomes linearly reduced with increasing sTR. In conclusion, the conventional resting state metrics (e.g., FC, ICA) were not markedly affected by different TRs (0.1–3 s). However, cardiorespiratory signals showed strongest aliasing in central brain regions in sTR 1–2 s. Pulsatile QPP and other dynamic analyses benefit linearly from short TR scanning.
Introduction
In 1995 Biswal and co-workers discovered functional connectivity (FC) in resting state blood oxygen level dependent (BOLD) signal of motor cortices by their continuous very low frequency fluctuations (VLF < 0.1 Hz) (Biswal et al., 1995). Initially, the VLF phenomena were linked to existing literature on low frequency physiological phenomena like vasomotor waves (Kiviniemi et al., 2000). The current view is that spontaneous neuronal activity avalanches synchronize brain activity in functionally connected areas and become visible in hemodynamic signals after a delay of few seconds (Liu and Duyn, 2013; Palva et al., 2013; Keilholz, 2014; Ma et al., 2016; Liu et al., 2018). The spread of such avalanches can be depicted by novel techniques such as inverse imaging (INI) and magnetic resonance encephalography (MREG), that sample functional magnetic resonance imaging (fMRI) data with short repetition times (TR) (Lin et al., 2012, 2018; Assländer et al., 2013; Rajna et al., 2015).
In addition to neuronal activity that is coupled to hemodynamics, the BOLD signal reflects multiple other signal sources, such as motion, physiological pulsations and technical artifacts. The physiological factors like cardiorespiratory changes modulate the BOLD signal and can mask the neuronally driven VLF activity in resting state BOLD signals (Wise et al., 2004; Birn et al., 2006). As the BOLD signal reflects blood oxygen level, also direct effects of the cardiorespiratory pulses themselves can be detected in fMRI data (Shmueli et al., 2007; Chang and Glover, 2009, 2010).
Previously, the research of cardiorespiratory brain pulsations has not gained as much interest because they have been deemed as noise. However, there is increasing evidence showing that these physiological signals or the “noise” it produces in TR BOLD data with long TR, can be used to measure disease-specific changes in patient groups (Makedonov et al., 2013; Tuovinen et al., 2017). This is strongly supported by the recent discovery of the glymphatic brain tissue clearance mechanism where the cardiovascular pulsations have been shown to drive the glymphatic brain clearance (Iliff et al., 2012; Nedergaard, 2013; Jessen et al., 2015). The short TR in 3D MREG can critically sample the spread of cardiovascular ~1 Hz and respiratory ~0.3 Hz pulsations and separate them from VLF (<0.1 Hz) quasi-periodic patterns (QPPs) (Kiviniemi et al., 2016).
The extent to which faster physiological pulsations alias over VLF BOLD signal, has been a prevailing uncertainty in BOLD fMRI. Since cardiac frequencies can be faster than 2 Hz, especially in animals and children, the critical sampling rate should be >4 Hz according to the Nyquist theorem, i.e., TR <0.25 s. However, most often the fMRI TRs are >0.4 s and therefore the data cannot critically sample faster cardiac signals (Liu, 2016). Consequently, aliasing between cardiac and the VLFs occur and may alter FC metrics. Additionally, the cardiorespiratory rhythms and their pressure modulations and physiological autonomic nervous system mediated counter-regulations, local vasomotor waves induce heart rate variability, which differ between patients and controls (van der Kooy et al., 2006; Thayer et al., 2010) causing yet another confounding factor in measures of FC that may require faster sampling.
Early literature on the sampling rate on FC measures usually utilizes single slice data that suffers from out of plane motion and other registration problems (Purdon and Weisskoff, 1998; Peltier et al., 2003; Kiviniemi et al., 2005). Recent studies on the fMRI sampling rate effects on resting-state FC has shown surprisingly small effects (Golestani et al., 2017; Demetriou et al., 2018). However, that and some other recent studies have usually been limited to < 4 Hz sampling rates for whole brain coverage (Cordes et al., 2014; Liu, 2016; Golestani et al., 2017; Chen et al., 2019) with different type of signal simulations extending below the critical 4 Hz limit.
However, more and more critically sampled 3D whole brain fMRI data has started to emerge, such as 0.136 sec TR VEPI (Posse et al., 2013), 100 ms MREG (Assländer et al., 2013; Lee et al., 2013; Kiviniemi et al., 2016; Raitamaa et al., 2018), 50 ms GIN (Boyacioglu and Barth, 2013) and currently leading 25 ms 3D whole brain scan INI (Chang et al., 2013b). This critically sampled data shows robust novel phenomena of the human brain physiology, such as propagating cardiorespiratory pulsations that both interact and modulate each other, depending on their anatomical proximity to pulsation sources (Posse et al., 2013; Kiviniemi et al., 2016; Raitamaa et al., 2018). Based on our observations, these novel signal changes cannot be comprehensively simulated due to their complex spatiotemporal pattern that is dynamically changing.
Therefore, in this study we used real critically sampled 0.1 s TR 3D single shot MREG data to explore how different sampling rates affect the results of the most commonly used resting state fMRI analysis tools. The 0.1 s TR MREG signal was downsampled to higher TRs ranging from 0.3 to 3 s. The use of subsampled TR (sTR) removes the confounding factors of imaging different TR values in separate scans and/or individuals, enabling identical physiological status and technological noise structure for comparing different sTRs. The single shot 3D data sampling scheme was further compared with interleaved slice sampling (INT) variant of the MREG data, and, with conventional interleaved gradient recalled echo-planar BOLD imaging (EPI BOLD).
The hypothesis was that the faster sTR produces both spatially more accurate brain maps and more accurate time series without aliasing of physiological pulsations. The quantified connectivity metrics included stationary spatial, time, and frequency domains, as well as dynamic analyses. In addition, the effects of aliasing were evaluated.
Materials and Methods
Participants
Ten healthy subjects (8 males, 23.8 ± 2.1 years old) were placed in the MRI scanner and asked to lay still and keep their eyes open and fixated on a cross on the screen while thinking of nothing particular (eyes open, resting state). Ear plugs were used to reduce scanner noise. Cushions were placed beside ears to restrict movement and to further reduce scanner noise. MREG (5 min) and EPI BOLD (5 min) sequences were scanned in said order. Written informed consent was obtained from each subject prior to scanning, in accordance with the Helsinki declaration. The study protocol was approved by the regional Ethical committee of Northern Ostrobothnia Hospital District in Oulu University Hospital.
Data Acquisition and Preprocessing
Subjects were scanned using Siemens 3T SKYRA scanner with 32-channel head coil. Additional cardiorespiratory data were collected using MRI-compatible multimodal neuroimaging system (Korhonen et al., 2014). MREG is a 3D single shot stack of spirals (SOS) sequence that under-samples k-space to reach a sampling rate of 10 Hz allowing critical imaging of physiological pulsations (Assländer et al., 2013). The SOS gathers k-space in 60 ms with spiral in/out repeating in every other turn continuously in positive z-direction to minimize the air-sinus off-resonance artifact (for more details, c.f. Assländer et al., 2013). The point spread function of the SOS-sequence is 3 mm with minimized off-resonance effects compared to other k-space undersampling strategies like concentric shells and spokes (Zahneisen et al., 2012; Assländer et al., 2013). Scanning parameters (TR = 100 ms, TE = 36 ms, flip angle = 5°, 3D matrix = 643, FOV = 192 mm) enabled scanning of the whole brain in 10 Hz with voxel size of 3 × 3 × 3 mm3. Conventional EPI BOLD scans were collected from the same subjects (TR = 2,150 ms, TE = 28 ms, flip angle = 15°, voxel size = 3 × 3 × 3 mm, matrix size 64*64, 45 slices = 47 ms in plane readout). In both methods, we used relatively low flip angles to minimize specific absorption rate (SAR), spin history effects, physiological pulsations and radio frequency (RF) artifacts in EEG in comparison to default flip angles (Gonzalez-Castillo et al., 2011; Assländer et al., 2013).
A reference image for MREG was acquired with a multi slice double gradient echo sequence with TR = 593 ms, TE = 2.46/4.92 ms, flip angle = 50°, dwell time = 4.9 us, FOV = 192 mm. The reference and raw data from the MREG sequence were transferred offline to a computing grid and reconstructed using the MATLAB tool provided with the sequence. The tool allows for a choice between several parameters for regularized reconstruction (Hugger et al., 2011); we selected L2-norm with finite difference operator (called “Total Variation” in the tool) and the regularization parameter was reduced to lambda = 0.15 from default 0.2 in order to obtain higher signal-to-noise ratio (SNR) images. Conjugate gradient optimization was also performed for 35 iterations for more robust convergence of images, c.f. Figure 1. Coil sensitivities were estimated from the reference image with the adaptive method, and dynamic off-resonance correction in k-space was used to minimize respiration and other motion related off-resonance artifacts from the data (Pfeuffer et al., 2002; Zahneisen et al., 2014). Anatomical 3D MPRAGE (TR = 1,900 ms, TE = 2.49 ms, TI = 900 ms, flip angle = 9°, FOV = 240 mm, 0.9 mm cubic voxel) images were used to register both MREG and EPI BOLD data into Montreal Neurological Institute (MNI) space.
Figure 1. Raw magnetic resonance encephalography (MREG) data with repetition time (TR = 0.1 s) and the subsampling scheme using an example arterial signal ([−2 30 16] mm in Montreal Neurological Institutional (MNI) space). (A) Raw spatial MREG image. (B) Raw timeseries of the arterial signal. (C) Processed and subsampled arterial signal and corresponding frequency amplitude spectra. Please note the aliased cardiac peak in the sTR 1–3 s in the amplitude spectra at ~0.2 Hz.
Both, EPI BOLD and MREG data were preprocessed with a typical FSL pipeline (Jenkinson et al., 2012). The data were high-pass filtered with cut-off frequency of 0.008 Hz (125 s). T1-relaxation effects were minimized using dummy scans (8 s). Motion correction was performed using FSL MCFLIRT. FSL BET was used for brain extraction (fractional intensity = 0.25, threshold gradient = 0.22, neck and bias-field correction). Images were spatially smoothed with 5 mm FWHM Gaussian kernel using fslmaths. Three-dimensional 3DMPRAGE anatomical images were used to register both the EPI BOLD and the MREG data into 4 mm MNI space prior to group ICA using FLIRT (Jenkinson and Smith, 2001; Jenkinson et al., 2002). FILM pre-whitening and data smoothness estimation was produced automatically by FSL MELODIC (Woolrich et al., 2001; Beckmann et al., 2006). FILM implements a robust method of correcting for auto-correlation in fMRI time-series which is (theoretically) independent of TR duration (Woolrich et al., 2001; Smith et al., 2004; Demetriou et al., 2018).
The MREG (TR = 0.1 s) data was downsampled to different sTR settings from TR 0.1 s to 0.3, 0.5, 1.0, 1.5, 1.8, 2.2, and 3.0 s, respectively by taking every nth [n = 3, 5, 10, 15, 18, 22, 30] sample from every voxel timeseries. In addition, an interleaved variant of the MREG data were computed by taking every 1,3,5…21; 2,4,6…22th axial slice (hence the INT sTR = 2.2 s) to emulate the interleaved EPI data gathering for comparing single shot MREG trajectory downsampled at 2.2 s sTR data and conventional gradient recalled EPI (TR = 2.15 s). MATLAB was used for MREG data downsampling and interleaving. Total of 100 datasets were obtained after processing (10 sTR settings for 10 subjects).
As one of the focuses was to assess physiological signal aliasing effects, we retained the physiological pulsations in the data as much as possible. Therefore, cerebrospinal fluid (CSF) and white matter were not regressed out from the datasets like they often are in functional connectivity analyses. Global signal was analyzed but not regressed, as the benefit of its regression is still under debate (Murphy and Fox, 2017). Furthermore, the datasets were not de-spiked, since there is no clear consensus yet what kind of de-spiking is advisable to apply to ultra-short TR fMRI data, especially since most aggressive de-spiking (AFNI 3dDespike -NEW25) removes some of the physiological pulsations from the data (Raitamaa et al., 2018).
Time and Spatial Domain Analysis (ICA, CV/tSNR, DPARSF)
Group PICA was computed for all 10 sTR settings. For every group ICA run, 70 independent components were calculated using FSL MELODIC in default setting (Kiviniemi et al., 2009). FSL function fslcc was used to calculate correlation values between different PICA components calculated for every sTR and 42 resting state network templates defined earlier (Kiviniemi et al., 2009; Abou-Elseoud et al., 2010). Default mode network (DMN) posterior cingulate cortex (PCC), cerebral artery, visual, auditory, motor and ventral attention network components were selected as interesting reference components that were visualized to show the effect of changing sTR (Beckmann et al., 2005; Kiviniemi et al., 2009; Smith et al., 2009).
CV is a standardized measure used in e.g., engineering and physics, which describes the variability of a dataset compared to its mean. CV was used as a metric for the variation of physiological fluctuations in the signal. Recently, the CV of BOLD signal (CVBOLD) has been shown to be altered by the disease processes (Makedonov et al., 2016; Tuovinen et al., 2017; Kananen et al., 2018). CV was calculated for each subject and the calculation was carried out for every voxel timeseries:
where μ is the standard deviation, and σ is the mean of the voxel timeseries. Mean images of the resulting CV maps were computed. For statistical analyses, PCC region of interest (ROI) values were computed for every sTR individually and compared to mean reference (TR = 0.1 s) map. In addition, temporal signal-to-noise ratio (tSNR) values were computed for each sTR using white matter (WM) and gray matter (GM) ROIs. The ROI areas were obtained via FSL atlas tools using a probabilistic threshold value of 50. A mean value from both WM and GM were calculated from each subject.
DPARSF (FC, ReHo, ALFF/fALFF)
For each FSL pre-processed dataset, seed-based functional connectivity (FC), regional homogeneity (ReHo), and amplitude and fractional amplitude of low frequency fluctuations (ALFF/fALFF) were calculated using Data Processing Assistant for Resting-State fMRI (DPARSF) V4.3_171210 MATLAB software package (Chao-Gan and Yu-Feng, 2010). The processing steps in DPARSF are described for each analysis method separately below. Global signal regression or white matter/CSF signal regressions were not used.
FC calculates seed-based correlation values between mean timeseries of a selected ROI area and the rest of the brain voxel timeseries. Six seed areas (5 mm spherical ROI) were selected for the analyses: DMN PCC (0 −53 26 mm in MNI), cerebral artery (0 32 16 mm in MNI), visual (2 −86 16 mm in MNI), auditory (50 −2 −8 mm in MNI), motor (2 −18 60 mm in MNI) and ventral attention network (-46 18 32 mm in MNI). Mean FC maps for each sTR and ROI were computed to see whether the changing of sTR influence the connectivity. For statistical analyses, correlation coefficient values were calculated between the mean ROIs (TR = 0.1 s) and other sTR settings individually using FSL function fslcc. Additionally, sTR 2.2 s (INT) and sTR 2.2 s were compared to test the effect of interleaving in MREG data.
ReHo measures the degree of regional synchronization of neighboring areas by calculating Kendall's coefficient of concordance (KCC) from the timeseries of every voxel and compares the neighboring voxels (Kendall and Gibbons, 1990; Zang et al., 2004). Cluster size of 27 voxels was used. Spatial smoothing (FWHM: kernel size [4 4 4]) was applied after ReHo calculations. After individual ReHo computations, mean spatial maps were calculated. Furthermore, ReHo KCC values were computed from the mean PCC ROI for all sTR settings.
Frequency Domain Analysis
Fast Fourier transformation (FFT) amplitude spectra from a single subject arterial region ([−2 30 16] mm in MNI) (Figure 1) and global image signal, venous (high respiratory power) and arterial (high cardiac power) amplitude spectra were analyzed to observe the impact of sTR in different physiological frequency bands. Additionally, normal distributions curves (bin size 5) of histogram of demeaned global signal in every sTR were computed.
To examine, whether the change in sTR affects the VLF content, ALFF and fALFF were computed. ALFF calculates the sum of spontaneous low frequency activations of a selected frequency band (Zang et al., 2007). FALFF was used to estimate the ratio of amplitude spectrum of VLF in comparison to the whole FFT spectrum (Zou et al., 2008). Frequency band of 0.01–0.1 Hz was selected for the computations.
Furthermore, we studied how the mean FFT power of VLF (0.01–0.1 Hz), respiratory (0.12–0.35 Hz) and cardiac (0.9–1.3 Hz) pulsations occur in different sTR settings. The respiratory and cardiac frequency bands were chosen based on group-level minimum and maximum frequencies observed from the physiological recordings and datasets so that the band is as limited narrow as possible while still encompassing the cardiorespiratory frequencies of each subject. The same frequency bands were used for each subject. Each band was separated by a margin to minimize overlap.
Spectral maps of the frequency bands were calculated for each dataset using AFNI 3dPeriodogram. The function outputs a power spectrum for each voxel timeseries. The frequency bins corresponding to VLF, respiratory and cardiac bands were collected from the periodogram datasets based on individual physiological monitoring. Furthermore, the collected bins were summed and finally mean of the sums was computed for each sTR. The cardiac bins were not calculated for sTRs > 0.5 s, as the slower sampling rates cannot extract the cardiac frequencies.
The effect of signal aliasing was estimated in two ways: (i) spatially, by correlating the mean cardiac and respiratory FFT power maps to individual respiratory and VLF maps for each sTR setting, respectively, (ii) power spectrally, by calculating the FFT power change in the respiratory and VLF bands in the right medial artery ROI ([46 −2 −8] mm in MNI) with different sTRs. Mean 0.1 s TR MREG datasets were used as baseline maps and compared how much the physiological power starts to overlap per given sTR.
QPP Analysis
All datasets were bandpass filtered to VLF (0.01–0.1 Hz) band using AFNI 3dTproject. We used a modified pattern finding algorithm to obtain quasi-periodic patterns (QPPs) and evaluated their changes in signal intensities and pulse propagations (Kiviniemi et al., 2016; Raitamaa et al., 2018). Estimation of timing and length of VLF pulse (length: 105–146 time points) for every subject was obtained from the VLF filtered global signal from TR 0.1 s. Length and timing were adjusted to the subsampled datasets. Subject-specific 4D QPP maps were created for every sTR. For the analyses, MATLAB circshift was used to ensure that QPP maps were in same phase.
QPP strength (i.e., how closely it resembled the template) of VLF wave for each sTR was quantified by correlation between the average QPP pulse and the VLF pulse at the last iteration of the QPP algorithm (Thompson et al., 2014). The mean correlation coefficients of VLF pulses from the last iteration were extracted from each subject and were compared as a function of sTR.
To measure the repeatability of the detected QPP maps, the spatial correlation of subject-specific average QPP maps of sTR > 0.1 s were compared to TR 0.1 s average QPP map using MATLAB corrcoef. For this, each sTR QPP map was interpolated to match the reference TR 0.1 s map.
For group mean images and videos, TR 0.1 s QPP maps were set to length of 150 time points and other sTRs lengths were adjusted accordingly. For every sTR, average group QPP maps were created.
Statistical Tests
Paired sample t-tests (MATLAB ttest) were used for all statistical testing to compare the reference MREG data (TR = 0.1 s) with sTR results (Figures 3, 6, 7, 9A). Furthermore, sTR 2.2 s (INT) and sTR 2.2 s were compared to test the effect of interleaving in MREG data. In QPP repeatability analysis, linear fitting was used as a measure of statistical significance. Statistical significances are indicated as *p < 0.05, **p < 0.01, and ***p < 0.001.
Results
As an example of the image quality of MREG 0.1 s TR data, Figure 1A illustrates a raw spatial image of one subject in three planes. An example of cardiovascular pulse timeseries (ROI MNI: [-2 30 16] in mm) from the anterior cerebral artery are presented in Figure 1B. Respective subsampling scheme is presented in Figure 1C, in which six examples of the sTR time series are shown to highlight the effect of changing sTR.
Commonly Used Resting State Metrics (FC, ReHo, CV/tSNR, ICA)
Seed voxel correlation of FC analysis, ReHo and CV produced nearly identical results as a function of sampling rate. Visually, only small differences could be detected in the spatial functional connectivity measures in group level maps as a function of sampling rate in Figure 2.
Figure 2. Different group-level resting state metrics with changing subsampled repetition times (sTR) including functional connectivity (FC), regional homogeneity (ReHo), coefficient of variation (CV) and independent component analysis (ICA). FC and ICA analyses comprise default mode network (DMN), artery (Art), visual (Vis), auditory (Aud), motor (Mot) and ventral attention network (VAN). For comparison, echo-planar imaging (EPI) maps below tends to show spatially more widespread FC maps and relatively similar ICA maps depending in the RSN.
Functional connectivity maps (Figure 2) showed nearly identical spatial results, which were also highly comparable to conventional 2.15 s TR EPI BOLD data or with interleaved MREG 2.2 s sTR data. The quantified connectivity measures presented no significant differences excluding INT (sTR = 2.2 s) which had significantly (p < 0.01) lower correlation coefficients compared to TR 0.1 s and compared to sTR 2.2 s (Figure 3A and Supplementary Figure 1). In addition, EPI scans had marginally elevated correlation coefficients compared to the MREG (TR = 0.1 s).
Figure 3. Resting state metrics in different subsampled repetition times (sTR) compared to reference maps (TR = 0.1 s). (A) Individual functional connectivity (FC) values from default mode network (DMN), posterior cingulate cortex (PCC ROI) correlated to reference FC map. (B) Regional homogeneity (ReHo) Kendall's coefficient of concordance (KCC) values from PCC ROI compared to mean reference map (TR = 0.1 s). (C) Coefficient of variation (CV) values from PCC ROI compared to mean reference map. (D) Correlation values between group probabilistic independent component analysis (PICA) components in different sTR settings and previously acquired 42 resting state network ICA templates. Significant differences between TR 0.1 and other sTR settings are marked with *p < 0.05, **p < 0.01, and ***p < 0.001. Outliers are marked with ♦.
Similarly, ReHo were visually identical at different sampling rates (Figure 2). However, with interleaved sampling of 1,3,5…,2,4,6 of axial MREG slices were collected with a final 2.2 s sampling rate, the ReHo values decreased significantly (p < 0.001), when compared to both single shot k-space trajectory with TR 0.1 s and sTR 2.2 s (Figure 3B). Furthermore, conventional EPI showed significantly (p < 0.001) higher and sTR 2.2 lower (p < 0.05) values.
CV of the image data were nearly identical with respect to the altering sTR (Figure 2). The EPI BOLD data had significantly (p < 0.05) lower CVs compared to MREG (TR = 0.1 s), while downsampling or interleaved data gathering had no effect on the CV values. The same effect can be seen from the mean tSNR values from WM and GM ROIs where downsampling did not show any significant change but EPI BOLD values were significantly (p < 0.001) higher compared to MREG (Supplementary Figure 2).
ICA detected RSNs in all sTR settings, but there was some variability between the results (ICA in Figure 2). The conventional EPI BOLD and interleaved data showed similar results as well. In spatial consistency analysis of the 42 RSN template correlations, all sTR settings produced very similar results, except for EPI, which had significantly higher correlation values (p < 0.01).
Mean Spectral Metrics
The group mean FFT amplitude spectra of global image signal, respiratory and cardiac signals are presented in Figure 4. In global signal spectra, the amplitude peaks were relatively higher in VLF band with increasing sTR. This might occur due to cardiorespiratory aliasing. The respiratory amplitude peaks from sagittal sinus and cardiac amplitude peaks from right middle cerebral artery presents multiple peaks in the spectra due to normal variations in respiratory and cardiac rates. Please notice, these peaks were lost as the sTR was increased >0.5 s. The same effect was observed in the arterial signal on a single subject level (Figure 1). EPI BOLD presented the lowest VLF amplitudes in all images.
Figure 4. Group-level fast Fourier transform (FFT) amplitude and histogram distribution analyses in different subsampled repetition times (sTRs). (A) FFT amplitude spectra of global image signal where the 0.1 TR data is cut to 2 Hz from original 5 Hz. (B) FFT amplitude spectra of respiratory region of interest (sinus sagittal, ROI, [2 −98 4] mm in Montreal Neurological Institute (MNI) space). (C) FFT amplitude spectra of cardiac ROI (right medial artery) [46 −2 −8] mm in MNI. Due to differences in individual cardiorespiratory rates, several frequency peaks can be detected at the group level. (D) Mean normal distribution curves from demeaned global signal histograms show lowering and widening of the distributions as a function of sTR.
The global signal distribution curves illustrated the differences in statistical power between sTR values and revealed an exponential decay in histogram counts and widening of the distributions per each sTR. Interestingly the 2.15 EPI BOLD distributions are highly similar in shape with MREG 0.1 s data but >5 times smaller.
BOLD Signal Frequency Amplitude and Power Mapping
ALFF results indicated that the increasing sTR raised the VLF power of the images (Figure 5). Conventional EPI data showed lower ALFF values than MREG data, but the interleaving had no effect. FALFF results were most clearly affected by the increasing sTR (Figure 5). FALFF showed an increase over the sTR values which is due to the proportional increase of the lowest frequencies due to the reduction of the spectral coverage in higher sTR values.
Figure 5. Mean fast Fourier transform (FFT) amplitude (ALFF/fALFF) and power (VLF, Resp, Card) encoding maps for each subsampled repetition times (sTRs). The cardiac power can be detected until 0.3 s and respiratory until 2.2 s sTR. However, the respiratory and cardiac power images started to overlap in sTR > 0.5 s, increasing to sTR 1–2 s. Very low frequency (VLF) power showed a subtle but steady increase with increasing sTR.
Group mean FFT amplitude and power encoding maps (Figure 5) revealed how the cardiac, respiratory and VLF frequency intensities changed in different sTR settings. The cardiac power started to fade away in sTR >0.3 s and respiratory power above 2.2 s sTR. Furthermore, the cardiac power distributions started to overlap on top of respiratory frequency maps as the sTR is increased >0.5 s indicating aliasing, c.f. Figures 5–7. The conventional 2.15 s TR EPI measurement had notably lower power (please notice different scaling for EPI BOLD in Figure 5). Interleaved 2.2 s sTR power was also lower compared to 2.2 s sTR single shot trajectory in VLF power images, which agrees with ReHo and mean global signal amplitude spectral changes.
The aliasing was quantified both spatially and in frequency power analyses. This was evaluated by measuring how much different frequency maps start to resemble each other spatially as a function of sTR (Figure 6). In critically sampled 0.1 s TR MREG data the cardiac, respiratory and VLF frequencies had a low ~0.4 mean spatial correlation. For example, the major arteries seen in cardiac frequency maps were lacking from the respiratory maps, while the respiratory band dominated in posterior CSF spaces (Figure 5). However, the spatial correlation between cardiac and respiratory power maps became increasingly more similar (p < 0.001) and started to increase as a function of the sTR until 1.0 s, after which the similarity plateaus at 0.65 level until sTR 3.0 s (Figure 6A). Furthermore, the spatial correlation between cardiac (p < 0.001) and respiratory (p < 0.01) vs. VLF increased steadily until sTR 1.8 s.
Figure 6. Spatial aliasing measure using spatial correlation values. (A) The mean MREG cardiac fast Fourier transform (FFT) power image of repetition time (TR = 0.1 s) correlated to individual respiratory power images as a function of sampled TR (sTR) and echo-planar imaging (EPI) scan. (B) Mean cardiac FFT image correlated to individual very low frequency (VLF) images. (C) Mean respiratory image correlated to individual VLF images. Significant differences between TR 0.1 s and other sTR settings are marked with **p < 0.01 and ***p < 0.001. Outliers are marked with ♦.
In frequency power analysis of a right medial artery ROI ([46 −2 −8] mm in MNI) near insula, the image signal power from the arterial ROI started to increase significantly (p < 0.001) both in the respiratory and VLF power band in sTR >0.3 s as a sign of aliasing (Figure 7A). The results indicate that once the cardiac pulsation is no longer critically sampled, the power starts to become aliased as a “respiratory” and “VLF” power as a function of sTR increase. When correlating the mean respiratory maps to individual VLF maps, the aliasing effect could also be seen in VLF range as the correlation values increased significantly (p < 0.01) after 0.1 s TR and almost linearly as the sTR increased. However, the influence of the sTR had a clearly less steep effect by the increasing sTR in the VLF. The conventional 2.15 s BOLD images showed low spatial correlation since the cardiac power is very low to begin with and so the spatial correlation is also low.
Figure 7. Cardiac aliasing measure using fast Fourier transform (FFT) power intensity values of arterial region of interest (ROI, [46 −2 −8] mm in Montreal Neurological Institute (MNI) space). (A) Intensity values of individual respiratory power images. (B) Intensity values of individual very low frequency (VLF) power images. Significant differences between repetition time (TR = 0.1 s) and other sampled TR (sTR) settings are marked with *p < 0.05, **p < 0.01, and ***p < 0.001. Outliers are marked with ♦.
QPP Analysis of Spreading BOLD Waves
The power of cardiac pulsation could not be critically detected in sampling rates >0.3 s (Figures 3–5) and therefore we did not quantify differences in the detection of cardiac pulse propagation like we did earlier (Raitamaa et al., 2018). Also, the respiratory power suffered from marked aliasing in sTRs above 0.5 s (Figures 5–7). Therefore, we quantified the performance of different sTR in detecting VLF (0.01–0.1 Hz) propagating QPP BOLD waves.
Figure 8 illustrates 15 s QPP waves captured with different sTR and how the spatiotemporal illustration of each wave changed as a function of sTR. The lowest sTR 0.1–0.5 s were down-sampled from 150, 50 and 30 to 15 images for illustration of the QPP wave spreading in axial plane. Videos of QPPs in TR of 0.1 s and 0.5, 1.0, 2.2 and 3.0 s are shown in Supplementary Material. As expected, much of the dynamics were lost as a function of increasing sTR.
Figure 8. Group mean very low frequency (VLF) waves as a function of subsampled repetition time (sTR). Please notice that the first three rows are downsampled from 150, 50, and 30 to 15 images, respectively to fit the picture.
We also quantified the QPP strength within each subject of the QPP VLF waves as a function of sTR (Figure 9A). The analysis revealed that 0.1 s TR data had highest spatial correlation of the detected VLF waves on average despite the largest number of brain volumes (i.e., 150 vs. 5 volumes between highest and lowest TR). There is a linear trend where the intra-individual detection accuracy of the detected VLF QPP waves gave lower (0.01 < p < 0.05) values as a function in sTR > 0.3 s excluding EPI.
Figure 9. Very low frequency (VLF) quasi-periodic pattern (QPP) stability measurements. (A) Intra-individual QPP strength as a function of subsampled repetition time (sTR). Significant differences between repetition time TR 0.1 s and other sTR settings are marked with *p < 0.05 and **p < 0.01. (B) Inter-individual correlation coefficients of QPP wave repeatability as a function of sTR. Correlation coefficients decreased statistically significantly (p < 0.001) by decreasing slope of sTR. Outliers are marked with ♦.
When comparing how similar QPP waves were detected between subjects, the 0.1 s MREG data was used as a reference due to highest VLF wave detection accuracy. Each sTR QPP map was interpolated to correspond the 0.1 s MREG time points. On average, the detected QPP waves became linearly less correlated as a function of sTR (Figure 9B). Conventional EPI had less robustly detected QPP waves. Interleaved timing of data acquisition showed no effect on the wave detection accuracy nor on repeatability of the detected QPPs.
Discussion
The study analyzed the effect of 3D image sampling rate on most popular fMRI metrics used. The time and spatial resting state fMRI metrics (FC, CV/tSNR, spatial ICA) were not markedly affected by sTR. In frequency domain analysis, the aliasing of the cardiorespiratory power seemed to increase signal power as a function of the sTR (sTR > 0.3 s). Importantly, aliasing effects occur mostly between cardiac and respiratory power. The VLF power increased also significantly as a function of sTR due to aliasing but the power of aliasing was smaller. In dynamic QPP analyses, shorter sTRs seemed to produce more stable results.
Also, to our surprise, the effect of sampling rate on non-dynamic rs-fMRI measures was not as strong as hypothesized. Most of the measures stayed stable within the range of sTR. Golestani and co-workers found surprisingly minimal effect on sampling rate on ReHo, ALFF and FC results that are in full agreement with our results (Golestani et al., 2017). Several researchers have found that higher sampling rate appears to be beneficial for resting state fMRI measures, especially for ICA/dual regression and our results also agree on this (Smith et al., 2013; Preibisch et al., 2015; Golestani et al., 2017; Demetriou et al., 2018.
Compared to MREG in general, the interleaved EPI BOLD had statistically significantly (p < 0.001) increased ReHo and tSNR. However, the f/ALFF, global signal amplitude spectrum baseline and signal variance as measured with CV are reduced in EPI compared to MREG. Interleaved MREG variant (INT sTR 2.2 s) also induced significantly (p < 0.001) lower ReHo compared to both TR 0.1 and to single shot variant (sTR 2.2 s). The conventional interleaving introduces TR/2 delay and interleaved multislab scanning introduces TR/2xN delay between subsequent scans, where N is the number of slabs. This introduces delays between neighboring voxels in longitudinal z-direction and alters regional similarity and frequency measures due to discontinuous sampling of the propagating cardiorespiratory pulses that traverse the brain repeatedly in addition to VLF pulses (Kiviniemi et al., 2016).
Aliasing in Cortical Connectivity Analysis
The current data illustrates that while the aliasing of cardiorespiratory pulsations is a significant factor in frequency domain, the currently used time and spatial domain analyses tools perform well in detecting robust resting state connectivity. In summary, the stationary spatial and temporal connectivity measures, whether local or long distance, were not significantly affected by changes in sTR. However, interleaved slice acquisition seems to affect some of these measures.
This study confirms the information from reduced data length analytics (Bright and Murphy, 2015) and is consistent with coactivation pattern (CAP) analysis (Liu and Duyn, 2013), where the information on spatial functional connectivity of regions can be depicted even in one single brain volume. However, the CAP analysis for instance, requires individual voxel level thresholding with the time domain signal mean/std.
ICA performs best with large data distributions and therefore conventionally spatial ICA is preferred over temporal ICA, since it offers larger distributions (Calhoun et al., 2001; Kiviniemi et al., 2003; Beckmann et al., 2005). Spatial ICA also uses BOLD signal's temporal variance which induces the non-Gaussian changes in the signal distributions, by which the statistical independence is inferred (Calhoun et al., 2001; Kiviniemi et al., 2003; Beckmann et al., 2005). The clear advantage of the short TR can be seen in Figure 4, where 0.1 s TR overall signal distribution histogram has almost 5 times higher and substantially larger distribution compared to other distributions. This is the basis of the statistical power advantage already depicted by several groups for using short TR measurements (Smith et al., 2013; Preibisch et al., 2015; Golestani et al., 2017; Demetriou et al., 2018). Furthermore, in the case of combined spatiotemporal ICA or temporal ICA alone approaches, the 0.1 s data most likely outperforms the slower sTRs simply due to statistically more valid signal histograms (Figure 4D).
The physiological pulses themselves modulate the BOLD signal (Birn et al., 2006, 2008; Chang and Glover, 2010; Chang et al., 2013a). Recently, ultra-fast MREG data was able to illustrate how the cardiac and respiratory pulsations propagate in repeated waves over the whole brain, giving rise to physiologically driven and modulated variance in the data (Kiviniemi et al., 2016; Raitamaa et al., 2018). However, most of the cardiorespiratory power is centered within central areas near CSF ventricles (respiratory) and major arteries and veins (cardiac). The VLF power dominance in the brain cortex seems to enable robust connectivity measurement despite imminent cardiorespiratory aliasing especially in sTRs 1–2 s.
The physiological signal sources need to be separated from the functional signals before any realistic interpretations of the neurovascular task or connectivity data can be drawn. Through the years, source separation tools such as the ICA have been developed to offer robust detection of functionally connected regions or constellations of regions, such as resting state networks (Calhoun et al., 2001; Kiviniemi et al., 2003; Beckmann et al., 2005; Griffanti et al., 2015; Vidaurre et al., 2017). However, in this respect also other measures, such as FC and CV are highly reproducible over a range of sTRs. Therefore, even though there is significant aliasing, especially in areas near medial cerebral artery and sagittal sinus neighborhood, the dominance in VLF induces robust functional connectivity in the cortex.
BOLD signal stability is often measured as tSNR or CV which is the inverse of tSNR. In this study, we calculated both. In terms of signal observations of hemodynamically convolved neuronal activation response, in our sampling scheme the CV or tSNR do not change as we only downsample the MREG data instead of scanning using different TRs. However, compared to EPI BOLD the CV values were significantly higher (Figures 2, 3) and tSNR values significantly lower (Supplementary Figure 2) which could be partially caused by higher flip angle and slightly shorter TE in EPI BOLD. In spite of all, both sequences have relatively long TEs, which makes them T2* -weighted.
Recent studies indicate how the CV reflects also physiological pulsation changes, rather than neuronally driven alterations. Furthermore, the CV has recently been shown to have high sensitivity to pathological condition even at individual level (Makedonov et al., 2016; Tuovinen et al., 2017, 2018). Our preliminary experience on disease related BOLD signal noise metrics suggests that at 1.8 s TR the BOLD signal CV can also be highly sensitive to pathological conditions (Tuovinen et al., 2018). Data from 0.1 s MREG CV measures seem to be at least as sensitive in intractable epileptic patient data (Kananen et al., 2018). This is in line with the CV results that do not alter as a function of TR but can be sensitive to differences in sequence parameters and k-space trajectory and other technical issues, c.f. Figure 2.
Aliasing of Physiological Pulsations in Central Areas
The FFT frequency power maps illustrate that the cardiac 1 Hz power is mainly detectable in sTR < 0.3 s. In the most critically sampled data, the 1 Hz cardiovascular signal pulses most prominently in the paravascular space in areas near major cerebral arteries and somewhat in the major venous sinuses. From the periarterial areas the cardiovascular pulses are convected into the CSF ventricles, centerline parasagittal CSF spaces, c.f. Figure 5. The respiratory power is nearly absent from the periarterial spaces and dominates more on the posterior and cortical perivenous structures and is relatively stronger in the posterior central CSF spaces. The VLF power tends to align mostly along the cortical gray matter without strong overlap with the midline parasagittal area.
Above 0.3 s sTR, the cardiac power became mostly aliased over respiratory power in central brain bordering CSF spaces. Thus, all critical analysis of respiratory power changes needs to be performed with data < 0.5 s TR. The most prominent aliasing occurred in sTRs between 1 and 2 s, where the respiration is still sampled critically but the cardiac is aliased over it; they become strongly mixed into sTR 1–2 s signal. This type of cardiorespiratory aliasing however seems to affect the measuring of FC on cortical structures minimally (Figures 2–3), since the connectivity occurs dominantly in VLF frequencies (< 0.1 Hz) that dominate in the cortex to begin with. Also, the effects of aliasing as a function of sTR on the VLF FFT power are less severe around sTR 1–2 s but increase further in sTRs >2.0 s (Figures 2, 5, 6).
The power of FFT amplitude spectra increased in VLF range as a function of sTR, which is a sign of cardiorespiratory aliasing (Kiviniemi et al., 2005; De Luca et al., 2006). The effect was highest in sTR 3.0 s. This reduces sensitivity to changes in physiological pulsations and furthermore does no longer have the capability to differentiate neither cardiorespiratory nor even different VLF fluctuation peaks as different sources. Furthermore, physiological pulsations seen in 0.1 s TR could not be detected on a group level global or individual voxel signal above 0.3 s sampling rate. Also, VLF power peak features became undetectable with vanishing power of the increased sTR.
Importantly, the f −α FFT amplitude spectrum curve became affected by the highest sTRs (Figure 4). This has an inevitable effect on measures of signal stationary metrics like Hurst exponent (H) (Bullmore et al., 2001; Wink et al., 2008), fractal dimension (Df) (Kiviniemi et al., 2005; Kiviniemi, 2008) since the power spectral intensity f(I) = f −α, where Df = (3-α)/2 and Df = 2-H. The f/ALFF as well as f −α metrics may then also be sensitive to physiological pulsations. In other words, results based on comparing patients with controls may suffer from cardiorespiratory differences between groups. To avoid these factors, physiological signals should be measured and FFT power spectral analyses should be performed on as short TR data as possible. This is our recommendation for future studies.
Dynamic Connectivity and Effect of TR
The recent discoveries of dynamic functional connectivity analytics of fMRI data show that there are marked changes over time (Hutchinson et al., 2012). Wavelet analyses and time windowed approaches have suggested that the connectivity of regions varies markedly over time (Chang and Glover, 2010; Kiviniemi et al., 2011; Smith et al., 2012; Liu and Duyn, 2013). Targeted averaging algorithms can detect quasi-periodic very low frequency BOLD signal pulses that travel over the RSN and connectivity gradient patterns (Majeed et al., 2011; Pan et al., 2013; Keilholz, 2014; Thompson et al., 2014).
In this study we evaluated the VLF QPP maps to quantify effects of sTR on dynamic BOLD connectivity metrics. First, the FFT analysis indicated that the faster cardiorespiratory pulsations could not be even evaluated as a function of sTR: the cardiac power is not visible above 0.5 s sTR and furthermore the cardiac power aliases over respiratory power. Secondly, the VLF pulse analysis indicated that compared to the most critically sampled 0.1 s TR data, both the QPP strength (i.e., spatial similarity of waves) and inter-individual repeatability of the QPPs become linearly reduced with increasing sTR. In other words, the QPP analysis significantly lost its accuracy as a function of sTR. Taken together, the analysis of dynamic and physiological pulsations benefits when performed on data with TRs < 0.3 s.
As with most of the dynamic connectivity metrics, the problem has been the reducing statistical power of the analysis when shorter and shorter connectivity epochs have been attempted to be analyzed. Reducing dynamic analysis window length reduces also the number of samples and degrees of freedom for statistical inferences. Our results support the use of short TR in the evaluation of dynamic connectivity metrics as (a) for mathematical procedures, short TR gives more time points per analyzed time window, (b) it enables critical sampling and differentiation of physiological pulsations, (c) it avoids aliasing of the physiological signal over the targeted phenomena, (d) enables mapping of modulations of the pulsations, and e) it offers markedly more accurate signal distributions for statistical inferences, c.f. Figure 4D.
Limitations and Strengths of the Study
The results may not be directly comparable due to multiple technical differences between interleaved EPI vs. MREG. However, despite the technical differences, the results in stationary connectivity are highly similar in a wide range of frequencies and RSNs. This study aimed to minimize confounding factors like intra-subject status changes (cardiorespiratory status, vigilance, mood, etc.) as well as technical factors (SNR, sequence parameters, sampling trajectory, scanner model, field strength) issues by using data from the same subjects scanned once with fast TR that was downsampled to higher sTRs. This enables the direct comparison of the effects of sampling rate only. EPI data was scanned for comparing the results to a more conventional technique.
However, the interleaved vs. single shot k-space sampling technique is not the only difference when comparing MREG and EPI data. As mentioned accurately by Golestani and co-workers, comparisons of long vs. short TR measurements need to be taken with a grain of physics, i.e., the EPI vs. MREG also have several other technical differences. In order to be fast, TR is very short in MREG which sensitizes the signal more to T1 inflow and steady state precession effects compared to conventional EPI BOLD (Liu, 2016). On the other hand, TE is relatively long in MREG which makes it T2* -weighted as is EPI BOLD. The low flip angles 5° vs. 15° in MREG and EPI BOLD, respectively, reduces the sensitivity to physiological pulsations on both methods as well (Gonzalez-Castillo et al., 2011). Furthermore, the system noise changes linearly or even quadratically with the TR due to imperfections as the sequence is repeatedly run. All these effects need to be considered while discerning the differences between EPI vs. MREG results and again the direct comparison between them is rather difficult.
According to (Glover, 2012), the spiral readout has reduced sensitivity to motion, shortened readout times, improved signal recovery in most frontal and parietal brain regions, and exhibited blurring artifacts instead of ghosts or geometric distortion. MREG combines Spiral-in/out trajectory which further has diminished susceptibility-induced signal dropout and increased BOLD signal. The EPI readout trajectory is subject to ghosts from off-resonance and gradient imperfections and is intrinsically sensitive to cardiac-induced pulsatile motion from substantial first- and higher order moments of the gradient waveform near the k-space origin (Glover, 2012). So, the artifact and BOLD sensitivity profiles are also different between EPI and MREG. In summary, looking at the similarity of the resting state connectivity measures with different techniques, it seems that human brain connectivity is quite a stable phenomenon that can be measured robustly with technically quite different scanning approaches.
Earlier, the use of two different sampling rates scanned at subsequent scanning sessions have been used to exclude aliasing as a source of very low frequency fluctuations as a source of resting state functional connectivity (Beckmann et al., 2005; Kiviniemi et al., 2005). Theoretically, the previous data that used different scan sessions and so the physiological pulses were not identical between the TRs and therefore not directly comparable. However, in studies with different sampling rate in different scan, also there the results seem to be highly similar to our data (Smith et al., 2013; Golestani et al., 2017).
While attempting to quantify aliasing and dynamic metrics, we used the original 0.1 s TR MREG data as a reference. This may in theory involve a bias since the original MREG data can also detect modulations of cardiorespiratory pulse amplitudes and timing variations, as seen in Figure 1. The slower sTRs have no way of depicting these modulations and therefore a comparison of spatial connectivity and signal stability may be somewhat biased toward undermining the accuracy of the fast data, since it is also sensitivity to modulation effects.
Conclusion
The overall the effect of sampling rate on most commonly used stationary rs-fMRI metrics is minimal. The dominance of the VLF power in gray matter overpowers aliasing effects and enables highly reproducible stationary connectivity results. The aliasing is most dominant between cardiac and respiratory pulsations in central structures near or within CSF spaces. Different technical imaging approaches (e.g., interleaved EPI vs. SOS MREG) yield differential connectivity metrics stemming from multiple spin acquisition differences. Despite these differences, the results from different techniques give fairly good spatial agreement of human brain connectivity. Interleaved scanning of slices seems to introduce inaccuracies in some analyses due to discontinuous sampling of physiological signals. The analysis of dynamic connectivity and frequency based physiological pulsation benefits most from faster scanning.
Data Availability
The datasets generated for this study are available on request to the corresponding author.
Author Contributions
VKi and VKo: designed the study; NH, VKi, VKo, and VR: collected the data; NH, LR, VKi, VKo: analyzed the data; NH, VKi, VKo, JaK, LR, HH, VB, TT, AR and JuK: wrote the manuscript.
Funding
This work was supported by Jane and Aatos Erkko Foundation grant (VKi), Academy of Finland and Aivosäätiö TERVA grant 314497 (VKi), Academy of Finland Grants 123772 & 275342 (VKi), The SalWe Research Program for Mind and Body (Tekes—the Finnish Funding Agency for Technology and Innovation, Grant No. 1104/10) (VKi), Finnish Medical Foundation (VKi, TT), Finnish Neurological Foundation, KEVO grants from Oulu University Hospital (VKi), Epilepsy Research Foundation (JaK), Finnish Cultural Foundation, North Ostrobothnia Regional Fund (JaK), Orion Research Foundation sr (JaK, TT), Tauno Tönning Foundation (JaK) and Health and Biosciences Doctoral programme -grant (HH).
Conflict of Interest Statement
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.
Acknowledgments
Thanks to all study subjects for participation and thanks to CSC (Tieteen Tietotekniikan Keskus) for supporting the analysis tools.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins.2019.00279/full#supplementary-material
Video 1. Mean very low frequency quasi-periodic patterns (VLF QPPs) for repetition time (TR = 0.1 s).
Video 2. Mean very low frequency quasi-periodic patterns (VLF QPPs) for subsampled repetition time (sTR=0.3 s).
Video 3. Mean very low frequency quasi-periodic patterns (VLF QPPs) for subsampled repetition time (sTR = 0.5 s).
Video 4. Mean very low frequency quasi-periodic patterns (VLF QPPs) for subsampled repetition time (sTR=1.0 s).
Video 5. Mean very low frequency quasi-periodic patterns (VLF QPPs) for subsampled repetition time (sTR=1.5 s).
Video 6. Mean very low frequency quasi-periodic patterns (VLF QPPs) for subsampled repetition time (sTR=1.8 s).
Video 7. Mean very low frequency quasi-periodic patterns (VLF QPPs) for subsampled repetition time (sTR=2.2 s).
Video 8. Mean very low frequency quasi-periodic patterns (VLF QPPs) for subsampled repetition time (sTR=3.0 s).
Video 9. Mean very low frequency quasi-periodic patterns (VLF QPPs) for conventional echo-planar imaging (EPI) acquisition with repetition time (TR=2.15 s).
Video 10. Mean very low frequency quasi-periodic patterns (VLF QPPs) for interleaved variant with subsampled repetition time (sTR=2.2 s).
References
Abou-Elseoud, A., Starck, T., Remes, J., Nikkinen, J., Tervonen, O., and Kiviniemi, V. (2010). The effect of model order selection in group PICA. Hum. Brain Mapp. 31, 1207–1216. doi: 10.1002/hbm.20929
Assländer, J., Zahneisen, B., Hugger, T., Reisert, M., Lee, H. L., LeVan, P., et al. (2013). Single shot whole brain imaging using spherical stack of spirals trajectories. Neuroimage 73, 59–70. doi: 10.1016/j.neuroimage.2013.01.065
Beckmann, C. F., DeLuca, M., Devlin, J. T., and Smith, S. M. (2005). Investigations into resting-state connectivity using independent component analysis. Philos. Trans. R. Soc. Lond. B Biol. Sci. 360, 1001–1013. doi: 10.1098/rstb.2005.1634
Beckmann, C. F., Jenkinson, M., Woolrich, M. W., Behrens, T. E., Flitney, D. E., Devlin, J. T., et al. (2006). Applying FSL to the FIAC data: model-based and model-free analysis of voice and sentence repetition priming. Hum. Brain Mapp. 27, 380–391. doi: 10.1002/hbm.20246
Birn, R. M., Diamond, J. B., Smith, M. A., and Bandettini, P. A. (2006). Separating respiratory-variation-related fluctuations from neuronal-activity-related fluctuations in fMRI. Neuroimage 31, 1536–1548. doi: 10.1016/j.neuroimage.2006.02.048
Birn, R. M., Smith, M. A., Jones, T. B., and Bandettini, P. A. (2008). The respiration response function: the temporal dynamics of fMRI signal fluctuations related to changes in respiration. Neuroimage 40, 644–654. doi: 10.1016/j.neuroimage.2007.11.059
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., and Barth, M. (2013). Generalized iNverse imaging (GIN): ultrafast fMRI with physiological noise correction. Magn. Reson. Med. 70, 962–971. doi: 10.1002/mrm.24528
Bright, M. G., and Murphy, K. (2015). Is fMRI “noise” really noise? resting state nuisance regressors remove variance with network structure. Neuroimage 114, 158–169. doi: 10.1016/j.neuroimage.2015.03.070
Bullmore, E., Long, C., Suckling, J., Fadili, J., Calvert, G., Zelaya, F., et al. (2001). Colored noise and computational inference in neurophysiological (fMRI) time series analysis: resampling methods in time and wavelet domains. Hum. Brain Mapp. 12, 61–78. doi: 10.1002/1097-0193(200102)12:2<61::AID-HBM1004>3.0.CO;2-W
Calhoun, V. D., Adali, T., Pearlson, G. D., and Pekar, J. J. (2001). A method for making group inferences from functional MRI data using independent component analysis. Hum. Brain Mapp. 14, 140–151. doi: 10.1002/hbm.1048
Chang, C., and Glover, G. H. (2009). Relationship between respiration, end-tidal CO2, and BOLD signals in resting-state fMRI. Neuroimage 47, 1381–1393. doi: 10.1016/j.neuroimage.2009.04.048
Chang, C., and Glover, G. H. (2010). Time–frequency dynamics of resting-state brain connectivity measured with fMRI. Neuroimage 50, 81–98. doi: 10.1016/j.neuroimage.2009.12.011
Chang, C., Metzger, C. D., Glover, G. H., Duyn, J. H., Heinze, H. J., and Walter, M. (2013a). Association between heart rate variability and fluctuations in resting-state functional connectivity. Neuroimage 68, 93–104. doi: 10.1016/j.neuroimage.2012.11.038
Chang, W. T., Nummenmaa, A., Witzel, T., Ahveninen, J., Huang, S., Tsai, K. W., et al. (2013b). Whole-head rapid fMRI acquisition using echo-shifted magnetic resonance inverse imaging. Neuroimage 78, 325–338. doi: 10.1016/j.neuroimage.2013.03.040
Chao-Gan, Y., and Yu-Feng, Z. (2010). DPARSF: a MATLAB toolbox for“ pipeline” data analysis of resting-state fMRI. Front. Syst. Neurosci. 4:13. doi: 10.3389/fnsys.2010.00013
Chen, J. E., Polimeni, J. R., Bollmann, S., and Glover, G. H. (2019). On the analysis of rapidly sampled fMRI data. Neuroimage. 188, 807–820. doi: 10.1016/j.neuroimage.2019.02.008
Cordes, D., Nandy, R. R., Schafer, S., and Wager, T. D. (2014). Characterization and reduction of cardiac-and respiratory-induced noise as a function of the sampling rate (TR) in fMRI. Neuroimage 89, 314–330. doi: 10.1016/j.neuroimage.2013.12.013
De Luca, M., Beckmann, C. F., De Stefano, N., Matthews, P. M., and Smith, S. M. (2006). fMRI resting state networks define distinct modes of long-distance interactions in the human brain. Neuroimage 29, 1359–1367. doi: 10.1016/j.neuroimage.2005.08.035
Demetriou, L., Kowalczyk, O. S., Tyson, G., Bello, T., Newbould, R. D., and Wall, M. B. (2018). A comprehensive evaluation of increasing temporal resolution with multiband-accelerated protocols and effects on statistical outcome measures in fMRI. Neuroimage 176, 404–416. doi: 10.1016/j.neuroimage.2018.05.011
Glover, G. H. (2012). Spiral imaging in fMRI. Neuroimage 62, 706–712. doi: 10.1016/j.neuroimage.2011.10.039
Golestani, A. M., Kwinta, J. B., Khatamian, Y. B., and Chen, J. J. (2017). The effect of low-frequency physiological correction on the reproducibility and specificity of resting-state fMRI metrics: functional connectivity, ALFF, and ReHo. Front. Neurosci. 11:546. doi: 10.3389/fnins.2017.00546
Gonzalez-Castillo, J., Roopchansingh, V., Bandettini, P. A., and Bodurka, J. (2011). Physiological noise effects on the flip angle selection in BOLD fMRI. Neuroimage 54, 2764–2778. doi: 10.1016/j.neuroimage.2010.11.020
Griffanti, L., Dipasquale, O., Laganà, M. M., Nemni, R., Clerici, M., Smith, S. M., et al. (2015). Effective artifact removal in resting state fMRI data improves detection of DMN functional connectivity alteration in Alzheimer's disease. Front. Hum. Neurosci. 9:449. doi: 10.3389/fnhum.2015.00449
Hugger, T., Zahneisen, B., LeVan, P., Lee, K. J., Lee, H. L., Zaitsev, M., et al. (2011). Fast undersampled functional magnetic resonance imaging using nonlinear regularized parallel image reconstruction. PLoS ONE 6:e28822. doi: 10.1371/journal.pone.0028822
Hutchinson, J. B., Uncapher, M. R., Weiner, K. S., Bressler, D. W., Silver, M. A., Preston, A. R., et al. (2012). Functional heterogeneity in posterior parietal cortex across attention and episodic memory retrieval. Cereb. Cortex 24, 49–66. doi: 10.1093/cercor/bhs278
Iliff, J. J., Wang, M., Liao, Y., Plogg, B. A., Peng, W., Gundersen, G. A., et al. (2012). A paravascular pathway facilitates CSF flow through the brain parenchyma and the clearance of interstitial solutes, including amyloid β. Sci. Transl. Med. 4:147ra111. doi: 10.1126/scitranslmed.3003748
Jenkinson, M., Bannister, P., Brady, M., and Smith, S. (2002). Improved optimization for the robust and accurate linear registration and motion correction of brain images. Neuroimage 17, 825–841. doi: 10.1006/nimg.2002.1132
Jenkinson, M., Beckmann, C. F., Behrens, T. E., Woolrich, M. W., and Smith, S. M. (2012). Fsl. Neuroimage 62, 782–790. doi: 10.1016/j.neuroimage.2011.09.015
Jenkinson, M., and Smith, S. (2001). A global optimisation method for robust affine registration of brain images. Med. Image Anal. 5, 143–156. doi: 10.1016/S1361-8415(01)00036-6
Jessen, N. A., Munk, A. S., Lundgaard, I., and Nedergaard, M. (2015). The glymphatic system: a beginner's guide. Neurochem. Res. 40, 2583–2599. doi: 10.1007/s11064-015-1581-6
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
Keilholz, S. D. (2014). The neural basis of time-varying resting-state functional connectivity. Brain Connect. 4, 769–779. doi: 10.1089/brain.2014.0250
Kiviniemi, V. (2008). Endogenous brain fluctuations and diagnostic imaging. Hum. Brain Mapp. 29, 810–817. doi: 10.1002/hbm.20582
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., Kantola, J. H., Jauhiainen, J., Hyvärinen, A., and Tervonen, O. (2003). Independent component analysis of nondeterministic fMRI signal sources. Neuroimage 19, 253–260. doi: 10.1016/S1053-8119(03)00097-1
Kiviniemi, V., Ruohonen, J., and Tervonen, O. (2005). Separation of physiological very low frequency fluctuation from aliasing by switched sampling interval fMRI scans. Magn. Reson. Imaging 23, 41–46. doi: 10.1016/j.mri.2004.09.005
Kiviniemi, V., Starck, T., Remes, J., Long, X., Nikkinen, J., Haapea, M., et al. (2009). Functional segmentation of the brain cortex using high model order group PICA. Hum. Brain Mapp. 30, 3865–3886. doi: 10.1002/hbm.20813
Kiviniemi, V., Vire, T., Remes, J., Elseoud, A. A., Starck, T., Tervonen, O., et al. (2011). A sliding time-window ICA reveals spatial variability of the default mode network in time. Brain Connect. 1, 339–347. doi: 10.1089/brain.2011.0036
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
Korhonen, V., Hiltunen, T., Myllylä, T., Wang, X., Kantola, J., Nikkinen, J., et al. (2014). Synchronous multiscale neuroimaging environment for critically sampled physiological analysis of brain function: hepta-scan concept. Brain Connect. 4, 677–689. doi: 10.1089/brain.2014.0258
Lee, H. L., Zahneisen, B., Hugger, T., LeVan, P., and Hennig, J. (2013). Tracking dynamic resting-state networks at higher frequencies using MR-encephalography. Neuroimage 65, 216–222. doi: 10.1016/j.neuroimage.2012.10.015
Lin, F. H., Polimeni, J. R., Lin, J. L., Tsai, K. W., Chu, Y. H., Wu, P. Y., et al. (2018). Relative latency and temporal variability of hemodynamic responses at the human primary visual cortex. Neuroimage 164, 194–201. doi: 10.1016/j.neuroimage.2017.01.041
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
Liu, T. T. (2016). Noise contributions to the fMRI signal: an overview. Neuroimage 143, 141–151. doi: 10.1016/j.neuroimage.2016.09.008
Liu, X., De Zwart, J. A., Schölvinck, M. L., Chang, C., Ye, F. Q., Leopold, D. A., et al. (2018). Subcortical evidence for a contribution of arousal to fMRI studies of brain activity. Nat. Commun. 9:395. doi: 10.1038/s41467-017-02815-3
Liu, X., and Duyn, J. H. (2013). Time-varying functional network information extracted from brief instances of spontaneous brain activity. Proc. Natl. Acad. Sci. U.S.A. 110, 4392–4397. doi: 10.1073/pnas.1216856110
Ma, Y., Shaik, M. A., Kozberg, M. G., Kim, S. H., Portes, J. P., Timerman, D., et al. (2016). Resting-state hemodynamics are spatiotemporally coupled to synchronized and symmetric neural activity in excitatory neurons. Proc. Natl. Acad. Sci. U.S.A. 113, E8463–E8471. doi: 10.1073/pnas.1525369113
Majeed, W., Magnuson, M., Hasenkamp, W., Schwarb, H., Schumacher, E. H., Barsalou, L., et al. (2011). Spatiotemporal dynamics of low frequency BOLD fluctuations in rats and humans. Neuroimage 54, 1140–1150. doi: 10.1016/j.neuroimage.2010.08.030
Makedonov, I., Black, S. E., and MacIntosh, B. J. (2013). BOLD fMRI in the white matter as a marker of aging and small vessel disease. PLoS ONE 8:e67652. doi: 10.1371/journal.pone.0067652
Makedonov, I., Chen, J. J., Masellis, M., MacIntosh, B. J., and Initiative, A. D. N. (2016). Physiological fluctuations in white matter are increased in Alzheimer's disease and correlate with neuroimaging and cognitive biomarkers. Neurobiol. Aging 37, 12–18. doi: 10.1016/j.neurobiolaging.2015.09.010
Murphy, K., and Fox, M. D. (2017). Towards a consensus regarding global signal regression for resting state functional connectivity MRI. Neuroimage 154, 169–173. doi: 10.1016/j.neuroimage.2016.11.052
Nedergaard, M. (2013). Garbage truck of the brain. Science 340, 1529–1530. doi: 10.1126/science.1240514
Palva, J. M., Zhigalov, A., Hirvonen, J., Korhonen, O., Linkenkaer-Hansen, K., and Palva, S. (2013). Neuronal long-range temporal correlations and avalanche dynamics are correlated with behavioral scaling laws. Proc. Natl. Acad. Sci. U.S.A. 110, 3585–3590. doi: 10.1073/pnas.1216855110
Pan, W. J., Thompson, G. J., Magnuson, M. E., Jaeger, D., and Keilholz, S. (2013). Infraslow LFP correlates to resting-state fMRI BOLD signals. Neuroimage 74, 288–297. doi: 10.1016/j.neuroimage.2013.02.035
Peltier, S. J., Polk, T. A., and Noll, D. C. (2003). Detecting low-frequency functional connectivity in fMRI using a self-organizing map (SOM) algorithm. Hum. Brain Mapp. 20, 220–226. doi: 10.1002/hbm.10144
Pfeuffer, J., Van de Moortele, P. F., Ugurbil, K., Hu, X., and Glover, G. H. (2002). Correction of physiologically induced global off-resonance effects in dynamic echo-planar and spiral functional imaging. Magn. Reson. Med. 47, 344–353. doi: 10.1002/mrm.10065
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
Preibisch, C., Castrillón G J, G., Bührer, M., and Riedl, V. (2015). Evaluation of multiband EPI acquisitions for resting state fMRI. PLoS ONE 10:e0136961. doi: 10.1371/journal.pone.0136961
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., Korhonen, V., Huotari, N., Raatikainen, V., Hautaniemi, T., Kananen, J., et al. (2018). Breath hold effect on cardiovascular brain pulsations–A multimodal magnetic resonance encephalography study. J. Cereb. Blood Flow Metab. 11:271678X18798441. doi: 10.1177/0271678X18798441
Rajna, Z., Kananen, J., Keskinarkaus, A., Seppänen, T., and Kiviniemi, V. (2015). Detection of short-term activity avalanches in human brain default mode network with ultrafast MR encephalography. Front. Hum. Neurosci. 9:448. doi: 10.3389/fnhum.2015.00448
Shmueli, K., van Gelderen, P., de Zwart, J. A., Horovitz, S. G., Fukunaga, M., Jansma, J. M., et al. (2007). Low-frequency fluctuations in the cardiac rate as a source of variance in the resting-state fMRI BOLD signal. Neuroimage 38, 306–320. doi: 10.1016/j.neuroimage.2007.07.037
Smith, S. M., Beckmann, C. F., Andersson, J., Auerbach, E. J., Bijsterbosch, J., Douaud, G., et al. (2013). Resting-state fMRI in the human connectome project. Neuroimage 80, 144–168. doi: 10.1016/j.neuroimage.2013.05.039
Smith, S. M., Fox, P. T., Miller, K. L., Glahn, D. C., Fox, P. M., Mackay, C. E., et al. (2009). Correspondence of the brain's functional architecture during activation and rest. Proc. Natl. Acad. Sci. U.S.A. 106, 13040–13045. doi: 10.1073/pnas.0905267106
Smith, S. M., Jenkinson, M., Woolrich, M. W., Beckmann, C. F., Behrens, T. E., Johansen-Berg, H., et al. (2004). Advances in functional and structural MR image analysis and implementation as FSL. Neuroimage 23, S208–S219. doi: 10.1016/j.neuroimage.2004.07.051
Smith, S. M., Miller, K. L., Moeller, S., Xu, J., Auerbach, E. J., Woolrich, M. W., et al. (2012). Temporally-independent functional modes of spontaneous brain activity. Proc. Natl. Acad. Sci. U.S.A. 109, 3131–3136. doi: 10.1073/pnas.1121329109
Thayer, J. F., Yamamoto, S. S., and Brosschot, J. F. (2010). The relationship of autonomic imbalance, heart rate variability and cardiovascular disease risk factors. Int. J. Cardiol. 141, 122–131. doi: 10.1016/j.ijcard.2009.09.543
Thompson, G. J., Pan, W. J., Magnuson, M. E., Jaeger, D., and Keilholz, S. D. (2014). Quasi-periodic patterns (QPP): large-scale dynamics in resting state fMRI that correlate with local infraslow electrical activity. Neuroimage 84, 1018–1031. doi: 10.1016/j.neuroimage.2013.09.029
Tuovinen, T., Kananen, J., Rytty, R., Moilanen, V., Elseoud, A. A., Remes, A. M., et al. (2018). Altered BOLD signal variation in Alzheimer's disease and frontotemporal dementia. BioRxiv [preprint] 455683, 1–47. doi: 10.1101/455683
Tuovinen, T., Rytty, R., Moilanen, V., Abou Elseoud, A., Veijola, J., Remes, A. M., et al. (2017). The effect of gray matter ICA and coefficient of variation mapping of BOLD data on the detection of functional connectivity changes in Alzheimer's disease and bvFTD. Front. Hum. Neurosci. 10:680. doi: 10.3389/fnhum.2016.00680
van der Kooy, K. G., van Hout, H. P., van Marwijk, H. W., de Haan, M., Stehouwer, C. D., and Beekman, A. T. (2006). Differences in heart rate variability between depressed and non-depressed elderly. Int. J. Geriatr. Psychiatry 21, 147–150. doi: 10.1002/gps.1439
Vidaurre, D., Smith, S. M., and Woolrich, M. W. (2017). Brain network dynamics are hierarchically organized in time. Proc. Natl. Acad. Sci. U.S.A. 114, 12827–12832. doi: 10.1073/pnas.1705120114
Wink, A. M., Bullmore, E., Barnes, A., Bernard, F., and Suckling, J. (2008). Monofractal and multifractal dynamics of low frequency endogenous brain oscillations in functional MRI. Hum. Brain Mapp. 29, 791–801. doi: 10.1002/hbm.20593
Wise, R. G., Ide, K., Poulin, M. J., and Tracey, I. (2004). Resting fluctuations in arterial carbon dioxide induce significant low frequency variations in BOLD signal. Neuroimage 21, 1652–1664. doi: 10.1016/j.neuroimage.2003.11.025
Woolrich, M. W., Ripley, B. D., Brady, M., and Smith, S. M. (2001). Temporal autocorrelation in univariate linear modeling of FMRI data. Neuroimage 14, 1370–1386. doi: 10.1006/nimg.2001.0931
Zahneisen, B., Assländer, J., LeVan, P., Hugger, T., Reisert, M., Ernst, T., et al. (2014). Quantification and correction of respiration induced dynamic field map changes in fMRI using 3D single shot techniques. Magn. Reson. Med. 71, 1093–1102. doi: 10.1002/mrm.24771
Zahneisen, B., Hugger, T., Lee, K. J., LeVan, P., Reisert, M., Lee, H. L., et al. (2012). Single shot concentric shells trajectories for ultra fast fMRI. Magn. Reson. Med. 68, 484–494. doi: 10.1002/mrm.23256
Zang, Y., Jiang, T., Lu, Y., He, Y., and Tian, L. (2004). Regional homogeneity approach to fMRI data analysis. Neuroimage 22, 394–400. doi: 10.1016/j.neuroimage.2003.12.030
Zang, Y. F., He, Y., Zhu, C. Z., Cao, Q. J., Sui, M. Q., Liang, M., et al. (2007). Altered baseline brain activity in children with ADHD revealed by resting-state functional MRI. Brain Dev. 29, 83–91. doi: 10.1016/j.braindev.2006.07.002
Keywords: resting state, magnetic resonance encephalography, aliasing, pulsations, quasi-periodic patterns
Citation: Huotari N, Raitamaa L, Helakari H, Kananen J, Raatikainen V, Rasila A, Tuovinen T, Kantola J, Borchardt V, Kiviniemi VJ and Korhonen VO (2019) Sampling Rate Effects on Resting State fMRI Metrics. Front. Neurosci. 13:279. doi: 10.3389/fnins.2019.00279
Received: 26 October 2018; Accepted: 08 March 2019;
Published: 02 April 2019.
Edited by:
Jean Chen, University of Toronto, CanadaReviewed by:
Fahmeed Hyder, Yale University, United StatesStefan Posse, University of New Mexico, United States
Copyright © 2019 Huotari, Raitamaa, Helakari, Kananen, Raatikainen, Rasila, Tuovinen, Kantola, Borchardt, Kiviniemi and Korhonen. 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: Niko Huotari, niko.huotari@oulu.fi