Search tips
Search criteria 


Neuroimage. 2017 November 15; 162: 93–105.
PMCID: PMC5711428

Cardiac cycle-induced EPI time series fluctuations in the brain: Their temporal shifts, inflow effects and T2* fluctuations


The cardiac-induced arterial pressure wave causes changes in cerebral blood flow velocities and volumes that affect the signals in echo-planar imaging (EPI). Using single-echo EPI time series data, acquired fast enough to unalias the cardiac frequency, we found that the cardiac cycle-induced signal fluctuations are delayed differentially in different brain regions. When referenced to the time series in larger arterial structures, the cortical voxels are only minimally shifted but significant shifts are observed in subcortical areas. Using double-echo EPI data we mapped the voxels’ “signal at zero echo time”, S0, and apparent T2 over the cardiac cycle. S0 pulsatility was maximised for voxels with a cardiac cycle-induced timing that was close to the arterial structures and is likely explained by enhanced inflow effects in the cortical areas compared to subcortical areas. Interestingly a consistent T2 waveform over the cardiac cycle was observed in all voxels with average amplitude ranges between 0.3-0.55% in grey matter and 0.15–0.22% in white matter. The timing of the T2 waveforms suggests a partial volume fluctuation where arteriolar blood volume changes are counterbalanced by changes in CSF volumes.

Keywords: Ultra-fast EPI, Multiband, Multi-echo, Cardiac pulsatility, Partial volume, Physiological noise, Cerebral compliance

1. Introduction

Echo-planar imaging (EPI) with sampling after multiband excitation achieves temporal resolution in the sub-second regime whilst maintaining a good slice coverage and spatial resolution (Feinberg and Setsompop, 2013, Feinberg et al., 2010, Moeller et al., 2010). In combination with improved signal-to-noise ratio (SNR) at ultra-high magnetic field, these fast acquisition techniques facilitate the study of spatio-temporal phenomena at finer scales and higher frequencies (Setsompop et al., 2016). The ability to depict higher frequency regimes without alias turns cardiac cycle-induced fluctuations, typically recorded as nuisance signal, from “physiological noise” into a valuable signal that may carry information about cerebrovascular mechanisms. So far, studies have addressed the spatial pattern of cardiac cycle-related fluctuations in EPI time series (Dagli et al., 1999, Tong and Frederick, 2014, Kiviniemi et al., 2016), but little is known about the temporal signatures and underpinning contrast mechanisms. Generically referred to as “pulsatility” these fluctuations are of interest for research into cerebral compliance in response to the arterial pressure wave and the impact of pulsatility changes on tissue (Robertson et al., 2010, Webb et al., 2012). A better understanding of cardiac cycle-induced EPI signals might further be of use to improve fMRI physiological noise cleaning.

This work aims to investigate the timing of cardiac cycle-induced EPI signal fluctuations in the brain, and their associated MR parameters. Here, the voxel's signal at zero-echo time, S0, and its apparent transverse decay constant, T2, underpin these fluctuations. Two types of data set were acquired:

  • 1.
    Single-echo EPI with a repetition time (TR) of 328 ms. These scans recorded the unaliased EPI time series up to the cardiac frequency regime. We used these data to calculate the temporal shifts in the cardiac cycle-induced signal fluctuations.
  • 2.
    Double-echo EPI. These scans were acquired to fit S0 and the apparent T2 values. Each measurement time point was phase-locked to an externally measured cardiac trigger signal to subsequently create voxel-wise S0 and T2 waveforms over the cardiac unit cycle.

The single-echo EPI data revealed temporal shifts in the cardiac cycle-related time series between and within tissue types and between brain areas, which we refer to as “cardiac phase shift”. We relate the S0 and T2 pulsatility over the cardiac cycle from the double-echo EPI data to these shifts. For the remainder of this text we will refer to “pulsatility” as the fluctuation of a variable around its mean value over the course of the cardiac period. We will further use the term “phase” exclusively to denote the temporal position within the cardiac cycle (not to be confused with complex k-space or image phase).

2. Theory

The arterial pressure wave propagation speed is on the order of 10 m/s (Asmar et al., 1995). It is therefore beyond the temporal resolution of EPI time series and can be regarded as being experienced by all parts of the brain instantaneously. However, the brain is confined by the skull and the pressure wave is tightly coupled to much slower changes in local cerebral blood volume (CBV) and flow velocity (CBFV) (Wagshul et al., 2011, Van De Vosse and Stergiopulos, 2011). These CBV and CBFV changes affect the MR signal and manifest in cardiac frequency fluctuations in the EPI time series.

2.1. EPI signal model

Let τc be the systole-to-systole time period that spans the cardiac unit cycle τ with τ(0,,τc]. We assume that cardiac cycle-induced fluctuations are independent of the absolute cycle duration, but nevertheless have a fixed phase relationship to the unit cycle (Hu et al., 1995). The voxel signal is a combination of intra- and extravascular signals (Marques and Bowtell, 2008) and can be expressed as a multi-exponential decay (Obata et al., 2004, Zhao et al., 2007, Uludag et al., 2009)


where TE denotes the echo time and i the compartments of grey matter, white matter, arterioles, capillaries, venules and CSF. vi, ρi and Mz,i are the respective partial volume fractions, proton densities and longitudinal magnetisations. We explicitly stated vi as a function of the cardiac cycle τ as blood volumes within a voxel may change from systole to diastole. The magnetisation Mz,i(τ) will vary with inflow effects (see discussion of S0 in the following section). Further, Mz,i could fluctuate with variations in the longitudinal relaxation constant T1 due to oxygenation changes. Eq. (1) assumes that T2,i for a given compartment does not vary over the cardiac cycle, as baseline oxygen extraction rates are not related to the cardiac cycle to our knowledge.

An experimental measurement of all components in the multi-exponential model is difficult and often a simplified mono-exponential decay is used instead (Zhao et al., 2007, Bianciardi et al., 2016)


The above model allows a measurement of the voxel's S0 and apparent T2 using double- or multi-echo data. An approximation SmonoSmulti will be limited to voxels with small partial volume contamination. If partial volume contamination is large, e.g. in a voxel consisting of half grey matter and half CSF, the mono-exponential decay becomes inaccurate because of the large difference in T2 values between the compartments. We use simulations to calculate the theoretical error of estimating apparent S0 and T2 values from the mono-exponential fit for different partial volume combinations and test to what extent an approximation SmonoSmulti is valid (see Methods).

Different physiological mechanisms cause fluctuations in the S0(τ) and T2(τ) parameters. Inflow effects and partial volume fluctuations will manifest in S0(τ). This can be written more explicitly as S0(T1,M0,TR,α,CBV(τ),CBFV(τ)) (Bianciardi et al., 2016), where M0 is the equilibrium magnetisation and α the radio-frequency (RF) flip angle. CBV(τ) and CBFV(τ) determine the amount of already RF-excited spins being replaced with fresh ones in the voxel. CBFV pulsatility will only affect S0 if the spins experience multiple RF-excitations. In a classical single-slice acquisition, this is determined by the minimum blood velocity to traverse the slice. This simple theory will not hold for multiband acquisitions, where slower spins can experience fewer RF-pulses than faster ones if their position is out of synchrony with the slice acquisition order. Complicated vascular pathways through the slices require modelling of the vascular tree which is beyond the scope of this paper. However, it must be acknowledged that inflow-related magnetization enhancement due to CBFV still occurs with multiband excitation. In our experiments we used TRs of 328 ms (single-echo data) and 405 ms (double-echo data). Blood spins in large arterial structures, such as the middle cerebral artery (MCA) travel at about half a metre per second (Meckel et al., 2013) and traverse the imaging slab within one TR. Hence, S0 pulsatility in voxels containing these structures should be induced by CBV and not CBFV pulsatility (see Bianciardi et al. for a detailed discussion (Bianciardi et al., 2016). Blood flow velocities are reduced in the smaller arterial branches, for example velocities in arteries in the basal ganglia are on the order of 3–6 cm/s (Bouvy et al., 2016) and here a mixed contribution of CBV and CBFV pulsatility to the inflow S0 effect is expected. Blood flow velocities at the arteriolar and capillary level are on the order of less than a millimetre per second and spins will be saturated, rendering inflow effects negligible.

Arterioles, capillaries and venules should not alter their oxygen extraction rate between systole and diastole, i.e. single compartment T2-values are assumed to be constant at each vascular level (although oxygen saturation will decrease down the vascular tree). However, we do expect the intravoxel microvascular and CSF partial volumes to fluctuate from systole to diastole, which should change the apparent voxel T2.

2.2. Temporal shifts in the cardiac cycle-related signals

Fig. 1 shows a schematic presentation of the cardiac cycle-induced shifts in EPI time courses. A shift at the cardiac frequency fc between two voxels can be obtained via the phase angle of their Fourier-transformed complex signals (more details will be presented in the Methods section). These phase shifts are likely related to underlying S0(τ) and T2(τ) variations that differ between tissues and brain areas, depending on levels of vascularisation and timing of cardiac cycle-induced variations in inflow and partial volume effects. Additionally, cardiac cycle-induced tissue displacement varies in magnitude and timing in different brain areas. Non-rigid brain motion is known to be maximal in the brain stem and propagates radially towards the outer parts of the brain (Poncelet et al., 1992, Zhong et al., 2009, Soellinger et al., 2009) and could cause timing variations in the EPI signals.

Fig. 1
The temporal shift in the cardiac cycle-related EPI time series between voxels r and k can be defined by the phase angle ϕ between their complex signals S~ (obtained via Fourier transform) at the cardiac frequency fc.

3. Methods

3.1. Single-echo EPI: Cardiac phase mapping

We calculated each voxel's phase ϕ at the cardiac frequency fc using the following steps:

  • 1.
    A high temporal resolution external cardiac trace (e.g. pulse-oximetry) is measured during fMRI acquisition. It is sub-sampled to the acquired TR and then Fourier-transformed to obtain the expected cardiac power spectrum. The cardiac frequency fc is then defined as the component with the maximum power amplitude.
  • 2.
    The EPI time series S(t) in each voxel is Fourier-transformed into the frequency domain to yield S~(f).
  • 3.
    An arbitrary voxel r is selected to serve as a fixed phase reference (absolute phase values are re-referenced to the dominant value in an arterial mask at a later stage, see Post-processing section).
  • 4.
    The phase ϕk at voxel k is calculated as the angle between the complex values using the four-quadrant inverse tangent:



where ϕk,r is the phase shift induced by the timing offset Δtk,r in the acquisition between the slices of voxel k and the reference voxel r.

All steps were performed in MATLAB, R2016a (The MathWorks, Inc., Natick, MA, USA) and the code is provided as PhaseMapping.m. We used simulated data first to test the phase mapping code, in particular to guarantee a correct re-shift of the slice timing-induced phase offsets ϕk,r. A detailed description of the code testing is provided in the Supplementary Material (S1).

3.2. Double-echo EPI: S0 and T2 mapping

3.2.1. S0 and T2 calculation from double-echo data

The mono-exponential signal model Smono(τ) can be discretized into a fixed number, n, of equidistant bins:


A trigger signal sent from the scanner at each RF-excitation allows the post-ordering of measurement points to their appropriate acquisition bin. S0 and T2 are obtained from a first-order polynomial fit of the double-echo data to the linear form of Eq. (4):


S0 values were assigned into the bin closest to the time of the RF trigger (“zero-echo time”). T2 values were assigned to the bin closest to the time point halfway between the first and second echo (TE1 + TE2) / 2. Again, the multiband-specific timing needs to be taken into account for data ordering. An example timing scenario for a double-echo multiband dataset is provided in Fig. 2. The flip angle for this measurement was set to 90 to maximise the sensitivity to inflow effects (see (Bianciardi et al., 2016)). With T2(τ) and S0(τ) fluctuation maps calculated, we quantified pulsatility in a voxel as the maximum fluctuation amplitude around the mean in percent

  • T2,p was defined similarly.
Fig. 2
Cardiac cycle segmentation and timing of a double-echo multiband sequence. Timings are shown for bands 1 and 3. The pulse-ox signal (orange) marks the systole-to-systole period τc of the cardiac unit cycle that is segmented into equidistant bins ...

3.2.2. Validity of mono-exponential vs. multi-exponential signal model

As described above we used the mono-exponential signal model Smono to fit S0 and T2 for each voxel from double-echo data. Ideally a multi-compartment model Smulti as in Eq. (1) is preferred, but more echoes would be needed yielding a trade-off in resolution for the same readout duration. Therefore we used simulations to test the limitations of the model approximation SmonoSmulti for different partial volume combinations in a grey matter voxel. Important here is the TE=0 intercept error which is the S0 estimate. Variations in T2 will introduce an apparent S0 fluctuation. Similarly, SNR limits the fitting accuracy of T2 and S0. To quantify the size of the coupling effect between T2 and S0 fitting and SNR we performed the following steps in MATLAB (the code is provided as MultiVsMono.m):

  • 1.
    Smulti decays were simulated (after Eq. (1), but with fixed parameters) for a grey matter voxel for different partial volume fractions vGM. Further, the voxel was assumed to contain microvascular volumes (see Table 1) for a typical cortical voxel and the remaining volume was taken up by CSF, so that ivi=1. This is the most extreme case of deviating T2 values and we expect the effect to be smaller in a voxel that consists of a mixture of grey and white matter, but no CSF. The parameters used to simulate Smulti were taken from the literature and are listed in Table 1.
    Table 1
    List of the parameters that were used in the simulations of the multi-exponential model Smulti (Eq. (1)).
  • 2.
    The simulated Smulti was then sampled at the echo times that were used in the experiments (TEs=[21.4,52.1] ms at a field strength of 7T), and Gaussian noise, scaled by the SNR value for each echo time, was added (the SNR for each echo was evaluated from the experimentally acquired double-echo data by taking the mean signal in the grey matter mask and dividing it by the standard deviation of a ghost-free area in the image corner). This was repeated 50 times, which was the average number of time points we had sampled for each cardiac cycle bin (total measurement points were 500, split into 10 bins, see MR-data and post-processing), to create a noisy array of simulated data points.
  • 3.
    We fitted S0 and T2 to the synthesised noisy data using the linearised mono-exponential model (Eq. (5)).
  • 4.
    The mono-exponential signal (Eq. (2)) was then calculated with the fitted S0 and T2 to assess the deviation from the noise-free multi-exponential signal.
  • 5.
    To test the coupling of T2 and S0 the above steps were repeated for a grey matter T2 variation of ±1%, which captures the maximum T2 pulsatility we observed in our results.

We further ran the simulation without the addition of noise. This case quantifies the coupling of T2 variations and the S0 fit in the case of noise-free data or “infinite” sampling and is as such the fundamental limit of the approximation SmonoSmulti. All the above steps were run for three partial volume scenarios vGM = 0.5,0.9 and 0.94. The first scenario represents the case of highest partial voluming, the second case matches the partial volume estimate threshold that we used in our tissue masking (see Post-processing section). The last case is at the limit of a pure grey matter voxel, assuming that the microvascular volume is 0.055 for a cortical voxel (see Table 1), leaving only a CSF volume fraction of 0.005.

4. MR data

We scanned 8 healthy subjects (age range: 23–34 ys) under institutional ethical approval. The scans were acquired on a 7T whole-body scanner (Magnetom 7T, Siemens, Erlangen, Germany) using a 1 channel Tx/32 channel Rx head coil (Nova Medical, Wilmington, MA, USA). Physiological signals were measured throughout the scans using a respiratory belt and a pulse-oximetry finger clip, recorded at 1 kHz (BioPac Systems, Goleta, CA, USA). The following scans were acquired during the session:

Tissue segmentation: T1-MPRAGE

T1-MPRAGE: TR = 2200 ms; inversion time, TI = 1050 ms; TE = 2.96 ms; flip angle, α=7; nominal resolution: 0.7 × 0.7 × 0.7 mm3; bandwidth, BW = 240 Hz/px; GRAPPA factor = 2; field of view, FOVread = 224 mm (read direction), FOVphase = 100% (phase direction), slices = 256, scan time 6 min.

Cardiac phase mapping: single-echo EPI

Multiband EPI (Moeller et al., 2010, Feinberg et al., 2010, Xu et al., 2013): TR = 334 ms; TE = 26 ms; inter-echo spacing in EPI echo train = 0.68 ms; α = 33°; multiband factor = 6; GRAPPA factor = 2; BW = 1852 Hz/px; nominal resolution: 2.2 × 2.2 × 2.2 mm3; FOVread = 200 mm; FOVphase = 100%; no partial Fourier acquisition; CAIPIRINHA FOV shifting applied; slices = 36; volumes = 2200; scan time 12 min. Subjects were instructed to stay awake with eyes open.

S0 and T2 mapping: double-echo EPI

Multiband EPI: TR = 405 ms; TE1 / TE2 = 21.4/52.1 ms; inter-echo spacing in EPI echo train = 0.65 ms; α = 90°; multiband factor = 3; GRAPPA factor = 2; BW = 2085 Hz/px; nominal resolution: 2.2 × 2.2 × 2.2 mm3; FOVread = 200 mm; FOVphase = 100%; no partial Fourier acquisition; CAIPIRINHA FOV shifting applied; slices = 15; volumes = 500; scan time 3.3 min. Subjects were instructed to stay awake with eyes open.

Additionally we acquired a GRE-based field map.

5. Post-processing

5.1. Tissue segmentation: T1-MPRAGE

The T1-weighted (T1w) image was run through the fsl_anat pipeline in FSL which performs skull stripping, bias field correction and calculates PVE maps (Jenkinson et al., 2012). An arterial mask was created from the T1w scan using its inherent bright signal of fast arterial blood in the larger arteries (Grinstead et al., 2010). A threshold was chosen on a subject-by-subject basis to create a binary arterial mask.

5.2. Cardiac phase mapping: single-echo EPI

The data were processed through FSL FEAT, including skull stripping, motion correction, field map correction, de-meaning and registration to the T1w structural and to MNI standard space. Appropriate transformations were applied to the PVE maps and the arterial masks for analysis in EPI space. The processed data were read into MATLAB to perform cardiac phase mapping as described in the Methods section. We created a “cardiac mask” which is a binary mask containing only those voxels that have EPI signals with sufficient cardiac modulation. This cardiac mask was derived from a general linear model fit between the voxel's EPI power and the pulse-ox power spectra (a detailed description is provided in the Supplementary Material S2), rejecting voxels with p-values above 0.05. PVE maps were binarised at a value of 0.9 to create tissue masks with minimal partial voluming. Finally the cardiac phase map was multiplied by the binary cardiac and tissue masks to produce tissue-specific cardiac phase distributions. The high PVE threshold for each mask and the targeting by the cardiac binary mask help ensure that we evaluate voxels that can be approximated with the mono-exponential signal model (see Results section 6.2.1) and that show sufficient signal modulation with the cardiac cycle. Exemplar tissue masks with sufficient cardiac signal modulation are shown in Fig. 3.

Fig. 3
Tissue-specific cardiac masks in a representative subject for a lower (a) and upper slice (b). The T1w-image is displayed in EPI space. The red and blue voxels mark grey and white matter voxels with a PVE of at least 0.9 (i.e. minimal partial voluming) ...

The absolute phase value in a subject is arbitrary and phase distributions between subjects will appear shifted. To synchronise phases between subjects we shifted each subject's phases such that the predominant phase value in the arterial mask distribution matched the value − π/2 (see arrow in Fig. 4a) in the Results section). For the remainder we will refer to this peak as the “dominant arterial phase peak”.

Fig. 4
a)c) Distributions of the cardiac frequency phase shift in the EPI time series in a single representative subject in the arterial, grey matter and white matter masks. d) is a whole brain control distribution at a non-significant frequency component, ...

We also created a control distribution to confirm that the phase mapping approach produces random phases for a low power frequency that only contributes to thermal noise (derived from a whole brain mask).

In addition to the tissue-specific analysis we evaluated phase values in different brain structures. To achieve this, we warped the phase maps and the cardiac binary masks into MNI space. Brain area masks from the Harvard-Oxford atlases (included in FSL) were thresholded at a probability of 0.5 and mean phases were calculated in those masks. The mean phase in an MNI voxel was calculated by converting the phase value to Cartesian coordinates, summing over each Cartesian component and reconverting to polar coordinates. In MATLAB this is implemented using the four-quadrant inverse tangent: ϕ¯=atan2ssinϕs,scosϕs, where s is the number of values per voxel to average over.

5.3. S0 and T2 mapping: double-echo EPI

Images were motion corrected first (McFLIRT, FSL) and then read into MATLAB to calculate S0(τj) and T2(τj) maps as described above. The first 10 TRs were omitted to account for non-steady state effects. We segmented each cardiac cycle into 10 bins. Registration to the single-echo EPI was performed and the transformation was applied to the S0(τj) and T2(τj) maps. To synchronize cycles' τj across subjects for group averaging, we shifted each subject's τj such that the minimum S0 in the arterial mask (min[S0(τj)]) occurred in the first bin (j = 1).

6. Results

6.1. Single-echo EPI: Cardiac phase mapping

A representative subject's cardiac phase distributions in the different tissue masks is displayed in Fig. 4, and the group average distributions are given in Fig. 5. The arterial mask distribution has a distinct peak, which we used for inter-subject phase synchronisation by aligning it to ϕ = −π/2 (Fig. 4, Fig. 5). This dominant arterial phase is present in grey as well as white matter, but is more distinct in the former. A second phase peak is observed in all masks, suggesting a bi-modal distribution of temporal shifts in the timing of cardiac cycle-induced signals. A control distribution, derived from a whole brain mask at a low power frequency component is plotted in Fig. 4, Fig. 5d). It demonstrates that the phase mapping algorithm produces a uniform distribution as expected for a frequency that only contributes to the noise floor.

Fig. 5
a)–c) Group averaged distributions of the cardiac frequency phase shift in the EPI time series in the arterial, grey and white matter masks. d) is a whole brain control distribution at a non-significant frequency component, to confirm that the ...

The tissue-specific distributions do not yield any information about spatial variations in the shifts. We therefore calculated the phase distributions in MNI space for the cortex and subcortical areas. Fig. 6 shows the tissue-specific grey matter phase distribution split into the area-specific distributions in the cortex, putamen, caudate, thalamus and pallidum (all histograms are displayed normalised). The distribution in the cortex (red) shows a single peak slightly shifted to the left of the dominant arterial phase value of – π/2 but the second bi-phasal peak from the overall grey matter distribution (light blue) is absent. However, the subcortical areas show distinct shifts away from the dominant arterial phase, indicating that the cardiac cycle-induced signals in these areas have offset timings to the signals in the arterial regions and the cortex. In summary Fig. 6 shows how the overall bi-modal grey matter phase distribution is the sum of area-specific shifted single-peak distributions.

Fig. 6
Area-specific distribution of cardiac phase values in MNI space. The light blue distribution on the top is the overall grey matter phase distribution that can be split into the cortical (red), putamen (pink), caudate (yellow), thalamus (green) and pallidum ...

6.2. Double-echo EPI: S0 and T2 mapping

6.2.1. Validity of mono-exponential vs. multi-exponential signal model

The simulated multi-exponential decay curves of a grey matter voxel with partial volume fractions of vGM = 0.5,0.9 and 0.94 are displayed in Fig. 7a) together with the mono-exponential signals (dashed lines) that result from a fitting of S0 and T2 using two echo times (TE = [21.4,52.1]ms). The deviation between Smulti and Smono (Fig. 7b)) increases with partial voluming as expected. This discrepancy causes an apparent S0 change if T2 fluctuates. Fig. 8a) shows the amount of S0 fit error that is introduced by a ±1% change in the grey matter T2, assuming noise-free signals. For low partial volume mixing, vGM = 0.9 (pink plot), the error is smaller than ±0.1% and for vGM = 0.94 it is negligible. This coupling effect is enhanced in the presence of noise. We found the grey matter mask SNR for the first echo to be 205 and 61 for the second echo, respectively. Fig. 8b) shows an exemplar plot of the error fluctuation when S0 is fit from 50 noisy double echo signals, taking the aforementioned SNR values into account. Besides the model error in Fig. 8a) there is an additional variation in the apparent S0 fluctuation from the noisy data. For low partial volume mixings of vGM>0.9 (pink and blue plots) the apparent S0 fluctuates below ±0.4%.

Fig. 7
a) Multi-exponential decays (solid lines) are plotted for a grey matter voxel with partial volumes of 0.5,0.9 and 0.94 (see Table 1 for all parameters). The dotted lines are the decay curves when a mono-exponential curve is fitted at echo times ...
Fig. 8
An apparent S0 fluctuation is introduced by a T2 fluctuation when the signal is approximated by a mono-exponential model. a) plots the S0 error for a variation in grey matter T2 of ±1% assuming noise-free data. b) plots the enhanced ...

6.2.2. Experimental S0 and T2 mapping

The fast single-echo EPI data revealed a variation in the timing of cardiac cycle-induced EPI signals. With the voxel-wise S0(τj) and T2(τj) waveforms from the double-echo EPI we can look into the specific timing of these parameters for voxels with different cardiac phase values ϕ. Fig. 9 shows the phase distribution from the single-echo data and the peak-to-peak amplitude fluctuations S0,p and T2,p plotted on top (derived from the double-echo data). S0,p is strongest at the dominant arterial phase value of – π/2 and becomes essentially unmeasurable at the second phase peak around + π/2. Interestingly the peak-to-peak T2,p fluctuation is persistent for all voxels and ranges between 0.3 and 0.55%. The detailed pulsatility waveforms (as demeaned percentage fluctuations) ΔS0(τj) and ΔT2(τj) at the dominant arterial phase value of – π/2 and the second peak value of +π2 are shown on the bottom of Fig. 9.

Fig. 9
a) The grey matter distribution of cardiac phase values (blue histogram, derived from the single-echo EPI data) with the peak-to-peak percentage fluctuation in S0 (green) and T2 (red) plotted on top. For the voxels at the dominant arterial phase ...

Three exemplar single-subject pulsatility ΔS0(τj) and ΔT2(τj) waveforms at the dominant phase value – π/2 are displayed in Fig. 10. From the simulation study of T2 and S0 coupling we anticipate an apparent fluctuation in S0 with T2 fluctuations. In the group averaged waveforms in Fig. 10b) and c), as well as the single-subject examples we see that S0 and T2 are neither fluctuating in phase, nor are they anti-cyclic. Generally we observed a range of shifts between ΔS0(τj) and ΔT2(τj) for different phase ranges, which can be seen best in the animated version of Fig. 9 which increments through the entire phase range (provided as in the Supplementary Materials). The single subject S0 fluctuation magnitude of ±1% further exceeds the estimated error that is introduced by a T2 fluctuation. However, the error bars on the S0 waveforms of 0.5% are on the order of the estimated apparent S0. Overall, the S0 fluctuation at the dominant arterial peak resembles characteristics of a typical arterial systole-to-diastole waveform (compare for example with the pulse-ox signal in Fig. 2), with a steep rise towards systole and a slower decrease towards diastole.

Fig. 10
Single-subject examples of the pulsatility ΔS0(τj) and ΔT2(τj) waveforms (plotted twice for better visualisation) from voxels that show a single-echo EPI fluctuation at the dominant arterial phase value (− ...

Supplementary videos related to this article can be found at

The following is the supplementary data related to this article:

Video 1:

Animated plot of the grey matter distribution of cardiac phase values (blue histogram, derived from the single-echo EPI data). The peak-to-peak percentage fluctuation in S0 (green) and T2 (red) is plotted in the same graph. The animation flicks through the difference phase ranges and displays the corresponding ΔS0 and ΔT2 fluctuation waveforms (derived from the double-echo data) of the voxels that had cardiac phase values in that range.

Click here to view.(1.6M, mp4)Video 1

White matter results are plotted in Fig. 11. Pulsatility in S0 and T2 are reduced compared to grey matter. Similarly the S0 pulsatility is maximised around – π/2 and T2 peak-to-peak fluctuations are around 0.2%. The detailed ΔS0(τj) and ΔT2(τj) waveforms for all phases can be clicked through in the animated version of Fig. 11 which is provided as in the Supplementary Materials).

Fig. 11
a) The white matter distribution of cardiac phase values (blue histogram, derived from the single-echo EPI data) with the peak-to-peak percentage fluctuation in S0 (green) and T2 (red) plotted on top. For the voxels at the dominant arterial phase ...

Supplementary video related to this article can be found at

The following is the supplementary data related to this article:

Video 2:

Animated plot of the white matter distribution of cardiac phase values (blue histogram, derived from the single-echo EPI data). The peak-to-peak percentage fluctuation in S0 (green) and T2 (red) is plotted in the same graph. The animation flicks through the difference phase ranges and displays the corresponding ΔS0 and ΔT2 fluctuation waveforms (derived from the double-echo data) of the voxels that had cardiac phase values in that range.

Click here to view.(1.6M, mp4)Video 2

The average T2 values in the grey and white matter masks, as well as the peak-to-peak fluctuations S0,p and T2,p at the dominant arterial peak and the second peak, are summarised in Table 2. The group averaged value of T2 within the grey matter mask was 29.2±1.2 ms and within the white matter mask was 24.1±0.8 ms. A strict comparison of the absolute T2 values with literature values is not straightforward due to multiple differences in the acquisition (resolution, pulse sequence etc.), however our values fall into the expected range from previously reported studies (see the list of values in Table 2).

Table 2
Overview of results. The left hand of the table summarises the group averaged T2 values in the grey and white matter masks (thresholded at a PVE > 0.9) from our double-echo data with literature values for comparison (where ...

7. Discussion

7.1. Cardiac phase distribution

Cardiac cycle-induced shifts in the EPI time series showed a bi-phasic distribution in all subjects. A dominant phase mode was identified and appears to originate in the large arterial structures. When disentangling the phase distributions into area-specific distributions we found that the subcortical areas are shifted with respect to the cortex and explain the overall distribution in grey matter. This lag phenomenon could have potentially similar origins as the spatio-temporal variations in the neuro-vascular haemodynamic response function that varies with cortical depth and has been associated with variations in microvascular volumes (Zhao et al., 2006, Puckett et al., 2016) and their timing of changes in blood flow velocities and volume. A better understanding of an equivalent cardiac haemodynamic response function could be important in physiological cleaning of fMRI data.

An alternative origin of the phase shifts is the rigid and non-rigid brain motion over the cardiac cycle. However, as we saw in the example masks of significantly cardiac cycle-modulated voxels (Fig. 3) we notice two things about the spatial distribution of these masks. Firstly, only a fraction of all brain voxels show a significant signal modulation and secondly, the spatial prevalence is close to larger arterial structures and perivascular CSF spaces. As such, we rule out that rigid brain motion was a major source of EPI signal fluctuations. However, non-rigid motion propagating from the brainstem outwards could be an alternative explanation. It has been measured to be delayed on average in the peripheral lobes by about 200–300 ms (Zhong et al., 2009). The peaks in the cortical and subcortical phase distributions are approximately spread apart by π/2 and even more in the thalamus which is closest to the brainstem (compare the plots in Fig. 6). This phase difference translates into about 250 ms for a standard heart beat of 1 s. On the contrary the cortical phases closely match the dominant arterial phase and the brainstem is known to exhibit the highest displacement at the time of bulk arterial blood inflow. Further cortical voxels show more significant cardiac cycle related EPI signal modulations than subcortical areas and white matter. We conclude that the cardiac cycle-induced EPI signal timing in the cortex and most subcortical areas is dominated by vascular changes and not brain motion.

7.2. S0 and T2 fluctuations over the cardiac cycle

7.2.1. Validity of mono-exponential vs. multi-exponential signal model

The voxel's EPI signal decay is the sum of multiple exponential decays from its tissue composition. A mono-exponential approximation allows the fitting of S0 and T2 parameters to study effects of inflow and partial volume fluctuations over the course of the cardiac cycle. If the voxel has minimal partial volume mixing (here we used vGM > 0.9 in simulations and a PVE > 0.9 to threshold tissue masks) the fitting error in S0 should be negligible for noise-free data. However, taking the SNR of the two echoes and the number of measurement points into account we estimated that apparent S0 fluctuations could reach levels of up to 0.4% for a grey matter T2 fluctuation of ±1%. This was in accordance with the size of the single-subject S0 error bars in a single cardiac cycle bin. Additional inaccuracies in T2 that subsequently couple inversely into S0 could arise from local variations in echo times that are due to small off-resonance effects, particularly around the orbitofrontal lobe (Deichmann et al., 2002).

7.2.2. S0 and T2 mapping

Voxels with significant cardiac modulation showed distinct S0 and T2 waveforms over the cardiac cycle, with the former resembling characteristic features of the arterial pressure wave with a steep increase to systole and a flatter decrease towards diastole. At the dominant arterial peak we found a slightly offset anti-cyclic behaviour between S0 and T2. Bianciardi et al. (Bianciardi et al., 2016) reported a similar timing relationship between the two parameters. They measured S0 and T2 waveforms in a region of interest in the carotid artery. Our scans did not cover as far down as the carotid, but it is a sensible assumption that the carotid area has a cardiac phase value at the dominant arterial phase and would show a comparable combination of S0 and T2 waveforms.

Generally, S0 peak-to-peak fluctuations vary strongly across the cardiac phase distribution, with the highest values being at the dominant arterial peak. S0 pulsatility is expected to be a proxy of CBV and CBFV pulsatility and its decline is consistent with the finding that subcortical areas showed shifted cardiac phase distributions. These areas are less vascularised, and have slower blood flow velocities and a reduced inflow effect is expected. Similar arguments apply to the reduced S0 pulsatility in white matter which is greatly reduced compared to grey matter.

A persistent T2 fluctuation was observed for all voxels that had a significantly cardiac cycle-modulated EPI time series. Oxygenation variations affect T2, however to our knowledge oxygen extraction does not vary over the course of the cardiac cycle as baseline metabolic extraction should be constant. An alternative explanation could be partial volume fluctuations. It is thought that the sum of all brain tissue and fluid compartments is constant at all times. Blood and CSF volume fluctuations have to be counterbalanced by an inverse fluctuation in one or both compartments (this is often referred to as the “Monro-Kellie doctrine” (Mokri, 2001). At the voxel level an increase in arteriolar blood volume during systole has to coincide with an outflow in either CSF or venous blood. Even small reductions in CSF volume during systole, and subsequent inflow towards diastole could explain the observed apparent T2 waveforms, as CSF's T2 greatly exceeds the T2 value of all other voxel compartments. Venous blood T2 is the shortest of all components in the voxel and a reduction in venous partial volume during a systolic increase in arteriolar blood volume should prolong the voxel's apparent T2. This contradicts the observed reduction in the voxel's apparent T2 during the systolic increase in S0. Recent glymphatic system research suggests that CSF enters the brain via the spaces surrounding the penetrating arteries (Iliff et al., 2012), moving down paravascular pathways driven by cardiac pulsation (Iliff et al., 2013). Changes in CSF partial volume during hypercapnia (Scouten and Constable, 2008) and neuronal activation (Scouten and Constable, 2007), both of which increase microvascular volumes, have further been reported in VASO-fMRI data. If CSF and blood volume fluctuations are the dominant cause of the T2 waveforms, the relative timing of S0 and T2 remains non-intuitive. The slight offset in timing between systolic S0 increase and T2 decrease suggests that arteriolar blood volume changes are delayed relative to the inflow effect. The T2 fluctuations were also observed in white matter, where we would not expect CSF, however interstitial fluid could similarly serve as an exchanging volume compartment.

Adams et al. (2017) (measured volumetric strain in grey and white matter at 7T over the course of the cardiac cycle using displacement encoding with stimulated echoes. They detected an increase in strain in both grey and white matter, with a larger effect in the former. Relating the effect to “microvascular swelling” their measurements potentially share a similar or the same phenomenon that underlies the grey and white matter T2 waveforms that we detected in our EPI data. Ultimately, more research is needed to investigate the various impacts of the cardiac cycle on MR signal parameters in the brain.

8. Conclusion

Cardiac cycle-induced EPI signals are temporally shifted from the cortex, that closely matches the timing of signals in arterial structures, towards subcortical areas that appear delayed. These shifts are associated with different inflow S0 and apparent T2 pulsatility over the course of the cardiac cycle. Variations in microvascular density and flow velocity are likely the cause of the variation in S0. T2 fluctuations are likely a result of partial volume fluctuations between blood and CSF or interstitial fluid volumes.


This work was supported by the Initial Training Network, HiMR funded by the FP7 Marie Curie Actions of the European Commission (FP7-PEOPLE-2012-ITN-316716). We thank the NIHR Oxford Biomedical Research Centre for support. The Wellcome Centre for Integrative Neuroimaging is supported by core funding from the Wellcome Trust (203139/Z/16/Z). We also thank the Dunhill Medical Trust for support. We thank the University of Minnesota Center for Magnetic Resonance Research for the provision of the multiband EPI sequence software.


Appendix ASupplementary data related to this article can be found at

Mathematical symbols

cardiac frequency
equilibrium magnetisation
longitudinal magnetisation
number of bins
error probability
reference voxel index
time-domain signal
frequency-domain signal
mono-exponential time-domain signal
multi-exponential time-domain signal
signal amplitude extrapolated to TE = 0
normalized peak signal amplitude
longitudinal relaxation time
effective transverse relaxation time
normalized peak effective transverse relaxation time
echo time
inversion time
repetition time
partial volume
flip angle
signal delay
proton spin density
cardiac unit cycle
systole-to-systole time period
cardiac cycle – related phase shift
mean value of x

Appendix A. Supplementary data

The following are the supplementary data related to this article:

Click here to view.(949 bytes, zip)


Adams A., Luijten P.R., Zwanenburg J.J.M. Proceedings of the International Society of Magnetic Resonance Imaging in Medicine. 2017. Measurements of cardiac related pulsatile volumetric strain in grey and white matter brain tissue with high resolution DENSE at 7T; p. 704.
Asmar R., Benetos A., Topouchian J., Laurent P., Pannier B., Brisac A.-M., Target R., Levy B.I. Assessment of arterial distensibility by automatic pulse wave velocity measurement validation and clinical application studies. Hypertension. 1995;26(3):485–490. [PubMed]
Bianciardi M., Toschi N., Polimeni J.R., Evans K.C., Bhat H., Keil B., Rosen B.R., Boas D.A., Wald L.L., Bianciardi M., Toschi N., Polimeni R., Evans K.C., Bhat H., Keil B., Rosen B.R., Boas D.A., Wald L.L. The pulsatility volume index: an indicator of cerebrovascular compliance based on fast magnetic resonance imaging of cardiac and respiratory pulsatility. Philos. Trans. Royal Soc. A: Math. Phys. Eng. Sci. 2016;374(2067):1–19. [PMC free article] [PubMed]
Boas D.A., Jones S.R., Devor A., Huppert T.J., Dale A.M. A vascular anatomical network model of the spatio-temporal response to brain activation. NeuroImage. 2008;40(3):1116–1129. [PubMed]
Bouvy W.H., Geurts L.J., Kuijf H.J., Luijten P.R., Kappelle L.J., Biessels G.J., Zwanenburg J.J.M. Assessment of blood flow velocity and pulsatility in cerebral perforating arteries with 7-T quantitative flow MRI. NMR Biomed. 2016;29(9):1295–1304. [PubMed]
Dagli M.S., Ingeholm J.E., Haxby J.V. Localization of cardiac-induced signal change in fMRI. NeuroImage. 1999;9(4):407–415. [PubMed]
Deichmann R., Josephs O., Hutton C., Corfield D., Turner R. Compensation of susceptibility-induced BOLD sensitivity losses in echo-planar fMRI imaging. NeuroImage. 2002;15(1):120–135. [PubMed]
Feinberg D., Setsompop K. Ultra-fast MRI of the human brain with simultaneous multi-slice imaging. J. Magn. Reson. 2013;229:90–100. [PubMed]
Feinberg D., Moeller S., Smith S.M., Auerbach E., Ramanna S., Glasser M.F., Miller K.L., Ugurbil K., Yacoub E. Multiplexed echo planar imaging for sub-second whole brain FMRI and fast diffusion imaging. PLoS One. 2010;5(12):e15710. [PubMed]
Grinstead J.W., Rooney W., Laub G. 2010. The origins of bright blood MPRAGE at 7 Tesla and a simultaneous method for T1 imaging and non-contrast MRA. (Proceedings of the International Society of Magnetic Resonance in Medicine). 1429.
Herscovitch P., Raichle M.E. What is the correct value for the brain–blood partition coefficient for water? J. Cereb. Blood Flow Metab. 1985;5(1):65–69. [PubMed]
Hu X., Le Tuong Huu, Parrish T., Erhard P. Retrospective estimation and correction of physiological fluctuation in functional MRI. Magn. Reson. Med. 1995;34(2):201–212. [PubMed]
Iliff J.J., Wang M., Liao Y., Plogg B.A., Peng W., Gundersen G.A., Benveniste H., Vates G.E., Deane R., Goldman S.A., Nagelhus E.A., Nedergaard M. A paravascular pathway facilitates CSF flow through the brain parenchyma and the clearance of interstitial solutes, including amyloid. Sci. Transl. Med. 2012;4(147) 147ra111–147ra111. [PMC free article] [PubMed]
Iliff J.J., Wang M., Zeppenfeld D.M., Venkataraman A., Plog B.A., Liao Y., Deane R., Nedergaard M. Cerebral arterial pulsation drives paravascular CSF-interstitial fluid exchange in the murine brain. J. Neurosci. 2013;33(46):18190–18199. [PubMed]
Ivanov D., Schäfer A., Deistung A., Streicher M.N., Kabisch S., Henseler I., Roggenhofer E., Jochimsen T.H., Reichenbach J., Uludag K., Turner R. In vivo estimation of the transverse relaxation time dependence of blood on oxygenation at 7 tesla. Proc. Int. Soc. Magn. Reson. Med. 2013;21:2472.
Jenkinson M., Beckmann C.F., Behrens T.E.J., Woolrich M.W., Smith S.M. Fsl. NeuroImage. 2012;62(2):782–790. [PubMed]
Kiviniemi V., Wang X., Korhonen V., Keinänen T., Tuovinen T., Autio J., LeVan P., Keilholz S., Zang Y.-F., Hennig J. Ultra-fast magnetic resonance encephalography of physiological brain activity–glymphatic pulsation mechanisms? J. Cereb. Blood Flow Metab. 2016;36(6):1033–1045. [PubMed]
Marques J.P., Bowtell R. Using forward calculations of the magnetic field perturbation due to a realistic vascular model to explore the BOLD effect. NMR Biomed. 2008;21(6):553–565. [PubMed]
Meckel S., Leitner L., Bonati L.H., Santini F., Schubert T., Stalder A.F., Lyrer P., Markl M., Wetzel S.G. Intracranial artery velocity measurement using 4D PC MRI at 3 T: Comparison with transcranial ultrasound techniques and 2D PC MRI. Neuroradiology. 2013;55(4):389–398. [PubMed]
Moeller S., Yacoub E., Olman C., Auerbach E., Strupp J., Harel N., Ugurbil K. 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. 2010;63(5):1144–1153. [PubMed]
Mokri B. The Monro–Kellie hypothesis applications in CSF volume depletion. Neurology. 2001;56(12):1746–1748. [PubMed]
Obata T., Liu T.T., Miller K.L., Luh W.M., Wong E.C., Frank L.R., Buxton R.B. Discrepancies between BOLD and flow dynamics in primary and supplementary motor areas: application of the balloon model to the interpretation of BOLD transients. NeuroImage. 2004;21(1):144–153. [PubMed]
Peters A.M., Brookes M.J., Hoogenraad F.G., Gowland P.A., Francis S.T., Morris P.G., Bowtell R. T2* measurements in human brain at 1.5, 3 and 7 T. Magn. Reson. Imaging. 2007;25(6):748–753. [PubMed]
Pohmann R., Speck O., Scheffler K. Signal-to-noise ratio and MR tissue parameters in human brain imaging at 3, 7, and 9.4 tesla using current receive coil arrays. Magn. Reson. Med. 2016;75(2):801–809. [PubMed]
Poncelet B.P., Wedeen V.J., Weisskoff R.M., Cohen M.S. Brain parenchyma motion: measurement with cine echo-planar MR imaging. Radiology. 1992;185:645–651. [PubMed]
Puckett A.M., Aquino K.M., Robinson P.A., Breakspear M., Schira M.M. The spatiotemporal hemodynamic response function for depth-dependent functional imaging of human cortex. NeuroImage. 2016;139:240–248. [PubMed]
Robertson A.D., Tessmer C.F., Hughson R.L. Association between arterial stiffness and cerebrovascular resistance in the elderly. J. Hum. Hypertens. 2010;24(3):190–196. [PubMed]
Rooney W.D., Johnson G., Li X., Cohen E.R., Kim S.-G., Ugurbil K., Springer C.S. Magnetic field and tissue dependencies of human brain longitudinal 1H2O relaxation in vivo. Magn. Reson. Med. 2007;57(2):308–318. [PubMed]
Rostrup E., Knudsen G.M., Law I., Holm S., Larsson H.B.W., Paulson O.B. The relationship between cerebral blood flow and volume in humans. NeuroImage. 2005;24(1):1–11. [PubMed]
Scouten A., Constable R. Applications and limitations of whole-brain MAGIC VASO functional imaging. Magn. Reson. Med. 2007;58(2):306–315. [PubMed]
Scouten A., Constable R.T. VASO-based calculations of CBV change: accounting for the dynamic CSF volume. Magn. Reson. Med. 2008;59(2):308–315. [PubMed]
Setsompop K., Feinberg D.A., Polimeni J.R. Rapid brain MRI acquisition techniques at ultra-high fields. NMR Biomed. 2016;29(9):1198–1221. [PubMed]
Soellinger M., Rutz A.K., Kozerke S., Boesiger P. 3D cine displacement-encoded MRI of pulsatile brain motion. Magn. Reson. Med. 2009;61(1):153–162. [PubMed]
Tong Y., Frederick B. deB. Studying the spatial distribution of physiological effects on BOLD signals using ultrafast fMRI. Front. Hum. Neurosci. 2014;8:1–8. [PubMed]
Uludag K., Müller-Bierl B., Ugurbil K. An integrative model for neuronal activity-induced signal changes for gradient and spin echo functional imaging. NeuroImage. 2009;48(1):150–165. [PubMed]
Van De Vosse F.N., Stergiopulos N. Pulse wave propagation in the arterial tree. Annu. Rev. Fluid Mech. 2011;43:467–499.
Wagshul M.E., Eide P.K., Madsen J.R. The pulsating brain: a review of experimental and clinical studies of intracranial pulsatility. Fluids Barriers CNS. 2011;8(1):5. [PubMed]
Webb A.J.S., Simoni M., Mazzucco S., Kuker W., Schulz U., Rothwell P.M. Increased cerebral arterial pulsatility in patients with leukoaraiosis: arterial stiffness enhances transmission of aortic pulsatility. Stroke. 2012;43(10):2631–2636. [PubMed]
Xu J., Moeller S., Auerbach E.J., Strupp J., Smith S.M., Feinberg D.A., Yacoub E., Uğurbil K. Evaluation of slice accelerations using multiband echo planar imaging at 3T. Neuroimage. 2013;83:991–1001. [PubMed]
Zhao F., Wang P., Hendrich K., Ugurbil K., Kim S.G. Cortical layer-dependent BOLD and CBV responses measured by spin-echo and gradient-echo fMRI: insights into hemodynamic regulation. NeuroImage. 2006;30(4):1149–1160. [PubMed]
Zhao J.M., Clingman C.S., Närväinen M.J., Kauppinen R.A., van Zijl P. Oxygenation and hematocrit dependence of transverse relaxation rates of blood at 3T. Magn. Reson. Med. 2007;58(3):592–597. [PubMed]
Zhong X., Meyer C.H., Schlesinger D.J., Sheehan J.P., Epstein F.H., Larner J.M., Benedict S.H., Read P.W., Sheng K., Cai J. Tracking brain motion during the cardiac cycle using spiral cine-DENSE MRI. Med. Phys. 2009;36(8):3413–3419. [PubMed]